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

-------