-------
acting on the surface of a variable depth basin of the form:
||= 2.0(1 - 2x)(y - y2)
- 2.0(1 - 2y)(x - x2)
0 <_ X
0 £ y
The wind stress creates a maximum wind stress curl (9r /3y) of -0.5
X
at x = 0.5 and y = 0.0 and 1.0.
The thermocline depth, £, contours for these conditions are shown in
1 /?
Fig. 24 for 6 = 1/24 where 6 = (A,/A ) ' and A^ is the vertical eddy
diffusivity in the hypolimnion. The £ solution contains boundary layers
along y = 0, y • 1.0 and x = 0 boundaries. These boundary layer regions
are more clearly shown in Fig. 25 which is a plot of the thermocline
depth along the x = 48.3 km (x = 0.5 in non-dimensional form) station.
Figure 25 shows the thermocline to have a large shape variation in the
cross wind direction.
The horizontal velocities at the surface and at 20 and 35 meter depths
are shown in Figs. 26, 27, and 28. The figures show the effect of the
shape of the thermocline on the velocities. The thermocline shape
produces the clockwise gyre shown in Fig. 27. Below the thermocline
at a depth of 35 m, geostrophic velocities of order 1.0-5.0 cm/sec
occur. The high geostrophic velocities near the boundaries are created
by the r\ boundary layers, where n = 5 + £.
The unit order wind stress gradient solution is very dependent on the
value of the hypolimnion eddy diffusivity. Figure 25 shows the shape
of the thermocline at x = 48.3 km for the same applied wind stress for
the cases where 6 = 1/24, 1/12 and °°. The 6 = °° case corresponds to the
solution where the velocities in the entire hypolimnion are zero,
which is the quasi-compensation assumption referred to above.
Case 2; Uniform Wind with Variable Depth
To demonstrate the importance of wind stress curl, the thermocline
63
-------
DISTANCE
CURRENT
MAGNITUDE
QE=
15 miles
20km
,lft/sec
20 cm/sec
-27.0 -23.0 -19,0
-. -. -, lRn
31.0 \-25.0 \-2I.O ! -.] 5'° Tll-0
\\i °
-7.0
Figure 24. Thermocline depth contours with wind
stress curl of order one (6 = 1/24).
64
-------
-lOi-
= 00
(ZERO HYPOLIMNION
VELOCITY APPROXIMATION)
20 40 60 80
DISTANCE, km
100
Figure 25. Thermocline profile at x = 48.3 km with wind
stress curl of order one.
DISTANCE
CURRENT
MAGNITUDE
OE
01
, (Smiles
20km
.Ift/sec
20 cm/sec
Figure 26. Horizontal surface velocities (6 = 1/24)
65
-------
DISTANCE
CURRENT
MAGNITUDE
.(Smiles
20km
.1 ft/sec
20 cm/sec
' -x S
' - v x
1 - v N
V • • ,
\ \
*"" """' S
•— -•
1 ^"~ — — — ~ f^*"
I
\. ...
x
\
\
\
\
1
^
.
\
1
1 •
\
\
\\
\
1
1 t
/
S
,
.
4 t
\
\' *
i '
I. '.
• f
% »
*
*
4
'
*
*
k
P
Figure 27. Horizontal velocities at 20 m depth.
f
T
i
\
\
t
l
l
V
\
•
\
V
\
\
\
\
*
•
\
I 4 < • • • •
j I » - • - *
I \ ' * * * *
v „,,*,.
\ ^ •*'**••
\ •.»•<-••
\ -.*»•»-•
X H ^ - * * -
- --• ---«.>,
*"••*"-•*->.
-""•""•"•x
"•-••••--'•>%\
• • — - - ^ •» J
~ • ' * i
^
V
\
t
•
*
x
\
\
V
\
\
\
\
^
\
t
r
.
*
'
•
Figure 28. Horizontal velocities at 35 m depth.
66
-------
contours for the same conditions as Case 1, except with a uniform wind,
were calculated. It can be shown that, without the wind stress curl,
there is very little cross-wind variation in the thermocline depth.
The zeroth order solution for the uniform wind case also yields a zero
hypolimnion pressure gradient. Therefore the uniform wind does not
create any zeroth order geostrophic inviscid currents in the hypolimnion
as in Case 1.
67
-------
SECTION VII
TIME-DEPENDENT, CONSTANT DENSITY MODELS
Except that the flow is now varying with time, the models discussed
in this section incorporate the same assumptions as those made in the
previous section, i.e., the water density is constant, the pressure
is hydrostatic, nonlinear acceleration and horizontal diffusion terms
can be neglected, the vertical eddy diffusivity is constant, and
surface displacements are small by comparison with the depth.
The basic model is presented and discussed in Section 7.1. Representa-
tive results are shown for (1) a spatially uniform wind varying as a
step-function with time, and (2) a spatially and time-varying wind with
the wind velocity determined from observations made during Storm Agnes
in June, 1972. Results from the latter case are compared with
observed surface displacements. In Section 7.2, the extension of this
model to take into account a finer grid near shore and coupling with
the coarser off-shore grid is discussed and results of the calculations
are given. With this basic model, some analytic investigations of
the important time scales governing the flow were made. In addition,
the rigid-lid and free surface models were investigated and compared.
These latter investigations are presented in Section 7.3.
7.1 AN OVERALL LAKE ERIE CIRCULATION MODEL
With the approximations discussed above, the equations of motion reduce
to
_ fv = _ + A <32)
at p ax v az2
68
-------
8V -L. f * J2 J. A 9 V /--3-3A
— + fu = T*- + A ——o (33)
3t p 9y v 9z2
^- = - Pg (34)
dz
The boundary conditions are
z = 0: PA ~ = T , pA — = T (35)
v3z x v9z y
w » -p- (36)
z = -h(x,y): u=0, v=0, w=0 (37)
For simplicity, it is also generally assumed in our calculations
(although not necessary for the calculations) that the water is initially
at rest and therefore
t = 0: u = 0, v = 0, w=0, C = 0 (38)
The hydrostatic equation, Eq. (34), can be integrated to give
p = Pa + pg(s-z) (39)
where p is the atmospheric pressure and may be a function of x, y,
3.
and t. From this and neglecting the variations in the atmospheric
pressure, it follows that
^P _ „„ <*£ 9? .3£ ftn\
— PS — > — = Pg \H(J}
9x 3x 3y 3y
and therefore p can be eliminated from Eqs. (32) and (33). It is
convenient to integrate the continuity equation over the depth. The
69
-------
result is
, .
8t 9x 3y U (41)
In many studies of lake and ocean hydrodynamics, the rigid-lid approxima-
tion is invoked. This assumption implies that no vertical motions are
allowed at the surface z = 0. The upper boundary condition, Eq. (36),
then reduces to
z = 0: w = 0 (42)
and the integrated continuity equation, Eq. (41), reduces to
f + U - •
Stretched Coordinate System
It is convenient to transform the equations from an x,y,z system to an
x»y,o system where 0 <^ a <^ 1 and 0 = 0 corresponds to the bottom of
the lake while a = 1 corresponds to the surface of the lake (Freeman,
Hale, and Danard, 1972). With this transformation and a uniform a grid
spacing, the spacing in the x,y,z system is smaller in shallow regions
of the lake and increases in proportion to the depth. This ensures that
in the shallow areas there is no loss of accuracy in the computations
due to lack of vertical resolution. The transformation also makes it
easier to incorporate the bottom topography into the numerical program
and to apply the boundary conditions.
The a coordinate is defined by
a = 1 + (44)
70
-------
where h = h(x,y) is the depth at the location (x,y). With this trans-
formation and the modifications discussed above, the governing
equations reduce to
3u ,. 3? . v 32u //r\
_ _ fv . _ pg_ + ^ _,, {45)
9v . . 3£ . Av 32v (46)
- + fu = - pg + p- ^
IS. + M + II = n (47)
3t 3x 3y
01 01
where U = / udz = h / uda and V = / vdz = h / vda. The boundary
,. . -h 0 -h 0
conditions are
. . 8u . . 3v . (48)
a =1: pA -r— = hT , pA — = hx
v 3a x v3a y
a = 0: u = 0, v = 0 (49)
At the coast, it is assumed that there is no net normal mass flux.
Spatial Grid and Time Differencing
To set up the finite difference scheme, the spatial grid shown in
Fig. 29 is used. The indices i,j,k refer to grid points in the x, y, a
directions respectively. Constant grid spacings Ax, Ay,Aa are used
in the x,y,a directions and, generally, Ax = Ay. The variables are
defined in such a way that the only variable occurring on boundary
segments parallel to the y-axis is U and the variable occurring on
boundary segments parallel to the x-axis is V. Such an arrangement
is convenient since the boundary condition of no normal flow through
the coast involves U or V on appropriate boundary segments.
71
-------
Scc.lt
18
If
-no ,
:LAJ
r 1
--J-*
ii
T
I
*f
k
T
ft'
Y1
T'
4,
V
" J"
A
j
t- -* — :
•4-
*
t — J —
t
L— O— .
— a— i
~t-
J
V
T
~t~-
1
T
"t"
A
,4-:
T
£
1
-r
i.
T
"I"
k — *— -;
1
4-
*
~j '
T
C-f-i
LJH
"t"1
?
I
4-
-A-
l
1
,-^_
14
' j,
1
L — 0— ;
y^
:-
•fl
'1~;
T
r - JL -
'~l
T
^_^_
*U^
AV,i
"5
4-
t- -•--«
— -^-.
rfj
>-i-.
T
i
i-
L
r
V
~J
Y
— A—
-A
r
A -
5
A _
1
d
T
t]
T
^
1,
Figure 29. Horizontal grid pattern in part of Lake Erie.
72
-------
The staggering of the variables in the horizontal plane makes it possible
to use centered differences in space to evaluate derivatives with
respect to x and y. Values for variables are sometimes needed at
points at which they are not defined. These values are obtained by
linearly interpolating from the values immediately adjacent to the point
under consideration. Forward differences in time are used in the
numerical calculations. The vertical diffusion terms in Eqs. (45)
and (46) are evaluated by the Du-Fort Frenkel scheme [Smith, 1965]
to ensure numerical stability.
The numerical procedure can be summarized as follows:
(1) Assume that the values at time t are known (initially at t., ,
n i
u = v = £ = 0).
(2) Use Eq. (45) to calculate un+1.
(3) Use Eq. (46) to calculate vn+1.
(4) Evaluate Un and v by numerically integrating u and v over the
depth.
(5) Set U and V equal to zero on appropriate boundary segments.
(6) Evaluate ?n+1 from Eq. (47).
This completes the calculation at time t _ and the computation for
the next time step can now be started with step 1.
Numerical Stability Considerations
The numerical model consists of a system of 3 differential equations
together with 2 numerical integrations. The determination of a
stability criterion for the complete system is quite difficult and,
in general, simple closed form results cannot be obtained for the
allowable time step. However, by examining the various physical
processes individually, some estimates of the stability limits can
be obtained. This procedure has been found to provide reasonable
guidelines for the choice of the allowable time step At [Crowley,
1968].
The physical processes that place the most stringent limitations on
At are vertical diffusion (particularly in shallow regions of the lake)
and surface gravity waves. The use of a Du-Fort Frenkel scheme
73
-------
makes it possible to overcome the limitations due to viscous diffusion
[Roache, 1972]. The scheme used in the present model is a modified form
of the Du-Fort Frenkel scheme in which forward time differencing is
used instead of a leap-frog scheme. However, even this modified
scheme is found to be stable in practice. The time step is therefore
limited only by surface gravity waves.
The stability criterion introduced due to the surface waves is that
the time step be less than the time it takes for the waves to move
across one grid spacing, the Courant-Friedrichs-Lewy condition
[Richtmyer and Morton, 1967]. A simple stability analysis for these
waves shows that for stability At < Ax/ v'gh, where At, Ax and h are
dimensional quantities [Richtmyer and Morton, 1967]. In practice, for
the full set of equations used in the present model, it is found that
the numerical procedure is stable as long as
At 1 — , (50)
where h is the maximum depth in the lake.
Determination of the Wind Stress
The relationship between the wind stress and the wind velocity is taken
to be that discussed in Section 4.2, i.e.,
T = p C, W W (51)
— a d a —a
In general, the winds are not uniform over the lake. Furthermore,
the wind velocities are only measured at a few points on the shores
of the lake. Therefore the wind stress over the remainder of the lake
has to be determined by interpolating from the known values at a few
points. The method of interpolation used in this study is similar
to that used by Platzman (1963). The stress at any point (x,y) in
74
-------
the lake is assumed to be given by
T(x,y,t) = Z Xm(x,y,Hm(t) (52)
where m is an index identifying a particular station, T (t) is the wind
m
stress at station m at the time t and X (x,y) is a dimensionless
m
interpolating function. T can be determined from Eq. (51) from the
known wind velocity at station m.
To specify the intei
point (x,y), define
To specify the interpolating function X (P), where P denotes an arbitrary
m
am(P) - (x-xm) + (^m> (53)
the square of the distance between P and P . Further, let
m
a (P) (54)
n
then
f (P)
zf~(P) (55)
m m
Instead of defining a (P) as in Eq. (53), any monotone function of the
m
distance between P and P could be used.
m
For each data point P , contours of the interpolating function X (P)
m m
can be determined. Results for the interpolating function obtained
from six data points around the lake are shown in Fig. 30. It can be
seen that the influence of any data point diminishes as the distance
from it increases. In particular, the influence of point 5 on the
stresses in the lake is small, since point 5 is at a larger distance
from the lake compared to the other points. From Eq. (55), it is
evident that the sum of the influences of all the data points at P
is unity, i.e.,
Z X (P) = 1. (56)
m m
75
-------
uJ
Q)
J
-------
Numerical Grid Pattern and Choice of Time-Step
To set up the numerical grid in the lake, the boundaries of the lake
are first approximated by straight line segments parallel to the x and
y axis (see Fig. 30). A uniform horizontal grid with Ax = Ay = 6.4 km
(4 miles) was used. In the vertical direction, the region from the
bottom to the top (a=0 to a=l) is divided into five intervals of length
Aa=0.2.
To determine the allowable time step for the calculations in Lake Erie,
the maximum depth of 67 meters is used. Then, with Ax=6.4 km, Eq.
(50) requires that At be less than or equal to about 50 seconds. With
this time step, the calculations showed no evidence of instability.
Results for a Uniform Wind
As the first example of the time dependent flow in Lake Erie, consider
a uniform wind blowing along the long axis of the lake with T = 1 dyne/
2 x
cm and T =0. According to Wilson's formula for the determination of
the wind stress from the wind speed, this corresponds to a wind
speed of 5.5 meters per sec. The wind direction is W32°S. Calculations
2 2
were made for eddy diffusivities of 25 cm /sec and 40 cm /sec. These
values are approximately the same as the values used by Gedney and
Lick (1972) for comparable winds (but steady state) and for which the
2
numerical results compared well with measured data. For A =25 cm /sec,
v 2
the calculations were carried out to about 57 hours, while for A =40 cm /
v
sec the calculations were terminated at about 21 hours.
2
A representative sample of the numerical results for A =25 cm /sec
is shown in Figs. 31 to 34. From Fig. 31, it can be seen that there
is a dominant period of about 14 hours. Henry [1902] has examined a
series of limnograms at the western and eastern ends of Lake Erie.
He found that the uninodal oscillation of the lake was most conspic-
uous for several days and it had a period of 14 hours. Evidently the
period observed in the numerical results corresponds to these uninodal
oscillations along the long axis, since the wind is in the direction
of this axis. The transverse oscillations and the effects of the
77
-------
o w
co en
z>
o
X
o
in
o
ro
o
OJ
CVI
O ~
o
CVJ
I
(D
I
4J
cfl
a)
a
o
•H
4J
I
cn
cfl
0
a)
P.
ca
•H
a)
o
CO
0)
60
•rl
78
-------
u
-------
1
fi
0)
g
f
o
•H
O
O
O
N
•H
M
s
CO
ro
0)
)-i
3
00
80
-------
CO
tfl
PQ
CU
4-1
CO
rt
w
4-1
o
g
a
a
o
o
o
4J
•H
a
o
rH
(1)
c
o
N
•H
M
O
33
-a-
ro
cu
M
(WDJHidHO
81
-------
irregular boundaries can be seen in some of the smaller oscillations
that are superimposed on the main oscillation.
Contours of the surface displacement at t = 9.5 hours are shown in
Fig. 32. The lowest point is at the western end of the lake and the
highest point is near Buffalo at the eastern end of the lake. The
peak set-up between Buffalo and Toledo (Buffalo lake level-Toledo
lake level) is about 40 cm. In general, except near the northern and
southern coasts and at small times, lines of constant £ are almost
normal to the wind direction. The effects of Coriolis forces and the
bottom topography prevent these contours from being exactly normal
to the wind.
Velocity profiles at two points in the lake are shown in Figs. 33 and
34. The difference in profiles in the shallow (Fig. 33) and deep
(Fig. 34) regions of the lake is evident. In the deeper region, the
top and bottom Ekman layers can be identified. Furthermore, in the
deeper regions, Coriolis forces are large compared to viscous forces,
and therefore the v velocities are large. In contrast to this, v is
much smaller in the Western Basin where the water is shallow and Coriolis
forces are not important. It can be seen (Fig. 33) that in the shallow
region, the approach to steady state is rapid and much less than a day.
In the deeper region (Fig. 34), the approach to steady state is much
slower with a relaxation time on the order of a few days.
Results for an Actual Storm
To verify the numerical model, consider next the case of an actual
storm on Lake Erie for which wind and lake level data are available.
Tropical storm Agnes stalled near Lake Erie and produced northerly
winds over the lake during June 21 to June 25, 1972. Wind data during
the storm is available at six stations around the lake. The wind
magnitude and direction at two different times are shown in Figs. 35 and
36. All of the observed data for the storm has been taken from
Paskausky (1973).
82
-------
en
M
o
,e
o
o
-*
cs
ON
CN
CN
at
c
c
o
CO
-------
UJ
UJ
«J.2
\,
o
ft
o
o
CM
rt
CN
CN
-------
From the wind velocities, the wind stress at the six locations can be
determined. However, the velocities over water are generally greater
than those over land. A factor of 1.0 to 1.5 is usually taken (Reid
and Bodine 1968, Gedney and Lick 1972). In the present case, a factor
of 1.25 gave reasonable agreement with observational results, as will
be seen below. Once the wind stress at the six locations is known,
the stress over the remainder of the lake can be determined by the
interpolation procedure discussed above.
Wind data for the storm is given at hourly intervals starting at 1000
hours on 6/22/72. Examination of the measured surface displacement
curves shows that by this time the lake level had been sufficiently
disturbed and the initial condition of zero surface displacement
everywhere in the lake would be inappropriate if calculations were
started at this time. However, about 10 hours earlier, i.e., at 2400
hours on 6/21/72, the lake level appeared to be relatively undisturbed
and this time was chosen to be the starting point for the numerical
computations. The wind velocities during the period from 2400 hours
on 6/21/72 to 1000 hours on 6/22/72 were determined by linearly
interpolating from zero to the observed value at 1000 hours on 6/22/72.
Furthermore, since only hourly values of the wind stress can be
determined from the observed wind data, any intermediate stress values
were obtained by linear interpolation between two hourly values. The
wind stress as a function of time at Cleveland and Toledo is shown in
Fig. 37. The rapid changes in the magnitude of the wind stress are
evident.
The numerical calculations were carried out with a time step of 50 sec
2 2
and for eddy diffusivities of 25 cm /sec and 40 cm /sec. For
2
A =25 cm /sec, numerical results for the surface displacement as a
function of time at Cleveland and at Port Stanley are shown in Figs.
38 and 39. The results compare favorably with the measured values of
the surface displacement. In particular, the numerical and observed
data have a similar variation.
2
Calculations were also made for an eddy diffusivity of 40 cm /sec.
The results for the surface displacement were practically indisting-
85
-------
CvJ
CO
UJ
2:
Q
CO
CO
LU
CO
O
z:
6 -
H -
o
2HOOHRS 1200 2MOO 1200
6/22/72 6/23/72
Figure 37. Magnitude of wind stress as a function of time during
storm Agnes.
86
-------
CLEVELAND, OHIO
MEASURED
CALCULATED
2400
1200
6/22/72
2400
1200
6/23/72
TIME
(HOURS)
Figure 38. Calculated and measured water level fluctuations at
Cleveland during storm Agnes.
87
-------
PORT STANLEY, ONTARIO
6/22/72
6/23/72
2HOO HRS 1200
2100
I2QO
10 r
0
-10
-20
•TIME
(HOURS)
MEASURED
CALCULATED
Figure 39. Calculated and measured water level fluctuations
for Port Stanley during storm Agnes.
88
-------
2
mshable from those for A =25 cm /sec. This shows that the magnitude
and direction of the wind velocity dominate the flow. Of course, after
the storm, as the surface oscillates and decays to an equilibrium
2
position, the differences between the results for A = 25 cm /sec and
2 v
A =40 cm /sec would be more significant, as indicated by the constant
wind case.
It was mentioned previously that a linear interpolation for the wind
velocity was used for the first ten hours. Varying the magnitude of
this initial wind set-up of course changed the surface displacement
during the first ten hours. However, after this time, the results
were practially independent of what wind velocity was assumed for the
first ten hours.
Further verification of the numerical model was obtained by comparing
the numerical results with the observed water level at four points
around the lake at the time of approximately peak set-up at Cleveland.
This is shown in Fig. 40 where the surface displacement contours
are shown at 0400 hours on June 23. It can be seen that the northern
end of the lake is at a much lower level than the southern end. The
measured values of £ at Port Stanley, Marblehead, Cleveland and
Barcelona are also shown in the figure. The calculated values at these
four points are -12, 5, 22 and 9 cm respectively and these are in
reasonable agreement with the measured values of -15, 6.4, 26, and 12
cm. Thus the spatial variation of £ obtained from the numerical
model is also in accord with observed data.
Although the above numerical procedure is correct, it is not very
efficient. Two alternate and useful procedures have been developed
for modeling of three-dimensional, time-dependent lake currents. Both
are more efficient than the method indicated above and both have been
used previously in atmospheric and oceanic modeling. One method is
that developed for lakes by Simons (1971, 1972, 1974, 1975). In this
method, the external (vertically-averaged) and internal modes of the fl'
field are treated separately. The time step used in the calculM :
of the external mode is limited by the stability requirement.fa indicated
89
-------
o
_J
LU
>
LU
_J
a:
LU
i
0
LU
CC.
3
in
<
LU
CO
_
LU
o
cc
<
CD
c
o
CO
M
O
O
O
4-1
n)
to
0)
C
oc
CD
OJ
Q
2
<
_J
LU
LU
-1
O
O
M-l
c
o
•H
4-1
td
OJ
OJ
u
cfl
M-l
M
P •
to eg
to «
i-i m
3 og
O
•u a)
C C
O 3
0)
M-
3
M
•H
90
-------
above. However, the calculation of the internal structure of the flow
is not limited by these requirements and a much larger time step can
be taken for these calculations. This procedure has been incorporated
into the above model and further computations are presently being made.
The second method (Paul and Lick 1973, 1974) also involves separating
the calculations into external and internal modes. However, the
external mode, is treated by using a rigid-lid assumption, an assumption
which effectively eliminates free surface gravity waves. More details
on this method are given in Section VIII.
7.2 THE NEAR SHORE
The time-dependent model described in the previous section has been
modified so as to describe in more detail the flow in a near-shore area
near Cleveland (Sheng and Lick, 1975). This area is the same as that
described in Section 6.3, where the flow was calculated by means of a
steady-state model. In this near-shore region, a one mile or 1/2
mile grid is necessary for adequate detail. This reduction in grid
size not only increases the number of grid points and therefore the
number of calculations per time step but also decreases the maximum
allowable size of the time step, At , and therefore increases the
max
number of time steps necessary for calculations over a specified time
interval.
It is unrealistic to describe the entire lake by means of a one mile,
or finer, grid. By comparison with the four mile grid used in the
overall lake calculation, a one mile grid would increase the number of
grid points by a factor of 16 and, in addition, reduce At by a
max
factor of 4.
In order to conserve computer time while still retaining the necessary
grid size in the near shore, the lake has been divided into regions,
each with its own grid size and, correspondingly, its own maximum
allowable size of time step. These regions are: (1) The near shore.
The grid size is one mile. The maximum depth is 22 m. This yields
a At of 22 sec. (2) The Eastern Basin. In this basin, an eight
IHclX
91
-------
mile grid is used. The maximum depth is 65m and At = 100 sec,
max
Although the grid is fairly crude, this will have comparatively little
influence on the motion in the near-shore region near Cleveland.
(3) The Central and Western Basins. A four mile grid is used in the
rest of the lake. The maximum depth is 26 m which gives a At of
max
82 sec.
It was decided to use a time step of 20 sec. in the near shore and
80 sec. in the rest of the lake. This satisfies the above require-
ments and, in addition, four time steps in the near shore are then
exactly equivalent to one time step in the rest of the lake, a
computational convenience.
The coupling between regions with different grid sizes and different
time steps can of course create problems. See a discussion of these
problems in Section V. To investigate these problems, the coupling
process was first studied for one and two dimensional systems. From
this experience, the present computational procedure evolved. In this
procedure, it can be shown that: (1) Mass is conserved exactly as
fluid moves from one region to another. (2) Pressure forces (due
to surface displacements) are transmitted smoothly across the interface
between two regions. (3) The results agree with the results obtained
with a single grid size and a single time step throughout the entire
domain, and therefore the results vary smoothly throughout the domain.
The coupling process can be illustrated by a discussion of the coupling
between the Eastern and Central Basins. The grid pattern near the
boundary between the two regions is shown in Fig. 41. The two grids
overlap over a region of four miles. Let all variables in the Eastern
Basin be denoted by the subscript e and all variables in the Central
Basin by the subscript c. The solution procedure is: (1) Along the
vertical line I =1 and I =40, calculate the u 's from the finite
ue uc e
difference equations. The u 's along this same line are then interpolated
c
from the u 's. (2) Along I =1 and I = 40, calculate the v 's
e ve vc c
and the £ 's from the finite difference equations. The v 's and the
£ 's along this line are then interpolated from the v 's and the £ 's.
92
-------
Dve,ve
Oue, Ve
* Cc
•0 vc,Vc
Ju £ Jv
uele ve
\x\\\\v
Iu 39 40
. f 39 40
Figure 41. Grid pattern at boundary between Eastern
and Central Basins.
93
-------
The coupling between the near-shore region and the Central Basin is
treated in a similar manner with obvious changes when the orientation
of the boundary is different.
Calculations are made in a similar manner to that described in the
previous section. Calculations were made initially to verify the
procedure. Representative results are shown in Fig. 42, the surface
displacement as a function of time and as a function of distance along
a line across the lake approximately perpendicular to the shoreline
at Cleveland. For these calculations, the wind stress was taken to
2 2
be T = 1.55 dynes/cm and T =1.1 dynes/cm with the wind starting
impulsively from rest at time zero. The smoothness of the results is
evident. Also noticeable is the 14 hour longitudinal seiche on which
is superposed the 2 1/2 hour transverse seiche. Preliminary results on
the effects of an island jetport have also been calculated (Sheng and
Lick, 1975)
7.3 ANALYTIC RESULTS
General characteristics of the time-dependent flow in a lake were
investigated (Haq, Lick, and Sheng 1974, Haq and Lick 1975) using
primarily analytic methods which were then verified and extended by
numerical methods. The time-dependent flow in a lake or ocean can
be thought of as the flow in a rotating container with a free upper
surface. Many different time scales are present and their relative
importance is determined by the physical parameters appearing in the
problem. In the investigation which is summarized here, time scales
of interest were explicitly identified and simple analytic formulas
for these time scales were derived. In addition, the rigid-lid
and free-surface models were investigated and compared. The analytic
results were verified by comparison with numerical calculations.
For this investigation, a relatively simple model was investigated,
i.e., a constant density, time-dependent, wind-driven flox^ in a basin
of constant depth. The emphasis is on values of parameters (especially
Ekman number) corresponding closely to those in the Great Lakes and in
particular to those in Lake Erie, the shallowest of the Great Lakes.
94
-------
UJ
LU
O
Hco
0)
c
•H
60
(3
O
iH
tfl
cfl QJ
rH
4J U
c
0) 4-1
S o)
n) ^J
tH cfl
CX -H
W
•H Q)
T) &
4J
QJ
u tn
to M
M-( O
M K
3 O
I
95
-------
The basic equations are those presented in Section 7.1. It is assumed
that the wind is spatially uniform and varies as a step function with
time. By means of Laplace transforms, an analytic (but quite complex)
solution can be obtained. The method of solution and the solution
itself can be simplified by using a substitute kernel technique, an
approximation shown to be generally valid for the range of parameters
representative of large lakes.
Results are discussed here for two simple cases, (1) an infinitely
long lake of width L and depth h, and (2) a square basin of length
L and depth h. For an infinitely long lake, it can be shown, for the
rigid-lid model, that the surface displacement decays with a time
scale t = 1/fy2, where y2 is given by
/E
y v 2
2irh
sinh , sin
cl
2irh
cosh .. 4- cos
d
2lTh
d
2trh /E* . ,_ 2trh ^ . 2irh
d y 2 Sinh d + Sin d
(57)
where d is the Ekman layer thickness. For a shallow lake, where
h/d«l, the time scale t is given by t = 4ir2h2/5fd2 and is
22
proportional to h /d or 1/E , i.e., a viscous diffusion time (Greenspan,
1968). For a deep lake, where h/d»l, the decay time is given by
t, = 2irh/fd and is proportional to l/i/E~, i.e., the usual spin-up
time for a rotating container with a rigid lid (Greenspan, 1968). »
For an infinitely long lake with a free surface, the motion is more
complex. Various asymptotic time scales can be identified and depend
2 2
on the parameter $, where g = gh/L f , as well as on the Ekman number.
In the limit as g-*» (the rigid-lid approximation), these times are
t , a diffusion or spin-up time depending on the value of h/d, and
t? = 2L/i/gh", the period of oscillation of surface gravity waves.
In this limit, the amplitude of the solution corresponding to gravity
waves goes to zero so that these solutions are absent from the rigid-
96
-------
lid results. As 3-K), the dominant times are t = 2-rr/f, an
222
inertial period (Proudman, 1952), and t = fL /IT y gh, a flux time scale,
This latter time scale arises from the fact that the lake surface is
displaced by the wind stress and it takes a certain amount of time,
given by t,, to move fluid from one side of the lake to the other under
the action of a wind stress or induced pressure gradient. In a deep
lake, this mass motion occurs through the Ekman layers.
For intermediate values of 3 corresponding to real lakes, the time
scales can be calculated from simple analytic formulas but cannot be
easily interpreted, being combinations of the above times. In general,
the solution consists of (a) a decaying part which has a decay time
varying from t to t, as 3 goes from °° to 0, and (b) an imposed
decaying oscillatory part with a decay time of the order of t., and a
period of oscillation varying from t9 to t. as 6 goes from °° to 0.
As g-*», the free-surface model approaches the rigid-lid model. For 3-*°°
and E «1, the results correspond to those for the spin-up of a
rotating container with a rigid-lid [Greenspan, 1968]. The rigid-lid
and free-surface results for E =1.0 and E = .025 are shown in Figs.
v
43 and 44 respectively. These results are for a lake with width L =
2
300 km, h = 10m, and T = 1 dyne/cm . Analytical results for the rigid-
lid case show that the surface displacement is given by
x 1] . -2.4ft
U 6 } (58)
for E =1.0 and
v
-e-'13ft)
(59)
for E = .025, where U is a reference velocity taken to be 10 cm/sec.
V R
These results are shown in Figs. 43 and 44 and the agreement with the
numerical results is seen to be fairly good.
97
-------
TIME (HOURS)
u
ANALYTICAL RESULT
NUMERICAL RESULT
-10--
RIGID LID
Figure 43. Surface displacement at x = 0 as a function of time
for an infinitely long lake.with L = 300 km, H = 10 m,
E = 1.0, T = 1 dyne/cm2.
98
-------
a)
4-1
•H 03
C 0)
G ct)
•H >
C 4-J
cfl C
-,
a t3
U-l
rH
n)
TD
C
•u CM
tt) o
0)
e
0) w
a
n) ^
iH 4-1
CX-H
cn >
•H
Qj (8
0 rH
fl
M-l M
3
CO
3
60
•H
O "4-4
H o
99
-------
From Figure 43 (for E =1.0), it can also be seen that the decay
time for this case as predicted by the rigidr-.lid model is much shorter
(by about a factor of 4) than that predicted by the free surface model.
However, no oscillations of the free-surface are evident due to the
large viscous damping in the lake. The dominant time scale for the free-
surface case can be shown to be of the order of t - 1.6. In contrast
with this, the decay time of the rigid-lid solution obtained from Eq.
(70) is 0.42.
Fig. 44 shows the results for a lower Ekman number (E = .025) and for
two different values of 3. Consider first the results for 6=0.11,
2
corresponding to a lake with L = 200 km, h = 10 m, and T=! dyne/cm . It
is evident that the decay time predicted by the rigid-lid solution is
once again shorter than the time obtained from the free surface model.
The oscillations of the free surface can be clearly identified. The
analytical results predict a decay time and period of oscillation of
11.5 and 4.5 respectively and these are in satisfactory agreement with
the numerical results. The case of (3=.025 corresponds to a lake with
2
L=630 km, h=10 m, and T=0.475 dynes/cm . The choice of these parameters
ensures that the rigid-lid results remain unchanged. For this small
value of g, the flux time scale t, becomes the dominant decay time.
The decay time and period of oscillation for this case are 31.4 and
5.7 respectively and these are in agreement with the numerical results.
For the time-dependent flow in a square basin, analytic results were
found in the rigid-lid approximation, but will not be discussed here.
Representative and interesting numerical results for free-surface
flows in a square basin are shown in Fig. 45, which gives the surface
displacement as a function of time for values of 3 of 0.98 and 3.92.
2
The wind stress is T =1 dyne/cm , T = 0. The case of 3 = 3.92
x Y
corresponds to L = 50 km and h = 10 m, values representative of the
Western basin in Lake Erie, while 3 = 0.98 corresponds to L = 100 km
and h = 10 m. A beat effect is apparent for £ = 3.92. R>r the free
inviscid oscillations in a lake, it can be shown (Lamb 1932, Golds-
borough 1931, Corkan and Doodson 1952) that two sets of waves are present,
100
-------
cfl
)-l
O
0)
•H
a
o •
•H O
4-1
a ii
§
H-l H
03 *
CN
w e
rt o
10 a)
o c
10
O H
II 10
-------
one rotating in a clockwise direction and the other rotating in a
counter-clockwise direction. The frequencies of these waves depend on
the depth and horizontal dimensions of the basin. For the present case,
the two opposing waves have frequencies which differ by only a
small amount.
For these calculations, the Ekman number was taken to be 0.025 with
2
a corresponding eddy diffusivity of A =2.5 cm /sec. For these same
v 2
conditions but with an eddy diffusivity of 25 cm /sec, a value more
representative of Lake Erie, damping of the surface displacement
is more evident and the beat phenomena is not as apparent.
It can be seen, at least for values of 0 which are typical of the Great
Lakes, that the results for the surface elevation £ are quite different
from those of B-*00, the rigid-lid limit. In particular it can be seen
from the analytic results that, in the rigid-lid limit, both the t and
t time scales are eliminated, i.e., no surface gravity waves or piling
up of water is permitted. Furthermore, the rigid-lid results predict
a much faster rate of approach towards the steady-state compared to the
free-surface results.
102
-------
SECTION VIII
TIME-DEPENDENT, VARIABLE DENSITY MODELS
Various models which include the effects of time dependence and variable
density have been developed by us and some of them are discussed in this
section. The basic assumptions used in all of these models are: (a) The
Boussinesq approximation is valid. This assumes that density variations
are small and can be neglected except in the gravity term. The
coupling between the momentum and energy equations is retained. A
consequence of the Boussinesq approximation is that the energy equation
reduces to a balance between convection and turbulent diffusion.
(b) Eddy coefficients are used to account for the turbulent diffusion
effects in both the momentum and energy equations. The horizontal
eddy coefficients are assumed constant but the vertical eddy coefficients
may vary depending on the temperature gradient and other parameters.
(c) The pressure is assumed to vary hydrostatically. In the problems
discussed here, variations in the vertical velocity are small enough
that the neglected terms in the vertical momentum equation are small
compared to the gravity term. A consequence of this approximation is
that the order of the system of equations is reduced as is the computa-
tional effort required for solution. (d) The rigid-lid approximation
is valid, i.e., w(z=0)=0. This approximation is used to eliminate
surface gravity waves and the small time scales associated with them.
The exclusion of these smaller time scales greatly increases the
maximum time step possible in the numerical computations. When surface
gravity waves are permitted, simple stability theory indicates that
the allowable time step is limited by the Courant-Friedrichs-Levy
condition, i.e., At < hx/vgh, where Ax is the grid size, g is the gravita-
tional constant, and h is the depth. With the rigid-lid approximation,
the limiting time is no longer At but is given generally either by
At = Ax/u, where u is a typical fluid velocity, or by At =Ax/u.,
where u. is the speed of an internal wave given by u. = /Apgh/p
103
-------
and Ap is the density difference across the wave interface. In either
case, At is generally much greater than Atf with the ratio At /At.
usually being on the order of 10 to 100. In the rigid-lid approxima-
tion, only the high frequency surface variations associated with gravity
waves are neglected. Internal variations with time are approximately
correct while the steady-state results are calculated correctly and
are the same as for the free-surface case.
Variable bottom topography is accounted for in the present numerical
model by a stretching of the vertical coordinate proportional to the
local depth, as was done in Section VII. This transformation ensures
that in the shallow areas there is no loss of accuracy in the computa-
tions due to lack of vertical resolution, a significant factor in
near-shore calculations. However, the transformation is not conformal
and the transformed diffusion terms involve cross-derivatives of the
spatial coordinates. For simplicity, the further approximation is
generally made that the diffusion terms containing derivatives of the
depth are negligible with respect to those diffusion terms containing
only the depth. This approximation, which greatly simplifies the dif-
fusion terms, is used in meteorological problems when topographic
variations are included (Phillips 1957, Smagorinsky, et al 1965).
The main applications of these models have been to river discharges
and thermal plumes from power plants flowing into large lakes (Paul
adn Lick 1973 a,b, 1974). Representative results of our calculations
are shown in Section 8.1 for the particular case of the thermal discharge
from the Point Beach power plant on Lake Michigan including realistic
geometry, buoyancy effects, wind stresses, and cross flows in the lake.
The same basic model has been applied to lake circulation problems and
representative results of these calculations are shown in Section 8.2.
Further extensions and applications of these models are presently
being made.
8.1 RIVER DISCHARGES AND THERMAL PLUMES
A great deal of work has been done recently trying to model the hydro-
dynamics of a thermal plume from a power plant (see Policastro and Tokar
104
-------
(1972) for a recent summary and comparison of existing models and field
data) and the similar problem of the discharge from a river. Most of
these investigations are by means of integral methods or rather simple
numerical models. Work similar to ours, i.e., three-dimensional numerical
models, has been done by Waldrop and Farmer (1974 a,b). However, their
model has a free surface (and therefore small limiting time steps At- £
Ax//gh ), no vertical stretching of the coordinates, and inaccurate
representation of the time behavior (because of approximations made to
expedite convergence).
Calculations have been made by us of simple river-lake flows with para-
meters typical of the Cuyahoga River entering Lake Erie (Paul and Lick
1973 a,b, 1974) and of the discharge from the Point Beach power plant
on Lake Michigan. Representative results of these latter calculations
are discussed here. Two cases are presented. These are the discharge
into a quiescent lake and the discharge into a lake with a cross-flow
and a surface wind present. The relevant parameters, based on field
measurements taken at the site of the outfall (Frigo, Frye, and
3
Tokar, 1974), are as follows: flow rate, 24.7 m /sec; outfall width,
10.8 m; outfall depth, 4.2m; ambient lake temperature, 8.3°C; discharge
temperature, 17,5°C; maximum discharge velocity, 0.9 m/sec. The bottom
topography is shown in Fig. 46. The outfall extends into the lake and
the discharge forms a 60° angle with the shore.
In presenting the basic equations, it is convenient to introduce the
following dimensionless variables:
wb
u* = — , v* = — , w* = —r-
"o «b uoho
AHt
x*= f ' y*= f > z*= t , t* = -^
bo bo ho b02
P"P
p* - r *^: 2 p* = u ' T* =
PghQFr PQ »
105
-------
-
in
LU
£t
o
X
V)
c
n)
I
a
O
n)
-------
,
Vo
A
Pr = -I*
K
u.
Fr
, U = u(x=0, y=0, z=0)
w - a
u _3h v 3b\
3s 9y
dc
dt
where b and h are the width and depth respectively of the discharge,
p is the density of the lake water and T is the air temperature.
vJ Ct
In the following, only non-dimensional variables will be used unless
otherwise explicitly noted. Therefore, in order to simplify the nota-
tion, the asterisks on these non-dimensional variables will be dropped.
The resulting system of transformed equations is then:
1 3(hu) , 1 3(hv) _9fi
h 3x h 3y 80
(60)
Re
Fr^
-dh / APdo
o
1
h
3h
9x
3 (huv)
3y
_
Ap
h
h
4
U
U
1
h
f\
3(ftu/
' 3a
•
r /
3 h3u
3x 1 3x
- R v = - -
o
/ \]
3 h3u
3y 1 3y I
. \
1 3 IY 3u
1
,
iz 3a \ 3a
\ /
3x
(61)
107
-------
I 9 (huv)
h 3x
Re
3_
3y
h/
1 3(hv )
h 3y
3 a
«
a 3h An
u Ap
, 1
' h
3 1 h3v
3x\ 3x
3
3y
h3v
i 9y,
i J
(62)
Pr-2^
-i. 5 (huT)
3x
I
h
h
+
3x
3x
3(hvT)
3y
h3T
3y
(63)
1 j*£ Re , .
h-37 = 7-2 (1 + p)
(64)
P=P(T)
The rigid-lid condition is difficult to apply in a numerical solution
of the above system of equations. To alleviate this difficulty, an
additional equation for the pressure which contains the rigid-lid
condition can be derived. This is accomplished by taking the diver-
gence of the vertically integrated horizontal momentum equations and
using the vertically integrated continuity and hydrostatic pressure
equations. The resulting equation then has the general form
108
-------
3 WP
s
3x I 3X / 8y
h P
T-H = F(u,v,w,T) (65)
oy
The term P is the integration constant resulting from the vertical
S
integration of the hydrostatic pressure equation and is the surface
pressure, i.e., the pressure at z = 0. This surface pressure can be
interpreted as a pressure due to a height of water above or below the
surface z = 0. In this way, surface displacements (neglecting the
transient motion due to gravity waves) can be calculated.
In the above equation, the time derivative of the vertical velocity at
the surface appears. The rigid-lid condition implies that this term
should be zero. However, due to numerical errors, some non-zero values
of this quantity are generally obtained. It can be shown that this
indicates that the continuity equation is not satisfied exactly. These
errors can accumulate with time if not accounted for properly. A
corrective procedure due to Hirt and Harlow (1967) is used to eliminate
these accumulative errors.
The boundary conditions are as follows. The inlet velocity and tem-
perature profiles are specified. The inlet velocity profile is a smooth-
ed average of that measured at the outfall while the inlet temperature
was taken as the constant value measured. The bottom and the shore
are taken as no-slip, impermeable, insulated surfaces. A heat transfer
condition proportional to the temperature difference (surface temperature
minus equilibrium temperature) and a wind-dependent stress are imposed
at the surface. The surface heat transfer condition is calculated from
the work of Edinger and Geyer (1965) while the stress acting on the
water surface due to the wind is calculated by the formulae developed by
Wilson (1960). The equation of state is obtained by a least squares
curve fit of density vs. temperature data.
Normal derivative pressure boundary conditions are derived from the
appropriate vertically integrated momentum equation. The outer x and
y boundaries for the numerical calculations must be taken at some finite
distance. The boundary conditions used here are that the normal
109
-------
derivatives of the velocities and the temperature are zero.
The values for the eddy coefficients were chosen as follows. From out-
fall measurements, the distance is determined that it takes the velocity
to decay to one-half of the inlet value. A reasonable value for the
horizontal coefficient is picked using previous experience. The ratio
of the horizontal to vertical coefficient is then chosen so that it gives
the same velocity decay as for the two-dimensional boundary layer jet
with bottom friction. The vertical eddy coefficient is taken as
dependent on the local vertical temperature gradient (see Eqs. (15)
and (16)). The dependence on temperature is the same as that suggested
by Munk and Anderson (1948) and, for the present case, is given by
aT
Av - a-gf (66)
where a and 8 are constants depending on the local conditions. The
constant a is equal to the vertical eddy diffusivity A under neutral
stability conditions and was taken to be 50 cm /sec while 3 was assumed
3
to be 200 cm /°C sec.
The equations are put into appropriate finite difference form in both
space and time. The horizontal velocities are defined at integral nodal
points, the vertical velocity is defined at half-integral nodal points,
the temperature is defined at half-integral nodal points in the horizontal
and integral nodal points in the vertical, and the surface pressure
is defined at half-integral nodal points in the horizontal. In the
derivation of the finite difference equations, variables are sometimes
required at points where they are not defined. In these circumstances,
the undefined quantity is taken as the simple average of the neighbor-
ing values.
The finite difference approximations to the equations are derived by
integrating the equations over nodal cells. Either the mid-point or
trapezoidal integration rule is used to evaluate these integrals.
Part of the rationale behind the arrangement of the variables is to
provide cells which lend themselves to easy use of the integration rules.
110
-------
The Euler explicit time scheme, where the time derivative is written
as a simple forward time difference and the rest of the equation is
evaluated with the previously calculated values, is used exclusively
in the present model. The three time level explicit time scheme used
by Paul and Lick (1973b) has been found to be not as stable as the
Euler explicit scheme. Even the use of the filtering procedure suggest-
ed by Bryan (1969) for the three time level scheme did not provide
enough stability when inertia dominated flows were calculated.
At each time step, the Poisson equation in the surface pressure has
to be solved. For this purpose, the alternating-direction-implicit
(ADI) method (Varga 1962, Wachspress 1966) is used and was found to be
particularly efficient.
After the temperatures are calculated by the explicit time scheme,
they are checked to see if static stability is satisfied, i.e., if the
temperatures decrease monotonically downward (assuming that density
increases with increasing temperature). When a static instability is
encountered, infinite mixing is assumed which just averages the temperature
over any unstable region.
The above procedure is carried out until the desired solution is obtained.
In both cases discussed below, a variable grid was used with the grid
size varying from about 2 m to about 100 m.
For the first case of the discharge into a quiescent basin, representa-
tive results of the calculations are shown in Figs. 47 to 49. Figs. 47
and 48 show the horizontal velocities at the surface and at a depth
of 2.1 m respectively while Fig. 49 shows the surface temperature
distribution. As the flow is discharged from the outfall, it is forced
towards the surface by the decreasing depth (seen in Fig. 46) as well
as by buoyancy. The large Froude number (Fr = 4.2) indicates that at
least initially buoyancy effects are not dominant. The higher
surface velocities on the left of the outfall centerline are due to the
fact that the depth decreases more rapidly in that region. After about
111
-------
O)
O)
E
in
ID
OJ
LO
en
CO
o
M
o
o
c
C
•H
CO
Q)
•H
4J
•H
O
O
O
cfl
M-l
3
bO
•H
112
-------
CO
o:
LJ
H
LJ
O
in
in
CJ
L-o
I
iH
M-l
Cfl
CO
O
M
O
I
CN
M-l
O
j:
4J
ft
0)
4-1
rt
to
0)
•H
4-1
•H
O
O
H
-------
ro
ro
4-1
I
to
U)
O
M
O
O
I
c
o
•rl
V4
4J
«)
•H
13
3
4J
nl
H
0)
QJ
O
cfl
U-t
QJ
S-i
3
60
114
-------
75 m from the outfall, the depth begins to increase. After this point,
the velocities at the 2.1 m depth remain small, with most of the flow
occurring near the surface.
A second calculation was made for a case when a cross-flow (current of
9.1 cm/sec) and an off-shore wind (approximately 5 m/sec) were present.
Representative results of the calculations are shown in Figs. 50 and
51, Fig. 50 shows the surface velocities and it can be seen that the
discharge is swept in the direction of the cross-flow. This is even
more evident in Fig. 51 which shows the surface temperatures. It can
be seen (compare Fig. 49) that the temperatures are displaced in the
direction of the cross-flow and are also displaced towards shore.
A comparison of the calculated results with field observations is shown
in Figs. 52 to 55. Figs. 52 and 53 are for the first case above and
show the temperature decay along the centerline and the isotherm
area for various temperatures. It can be seen that there is good
agreement between the predicted results and field observations.
Similar results are shown for the second case discussed above in
Figs. 54 and 55. Again good agreement between the predicted and obser-
ved values is evident. The calculations are presently being extended
to include more of the flow field.
Additional field data is available (Frigo, Frye, and Tokar, 1974) from
which one can determine more details of the flow field. However, the
flow, due to its turbulent nature (see Csanady (1973) for a general
description of the turbulent diffusion and nature of plumes), is highly
variable both in space and time. Continuous field measurements over
a sufficiently long period of time to average out these variations
must be made before more general comparisons can be made between our
calculations and observations. This has not been done yet. However,
the above calculated results, although limited, seem more than adequate
at this point.
115
-------
o
V
o
in
in
'CVJ
Lo
g
A
t'
0 0
V \ \'
v V . \l
V vV
' 0 V \(
*$/,/ .••'
r ^ ^ /.
xv/
w
CO
o
5-1
O
G
•H
^3
4-1
•H
&
co
-------
h-
I
to
CO
O
C
cfl
TJ
C
•H
f
•W
•H
3
4-1
cfl
S-i
CU
I"
01
cu
o
cfl
<4-(
M
3
CO
3
Ml
•H
117
-------
s\
D
2
D 9.
< 0
o
UJ
cr
D CL
< d
o
0
5
D
i i i i
q co
-------
tO
O
CO
•z.
o
a
o
o
< UJ
i-i cc
O UJ
_J O
UJ O
Lu 5:
10
O
en
en
o
o
fi
UJ
TJ
C
•H
!s
0
D
ro
O
(Y- w
"- cd
UJ D
i a
O g
(/> CD
4-1
O
tn
•H
a)
o
cd
cvJ
O
O
(X)
tO
OJ
cu
S-i
bO
119
-------
§
o
o
03
CO
o
o
a
c
CO
t)
a
2
O
tr
LJ
O
-O
rO
o
0 ii-
0 a
O
-O
CJ
O
-O
a:
UJ
UJ
UJ
o
I
to
o
UJ
2
-------
ID
r-9
D
m
_ O
D
a
z
o
UJ
Q O
i
a
r
o
CM
- UJ
2 5
tr
UJ
x
H-
o
tn
(O
-O
OJ
CO
o
i
en
en
o
G
cd
•a
•H
en
cd
a)
t-i
cd
a)
J3
W
O
en
•H
a)
o
LO
LO
3
W)
<3 <
121
-------
8.2 LAKE CIRCULATION MODELS
The most extensive work on three-dimensional, variable density models
for lakes is that due to Simons (1974, 1975), who used a free-surface
model. Extensive numerical computations have been made and the results
compared in detail with observations for two physical events. The first
was a calculation of the time-dependent flow during tropical storm
Agnes in June, 1972 when Lake Ontario was only weakly stratified
(Simons, 1974) and the second was a calculation of the flow in August,
1972 when Lake Ontario was strongly stratified (Simons, 1975). For these
calculations, the horizontal grid spacing was 5 km while the vertical
resolution consisted of four layers separated by horizontal levels
at 10, 20, and 40 m below the surface to approximate the location of
current meters used in the field studies.
For numerical efficiency, the calculations for the external and internal
modes were treated differently. This was done as follows. For each of
the four layers, the equations of motion were written in conservative
form similar to Eqs. 60 to 64. By summing these equations over the
four layers (equivalent to vertically averaging over the depth), one
obtains the equations for the external mode. In this manner, it is found
that much larger time steps can be taken for the internal mode computations
than those for the external mode computations. For Lake Ontario, a
time step of 100 sec was used for the external mode and a time step
of 15 min was used for the internal mode.
In the first calculation, the flow was essentially at constant
density. The vertical eddy viscosity was taken to be proportional to the
wind stress and of the form A = A + KT/P where A and K are constants.
v vo vo £
Good agreement with measurements was obtained by taking A = 25 cm /sec
and K = 100 sec.
In the second calculation, the flow was strongly stratified and variable
density effects were appreciable. However, turbulent diffusion of
momentum was treated in the same manner as for the first (constant density)
calculation with no dependence of the eddy viscosity on temperature
gradient. Turbulent diffusion of heat was not included in the calcu-
122
-------
lations but was determined by a comparison of the predicted and observed
results. In general, the predictions of surface water levels and currents
were quite good. Stratification was shown to have an appreciable effect
on the temporal variations of currents. Phenomena such as internal
waves were not simulated in a satisfactory manner because of the choice
of eddy coefficients. Temperature predictions showed general agreement
with large scale averaged observations. However, because of the neglect
of turbulent heat transfer, the detailed temperature distribution was
not given adequately.
From the differences in the calculated and observed temperature
distributions, a lake-wide diffusive heat flux was calculated. From
this, an effective eddy conductivity was determined. For this event,
2
effective eddy conductivities were about 1 cm /sec. For an earlier
period during Storm Agnes, effective eddy conductivities were determined
2
to be approximately 2 to 4 cm /sec, indicating the effect of higher
wind stresses on producing turbulence.
The rigid-lid model described in Section 8.1 has been used by us to
calculate the circulation in simple basins and in Lake Erie, especially
in the near-shore including a description of the flow in Cleveland Harbor.
More extensive reports on these calculations are forthcoming and they
will not be discussed here. Numerical experiments are presently being
performed in an attempt to understand the formation and decay of a
thermocline in a large lake and the accompanying modification to the
dispersion of a contaminant because of this stratification. Some
representative results of the flow in the Central Basin of Lake Erie
are discussed below.
For this Central Basin calculation, for convenience, the geometry was
approximated and the Central Basin was taken to be 60 miles wide and
100 miles long. The actual bottom topography was used except near shore
where a constant 6 ft. depth was assumed. The temperature was assumed
initially to be 75° F in the top 9.15 m and 55°F in the layer below that.
A variable eddy viscosity was assumed of the form described by Eq. (66)
123
-------
2 3
with q = 16,8 cm /sec and & = 84.0 cm /?C sec. Horizontal mixing was
neglected since this is, only important in narrow layers (less than 2
km wide) near shore,. A constant wind of 5,4 tn/sec from the South was
assumed.
Representative results of the calculations after approximately 14 hours
are shown in Figs. 56 to 59, Fig. 56 shows the horizontal velocities
at the surface, a generally uniform flow to the right of the wind and
along the longitudinal axis of the lake as might be expected. Varia-
tions near shore are due to local topography and depth variations.
Fig. 57 shows the velocities of a vertical and longitudinal plane
approximately in the center of the basin. The circulation consists
of one large cell in contrast to the two counter-rotating cells obtained
when a two-layer model is assumed (see Section 6.4 and Figs. 24 to 28).
In the two-layer model, it is assumed that there is no flow between the
epilimnion and hypolimnion while in the present case there is considerable
flow vertically at all depths. The two-layer model is valid in the limit
as the turbulent mixing and convection between the two layers becomes
negligible.
This turbulent mixing depends nonlinearly on the temperature gradient. In the
present calculation, this temperature gradient is small in the upper half
of the basin but much larger in the lower half of the basin as can be seen in
Figs. 58 and 59, which show the temperatures in vertical planes through the
center of the basin. Table 2 shows the correspondence between the temperatures
and the letters designating constant temperature lines in Figs. 58 to 61.
The strong stratification is apparent and reduces vertical mixing
appreciably. This stratification depends of course on the initial
conditions assumed but is strongly affected by the variable eddy coef-
2
ficients. For a constant eddy viscosity and conductivity of 16.8 cm /sec
equal to A in neutral stability, the same calculation produces the
temperature variations shown in Figs. 60 and 61. The smaller temperature
gradients and stronger mixing in this case are evident. Although the
temperature variations are strongly affected, the velocities are not
greatly modified. Qualitatively, the velocity variations are the same
as for variable eddy coefficients. Smaller vertical velocities are
124
-------
a
0
•H
m
Ul
_i
I
O
ck
o
tu
IK
8
.
.0
-lo o
CM/SEC.
o
f
Jo
U
U >-UI
z zo
< tu =>
*- a: t-
- 55
\\
\ \
\ \
\
\
\
\
^
\
\\
\
\
I
\
>\\\
\
I
1
II
\\
\\
\\
\\\
» \
\ ,
1
1 l\x
I
\\
\ 1 1 \\\
•H 0)
t3 ^O
•r) C
ti o
•H -iH
o
•H
5-1
0)
O
n)
4-t •
j_i o
co en
cU S
X!
4-J CNl
•U uO
cfl
CO
-------
-
/
- x \ \ *
,
1
u
HI
5 o
£ £
u
1
u
z§
Sc —
§1
1 / 1 I *
f
1 . 1 1 I » »
' 1 1 1 I I •
1
' I 1 I I •
•»!»..
J
* I 1 ' 1 » '
.
'I*/'*.
\^xV , , •
c
n!
CO
•H
X
nj
0 C
•n -H
d co
6 rt
M
60
C H
O cfl
HI_i
H
0) 0)
c u
cB
tH (U
a js
4J
Cfl M-l
U 0
•H
•U ^J
M CU
C1J 4-1
> a
(U
rj M
rj
•H
ai 4-j
•H
•U rC
•H 60
0 3
O O
CU J2
r> 4-1
P>
m
01
i-r
l^f\
1-t
126
-------
0 0
Figure 58. Temperatures in vertical plane along minor axis
and through the center of the Central Basin.
127
-------
g
O
I
O
J-j
cn
•H
X
cs
60
O
C
(8
n)
o
•H •
•M C
Vt -H
0) CO
> tfl
pa
CO S-l
-------
K I G
L J H F E
D DEFGH
Figure 60. Temperatures in vertical plane along minor axis and
through the center of the Central Basin. Constant
eddy coefficients.
129
-------
60
o w
t-t 4->
J3 C
4-1 01
•H
T) O
a -H
Cfl M-l
IM
CO 01
•H O
X O
n)
o -a
f-l QJ
C
60 tfl
C -w
o co
H C
tfl O
O
CU
(0
iH Cfl
tfl fQ
o
•H rH
4-1 Cfl
)-l M
01 4-1
> C
Q)
c o
•H
01
CO ^2
O) 4-1
4-1 O
0)
M )-i
a) a)
B
o>
H
01
$-1
oo
•H
130
-------
Table 2. LABELS FOR CONSTANT TEMPERATURE LINES
Contour Label
A
B
C
D
E
F
G
H
I
J
K
L
M
N
0
P
Contour Value
°C
14.1
14.7
15.4
16.0
16.6
17.3
17.9
18.6
19.2
19.8
20.5
21.1
21.8
22.4
23.0
23.7
noticeable with the maximum decrease on the order of 20%.
Strong horizontal temperature stratification occurs, especially near shore,
due to the fact that surface near-shore waters only mix with near-surface
waters while off-shore the surface waters mix with cooler bottom waters.
The asymmetry in the temperature profiles is due to convection with
downwelling occurring on the North shore and upwelling occurring on the
South shore.
More work is presently being done to investigate the transition from a
well-mixed flow to a strongly stratified flow. In particular, the near-
shore areas where upwelling and downwelling are important are being
investigated in detail.
131
-------
SECTION IX
REFERENCES
Allender, J.H. 1975. Numerical Simulation of Circulation and Advection-
Diffusion Processes in Saginaw Bay, Michigan. Ph.D. Thesis. University
of Michigan, Ann Arbor,
Berdahl, P. 1968. Oceanic Rossby Waves, A Numerical Rigid-Lid Model.
Lawrence Radiation. Laboratory, University of California, Livermore.
Report ITD-4500, UC-34.
Bonham-Carter, G., J.H. Thomas and D. Lockner. 1973. A Numerical
Model of Steady Wind-Driven Currents in Lake Ontario and The Rochester
Embayment Based on Shallow-Lake Theory. University of Rochester,
Rochester, N.Y.
Bonham-Carter, G. and J.H. Thomas. 1973. Numerical Calculation of
Steady Wind-Driven Currents in Lake Ontario and the Rochester Embayment.
Proc. 16th Conference on Great Lakes Research. International Association
for Great Lakes Research, p. 640-662.
Boyce, F.M. 1974. Some Aspects of Great Lakes Physics of Importance
to Biological and Chemical Processes. J. Fish. Res. Bd. Can. 31:689-730.
Browning, G., H. Kreiss, and J. Oliger. 1973. Mesh Refinement. Math.
Comp. 27:29-39.
Bryan, K. 1969. A Numerical Method for the Study of the World Ocean.
Comp. Phys., Vol. 4.
Casey, D.J. 1965. Lake Ontario Environmental Summary. U.S. Environ-
mental Protection Agency, Region V, Chicago.
Chen, J.H. and K. Miyakoda. 1974. A Nested Grid Computation for the
Barotropic Free Surface Atmosphere. Mon. Weather Rev. 102:181-190.
Corkan, R.H. and A.T. Doodson. 1952. Free Tidal Oscillations in a
Rotating Square Sea. Proceedings Royal Society, Ser. A:215.
Crowley, W.P. 1968. A Global Numerical Ocean Model, Part I. J.
Comp. Phys., Vol.3.
Csanady, G.T. 1973. Turbulent Diffusion in the Environment. D, Reidel
Publishing Company, Boston.
132
-------
Donelan, MI.A,,, F.C, Elder and P,F. Hamblin. 1974. Wind Stress from
Water Set-up, Proceedings 17th Conference, on Great Lakes Research.
International Assocarion for Great Lakes Research. p, 778-788,
2222
Douglas, J. 1955, On the Numerical Integration of 3 u/9x + 3 u/3y =
3u/3t by Implicit Methods. J. Soc. Indus. App. Math, 3:42-65,
Edinger, J.E. and J.C. Geyer, 1965. Heat Exchange in the Environment.
Edison Electric Institute, New York. Publication No. 65-902.
Freeman, N,G., A.M. Hale and M.B. Danard. 1972. A Modified Sigma
Equations Approach to the Numerical Modeling of Great Lakes Hydrodynamics.
J. Geophys. Res, 77:1050-1060.
Frigo, A., D. Frye, and J.V. Tokar, 1974. Field Investigations of
Heated Discharges from Nuclear Power Plants on Lake Michigan. Argonne
National Laboratories, Argonne, Illinois. Report ANL/ES-32.
Galloway, F.M. and S.J. Vakil. 1975. Criteria for the Use of
Vertical Averaging in Great Lakes Dispersion Models. Abstract, Proceed-
ings 18th Conference on Great Lakes Research.
Gedney, R.T. and W. Lick. 1970. Numerical Calculations of the Steady
State, Wind Driven Currents in Lake Erie. Proceedings 13th Conference
on Great Lakes Research, International Association for Great Lakes
Research. p. 829-838.
Gedney, R.T. and W. Lick. 1971. Numerical Calculations of the Wind-Driven
Currents in Lake Erie. Case Western Reserve University, Cleveland,
Ohio.
Gedney, R. and W. Lick. 1972. Wind-Driven Currents in. Lake Erie. J.
Geophys. Res. 77:2714-2723.
Gedney, R.T., W. Lick, and F.B. Molls. 1972, Effect of Eddy Diffusivity
on Wind-Driven Currents in a Two-Layer Stratified Lake. National
Aeronautics and Space Administration. TN D-6841.
Gedney, R.T., W. Lick and F.B. Molls. 1973. Effect of Bottom Topography,
Eddy Diffusivity, and Wind Variation of Circulation in a Two-Layer
Stratified Lake, National Aeronautics and Space Administration, TND-
7235.
Gedney, R.T., W. Lick, and F.B. Molls. 1973. A Simplified Stratified
Lake Model for Determining Effects of Wind Variation and Eddy Diffusivity.
Proceedings 16th Conference on Great Lakes Research. International
Association for Great Lakes Research, p. 710-722.
Goldsborough, G.R. 1931, The Tidal Oscillations in Rectangular Basins,
Proceedings Royal Society, Ser. A:132.
133
-------
Greenspan, H.P. 1968. The Theory of Rotating Fluids, Cambridge
'' - 1, /er s i t y Press. London,
lUnblin, P,,F,. 1969, Hydraulic and Wind-Induced f'irc'il-Jtlon in a Model
of Great. Lake. Proceedings 12th Conference on '^rear Lake Research.
International Association for Great Lakes Research,
Hamblin, P.""'. 1971. An. Investigation of Horizontal Diffusion ixi Lake
Ontario, Proceedings 14th Conference Great Lakes Research. International
Asso'.iation for Great Lakes Research, p. 570-577.
Haq, A, and W. Lick. 1975. On the Time-Dependent Flow ia a Lake. J.
v'Jeophys. Res. 80:431-437.
Haq, A. W. Lick, and Y.P. Sheng. 1975. The Time-Dependent Flow in
Large Lakes with Application to Lake Erie. Case Western Reserve
University, Cleveland, Ohio.
Henry, A.J. 1902. Wind Velocity and Fluctuations of Water Level of
Lake Erie. Bulletin. U.S. Weather Bureau,
Hirt, C.W. and F.H. Harlow. 1967. A General Corrective Procedure for
the Numerical Solution of Initial-Value Problems. J, Comp. Phys,,
Vol. 2.
Janowits, G.S. 1970. The Coastal Boundary Layers of a Lake when
Horizontal and Vertical Ekntan Numbers are of Different Orders of Magnitude.
Tellus, Vol. 22.
Janowitz, G.S. 1972. The Effect of Finite Vertical Ekman Number on the
Coastal Boundary Layers of a Lake. Tellus, Vol. 24.
Koh, R.C.Y. and Y.P. Chang. 1973. Mathematical Model for Barged Ocean
Disposal of Wastes. U.S. Environmental Protection Agency, Corvallis,
Oregon. EPA-660/2-73-029.
Lamb, H. 1932. Hydrodynamics. Dover Publications.
Matsuno, T. 1966. Numerical Integrations of the Primitive Equations by
a Simulated Backward Difference Method. J. Met. Soc. (Japan) 2:76-84.
Munk, W. H. and E.P. Anderson. 1948. Notes on the Theory of the
Thermocline. J. Mar. Res. 1:276-295.
Okubo, A, 1971. Oceanic Diffusion Diagrams, Deep-Sea Res. 18;789-
802.
Olsen, F.C.W. 1950. The Currents of Western Lake Erie. Ph.D. Thesis,
Ohio State University, Columbus.
134
-------
Oort, A.M. and A. Taylon. 1969. On the Kinetic Energy Spectrum Near
the Ground. Hon.. Weather Rev. 97:523-636,
Orlob, C.T. 1959. Eddy Diffusion in Homogeneous Turbulence. J.
Hyd. Div. Proceedings ASCE, HY9:75-101.
Paskausky, D.F. 1971. Weather Circulation in Lake Ontario. Proceedings
14th Conference on Great Lakes Research. International Association for
Great Lakes Research, p. 593-606.
Paskausky, D.F. 1973. Two Dimensional Prediction of Storm Surge on
Lake Erie, Appendix 2g. In: Preliminary Safety Analysis Report,
Perry Nuclear Power Plant. Cleveland Electric Illuminating Company,
Perry, Ohio.
Paul, J.F. and J. Prahl. 1971. Investigations of a Constant Temperature
Rectangular Jet. Case Western Reserve University, Cleveland, Ohio.
Paul, J.F. and W. Lick. 1973a. A Numerical Model for a Three-Dimensional,
Variable-Density Jet. Proceedings 16th Conference on Great Lakes
Research. International Association Great Lakes Research, p. 818-830.
Paul, J.F. and W. Lick. 1973b. A Numerical Model for a Three-Dimensional,
Variable-Density Jet. Technical Report. Case Western Reserve University,
Cleveland, Ohio.
Paul, J.F. and W. Lick. 1974. A Numerical Model of Thermal Plumes and
River Discharges. Proceedings 17th Conference on Great Lakes Research.
International Association for Great Lakes Research, p. 445-455.
Paul, J.F. 1975. A Note on the Vertically-Averaged Equations of Motion.
Case Western Reserve University, Cleveland, Ohio.
Peaceman, D.W. and H.H. Rachford. 1955. The Numerical Solution of
Parabolic and Elliptic Differential Equations. J. Soc. Ind. App.
Math. 3:28-41.
Phillips, N.A. 1957. A Coordinate System Having Some Special Advantages
for Numerical Forcasting. J. Meteor., Vol. 14.
Platzman, G.W. 1963. The Dynamic Prediction of Wind Tides on Lake
Erie. Meteor. Mono. 4:26.
Policastro, A.J., and J.V. Tokar. 1972. Heated-Effluent Dispersion in
Large Lakes: State-of-the-Art of Analytical Modeling, Part I. Argonne
National Laboratory, Argonne, Illinois. Report ANL/ES-11.
Proudman, J. 1953. Dynamical Oceanography. John Wiley and Sons, Inc.
135
-------
Rao, D,B. and T.S. Murty, 1970. Calculations of tbe Steady State Wind-
Driven Circulation in Lake Ontario. Arch. Meteor. Oeophys. Bio-
chira- A 14
Reid, R..O, and 3.R. Bodine. 1968. Numerical Model for Storm Surge?
in Galvt-ston Bay. J. Waterways Harbors. (Div. ASCE) 94:33-c-7.
Richtmyet, R,D, and K.W. Morton. 1967, Difference Methods for Initial-
Value Problems. Interscience Publishers.
Roache, P,J, 1972, Computational "Fluid Dynamics. Hermosa Publishers*
Albuquerque.
Sheng, Y.P. and W. Lick. 1973. Wind-Driven Currents in a Partially
Ice Covered Lake. Proceedings 16th Conference on Great Lakes Research.
International Association for Great Lakes Research. p. 1001-1008.
Sheng, Y,p, and W. Lick. 1975. The Wind-Driven Currents and Contaminant
Dispersion in the Near-Shore of! Large Lakes. Case Western Reserve
University, Cleveland,
Simons, T.J. 1971. Development of Numerical Models of Lake Ontario.
Proceedings 14th Conference on Great Lakes Research. International
Association for Great Lakes Research, p. 655-669.
Simons, T.J, 1972. Development of Numerical Models of Lake Ontario.
Proceedings 15th Conference on Great Lakes Reserach. International
Association of Great Lakes Reserach. p. 655-672.
Simons, T.J. 1974. Verification of Numerical Models of Lake Ontario,
Part I: Circulation in Spring, Early Summer. J. Phys. Ocean. 4:507-
523.
Simons, T.J. 1975. Verification of Numerical Models of Lake Ontario,
Part II; Stratified Circulations and Temperature Changes. J. Phys.
Ocean. 5:98-110.
Simons, T.J. 1976. Continuous Dynamical Calculations of Water Transports
in Lake Erie in 1970. J, Fish. Res. Bd. Can. To be published January
1976.
Smagorinsky, J. , G. Manabe and J.L. Holloway, 1965. Numerical Results
from a Nine-Level General Circulation Model of the Atmosphere. Mon.
Weather Rev. 93:1965.
Smith, G.D. 1965. Numerical Solution of Partial Differential Equations.
Oxford University Press.
Stoimnel, H. 1949. Horizontal Diffusion Due to Oceanic Turbulence.
J, Mar, Res, 8:199-225.
136
-------
Sundaram, T.R, and II. G. Rehm. 1970. Formation and Maintenance of
Thermoclines in Stratified Lakes Including the Effects of Power-
Plant Thermal Discharges. AIM Paper No. 70-238.
Varga, R.S. 1962,. Matrix Iterative Analysis. Prentice-Hall, Inc.
Wachspress, E.L. 1966. Iterative Solution of Elliptic Systems,
Prenctice-Hall Inc.
Waldrop, W.R. and R.C. Farmer. 1974a. Three-Dimensional Computation of
Buoyant Plumes. J. Geophys. Res. 79:1269-1276-
Waldrop, W.R. and R.C, Farmer, 1974b, Thermal Plumes from Industrial
Cooling Water, Proc. 1974 Heat Transfer and Fluid Mechanics Institute.
Stanford University Press, Stanford.
Welander, P. 1957. Wind Action on a Shallow Sea. Tellus 9:47-52.
Welander, P. 1961. Numerical Prediction of Storm Surges. Advances in
Geophysics, Vol. 8. Academic Press, New York.
Welander, P. 1968. Wind Driven Circulation in One and Two Layer Oceans
of Variable Depth. Tellus 20:1.
Wilson, B.W. 1960. Note on Surface Wind Stress Over Water on Low and
High Speeds. J. Geophys. Res. 65:3377-3382.
Witten, A.J. and J.H. Thomas. 1975. Calculation of Steady Wind-Driven
Currents in Lake Ontario with a Spatially Variable Eddy Viscosity.
Abstracts. Proceedings 19th Conference on Great Lakes Research, Inter-
national Association for Great Lakes Research.
Wurtell, M.G., J. Paegle and A. Sielecki. 1971. The Use of Open
Boundary Conditions with the Storm-Surge Equations. Mon. Weather Rev.
99:537-544.
137
-------
SECTION X
PUBLICATIONS
1, Gedney, R. and W. Lick. 1970. Numerical Calculations, of the Steady-
State, Wind Driven Currents in Lake Erie. Proceedings 13th Conference
on Great Lakes Reserach. International Association for Great Lakes
Research, p. 829-838.
2. Gedney, R.T. and W. Lick. 1971. Numerical Calculations of the Wind-
Driven Currents in Lake Erie. Case Western Reserve University,
Cleveland.
3. Gedney, R.T. and W. Lick. 1972. Wind-Driven Currents in Lake Erie.
J. Geophys. Res. 77:2714-2723.
4. Gedney, R.T., W. Lick, and F.B. Molls. 1972. Effect of Eddy
Diffusivity on Wind-Driven Currents in a Two-Layer Stratified Lake.
National Aeronautics and Space Administration, TN-D 6841.
5. Gedney, R.T., W. Lick and F.B. Molls. 1973a. Effect of Bottom
Topography, Eddy Diffusivity, and Wind Variation on Circulation in
a Two-Layer Stratified Lake. National Aeronautics and Space Adminis-
tration. TND-7235.
6. Gedney, R.T., W. Lick and F.B. Molls. 1973b. A Simplified Strati-
fied Lake Model for Determining Effects of Wind Variation and Eddy
Diffusivity. Proceedings 16th Conference on Great Lakes Research.
International Association for Great Lakes Research, p. 710-722.
7. Haq, A, and W. Lick. 1975. On the Time-Dependent Flow in a Lake,
J. Geophys. Res. 80:431-437.
8. Haq, A., W. Lick, and Y.P, Sheng. 1975. The Time-Dependent Flow
in Large Lakes with Application to Lake Erie, Case Western Reserve
University, Cleveland.
9, Kuhlman, J.K. 1974. Laboratory Modeling of Surface Thermal Plumes.
Ph.D. Thesis, Case Western Reserve University, Cleveland.
138
-------
10, Lick, W. 1976, Numerical Modeling of Lake Currents. Ann. Key.
Earth and Planetary Sci. Vol. 4. To be published.
11. Lick, W., J. Paul, and Y.P. Sheng. 1975, The Dispersion of
Contaminants in the Near-Shore Region, In: Mathematical Modeling
of Biochemical Processes in Aquatic Ecosystems (R. Canale. Editor).
Ann Arbor Science, Ann Arbor.
12. Paul, J.F, and J. Prahl. 1971. Investigations of a Constant
Temperature Rectangular Jet. Case Western Reserve University,
Cleveland.
13. Paul, J.F. and W. Lick. 1973a. A Numerical Model for a Three-
Dimensional, Variable-Density Jet. Proceedings 16th Conference
on Great Lakes Research. International Association for Great
Lakes Research, p. 818-830.
14. Paul, J.F. and W. Lick. 1973b. A Numerical Model for a Three-
Dimensional, Variable-Density Jet. Case Western Reserve University,
Cleveland. Technical Report.
15. Paul, J.F. and W. Lick. 1974. A Numerical Model for Thermal
Plumes and River Discharges. Proceedings 17th Conference on
Great Lakes Research. International Association for Great Lakes
Research, p. 445-455.
16. Paul, J.F. 1975. A Note on the Vertically-Averaged Equations of
Motion. Case Western Reserve University, Cleveland.
17. Sheng, Y.P. and W. Lick. 1973. Wind-Driven Currents in a Partially
Ice Covered Lake. Proceedings 16th Conference on Great Lakes
Research. International Association for Great Lakes Research.
p. 1001-1008.
18. Sheng, Y.P. and W. Lick. 1975. The Wind-Driven Currents and Con-
taminant Dispersion in the Near-Shore of Large Lakes. Case Western
Reserve University, Cleveland.
139
-------
TECHNICAL REPORT DATA
(Ptt'asc read Imlnictions on the reverse bcfo'e completing)
ni COOT ko'.
__EPA-600/3->6-_C)2g_
4. TITLt AND SUBTITLE
NUMERICAL MODELS OF LAKE CURRENTS
5, REPORT DATE
6. PERFORMING ORGANIZATION CODE
3. RECIPIENT'S ACCESSiOWNO.
April 1976 (Issuing Date)
7. AUTHOR(S)
Wilbert Lick
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORG \NIZATION NAME AND ADDRESS
Department of Earth Sciences
Case Western Reserve University
2040 Adelbert Road
Cleveland, Ohio 44106
10. PROGRAM ELEMENT NO.
1BA026
11. CONTRACT/GRANT NO.
R802359
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-ORD
15. SUPPLEMENTARY NOTES
16. ABSTRACT
As part of a research effort sponsored by the U.S. Environmental Protection Agency
to study the dispersion of contaminants in near-shore areas of large lakes, we have
developed numerical models which are capable of realistically describing the currents
throughout large lakes and, in particular, in the near-shore regions of these lakes.
This report summarizes our work to date on these, hydrodynamic models.
In our work, the emphasis has been on the development and use of three-dimensional
models. Three basic models and applications of these models are presented here.
The models are: (1) a steady-state, constant-density model; (2) a time-dependent,
constant-density model; and (3) a time-dependent, variable-density model. Each
model has its own limitations and has certain advantages over the others. Applica-
tions of each model, especially to flows in near-shore regions, have been made and
are discussed. Vertically averaged models have also been used by us, usually in
parametric studies, and a brief summary of these models is also given.
A list of all publications by us resulting from or pertinent to this project is
given in Section X of this report.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
Hydrodynamics
Circulation
Mathematical Models
Lakes
b.IDENTIFIERS/OPEN ENDED TERMS
Wind Driven Circulation
Models
Thermal Plume Models
Lake Erie
Lake Michigan
c. COSATI I icId/Guup
08H
19. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (This Report)
Unclassified
21. NO. OF PAGES
152
20 SECURITY CLASS (This pa^ej
Unclassified
22. PRICE
EPA Form 2220-1 (9-73)
140
US GQ'/ERNMENT PRINTING OFF(Ct 1976— f>
-------