SERA
United States
Duluth MN 55804
EPA 600 3 80 047
M;IV 1980
Research and Development
A Two-Mode
Free-Surface
Numerical Model for
the Three-
Dimensional Time-
Dependent Currents in
Large Lakes
-------
RESEARCH REPORTING SERIES
Research reports.of the Office of Research and Development, U.S. Environmental
Protection Agency, have been grouped into nine series. These nine broad cate-
gories were established to facilitate further development and application of en-
vironmental technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The nine series are:
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
6. Scientific and Technical Assessment Reports (STAR)
7. Interagency Energy-Environment Research and Development
8. "Special" Reports
9. Miscellaneous Reports
This report has been assigned to the ECOLOGICAL RESEARCH series. This series
describes research on the effects of pollution on humans, plant and animal spe-
cies, and materials. Problems are assessed for their long- and short-term influ-
ences. Investigations include formation, transport, and pathway studies to deter-
mine the fate of pollutants and their effects. This work provides the technical basis
for setting standards to minimize undesirable changes in living organisms in the
aquatic, terrestrial, and atmospheric environments.
This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.
-------
EPA-600/3-80-047
May 1980
A TWO-MODE FREE-SURFACE NUMERICAL MODEL
FOR THE THREE-DIMENSIONAL TIME-DEPENDENT
CURRENTS IN LARGE LAKES
by
Y. Peter Sheng and Wilbert Lick
Case Western Reserve University
Cleveland, Ohio 44106
Grant No. R-803704
PROJECT OFFICER
David M. DoIan
Large Lakes Research Station
Environmental Research Laboratory-Duluth
Grosse lie, Michigan 48138
ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
DULUTH, MINNESOTA 55804
-------
DISCLAIMER.
This report has been reviewed by the Environmental Research Laboratory,
U.S. Environmental Protection Agency, and approved for publication. Approval
does not signify that the contents necessarily reflect the views and policies
of the U.S. Environmental Protection Agency, nor does mention of trade names
or commercial products constitute endorsement or recommendation for use.
ii
-------
FOREWORD
The present report is part of a series of reports on research
concerned with understanding more thoroughly and being able to predict
the transport and fate of contaminants in lakes. This transport and
fate is a complex matter which involves physical, chemical, and bio-
logical processes, all of which need to be considered before the fate
of contaminants can be predicted.
This work has emphasized (1) the transport of contaminants throughout
the lake, especially in the near-shore, and (2) sediment-water inter-
actions including the biochemical transformations in the sediments and
the physical processes of sediment resuspension and deposition. Earlier
work was concerned mainly with these processes in a thermally non-
stratified lake. The later work has concentrated more on the significant
processes in a stratified lake.
J. David Yount, Ph.D.
Deputy Director
Environmental Research Laboratory-Duluth
iii
-------
ABSTRACT
A two-mode, free^surface model based on vertically-stretched coordinates
and a vertically-implicit scheme has been developed and applied to Lake Erie
under non-stratified conditions. A brief description of the general equa-
tions and boundary conditions is first given. The detailed equations for a
two-mode, free-surface model are then described. Finite-difference proce-
dures including the finite^difference equations are listed in detail. The
model is first applied to idealized basins and then to Lake Erie for two
wind conditions; (1) a uniform wind suddenly imposed, and (2) an actual
wind variable in both space and time. The results of the computations are
compared with the results from a single-mode, free-surface model and from a
rigid-lid model.
iv
-------
CONTENTS
Foreword ill
Abstract.... iv
Figures vi
Acknowledgments viii
1. Introduction 1
2. Conclusions and Recommendations 3
3. Equations and Boundary Conditions 5
4. General Description of Numerical Model 12
5. Applications to Idealized Basins 20
6. Applications to Lake Erie 26
References 54
Appendices
A. Derivation of the vertically stretched equations 57
B. An Implicit Scheme 60
-------
FIGURES
Number Page
3.1 Cartesian coordinates located at the nominal lake surface 6
4.1 Basic numerical grid structure 14
4.2 Schematics of computational procedure 19
5.1 A two-dimensional basin. 21
5.2 Integrated current and near-bottom current, and vertical
profile of horizontal currents at point 1 of Fig. 5.1 22
5.3 Integrated current and near-bottom current, and vertical
profile of horizontal currents at point 2 of Fig. 5.1 23
5.4 Results in a constant-depth square basin 25
6.1 Lake Erie bottom topography 27
6.2 Numerical grid structure, locations of wind stations, water
level gauges, and bottom topography of Western Lake Erie 28
613 Surface elevation computed by free-surface model as a function
of time at three locations in Lake Erie 30
6.4 Horizontal velocities at a constant 0.1 and 0.5 depth
from surface after 3 hours 31
6.5 Horizontal velocities at a constant 0.1 and 0.5 depth
surface after 24 hours 33
6.6 Wind speed as a function of time - 4 weather stations - Western
Basin, wind stresses at Toledo and Point Pelee 34
6.7 Comparison of surface elevation computed by free surface
model with water level measurement at 5 gauges 35
6.8a Horizontal velocities at 0.1 and 0.5 depth as a function
of time at station 1 (near Toledo) 36
6.8b Horizontal velocities at 0.1 and 0.5 depth as function
of time at station 2 (near Point Pelee)., ^,t 37
vi
-------
6.9 Horizontal velocities at constant 0.1 and 0.5 depth from
surface at 0600 3/7/76 39
6.10 Horizontal velocities at constant 0.1 and 0.5 depth from
surface at 1200 3/7/76 40
6.11 Horizontal velocities at constant 0.1 and 0.5 depth from
surface at 0000 3/8/76 41
6.12 Horizontal velocities at constant 0.1 and 0.5 depth from
surface at 1200 3/8/76 42
6.13 Large-scale eddy observed by ERTS-1 (LANDSAT-1) off the
tip of Point Pelee on 3/14/73 43
6.14 Horizontal velocities at constant 0.1 and 0.5 depth from
surface at 1200 3/8/76 44
6.15 Horizontal velocities at constant 0.1 and 0.5 depth from
surface at 0000 3/9/76 45
6.16 Horizontal velocities at 0.1 and 0.5 depth as function of
time at station 1 (near Toledo) 46
6.17 Horizontal velocities at 0.1 and 0.5 depth as function of
time at station 2 (near Point Pelee) 47
6.18 Time-averaged currents over 3 14-hour periods and the 42-hour
period starting 0000 3/7/76 station 1 (near Toledo) 49
6.19 Time-averaged currents over 3 14-hour periods and the 42-hour
period starting 0000 3/7/76 station 2 (near Point Pelee) 50
6.20 Daily-averaged horizontal velocities at constant 0.1 depth
from surface face on 3/7/76 51
6.21 Daily-averaged horizontal velocities at constant 0.1 depth
from surface on 3/8/76 52
vii
-------
ACKNOWLEDGMENTS
The writers wish to acknowledge with appreciation the support of this
project by the Environmental Protection Agency. Mr. Dave Dolan served as
Grant Project Officer.
The National Aeronautics and Space Administration, through Dr. Richard
Gedney, has given us considerable computer time for which we are grateful.
Mr. Frank Molls of NASA gave us considerable computer programming assistance.
viLi
-------
SECTION 1
INTRODUCTION
The need to quantitatively understand and assess the environmental im-
pact of power plant discharges, sewage outfalls, dredged material disposals,
and other contaminants has recently resulted in an increased effort to under-
stand the hydrodynamics of large lakes. Although some of the most important
aspects of wind-driven circulation in a large lake can be understood from the
analytical models of Ekman (1923) and Welander (1957), a detailed investiga-
tion of the currents in realistic situations often dictates the need of
numerical modeling.
In order to mathematically model and hence predict the dispersion and
fate of substances introduced into lakes, it is essential to be able to
accurately model the time-dependent currents in the lake. For accurate pre-
diction of the entire velocity field, three-dimensional, time-dependent
numerical models are often needed and numerous models of this type have been
developed. Reviews by Cheng et al. (1976) and Lick (1976a,b) have discussed
several of these three-dimensional models.
The simplest kind of three-dimensional model is the single-mode, free-
surface model which solves the equations of motion and continuity equation
with a single numerical time step. Such models are quite general and have
been applied to large lakes by Freeman et al. (1972), Katz and Kizlauskaus
(1973), Haq et al. (1974), Sheng (1975), and Sheng and Lick (1976). These
models consume a large amount of computer time since the maximum allowable
numerical time step is limited by the time it takes a surface gravity wave to
travel between two adjacent horizontal grid points, a very severe limitation.
In order to improve the numerical efficiency of the single-mode, free-
surface model, rigid-lid models have been formulated and applied to large
lakes by Liggett (1969), Paul and Lick (1974, 1979), and Bennett (1977). In
this type of model, it is assumed that the vertical velocity at the undis-
turbed location of the free-surface is zero. This eliminates the motions and
time scales associated with surface gravity waves and hence allows much
larger numerical time steps and reduced computational times.
Berdahl (1971) and Crowley (1969, 1970) compared the rigid-lid model
and single-mode, free-surface models in idealized oceanic basins and found
the results for Rossby waves predicted by the two models were very different,
although the rigid-lid model required much less computational time. Sheng
et al. (1978) compared the rigid-lid and free-surface models for Lake Erie.
They found that, for relatively short time intervals and in shallow waters,
significant differences between the results of the two models occur. The
-------
long-term, time-averaged circulations computed by the two models agreed well
in periods of strong winds and active seiching.
Simons (1972, 1974, 1975) developed a two-mode, free-surface model for
Lake Ontario based on a four-layer lattice. In this type of model, the ex-
ternal variables, i.e., the free-surface elevation and vertically integrated
currents, are treated separately from the internal, three-dimensional flow
variables. The free-surface gravity waves no longer limit the internal-mode
calculations and much larger time steps may be used for the internal mode cal-
culations with greatly reduced overall computational times.
One shortcoming of such a multi-layered model is the lack of vertical
resolution in the shallow parts of the basin. In order to have the same
vertical resolution in both the shallow and deeper parts of the basin, verti-
cally stretched coordinates (Philips 1957, Freeman et al. 1972), which result
in an equal number of uniformly-spaced vertical grid points over the entire
basin, can be used.
Since Lake Erie is the shallowest of the Great Lakes, the effects of
vertical diffusion and bottom shear stress are more pronounced than in other
Great Lakes. The numerical stability criterion associated with the vertical
diffusion terms in the governing equations is also much more stringent.
Vertically implicit schemes can be used to resolve this problem.
In the present study, a two-mode, free-surface model based on vertical-
ly-stretched coordinates and a vertically-implicit scheme has been developed
and applied to Lake Erie under non-stratified conditions. In the following,
a brief description of the general equations and boundary conditions is
first given. The detailed equations for a two-mode, free-surface model in
non-stratified situations are then described. Finite-difference procedures
including the finite-difference equations are listed in detail. The model
is first applied to idealized basins and then to Lake Erie for two wind con-
ditions: (1) a uniform wind suddenly imposed, and (2) an actual wind vari-
able in both space and time. The results of the computations are compared
with the results from a single-mode, free-surface model and from a rigid-lid
model.
-------
SECTION 2
CONCLUSIONS AND RECOMMENDATIONS
In the present study, a two-mode, free-surface numerical model capable
of realistically predicting currents in lakes under non-stratified conditions
has been developed. Vertically stretched coordinates have been used so as to
have the same vertical resolution in all parts of the basin. A vertically
implicit scheme has been used to eliminate numerical instabilities because of
the small vertical distances between grid points. The model is operational
and calculations have been made for both idealized and realistic basins and
wind conditions with specific applications to Lake Erie. The present report
describes the general numerical procedure and presents representative results
of the calculations.
This two-mode, free-surface model consumes considerably less computer
time than does a single-mode, free-surface model. The present model requires
more computer time than does a rigid-lid model but is superior to the rigid-
lid model in that it does not eliminate currents associated with free-surface
oscillations.
In a previous report (Lick, 1976a), a survey of numerical models of
lake currents was given and several recommendations on further work were made.
These recommendations are presently being implemented and are still valid.
In particular, it is suggested that a great deal of numerical experimentation
be done using these models in order to: (1) Understand the general charac-
teristics of flows in lakes, especially in the near-shore under stratified
conditions, e.g., upwelling and mass flux through the thermocline; (2) Cal-
culate details of specific flows and verify these by means of field observa-
tions, again especially in the nearshore. A combination of field observa-
tions, numerical experiments, and laboratory work needs to be done to deter-
mine more accurately the magnitude and functional dependences of the vertical
and horizontal eddy diffusion coefficients, the wind stress-wind velocity
relation, the wind factor ratio, and the bottom stress.
Hydrodynamic models are becoming quite detailed. Sufficient informa-
tion on their use is becoming available so that confidence may be placed in
their prediction. In particular, these models are sufficiently accurate so
that they can be used as bases for models of sediment resuspension and
transport, contaminant dispersion during dredging, phytoplankton and zoo-
plankton growth, community succession, etc. and should be so used. The
present hydrodynamic model, along with a wave-hindcasting model, has been
used by the present authors in a study of the resuspension, deposition, and
transport of sediments in the Western Basin of Lake Erie (Sheng and Lick,
1977b).
-------
Hydrodynamic models generally have not considered the effects of waves,
effects which become increasingly important as the depth of water decreases.
For an adequate understanding of the dispersion of contaminants in these
shallow, near-shore regions, further work needs to be done on the effects of
waves. The direct effects of waves in causing contaminant transport and
sediment resuspension and transport as well as the indirect effects in causing
long-shore currents need to be investigated, made quantitative, applied to
specific cases, and verified by means of field observations.
The hydrodynamics and dispersion of contaminants in a lake, particularly
when the lake is thermally stratified, depend strongly on the turbulent pro-
duction, dissipation, and transport processes. In the present analysis and
in most analyses of the dynamics of a lake, the concept of an eddy coefficient
has been used to describe turbulence. Recently, more sophisticated descrip-
tions of turbulent transport processes have been developed (e.g. Kolmogorov,
1942; Harlow and Nakayama, 1967; Launder and Spalding, 1972; Donaldson, 1973;
Lewellen, Teske, and Sheng, 1978). These models attempt to describe the
dynamics of the turbulence more accurately by means of transport equations
for certain turbulent properties of the flow. A comparative study of eddy-
viscosity and these transport models as applied to lake dynamics would be
useful. However, a major difficulty with the models, especially with regard
to lake dynamics, is a deficiency in our knowledge of the sources of turbu-
lence and how to parameterize the magnitude of these sources. This is the
part of the investigation which needs to be emphasized.
-------
SECTION 3
EQUATIONS AND BOUNDARY CONDITIONS
The basic equations used in the modeling of lake currents are the
hydrodynamic equations of conservation of mass, momentum and energy plus an
equation of state:
3x
at
3t
3y ' 3z "
f\
13u . 3uv , 3uw
3x
3y
Juv +iv
3x 3y
3z
fu _ I J£ + JL
U p 3y 3x \
+ -
T 9y
jhr
3y
3y
3v
3z
3z
3t
= -Pg
13uT 3vT , 3wT
3x 3y 3z
3x
3z
3z
P = P(T)
(3.3)
(3.4;
(3.5)
(3.6)
where x and y are the horizontal coordinates, z is the vertical coordinate
pointing vertically upward to form a right-hand coordinate system with x and
y (Figure 3.1), u,v, and w are the three-dimensional velocities in the x,y,
and z directions, t is time, f is the Coriolis parameter equal to 10 sec~l,
g is the gravitational acceleration, p is the pressure, p is the density, T
is the temperature, AJJ and Kg are the horizontal eddy coefficients, and Ay
and Ky are the vertical eddy coefficients.
Several assumptions are implicit in these equations. These are; (a)
the pressure is assumed to vary hydros tatically in the z direction, (b) the
Boussinesq approximation is valid, i.e., the effects of density variations
are considered unimportant except in the bouyancy term, (c) eddy coefficients
are used to account for turbulent diffusion effects, and (d) heat sources and
sinks in the fluid are neglected.
At the free surface of the lake, the appropriate boundary conditions
are: (a) The wind stress is specified,
_w
3v
w
3z
(3.7)
-------
r DISPLACED
\ LAKE SURFACE
--LAKE BOTTOM
NOMINAL
LAKE SURFACE
Figure 3.1 Cartesian coordinates located at the nominal lake surface
-------
w w
where T and T are the wind stresses in the x and y directions respectively
and are functions of the wind velocity; (b) the.kinematic condition is sat-
isfied,
where £ is the elevation of the free-surface; (c) The dynamic condition is
satisfied,
P = Pa (3.9)
where p is the atmospheric pressure; and (d) The heat flux is specified,
c*
*V ~3z~ = qs = h (T"V (3.10)
where Te is the equilibrium air temperature at which the surface heat flux is
zero and fi is the surface heat transfer coefficient.
At the bottom of the lake, the boundary conditions are: (a) a no-slip
or zero-flow condition is assumed or it is assumed that:
- ~ pCd VB ^B (3.11)
o
where T_ is the bottom shear stress vector, p is the water density, C, is the
skin-fraction coefficient, VR is the magnitude of the bottom current while
Vg is the bottom current vector; and (b) the heat flux or temperature is
specified,
3T
^"a^" qB °r T = TB (3-12)
It is convenient to write the above system of equations and boundary
conditions in non-dimensional form by defining the following non-dimensional
quantities:
(u*, v*, w*) = (u, v, wL/H) / Ur
(x*, y*, z*) = (x, y, zL/H) / L
A «fc ^ W W
t* = tf , p* = P/PvUvfL
P* - P/Pr, T* = (T-Tr)/Tr,
V - V
-------
Suppressing the asterisk (*) for clarity, the non-dimensional governing
equations become:
+ * (3
(3-18)
2
where R_ = Rossby number = U /fL, £ =horizontal Ekman number = (A^) /fL ,
2 1/2
EV = vertical Edman number = (Ayr) /f H , Fr = Froude number = U /(gh) ' ,
Pe^j = horizontal Peclet number = U L/Cl^ ), and Pe = vertical Peclet number
It is convenient to transform the equations from an x, y, z coordinate
system to a vertically-stretched x, y, o coordinate system where
o = z/h(x,y) (3.19)
where h(x,y) is the local depth of the lake. With this transformation and
a uniform a grid spacing, the number of grid points in the shallow areas of
the lake is the same as in the deeper areas. Hence there is no loss of
accuracy in the computations due to lack of vertical resolution. Freeman
et al. (1972) used this method for investigating the circulation in Lake
Huron. The equations resulting from this coordinate stretching are similar
to the sigma equations derived by Phillips (1957).
The details of the derivation of the stretched equations are developed
in Appendix A. The final set of non-dimensional equations are:
Continuity
3(hu) + 3(hv) jco = (3f2Q)
3x 3y So-
x-momentum
i«- ^f 3(huv) 3Xhuv) , .3(&m)1 3c
^ 1^^~ I * *f - > .p*. f\, . T^^^^ j * ^r _ i - « g _|_
9t h 3 x ^ 3y 3o V 3x ax
-------
V 3
3v
v3o
y-momentum
3v = ^ [3(huv) 3(hv2) h3((ov)l
3t h [ 3x 9y aa J ~ U ~
__
3y
- a
IS
h
Hydrostatic
3p ,
^=-ph
Energy
-1JL(A J2V)
2 3a v 3cy
h
3x
?eHh
T* T a
. JD Li o
+ 9
Pe Hh/ 3a
v
Equation of state
P = P(T)
where a and a are defined as:
x y
fo
f° /
3p 3h
! 3x 3x
L >a \ J
r
'o
3p . 3h
ly" 0 ly
a
\
I .
£H Pda + ap
pda + ap
ya
and where
_ w a
w ~ h ~ h
3hN
(3.21)
(3.22)
(3.23)
(3.24)
(3.25)
(3.26)
(3.27)
(3.28)
The higher order terms (H.O.T.) in the horizontal diffusion terms con-
tain bottom slopes and/or their products and hence are generally small com-
pared to the listed leading terms when H/L « 1. It should be noted that in
deriving Eqs. (3.21) and (3.22), the vertically-integrated form of the hydro-
static equation was substituted into the x- and y- momentum equations
(Appendix A). Horizontal eddy coefficients were assumed to be independent
of space and time and hence Ag = 1 and KJJ = 1. In general, the vertical eddy
coefficient is a function of the wind, local depth, and vertical temperature
gradient. The following general form has been used by several authors
(Munk and Anderson, 1948; Sundaram and Rehm, 1970):
Ri)
m
(3.29)
-------
Kv = Kvo (1 + °2 Ri) (3.30)
where A and K are the eddy coefficients in the absence of any thermal
stratification, o^, a2, m, and n are empirically determined constants, and
Ri is the Richardson number which is the ratio of buoyancy and inertia
forces:
a
0.30
where u is the mean horizontal velocity. Several linear or quadratic forms
of the equation of state have been used by various authors (Paul and Lick,
1974; Simons, 1974).
The appropriate non-dimensional boundary conditions in the vertically-
stretched coordinate system are as follows. Since the free-surface elevation
is generally small compared to the water depth, it is generally valid to
apply the free-surface boundary conditions at the undisturbed lake surface
(a - 0).
Surface, a = 0:
~ ht .
3u _ x 3y
3a " A * 30
Bottom, o = -1:
u - v -' u = 0
1 3T _
h3^=qB °r T = TB
River inflow or outflow:
u = u(x,y,a)
v = v(x,y,a)
T « T(x,y,a)
Shore: u » v = 0
i£= n
3n °
Notice that the higher order terms which contain bottom slopes are neglected
in the bottom boundary condition.
During time periods when thermal stratification in the lake is insigni-
ficant and hence the density and temperature are approximately uniform,
10
-------
the energy equation and equation of state are decoupled from the momentum
equations. The terms a^ and ay in Eqs. (3.21) and (3.22) become zero. In
Lake Erie, this constant density condition is approximately valid from late
Fall to early Spring. For the shallow Western Basin of Lake Erie, this
assumption is approximately valid throughout the entire year.
11
-------
SECTION 4
GENERAL DESCRIPTION OF NUMERICAL MODEL
Basic Equations of the Two-Mode Model
As a first step in the development of a two-mode, free-surface model, a
non-stratified lake is considered. The general procedure of the two-mode
model remains basically unchanged when the effects of thermal stratification
are included.
In single-mode, free-surface models for a non-stratified lake, Eqs.
(3.20), (3.21), (3.22), (3.23) and associated boundary conditions are solved
directly by finite-difference techniques. Due to the presence of free-
surface gravity waves, the numerical time step is limited by the Courant-
Frederick-Levy condition which, in dimensional form, is
max
where Ax is the horizontal grid size, hmax is the maximum depth of the lake,
and 'rghmax is the maximum speed of propagation of surface gravity waves.
For Lake Erie, with a maximum depth of 61 m and Ax of 1km, the maximum allow-
able At is only 29 sec.
To alleviate this problem, a two-mode, free-surface model which treats
the external and internal modes separately can be used. The external mode
of the flow, i.e., the free-surface elevation and vertically-integrated
currents, are calculated from the following vertically integrated continuity
equation and equations of motion.
8f "B fhr N*"x h ' " Sv x""xv h 'I ' * l*3v ' liHl S 2 ' 3 2 I ' ^V^1 ""^ v^«3)
" 21 i 2 9
B (4>4)
-------
system and are given by
EH 2
B
x
B
and a , a and a are defined as
x y xy
9h1
3y
2-
?-i
<_9hT
3yl
3u
3a
3v
9a
a = -1
a = -1
a
a
and a = <^) (^
The internal mode is described by gradients in the vertical direction.
The governing equations for the internal flow variables thus consist of the
equations of motion and continuity equation in terms of differences between
adjacent grid points in the vertical direction:
3(hAu.)
k
3x
3(Au.)
K.
3t
+ Av,+
k
9y
R
D
h
_9x
H 1 9 °
h |_9x *"
T 11
r, . ,
[uA^XJ
( An \
Auv
IV
9x
3o
= 0
3y
h-
+ ^-i|A
, 2 9a v 9a
h
T^j]
(4.5)
(4.6)
9t
An -A-
AUk+ h
where u, ,
the k-tn
f 3
[ 9x
8(Av )\
Vi '
o3C j
i
+ ^
3(Au^
9y
+ EV 9
h2 8a
v
vr- and o)v are fluid velocities in the
iv tv
grid point in the vertical direction, k
3(Avk)
r dO
x, y,
ranges
1
and a
from
directions
1 near the
(4.7)
bottom to k = kffiax near the surface, and A^ = Fk+l ~ F is defined as the
difference in the variable F between the (k+ l)-th and the k-th grid points.
Higher order terms in the horizontal turbulence terms are listed in Appendix
A.
In order to circumvent the non-linearity of the equations, a time-
lagged approximation is introduced when integrating Eqs. (4.3), (4.4), (4.6)
and (4.7) with respect to time, i.e., the non-linear terms are evaluated by
means of values of u, v, and u> at the previous time interval.
Grid Structure
The basic grid structure is shown in Figure 4.1. The indices i, j, and
k correspond respectively to x, y, and o coordinates. In the staggered grid
in the x-y plane, the variables u and U are calculated at the mid-points of
the two boundaries parallel to the y-axis, v and V are calculated at the mid-
13
-------
i.j.k+l
Ht
ij.kf I
/
U
t
tl
Figure 4.1 Basic numerical grid structure.
-------
points of the two boundaries parallel to the x-axis, and w and C are calcu-
lated at the center of the grid. In the o direction, co is calculated at
full grid points and u and v are calculated at the mid-points. In the
numerical calculation for an actual lake, the shoreline is fitted with a
rectangular grid such that u equals zero or is determined by the river flow
along a shoreline parallel to the y-axis and v equals zero or is determined
by the river flow along a shoreline parallel to the x-axis.
The use of staggered grid avoids calculating all variables at all grid
points and also greatly simplifies the evaluation of gradient terms in the
finite difference equations. To evaluate variables at grid points where they
are not defined, simple averaging is performed:
4
(4.8)
where u(i,j) indicates the grid point at which u-j^ is evaluated. Details
of evaluating the variables at various types of grid points can be found in
Sheng (1975).
Finite Difference Equations
Various one and two-step methods were used to solve the external equa-
tions in a rectangular lake with constant depth or constant bottom slope.
Based on extensive numerical experiments with these schemes, a two-step
method similar to the Lax-Wendroff two-step method (Richtmeyer, 1967) was
chosen for the numerical integration of Eqs. (4.2), (4.3) and (4.4). This
scheme gives numerically-stable and accurate results and also does not have
the time-splitting problem generally associated with the usual three-time-
level schemes. It is convenient to write Eqs. (4.2), (4.3), and (4.4) as:
aiT
J£- E(U,V,C,u,...) + V (4.9)
sv
,C,v,...) -U
,...) (4.n)
where E and F represent the right-hand-sides of Eqs. (3.3) and (3.4) exclu-
ding the Coriolis term and G represents the right-hand-side of Eq. (3.2).
These equations are solved in the following two steps. First, evaluate the
variables at t + A te/2 from the results at t and the previous internal time
level, i.e.,
Ae
Second, evaluate the variables at t + At from the results at t, t + Ate/2,
15
-------
and the previous internal time level, i.e.,
= U* + Ate E* tJ*1/2.f*1/2,lMf2.u*....<) + At/1
= V11 + At F* (Un+1/V+1/2,?n+1/2,vm,...) - At1
e
e G*
where At indicates the external numerical time step (limited by Eq. (3.1),
n and n + 1 are indices indicating the previous and present external time
levels, m indicates the previous internal time level, and E*, F*, and G*
are the finite-difference forms of E, F«and G. Spatial derivatives are
central differenced. For example, for -g~ at the point u(i,j) where u is
evaluated (see Figure 4.1), we have
Ax
is evaluated by
A rigorous analysis of the numerical stability of the above system of
equations is extremely complicated, if not impossible. Instead, stability
criteria for individual terms can be established and used as guidelines in
choosing the proper Atfi. It can be shown that for a deep lake the stability
condition imposed by Eq. (A.I) is generally much more restrictive than those
imposed by the horizontal turbulence or non-linear inertia terms.
Internal Mode
The numerical stability criterion associated with the vertical turbu-
lent diffusion term can be studied by considering the one-dimensional heat
equation in dimensionless form:
3u _ Ev 32u ., .
a- -- T T (4.20)
3t h2 3o2
The stability criterion for the forward-time, central-space or FTCS method
for the finite-difference form of Eq. (4.20) is:
E
At< T^ - (4.21)
2ti (Aor
2
For conditions typical of Lake Erie, h = 10m, Ay = 25cm /sec, Aa = 1/6,
and At < 600 sec. However, in the near-shore region where h is only 2m,
16
-------
At < 24 sec. Hence the stability criterion (4.21) is much more restrictive
than those associated with the other terms in the equations of motion. To
remove the restriction (4.21), we tried three schemes: (1) Du-Fort Frankel,
(2) Hopskotch, and (3) Implicit (see Appendix B) . The implicit scheme was
found to be unconditionally stable and also did not give any non-physical
oscillations in the results as did the other two schemes. The implicit
scheme is also relatively easy to implement even for the three-dimensional
system and was therefore adopted for the present study. Applying the verti-
cally implicit scheme to Eqs. (4.6) and (4.7), the difference equations are
as follows :
(Auk)
m+1
(Auk)m + At^Coriolis) + (Horizontal Turbulence)
+ (Vertical Turbulence) m+1 + (Non-Linear Inertia)m]
(4.22)
m+1
(Avk)m = (Avk)m + Ati[(Coriolis)m + (Horizontal Turbulence)
+ (Vertical Turbulence)4"1 + (Non-Linear Inertia)]
(4.23)
where At^ indicates the internal numerical time step and m and m + 1 are
indices indicating the previous and present internal time levels (t and t
+ At±).
Although the horizontal eddy coefficient is much larger than the verti-
cal eddy coefficient, the horizontal grid size is generally several orders
of magnitude larger than the vertical grid size. Hence the stability limit
imposed by the horizontal turbulence terms is generally not as restrictive
as to warrant the use of an implicit scheme.
At every horizontal location i,j where u and v are calculated, there are
grid points in the a direction and hence lmax -1 equations for Au^ or
In order to solve for the horizontal velocities at all ^nax levels,
two additional equations are provided from the external mode equations with
the bottom stress terms treated implicitly:
(4.24)
(4.25)
m Dmfl-Ate/At1 +
V+1 = V+1-Ate/At
i +
,...) _ At
F* (V+1,...) + At
e " ' e
The vertical velocity u> is computed from
m+1
W.
m+1
rAa
Ax
m+1
,
m+1
m+1 , m+1
j)Vi,j,k " hv(i,j-D ^.J-l,
(4.26)
and the surface elevation £ can be calculated from
17
-------
rm+l-At /At .
W.27)
Numerical Procedure
To maintain numerical stability, the numerical time step for the exter-
nal mode calculation Ate is limited by Eq. (4.1), while the time step for
the internal mode At£ is several times larger.
The computational procedure is schematically shown in Figure 4.2. At a
typical time level, the three-dimensional velocities U,V,CD and the external
variables U, V, c are known (from the previous computations or the given ini-
tial condition) , and the bottom shear stress Tg can be evaluated from the
vertical profile of velocities. Subsequently, a relatively small time step
Ate is used to integrate the vertically-integrated equations by explicit
marching in time. Values of U, V, and £ are evaluated from Eqs. (4.12)
through (4.17) whereas the bottom stress in the momentum equations is not
updated until after N - 1 steps, where N = At^/Ate and is an integer. The
internal mode equations (4.22) and (4.23) are then integrated with a rela-
tively large time step At^ while simultaneously integrating the vertically
integrated momentum equation from (N - l)At_ to NAte. For consistency, both
the vertical diffusion terms in the internal equations and the bottom stress
terms in the external equations are treated implicitly. The vertical velo-
city 0), surface elevation 5, and bottom stress TB can then be calculated
from the appropriate equations. This cyclic procedure is repeated until
the computation is complete.
From formal analysis, the use of a vertically implicit scheme indicates
unconditional numerical stability and hence the internal numerical time
step is only limited by the convective terms. However, in a shallow lake,
the bottom stress is quite significant compared to the wind stress and
varies appreciably with time in transient situations. Hence the use of an
internal time step At^ larger than the time scale associated with the varia-
tion of bottom stress, i.e., H~. , would result in an oscillation of the
solution leading to numerical instability. This additional constraint on
the numerical time step is
where hj^ indicates the minimum depth of the basin.
(4.28)
Since Ate is proportional to Ax whereas At^ in a shallow lake is general-
ly limited by Eq. (4.28) and hence independent of Ax, the external mode cal-
culation generally takes proportionally less time in a coarser numerical
grid relative to the internal mode calculation.
18
-------
B
t EQUATIONS
FOR VERTICAL
SHEAR)
I
A te ( EXPLICIT)
Atj
( IMPLICIT
VERTICAL
TURBULENCE)
» «r
I A te
( EXPLICIT)
I
Ate ( IMPLICIT BOTTOM
FRICTION )
. v, w, TB, u, v,C
(VERTICALLY
> INTEGRATED
EQUATIONS)
Figure 4.2. Schematics of computational procedure.
-------
SECTION 5
APPLICATIONS TO IDEALIZED BASINS
A Two-Dimensional Basin
Extensive numerical experiments were carried out for a two-dimensional
lake to test the performance of the two-mode hydrodynamic model. Some
representative results are described here.
As shown in Figure 5.1, the two-dimensional lake represents a vertical
cross-section of a three-dimensional lake. The basin was assumed to be 2 m
deep at the shallow end with the depth increasing linearly to 30 m at the
deep end over a width of 100 km. A uniform wind stress of 1 dyne /cm was
assumed to be suddenly Imposed at the surface and the flow variables were
calculated for 48 hours. Representative results at two locations, one near
the shallow shore and one near the deep shore, are shown in Figures 5.2 and
5.3.
At point 1 near the shallow shore, as shown in Figure 5.2, the time
histories of the vertically integrated current and the horizontal current
near the bottom are very similar. Both showed significant oscillation, with
a distinct period of approximately 5 hours associated with the seiche oscil-
lation in the basin. A steady state is approached in about 24 hours. The
vertical profiles of horizontal velocities at several times are also shown
in Figure 5.2 The results indicate that, due to the shallow depth, the
effect of bottom friction is dominant during the transient stage. The use
of a At.^ greater than the vertical diffusion time of *L- = 4600 sec (based
on h = 4.8 m and Av = 25 cm/sec) resulted in oscillations in the solution.
The results computed with Ati = 3600 sec = 1 hr and Ate = 600 sec = 10 min
agree well with those computed with At^ = Ate = 10 min.
At point 2 near the deep end, as shown in Figure 5.3, the vertically
integrated current indicates a transient oscillation with a 5-hour period
lasting for about 24 hours. However, the bottom current slowly increased
with time with little oscillation. The effect of bottom friction is rela-
tively less important at point 2 than at point 1.
A Constant-Depth Square Basin
As another example, flow in a square basin 50 km by 50 km with a depth
of 10 m was studied. A uniform wind stress of 1 dyne/cm2 in the x-direction
was suddenly imposed. The vertical eddy viscosity was assumed to be 25 cm2/
sec, a typical value for Lake Erie under moderate wind conditions. The
results are computed with Ax = 5 km, Aa = 1/6, Ate = 5 min and Ati = 1 hr.
20
-------
N>
2M
I DYNE / CM
30 M
U
w
-*-
AX
o
AX = 4 KM
Ao-= 1/6
Figure 5.1 A Two-dimensional basin.
-------
to
to
hours
-11/12
-2«-
Figure 5.2.
(a) Integrated current and near-bottom current, and
(b) Vertical profile of horizontal currents at point 1 of Fig. 5.1
-------
co
24
32
40
h-
u
-------
The surface elevations at two locations in the basin, shown in Figure 5.4,
clearly indicate a period of oscillation of 2.5 hours corresponding to the
seiche oscillation in the basin. As can be seen, the surface elevations
compare well with those computed by a single-mode free-surface model with
At = 5 min.
24
-------
50 KM-
°l
C(CM)
4r
C,
TWO-MODE, Ate » 5min, A ti * I hr
SINGLE-MODE. At = 5 min
Figure 5.4 Results in a constant-depth square basin.
-------
SECTION 6
APPLICATIONS TO LAKE ERIE
Lake Erie (Figure 6.1) is composed of three basins, a shallow Western
Basin, a larger and relatively flat, but deeper Central Basin, and a cone-
shaped and deep Eastern Basin. The average depths of the Western, Central,
and Eastern Basin are 7, 18, and 24 meters respectively.
The three-dimensional free-surface model as described in previous sec-
tions has been applied to Lake Erie. The results are being used to compute
the transport and resuspension of sediments in the Western Basin of Lake
Erie. The output is being compared with data from coordinated aircraft/ship
surveys by the National Aeronautics and Space Administration. (Sheng and
Lick, 1977a; 1977b.)
For sufficient accuracy, a 1.6 km horizontal numerical grid is used in
the Western Basin and a portion of the Central Basin and is coupled with a
12.8 km grid in the rest of Lake Erie (Figure 6.2). Five grid points in the
vertical direction are used. The flows in the fine and coarse grids are
coupled dynamically and calculated simultaneously in time to allow full
interaction between the two grids (Sheng, 1975; Sheng and Lick, 1976). The
coupling scheme has been shown to produce smooth results across the fine and
coarse grids without appreciable computational reflections. In the present
analysis, three islands in the Western Basin are included as are river flows
from the Detroit River, the Maumee River, and the Sandusky River and outflow
to the Niagara River.
As a boundary condition, the wind stress at the lake surface needs to
be specified. For the relation between the wind stress and wind velocity,
it is assumed that
i" ' 'a°D "A (S-U
where t_w is the wind stress vector, pa is the air density, CD is a skin-
friction coefficient, Wa is the magnitude of the wind velocity while V^ is
the wind velocity vector. The coefficient CD is taken to be 0.00273 for
strong winds (Wa exceeds 6 m/sec) and 0.00166 for light winds, where W^ is
measured at 6 m above the lake surface (Wilson, 1960).
To obtain the wind stress over the entire lake under variable wind
conditions, an interpolation scheme has to be applied to the irregularly
spaced date recorded at weather stations around the lake. Interpolation
with weighted averages based on an inverse square of the distance from wea-
ther stations has been used (Platzman, 1963):
26
-------
I
0 10 tO 50 40
0
*0-
40-
fO-
0-
loo-l
ItO-l
uo-
too-
tto-
WCSTIMN §»sm
CINTML B»»IN
EASTERN BASIN
LAKE ERIE LONGITUDINAL
CROSS SECTION
Figure 6.1 Lake Erie bottom tooography.
-------
* DETROIT
+ WINDSOR
BAR PT.
KINGSVILLE
,PT. PELEE
r
/
/
+
LONDON
4 WIND STATIONS
A WATER LEVEL GAUGES
0 10 Ml.
+ SIMCOE
8-MILE GRID
4-CLEVELAND
CONTOUR
LABEL
A
B
C
D
E
F
S
K
P
CONTOUR -T^
VALUE XA,f
10 FTK, V
20 Tl
30 \\
40 C^A
50 \
60
BASS ISLANDS
KELLEY'S
ISLAND
PELEE ISLAND
Figure 6.2 (a) Numerical grid structure of Lake Erie and locations
of wind stations and water level gauges.
(b) Bottom topography of Western Lake Erie.
28
-------
rw(x,y,t) = I Xm (x,y) tw(t) (6.2)
m m
w
where ^ is the wind stress as derived from the wind velocity at station m,
where Xm is a dimensionless weighting function given by
An
(6.3)
and where rm - [(x-xm) + (y-ym) ]°"5 is the distance between the point (x,y)
and the station (x^y^ . To account for the difference between winds over
the water and over the land, the stress actually used in the calculations was
the stress from Eq. (6.2) multiplied by 1.5. In the interpolation, data from
nine weather stations (Figure 6.2) were used.
In the numerical computations, the external mode of the free-surface
model is limited by a time step Ate £ x/(2gh)°*5 = 81 sec, while the internal
mode is limited by At± < h|in/2Av = 1500 sec (based on h,^ = 2.7m = 9 ft).
The external time step was taken to be 75 sec and the internal step to be
900 sec. A vertical eddy coefficient Av of 25 on2/sec was assumed. When the
horizontal turbulence terms were included in the equations, a horizontal eddy
coefficient AJJ of 10°cm /sec was assumed.
Constant Wind
As a first example, the time-dependent currents caused by a suddenly
imposed constant wind of 5.2m/sec from W67S were calculated. It was assumed
that the lake was calm initially. Calculations were carried out for 24 hours.
Although a suddenly imposed, constant wind is an idealized situation, it
serves as a simple example. Some weather systems may be approximated by such
a wind.
The computed results indicate that initially mass is transported in
the direction of the wind from the western end of the lake to the eastern
end. About seven hours after the start of the wind, the mass flow in the
interior of the lake begins to reverse its direction while the mass flow in
the shallow coastal regions is still in the direction of the wind. The sur-
face elevation versus time at three locations in the western portion of the
lake is shown in Figure 6.3. Periods of seiche oscillation in the longi-
tudinal direction (14 hours) and in the transverse direction (2 to 3 hours)
are apparent. The surface elevation approaches a steady state in about 24
hours. In the Central and Eastern Basins, however, the steady state is
approached in about two to four days. The horizontal velocities at three
hours after the initiation of the flow are shown in Figure 6.4 for every 3.2
km interval in the x and y directions. Figure 6.4 shows the horizontal
velocities near the surface (o = 0.1). It can be seen that the near-surface
currents are from the Western into the Central Basin and are deflected some-
what to the right of the wind. A similar flow pattern is also apparent at
mid-depth (a = 0.5) as shown in Figure 6.4b. It can be shown that the com-
puted results continue to show currents from the Western to the Central
Basins at all depths until about 7 hours at which time strong opposing cur-
29
-------
-10
CM
8
12
16
20
24 HOURS
TIME
-20
Figure 6.3
Surface elevation computed by free-surface model as a function of time
at three locations in Lake Erie. The wind is constant with a velocity
of 5.2 m/sec from W 67 S.
-------
WIND
10
20 Ml.
DISTANCE
I FT / SEC
0 20 KM
0 I FT/SEC 0
CURRENT [£? CURRENT
020 CM/SEC 0 20 CM/SEC
(a)
(b)
Figure 6.4
Horizontal velocities at a constant (a) 0.1 depth and (b) 0.5 depth from
surface after 3 hours. The wind is constant with a velocity of 5.2 m/sec
from W 67 S
-------
rents develop. The currents continue to oscillate with time until about 24
hours when a steady state is approximately reached. The steady-state cur-
rents as shown in Figure 6.5 compare well with the results of a steady-state
model (Gedney and Lick, 1971).
Variable Wind
The time-dependent flow in Lake Erie caused by actual wind conditions
during a two-day period from March 7 through March 8, 1976 was also calcu-
lated. For this period, Figure 6.6 shows the wind speed as a function of
time at several weather stations around the lake. It can be seen that the
winds were generally quite strong on March 7, on the order of 20 knots, and
gradually decayed on March 8 to below 10 knots. To generate the required
stress field for the numerical computation, hourly wind data over the regular
numerical grid were obtained by interpolation of the irregularly spaced wind
data at 9 weather stations (Figure 6.2). The lake was assumed to be calm
initially at 00:00 on March 7. Calculations were made for 48 hours until
00:00 March 9.
As a verification of the model, the computed surface elevations were
compared with water level data recorded by several gauges on Lake Erie
operated by the U.S. Ocean Survey. The results for the surface elevations
at Toledo, Bar Point, Point Pelee, Cleveland, and Sturgeon Point are shown
in Figure 6.7. The observed curves were obtained by subtracting the 24-hour
running mean of the sum of all stations from the recorded lake level data.
The longitudinal seiche period of 14 hours, on which is superimposed the
transverse seiche period of 2 to 3 hours, is again apparent. Reasonable
agreement between computed and observed data exists at all stations except
during the initial 12 hours.
The largest discrepancy between computed and observed data occurs at
Bar Point. This is probably due to the incorrect specification in the
analysis of the Detroit River which would clearly affect the water level in
the vicinity of the river. Although the flow rate at the mouth of the De-
troit River is generally time-varying and should depend on wind conditions,
detailed measurements of the flow rate during the time period of this study
were not available and the flow was assumed to be constant.
Figure 6.8 shows the horizontal velocities versus time near the surface
(o = 0.1) and at mid-depth (a = 0.5) at two locations in the Western Basin
(indicated at points 1 and 2 in Figure 6.2). The longitudinal and transverse
seiches are clearly identifiable at both locations. The transverse seiche is
more pronounced at point 1.
By comparing the time histories of the velocities at points 1 and 2
with the time histories of the wind stresses and surface elevations at near-
by stations, it is apparent that the currents at point 1 correlate more
closely with the wind stresses whereas the currents at point 2 correlate
more closely with the surface elevations. This can be understood as follows.
At point 1 near Toledo where the depth is shallow and the winds are strong,
the currents at all depths are primarily driven by the wind, and the influ-
ence of the pressure gradient caused by tilting of the free surface is
32
-------
WIND
DISTANCE
0 I FT/SEC
10 20 MILES
0 10 20 KM
0 I FT/SEC
CURRENT
0 20 CM/SEC ° 20 CM /SEC
'
.,,
Figure 6 5. Horizontal velocities at a constant (a) 0.1 depth and (b) 0.5
depth from surface after 24 hours. The wind is constant with
a velocity of 5.2 m/sec from W 67 S.
-------
30-
KNOTS
D -
C -
T -
P -
DETROIT
CLEVELAND
TOLEDO
POINT PELEE
(a)
48 HOURS
-4-
Figure 6.6 (a) Wind speed as a function of time at four weather stations
around Western Basin of Lake Erie.
(b) Wind Stresses at Toledo.
(c) Wind stresses at Point Pelee.
-------
STURGEON POINT
48 HOURS
FREE SURFACE MODEL
WATER LEVEL GAUGE
Figure 6.7 Comparison of surface elevation computed by free
surface model with water level measurement at 5
gauges.
35
-------
CM /SEC
12
24
36
48
HOURS
FREE SURFACE MODEL
Figure 6.8a. Horizontal velocities at 0.1 depth and 0.5 depth as
function of time at station 1 (near Toledo).
36
-------
t
CM/SEC
12
1
24
i
36
1
48 HOURS
4O-
FREE SURFACE MODEL
Figure 6.8b. Horizontal velocities at 0.1 depth and 0.5 depth as
function of time at station 2 (near Point Pelee).
37
-------
relatively small. However, at point 2 near Ft. Pelee where the depth is
greater and the winds are weaker during the two-day period, the currents
at all depths are more influenced by the pressure gradient which in turn
is significantly affected by surface oscillations due to gravity waves.
The horizontal distribution of currents at several times are shown in
Figures 6.9 to 6.12. During the first 24 hours, due to the strong winds and
the start-up of the flow, the instantaneous currents are generally quite
strong. The currents (at a = -0.1 and -0.5) at 0600 3/7/76 and 1200 3/7/76
are shown in Figures 6.9 and 6.10 respectively. It is interesting to note
that the currents between the islands and in the Pelee Passage are generally
quite strong compared to the currents in other regions. The near-surface
currents and the mid-depth currents are primarily in the same direction.
At 0000 3/8/76, although the winds are predominantly from the North-
west, the results show strong opposing currents at all depths (Figure 6.11)
due to the mass flux from the Central Basin into the Western Basin associated
with the free-surface movement. Maximum currents near the surface and at
mid-depth are on the order of 100 cm/sec and 60 cm/sec, respectively.
At 1200 3/8/76, the winds have decreased further and are relatively
mild. As shown in Figure 6.12, the currents are still quite significant. It
is of particular interest to note the presence of a large eddy off Point
Pelee. From the calculation, it can be shown that this large eddy existed
throughout the water column for approximately 1.5 hours. The presence of
such an eddy off Point Pelee of similar scale has been observed by the
LANDSAT satellites quite frequently. Figure 6.13 shows the eddy observed on
April 14, 1973. Unfortunately, LANDSAT satellites were not over Lake Erie
on 3/8/76. The horizontal currents computed by a rigid-lid model indicate
very weak currents at 1200 3/8/76 (Figure 6.14). The fact that the large
eddy is not present in the rigid-lid results indicates that seiche oscilla-
tion is an important cause of the observed eddy.
At 0000 3/9/76, the winds have decreased further and the flow in the
northern portion of the Western Basin showed opposing currents at mid-depth
and near the surface, indicating a quasi steady-state (Figure 6.15).
Comparison With A Rigid-Lid Model
The above calculations of the time-dependent flow in Lake Erie were
repeated using the rigid-lid model of Paul and Lick (1974, 1979). The
results of the present two-mode free-surface model and the rigid-lid model
were compared by Sheng et al. (1978). For relatively short time intervals,
significant differences between the results of the two models occur as
expected since the seiche motion is eliminated in the rigid-lid model. The
horizontal velocities at point 1 near Toledo and point 2 near Pt. Pelee as
computed by the two models are shown in Figures 6.16 and 6.17. Based on
comparisons made at four locations in the Western Basin, the two results
differ most at point 2, the pressure-driven station, and differ least at
point 1, the wind-driven station.
The rigid-lid results tend to represent a filtered version (i.e., a
38
-------
10
2O MILES
DISTANCE E
0 K) 20 KM
0 I FT / SEC
CURRENT ff
MAGNITUDE 0 2Q CH/SEC
CURRENT
MAGNITUDE
FT/SEC
020CM/SEC
-/'--.
. *s ' .
- r " * * «
U"b ' ' * «
* v » » i '
* *. ». *"
to****»_^*" ** *" **
* * w v *
-------
10 20 MILES
DISTANCE E
0 10 20 KM
CURRENT °JFT/SEC
MAGNITUDE 0 20 CH/SEC
CURRENT
MAGNITUDE
0 I FT/SEC
0 20 CM/SEC
(b
Figure 6.10 Horizontal velocities at a constant (a) 0.1 depth and
(b) 0.5 depth from the surface at 1200 3/7/76.
-------
0 5 FT/SEC
&
0 I M/SEC
20 MILES
20 KM
CURRENT
MAGNITUDE
0 I FT/SEC
H3
0 20 CH / SEC
Figure 6.11 Horizontal velocities at a constant (a) 0.1 depth and
(b) 0.5 depth from the surface at 0000 3/8/76.
-------
CURRENT
0 I FT/SEC
EP
0 20 CM/SEC
20 MILES
20 KM
I FT/SEC
20 CM / SEC
Vv
i
11
'
l_* *
^4
1
(a.
i
i
i
"U ' .
1
1
i
i
_J
fl
1
1
1
1
;
i
i
%
V
If"
"
T* ,
Kj
n
i h^
»^~K
iv -l_
ft ft. w fc 1
* v v
III 1
1 v »
»». fc *
v \ *
\ V * » »
s \ * \ »
iiii ]G
' ' 1 ' ". '.
blXf:;
i
i
i
ii
*!,
n
^ ' ' '.
kj
-
P
--1
i
i
'
i ,
.
t
«.
i t
i i
i i
*
%
» *
.
\
\
t
i i
i i
i i
1 V
» I
v
* -*
"<>
' "<
;'////
T ' / / /;
s / t 1
f f I \
t I
I *
I '
»v
»v %
\ V
Figure 6.12
Horizontal velocities at a constant (a) 0.1 depth and
(b) 0.5 depth from the surface at 1200 3/8/76.
-------
SUSPENDED
SEDIMENT
LOAD
ARCHIMEDES
SPIRAL FLOW LINES
. . « w *
& .'
Figure 6.13.
Large-scale eddy observed by ERTS-1 (LANDSAT-1) off the
tip of Point Pelee on 3/14/73.
43
-------
CURRENT
0 I FT / SEC
E?
0 20 CM /SEC
20 MILES
I FT/SEC
20 KM
20 CM/SEC
i
Figure 6.14
Horizontal velocities at a constant (a) 0.1 depth and
(b) 0.5 depth from the surface at 1200 3/8/76.
-------
0 I FT/SEC Q 10 20 MILES CURRENT °-
CURRENT gp DISTANCE I . ' . MAGNITUDE
MAGNITUDE Q^O CM/SEC 0 10 20 KM
I FT/SEC
20 CM / SEC
Figure 6.15 Horizontal velocities at a constant (a) 0.1 depth and
(b) 0.5 depth from the surface at 0000 3/9/76.
-------
CM / SEC
24
36
48 HOURS
"o.,I0
Figure 6.16.
FREE SURFACE MODEL
RIGID-LID MODEL
Horizontal velocities at 0.1 depth and 0.5 depth as
function of time at station 1 (near Toledo).
46
-------
CM/SEC
48 HOURS
-20
Figure 6.17.
FREE SURFACE MODEL
RIGID-LID MODEL
Horizontal velocities at 0.1 depth and 0.5 depth as
function of time at station 2 (near Point Pelee).
47
-------
version in which the seiche oscillations are eliminated) of the free-surface
results. Hence it is of interest to compare the two models in terms of
time-averaged currents over periods of one or more dominant seiche cycles.
The instantaneous currents at points 1 and 2 have been averaged over three
14-hour periods (0-14, 14-28, 28-42) and also over the entire 42-hour period
and are shown in Figures 6.18 and 6.19. As can be seen in Figure 6.18, the
time-averaged currents at point 1 as computed by the two models agree well
except during the third seiche cycle. However, as shown in Figure 6.19, the
time-averaged free-surface and rigid-lid currents at point 2 are significant-
ly different. It is interesting to note that the free-surface and rigid-lid
currents are still quite different when averaged over three complete seiche
cycles (42 hours).
The instantaneous currents computed by the rigid-lid model are general-
ly weaker than the free-surface currents (see Figures 6.12 to 6.14). The
rigid-lid currents indicate a much faster response to the wind and generally
agree very well with the results of a steady^state model (Gedney and Lick,
1971).
Since the winds were generally strong during the first day and lighter
during the second day, it is of interest to compare the daily averaged free-^
surface and rigid-lid currents during the two-day period. Daily averaged
currents have also been used in studies of long-term, time-dependent disper-
sion of pollutants in Lake Erie (Lam and Jaquet, 1976). Thus it is important
to investigate the differences in long-term averaged currents between the
free-surface and the rigid-lid models, particularly in the study of the dis-
persion of pollutants in a lake. As shown in Figures 6.20 and 6.21, the
daily-averaged currents computed by the two models agreed quite well on
March 7 but differed considerably on March 8. This can be understood as
follows. Although the effects of seiche oscillations were important through-
out much of the two-day period, the winds were also quite strong on the
first day. Hence, the direct wind-driven effect on the currents was domi-
nant over the effects of seiche oscillations. On March 8, the winds de-
creased significantly and the effects of sieche oscillations became relative^-
ly more important. The rigid-lid currents thus differed considerably from
the free-surface currents even when averaged over 24 hours. Although not
presented here, it can be shown that the differences between the daily-'
averaged free-surface and rigid-lid currents was more pronounced at depths
than near the free surface.
Due to the observed differences in long-term time-
-------
-O.I
R\\F
vo
0 5
0 5
5 cm/sec
-------
-------
10 20 MILES
FREE-SURFACE
DISTANCE E
0 10 20 KM
CURRENT gJ FT/SEC
MAGNITUDE 0 20 CM / SEC
ill V V \\V
VI V
V
mi
RIGID-LID
Figure 6.20 Daily-averaged horizontal velocities at a constant 0.1 depth from surface face on 3/7/76.
-------
0
DISTANCE E
10 20 MILES
0 10 20 KM
CURRENT
MAGNITUDE
FT/SEC
CM / SEC
>:.
' <
FREE-SURFACE
RIGID-LID
Figure 6.21 Daily-averaged horizontal velocities at a constant 0.1 depth from
surface on 3/8/76.
-------
currents as predicted by the rigid-lid model are generally smaller than
those predicted by the free-surface model, the use of a rigid-lid model
would generally result in smaller bottom shear stresses and less sediment
resuspended from the bottom.
Effect of Non-Linear Inertia and Horizontal Turbulence
In Lake Erie under non-stratified conditions, Rg « 1 and Eg « 1.
Hence non-linear inertia and horizontal turbulence can generally be neglec-
ted compared to Coriolis and vertical turbulence, respectively. In a
numerical study, Sheng (1975) found that the thicknesses of the coastal
boundary layers in a non-stratified Lake Erie were on the order of 1 Km or
less. Hence the non-linear inertia and horizontal turbulence terms, which
are important within the coastal boundary layers, can be neglected if a
horizontal numerical grid larger than 1 Km is used. In the present study
with a 1.6-Km grid, the time-dependent currents in the Western Basin were
computed with and without the non-linear inertia and horizontal turbulence
terms in the governing equations and the results compared for up to 12 hours.
For both the constant wind and the variable wind cases, the horizontal cur-
rents computed with and without these terms were almost identical throughout
the entire basin. Hence the non-linear inertia and horizontal turbulence
terms were neglected during subsequent calculations, resulting in significant
saving of computer time. However, in a stratified lake, thermal convection
generally has a strong influence on the temperature field, which in turn
significantly influences the flow field, and hence the nonlinear inertia (or
convective) terms must be retained.
53
-------
REFERENCES
Bennett, J.R., 1977; A Three-Dimensional Model of Lake Ontario's Summer
Circulation, I. Comparison with Observations, J. Phys. Oceanogr., ]_,
591-601.
Berdahl, P., 1971: Oceanic Rossby Waves; A Numerical Rigid-Lid Model.
. Rept. ITD - 4500, UC-34, Lasrence Radiation Lab., U. of Calif.,
Livermore.
Cheng, R.T., T.M. Powell and T.M, Dillon, 1976; Numerical Models of Wind-
Driven Circulation in Lakes, Appl. Math. Modelling, !_, 141-159.
Crowley, W.P., 1969; A Numerical Model for Viscous Free Surface Barotropic
Wind-Driven Ocean Circulation, J. Comp. Phys., 5^.
Crowley, W.P., 1970; A Numerical Model for Viscous Non~Divergent, Barotropic,
Wind-Driven, Ocean Circulation, J. Comp. Phys., 6^.
Donaldson, D.duP., 1973: Atmospheric Turbulence and the Dispersal of Atmo-
spheric Pollutants, AMS Workshop of Micrometeorology, Boston, 313-390,
Also EPA-R4-73-016a.
Freeman, M.J., A.M. Hale and M.B. Danard, 1972: A Modified Sigma Equations
Approach to the Numerical Modelling of the Great Lakes Hydrodynamics,
J. Geo. Res., 77. 1050-1060.
Gedney, R.T. and W. Lick, 1971: Wind-Driven Currents in Lake Eric, J. Geo.
Res., 77, No. 15.
Haq, A., W. Lick, and Y.P. Sheng, 1974: The Time-Dependent Flow in Large
Lakes with Application to Lake Erie, Tech. Rept., Earth Sci. Dept.,
Case Western Reserve Univ., Cleveland, Ohio.
Harlow, F.H., and P.I. tfakayama, 1967: Turbulence Transport Equations, The
Physics of Fluids, Vol. 10, No. 11.
Kizlauskaus, A.G. and P.L. Katz, 1973: A Numerical Model for Summer Flows in
Lake Michigan, Arch, Meteor. Geophys. BiokJLim..
KoLnogorov, A.N., 1942: Equations of Turbulent Motion of an Incompressible
Turbulent Fluid, Izv. Akad. Nauk. SSSR Ser. Phys. VI, No. 1-2.
Lam, D.C.L. and J.-M. Jaquet, 1976: Computations of Physical Transport and
Regeneration and Phosphorus in Lake Erie, Fall 1970, J. Fish. Res. Bd.
Can.. 33, 550-563. '
54
-------
Launder, B.E. and D.B. Spalding, 1972: Mathematical Models of Turbulence.
Academic Press, New York, N.Y.
Lewellen, W.S., M. Teske, and Y.P. Sheng, 1978: Progress Report on Modeling
Tornado Dynamics (An Assessment of the Sensitivity of Model Results to
the Turbulence Model), A.R.A,P. Tech. Memo. 78-5.
Lick, W., 1976a: Numerical Models of Lake Currents, U.S. Environmental Pro-
tection Agency, Washington, B.C.
Lick, W., 1976b: Numerical Modeling of Lake Currents, Annual Review of Earth
and Planetary Sciences, Vol. 4.
Liggett, J.A., 1969: Unsteady Circulation in Shallow, Homogeneous Lakes, JN.
Hydr. DIV., Proc. ASCE, J35 (HY4), 1273-1288.
Munk, W.H. and E.R. Anderson, 1948: Notes on A Theory of the Thermodine, J.
of Mar. Res., 7.
Paul, J. and W. Lick, 1974; A Numerical Model of Thermal Plumes and River
Discharges. Proc. 17th Conf. GLR, IAGLR.
Paul, J. and W. Lick, 1979; A Rigid-Lid Lake Circulation Model. U.S. EPA,
Washington, D.C.
Phillips, N.A., 1957: A Co-ordinate System Having Some Special Advantages
for Numerical Forecasting. J. Meteorol., 14.
Platzman, G.W., 1963: The Dynamic Prediction of Wind Tides on Lake Erie,
Meteorology Monograph, 4.
Richtmeyer, R.D. and K.W, Morton, 1967: 'Difference Methods for Initial Value
Problems', John Wiley, N.Y.
Roache, P.T., 1972: 'Computational Fluid Dynamics', Hermosa Pub., Albur-
qureque, N.M.
Sheng, Y.P., 1975: Wind-Driven Currents and Dispersion of Contaminants in
the Near-Shore of Large Lakes. Ph.D. thesis, School of Engrg., Case
Western Reserve Univ., Also as Rept. H-75-1, U.S. Army Engrg. Water-
ways Exp. Stn., Vicksburg, Miss. (NTIS AD-A017694).
Sheng, Y.P. and W. Lick, 1976; Currents and Contaminant Dispersion in the
Near-Shore Regions and Modification by Jetport. J. Great Lakes Res.
2^ 402-414.
Sheng, Y.P. and W. Lick, 1977a: Numerical Modelling of the Time-Dependent
Flow and Dispersion of Sediments in a Shallow Lake, presented at the
1977 AGU Midwest Meeting, Purdue Univ., Sept. 26-28.
55
-------
Sheng, Y.P. and W. Lick, 1977b: A Three-Dimensional Numerical Model for
the Time-Dependent Flow and Dispersion of Sediments in the Western
Basin of Lake Erie, presented at the 1977 AGU Fall Meeting, San
Francisco, Dec. 5-9, Transaction AGU, 58, 12.
Sheng, Y.P., W. Lick, R. Gedney and F. Molls, 1978: A Comparison of a
Free-Surface and a Rigid-Lid Model, J. Phys. Oceanogr., 8^.
Simons, T.J., 1972: Development of Numerical Models for Lake Ontario.
Proc. 15th Conf. GLR, IAGLR, 655-672.
Simons, T.J., 1974: Verification of Numerical Models of Lake Ontario
Part I: Circulation in Spring, Early Summer, J. Phys. Oceanogr.,
5, 98-110.
Sundaram, T.R. and R.G. Rehm, 1970: Formation and Maintenance of Thermo-
clinics in Stratified Lakes Including the Effects of Power Plant
Thermal Discharges, AIAA Paper No. 70-238.
Wilson, B.W., 1960: Note on Surface Wind Stress Over Water at Low and
High Speeds, J. Geo. Res., 65, 3377-3382.
56
-------
APPENDIX A
DERIVATION OF THE VERTICALLY STRETCHED EQUATIONS
In order to incorporate variable bottom topography into a general lake
circulation model, vertically stretched coordinates can be used. This pro-
cedure maps the variable depth basin to a constant depth one (with modified
equations) and results in an equal number of evenly spaced vertical grid
points in both the shallow and deeper parts of the basin.
The coordinate system is transformed from x,y,z, to a,6,cr, where a=x,
3 = y, and o= z/h(x,y). Hence it can be shown that (A-l)
(A-2)
3a
3x '
32a
3x2
3a| o
3y| h
32o
* 2
3h
3ct '
a 32h
h , 2
3x
3h
33
1 2o
'h2
/3h
I 3x
£
h
2 a 3h
(A-3)
_
3x
_3 o 3h 3
3d h 3a 3a
(A-4)
3 3 o_ 3h 3
liy ~ 33 h 33 3o
(A-5)
JL_ = I _JL
3z h 3a
(A-6)
and
3x
3a
2 ^2 ..
3a 3 4. 3 o _o
"3x" 3a3a "^T 3a
22
|3a| 3
13x I 2
(A-7)
3y
3a
3y
3y
2 30
37
2_i!
302
(A-8)
57
-------
2 2
T ~ ~? 2" (A-9)
3zZ ti 3a
Define
_ _L d* - L. ** (
U ~ R_ dt ~ 1L. dt
D D
1 dy 1 dB /.
v = *- = -- i A
V IL, dt IL dt ^
D D
1 dq
where u is related to the vertical velocity w by:
The expression
M + R I 3uA + 9vA . 3wA
3t ^ 3x 3y 3z
(A-14)
where A is the appropriate dependent variable (1 for the continuity equation
u, v, w for the three momentum equations, and T for the energy equation), can
be transformed into the following form:
3A , _ I 3 / »\ o 3h 3 / »\ i 3 / A \ 0" 3h 3 / , \ .
-^ + RT, Kr («A) - T- -r- £- (uA) + vs- (vA) - T- -TT- ^- (vA) +
k) (A-15)
11 OU 1
Substituting (A-13) into (A-15), we get
The horizontal turbulence terms
where A is u, v, or T, can be transformed into the following form:
58
-------
3ct
3d
Jh _3A
3d 3a
j}_
33 x"33
3h _3A
33 3o
_3h ^_
3d 3a
3A,
lh
33
3h
3B
3A
3a
(A-18)
The first two terms in (A-18) are generally much larger than the other
terms which contain bottom slopes or their products and hence are of higher
order in a basin with H/L«1. However, these higher order terms may be com-
parable to the leading terms in regions of sharp bottom slopes.
The vertical turbulence term
(A-19)
can be transformed into
Y. _ (A
,2 3a
(A-20)
The horizontal turbulence terms in the vertically integrated equations
in an x,y,z system are
-h
where A is u or v.
ay
dz
(A-21)
The corresponding expression in the transformed a,$,a system can be ob-
tained by vertically integrating Eq. (A-18) from a = -1 to a = 0 and is:
3d
33
.
h I
«
33
(A-22)
a=-l
where A is defined as
= h Ada
-1
The transformed equations are as shown in Eqs. (2.20) through (2.24).
(A-23)
59
-------
APPENDIX B
AN IMPLICIT SCHEME
The heat equation is here considered in non-dimensional form:
ia-i ifs. CR n
at "2 ,2 (B'1)
h 3a
The forward-time central-space (FTCS) scheme solves Eq, (B-l) by means of the
following finite-difference equation:
- " (B-2)
2
where a = Ev/(hAo) , n and n+1 are indices representing previous and present
time levels, K indicates the vertical level, and t is the non-dimensional
time step. Numerical stability requires the following criterion to be satis-
fied:
. At < Io (B-3)
This condition is very restrictive for a shallow lake such as Lake Erie.
To alleviate this problem, several schemes can be used. The Du-Fort Frankel
scheme uses the following form:
n+1 n-1 , , n n+1 n-1 . n v A. , ,%
"K = UK + a(uK+i - "K ~ "K + Vi} At (B~4)
Hence
" + u) At (B-5)
It can be shown that the Du->Fort Frankel scheme is unconditionally stable
but gives oscillatory results due to time-splitting (Roache, 1972). The
Hopskotch scheme is also unconditionally stable and uses the following forms:
Odd n, even K and even n, odd K;
n+1 n . . / n n . n N /TI/:\
"K " "K + aAt (UK+I ^ 2uK + "K^ (B^6)
60
-------
Odd n, odd K and even n, even K:
n+1 n , n+1 n+1 x ,_ ....
UR = UR + aAt (u^ - 2uK + u^) (B-7)
Both the Du-Fort Frankel and the Hopskotch schemes were used to solve Eq.
(B-l) and the results contained non-physical oscillations during the trans^
ient stage, although both schemes converged to the same steady state.
Sheng (1975) modified the Du-Fort Frankel scheme and applied the fol-
lowing form to compute the dispersion of contaminants in the near-shore
regions of Lake Erie:
n+1 n , n n+1 n-1 n , , ,
U = U + aAt (u - U - U + ) (B-8)
Relatively large numerical time steps up to 1 hour were used in the calcula-
tion with a horizontal spatial grid of 0.4 Km.
In order to have a numerically stable scheme and also alleviate the
physical oscillations, a completely implicit scheme can be used:
n+1 n , n+1 0 n+1 n+1.
U = U + aAt
The following boundary conditions were used:
u = 0 at a = -1 (B-10)
w
11 - -- at a - 0 (B-ll)
61
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1. REPORT NO.
EPA-600/3-80-047
3. RECIPIENT'S ACCESSION NO.
4. TITLE AND SUBTITLE
A Two-Mode Free-Surface Numerical Model for the Three-
Dimensional Time-Dependent Currents in Large Lakes
5. REPORT DATE
May 1980 issuing date
6. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
Y. Peter Sheng and Wilbert Lick
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Department of Earth Sciences
Case Western Reserve University
2040 Adelbert Road
Cleveland, Ohio 44106
10. PROGRAM ELEMENT NO.
1BA769
11. CONTRACT/GRANT NO.
R803704
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Duluth, Minnesota 55804
13. TYPE OF REPORT AND PERIOD COVERED
Final
14. SPONSORING AGENCY CODE
EPA/600/03
15. SUPPLEMENTARY NOTES
16. ABSTRACT
A two-mode, free-surface model based on vertically-stretched coordinates and
a vertically-implicit scheme has been developed and applied to Lake Erie under
non-stratified conditions. A brief description of the general equations and
boundary conditions is first given. The detailed equations for a two-mode, free-
surface model are then described. Finite-difference procedures including the
finite - difference equations are listed in detail. The model is first applied to
idealized basins and then to Lake Erie for two wind conditions; (1) a uniform
wind suddenly imposed, and (2) an actual wind variable in both space and time.
The results of the computations are compared with the results from a single-mode,
free-surface model and from a rigid-lid model.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lOENTIFIERS/OPEN ENDED TERMS C. COSATI Field/Group
Hydrodynamics
Circulation
Mathematical Models
Lakes
Lake Current Models
Lake Erie
08H
18. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (This Report)
UNCLASSIFIED
21. NO. OF PAGES
70
20. SECURITY CLASS (Thispage)
UNCLASSIFIED
22. PRICE
EPA Form 2220-1 (Rex. 4-77) PREVIOUS EDITION is OBSOLETE
62
mmms OFFICE. isw-657-M4/5668
------- |