United States
Environmental Protection
Agency
Environmental Research
Laboratory
Corvallis OR 97330
EPA-600/3-78-073
July 1978
Research and Development
&EPA
Enhanced Hydrodynamical
Numerical Model for
Near-Shore Processes
-------
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 nme series are:
1 Environmental Health Effects Research
2. Environmental Protection Technology
3 Ecological Research
4 Environmental Monitoring
5 Socioeconormc 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-78-073
July 1978
ENHANCED HYDRODYNAMICAL-NUMERICAL MODEL
FOR NEAR-SHORE PROCESSES
by
R. A. Bauer
and
A. D. Stroud
Compass Systems, Inc.
San Diego, California 92109
Contract Number 68-03-2225
Project Officer
R. J. Callaway
Marine and Freshwater Ecology Branch
Corvallis Environmental Research Laboratory
Corvallis, Oregon 97330
Corvallis Environmental Research Laboratory
Office of Research and Development
U. S. Environmental Protection Agency
Corvallis, Oregon 97330
-------
DISCLAIMER
This report has been reviewed by the Corvallis Environmental Research
Laboratory, U. S. Environmental Protection Agency, and approved for publica-
tion. 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.
11
-------
FOREWORD
Effective regulatory and enforcement actions by the Environmental Protection
Agency would be virtually impossible without sound scientific data on pollu-
tants and their impact on environmental stability and human health. Respon-
sibility for building this data base has been assigned to EPA's Office of
Research and Development and its 15 major field installations, one of which
is the Corvallis Environmental Research Laboratory (CERL).
The primary mission of the Corvallis Laboratory is research on the effects
of environmental pollutants on terrestrial, freshwater, and marine ecosystems;
the behavior, effects and control of pollutants in lake systems; and the de-
velopment of predictive models on the movement of pollutants in the biosphere.
A. F. Bartsch
Director, CERL
111
-------
PREFACE
The increasing concern with local environmental impacts created by
various discharges into coastal waters has intensified the requirement for
information regarding near-shore currents and exchange processes. A great
deal of insight into near-shore dynamics can be gained through the use of
sophisticated numerical tools. Simulation of processes in a numerically
analogous regime can preclude learning about coastal current dynamics as the
drama of a Santa Barbara oil spill engulfs beaches and birds. In the long
run there is no substitute for a well-planned comprehensive direct measure-
ment program, but planning and executing such a measurement program without
the insights gained from numerical tools can prove very costly indeed.
Here the primary goal has been to extend the applications range by en-
hancing the capabilities of an existing numerical tool, the Hydrodynamical-
Numerical (HN) model, developed by W. Hansen and T. Laevastu. As the em-
phasis of much earlier work was not aimed at near-shore dynamics, certain
enhancements seemed in order. A secondary goal has been to more fully
document the optimized HN model version, developed by R. Bauer, which pro-
vides the basis of the effort.
This report provides a more comprehensive documentation of the HN model
and documents the experimental enhancements. A follow-on Users Guide will
supplement this report in several important respects, namely, logic and
coding details, implementation guidance and operational guidance.
IV
-------
ABSTRACT
An optimized version of a multilayer Hansen type Hydrodynamical-
Numerical (HN) model is presented and discussed here as the basis for the
following experimental extensions and enhancements developed to more ap-
propriately handle near-shore processes:
•Non-linear term extension to facilitate small-mesh studies of near-
shore, including river inflow dynamics;
•Layer disappearance extension to enable appropriate procedures in
tidal flat and marshy regions, as well as some down/upwelling cases;
•Thermal advection enhancement for treatment of thermal pollution
cases by method of moments coupled with heat budget procedures for
dynamic plume development experiments;
•Monte Carlo diffusion enhancement to deal with dispersion via sta-
tistical methods and comparison to the method of moments experiments.
Extensive efforts were invested in determining reasonable and appropriate
boundary conditions for both the basic model and the extended versions
presented here.
This report was submitted in fulfillment of contract No. 68-03-2225 by
Compass Systems, Inc. under the sponsorship of the U. S. Environmental Pro-
tection Agency. This work covers the period from June 1975 to March 1977
and was completed March 31, 1977.
-------
CONTENTS
Foreword
Preface i-v
Abstract v
Figures viii
Tables i-x
Abbreviations and Symbols x
Acknowledgment xiii
1. Introduction 1
2. Conclusions and Recommendations 2
3. Numerical Methods and Techniques 4
4. Experimental Extensions of the HN Model 20
5. Results and Discussion 31
References 46
Appendix
A. Verification of the Model 49
VII
-------
FIGURES
Number Page
1 Scheme of the numerical grid 6
2 Orientation of the numerical grid 7
* — *
3 U, £ (other than C), and V computational stencils at closed
boundaries 12
4 U, V computational stencils at open boundaries 12
5 £ and £ computational stencils at open boundaries 13
* *
6 U, V computational stencils at open boundaries 13
7 U, V computational stencils at open boundaries 14
8 £\ tJ, V computational stencils at closed boundaries 15
9 Fractional z index notation 20
10 U, V non-linear computational stencils at closed boundaries . . 22
11 U, V non-linear computational stencils at open boundaries ... 23
12 The (n,m) cell 24
13 h depth profile for coastal grid 31
14 Tidal flat grid 32
15 Tidal flat grid bathymetry (depths in cm) 33
16 Closure configuration TFC: t=17700 sec 40
17 Closure configuration TFC: t=18090 sec 40
18 Closure configuration TFC: t=18240 sec 41
19 Closure configuration TFC: t=18390 sec 41
20 Closure analysis for TFC cell (6, 13) 42
Vlll
-------
Number Page
21 Post-reopening analysis for TFC cell (6, 13) 42
22 Maximum seaward advance of tidal flat F234 43
23 TFC continuous source "plumes" 45
A-l San Onofre outfall area grid bathymetry 50
A-2 San Onofre outfall tidal and wind prescription August 2, 1972 . 51
A-3 Infrared flight 0958-1107 compared to values after 15 source
hours 53
A-4 Infrared flight 1226-1344 compared to values after 17 source
hours 53
A-5 Infrared flight 1503-1617 compared to values after 19 source
hours 54
TABLES
Number Page
1 Boundary Condition Confusion 11
2 Typical Tidal Flat Grid and Run Specifications 34
3 Flat F3 Formation; Non-linear vs Linear 36
4 Flat F4 Formation; Non-linear vs Linear 37
5 Flat Fl Formation; Non-linear vs Linear 38
6 Flat F2 Formation; Non-linear vs Linear 39
IX
-------
LIST OF ABBREVIATIONS AND SYMBOLS
ABBREVIATIONS
F1-F4
F1234
F234
F34
HN
I
m,M
MM
MP
n,N
NE
NM •
NMMP
NP
NPMM
TFC
—bathymetric features identified by notes 1-4 on Figure 14;
also tidal flat segments associated with bathymetric features
—tidal flat associated with combined features F1-F4
--tidal flat associated with combined features F2-F4
—tidal flat associated with combined features F3 and F4
—Hydrodynamical-Numerical
—the array index (n,m)=(m-l)NE+n
--the column index (m=l,2,...,ME-1,ME)
—the array index (n,m-l)=I-NE(m minus)
--the array index (n,m+l)=I+NE(m plus)
—the row index (n=l,2,...,NE-1,NE)
--number of rows for array variables
—the array index (n-l,m)=I-l(n minus)
--the array index (n-l,m+l)=l+NE-l(n minus, m plus)
—the array index (n+l,m)=1+1(n plus)
—the array index (n+l,m-l)=I-NE+l(n minus, m plus)
—Tidal Flat Convective (nonlinear tidal flat case)
SYMBOLS
A
An
Aw
c
c
C
d
Ea
Ew
f
F
Fx,Fy
9
h
H
H
i
j
k
K
-area variable
-noon altitude of the sun in degrees
-ambient water temperature
-cloudiness factor in tenths of the sky
-mean specific heat
-conservative quantity (concentration or volume)
-depth (thickness) of the heated fluid
-moist vapor pressure of air
-saturation vapor pressure of air
-coriolis parameter
-center of concentration (one dimensional case)
-x,y coordinates of center of conservative quantity
-acceleration of gravity
-initial layer thickness
-dynamic layer thickness
-mean thickness at z point
-layer indicator (i=l,k)
-general index
-number of layers in the HN model
-new heat quantity (added heat)
-------
M —zeroth order moment of conservative quantity
M —first order moment of conservative quantity
M —second order moment of conservative quantity
N —a unit vector normal to the indicated boundary
p —stochastic variate
P —population of stochastic variates
qL —effect of external sources on layer surface deviation
Q —net heat flux
Qb —net atmospheric-oceanic (long wave ) radiative heat flux
Qe —evaporative/condensative heat flux
Qh —net atmospheric-oceanic sensible heat flux (convective trans-
fer or conduction)
Qm —heat flux resulting from human efforts
Qr —reflected (short wave ) radiative heat flux
Qs —total incident
O —total clear sky incident (short wave) radiative heat flux
*os
r —friction coefficient
R —width of a conservative quantity or a random value between zero
and one
Rx,Ry —x,y widths of a conservative quantity
s —distance variable
S —stochastic coefficient
t —time variable
t, —length of a day (sunrise to sunset in minutes)
d
t —initial time (start of a computation)
T —temperature variable
T —unit vector in the tangential direction
Ta —surface air temperature
Tw —mean surface water temperature for a cell
u,v —indicators of points of computation for U and V, respectively
u[P],v[P] —u,v turbulence component assigned the value of p
U,V —vertically integrated velocity components for the x,y direc-
tion, respectively
IJ,\7 —the "smoothed" U and V components, respectively (result of
^ ^ eddy viscosity)
U,V —linearly interpolated u,v velocity components at the u,v
points, respectively
U —relative humidity (%)
W -.-wind speed (m/sec)
x,y —horizontal space coordinate variables
X,y —horizontal components in the x,y direction, respectively, of
external forces including the air-sea interface
z —indicator of points of computation for £
XI
-------
a — smoothing parameter
a --the mean altitude of the sun for a short period
6 — directional parameter
e --directional parameter
£ — half grid length
T — half time step
p — density of fluid (assumed uniform through layer)
C — deviation of the top of the layer from an arbitrary level
_ surface
£ — smoothed value of C (result of eddy diffusivity)
v — vector velocity (resultant of U and V)
— vector "dot" product
"•"
" "
— superscript dot represents the operator 9/8t
I — interpreted as "evaluated on"
9N — an element of the normal to the indicated boundary
( )x — represents the operator 9/9x
( )y — represents the operator 9/9y
2 2222
V — represents the operator 9 /9x +9 /9y
A — the difference analog to 3
- — "approximately equal to" _ _
V — represents the vector operator i9/9x+j9/9y (where i and j
represent unit normal vectors in the x,y direction, respec-
tively)
XII
-------
ACKNOWLEDGMENTS
We wish to thank the Naval Environmental Prediction Research Facility
for allowing use of their facilities; S. Larson, K. Rabe and J. Harding for
their many constructive suggesions; and the NPERF computer support per-
sonnel for their assistance.
We wish to thank Dr. T. Laevastu for his guidance and to acknowledge
all of his innovative enhancements to the original Hansen model, which form-
ed the basis for the optimized model.
We wish to express our appreciation to R. Callaway, Project Officer,
for providing the proper environment to perform theoretical studies of the
model and wish to note the thoughtful user oriented feedback from C. Koblin-
ski, also of the Corvallis laboratory.
Finally, we wish to express our thanks to S. Raines for preparing the
illustrations and tables and to D. Bauer for typing the report.
Xlll
-------
SECTION 1
INTRODUCTION
The Hansen type optimized multilayer Hydrodynamical-Numerical (HN)
model described by Bauer (1) is a practical tool which has been used
successfully to study the dynamics of numerous coastal areas. The optim-
ized version of the HN model combines the vertically integrated single
layer HN model originally developed by Professor W. Hansen (2), (3), Uni-
versity of Hamburg, Germany, and the multilayer multiple-open boundary HN
model proposed by Hansen and developed by Dr. T. Laevastu (see Laevastu,
et al., (4)-(9)).
The thrust of much prior development of the HN model was toward large
coastal areas, semi-enclosed seas and even open ocean areas. In order to
adequately focus the model on near-shore problems certain enhancements were
required. First, non-linear terms of the dynamical equations not essential
for much earlier work are required for near-shore scales of motion. Second,
layer disappearance problems, tidal flats and marshes, for example,' handled
inadequately by prior versions of the model, required proper treatment in
near-shore cases. Additionally, evaluation and extension to thermal advec-
tion problems, as well as Monte Carlo dispersion procedures, were deemed
essential to show the utility of the numerical model procedure in practical
applications.
-------
SECTION 2
CONCLUSIONS AND RECOMMENDATIONS
The non-linear (convective) term modifications have been successfully
implemented in the single layer optimized version of the HN model and should
be incorporated whenever an HN application is expected to have significant
current gradients. Experience indicates that the non-linear terms do not
play a significant role when the grid length in the difference scheme exceeds
one kilometer. When the grid length is less than one kilometer, however,
current gradients tend to make the non-linear terms significant. The non-
linear modifications as developed were designed for multiple layer applica-
tions (from one to three layers in the current version) and should present no
difficulties when applied to multilayer cases. No testing has been per-
formed for any multilayer case, however.
The "tidal flat" modifications have been successfully implemented in the
surface layer and are practicable for inclusion in any HN application where
layer disappearance phenomena are an essential feature of the dynamics. Use
of these modifications would generally imply use of the non-linear terms as
well. The multilayer case was designed to handle multilayer disappearance
problems, e.g., down/upwelling through an initially extant layer. Again, no
extensive testing of the multilayer feature has been performed and some im-
portant upwelling cases are not handled adequately, e.g. extension of a lower
layer horizontally into a region where the layer was not initially extant.
These latter cases were not covered due to basic model design change require-
ments .
Coupling of the heat budget formulation with the conservative duplex ad-
vection of "added heat" and volume of heated fluid has been successful within
the assumptions made. Verification with a real case would be desirable to
validate the procedural verification. The assumption of advection without
diffusion is somewhat restrictive, thus experiments toward inclusion of a
controlled diffusive scheme would be desirable. The method solves the prac-
tical problem in its present form, but could be improved as to display of
temperature anomaly fields.
The Monte Carlo technique compares quite favorably to the conservative
advection by method of moments in terms of quantitative results. For test
cases, the cost of the Monte Carlo procedure was roughly an order of mag-
nitude less than the method of moments for approximately comparable quanti-
tative resolution with the same temporal resolution. The cost of any Monte
Carlo method is proportional to the statistical sample size, where the cost
of the method of moments is proportional to the number of affected grid
cells.
The present work included only procedural verification of the extensions
and enhancements. No attempt has been made to simulate any particular real
-------
case for comparison with real data sources. Verification of this latter
type should be accomplished prior to any attempt at predicting a real case
using these experimental procedures.
Various subtle aspects of the HN model and its associated boundary con-
ditions should be investigated further. For example, a tangential free
slip condition may enhance computations near unprescribed open boundaries
by reducing the transverse frictional effects currently associated with
these artificial grid bounds.
-------
SECTION 3
NUMERICAL METHODS AND TECHNIQUES
THE HYDRODYNAMICAL-NUMERICAL MODEL
The Hydrodynamical-Numerical (HN) model is an explicit numerical dif-
ference scheme based on leap-frog integration of the two dimensional Eu-
lerian form of the hydrodynamical equations through time to obtain a dynam-
ical boundary-value solution of tidal order.
The Hydrodynamical Equations
The following differential form of the hydrodynamic equations was de-
rived by integration of both x and y velocity components through a layer of
assumed uniform density to achieve vertically averaged mean components.
Terms to account for rotation of the earth, wind and tidal forcing, and dis-
sipative bottom (and upper layer) frictional effects have been included. In
multilayer cases the dissipative frictional effects are grossly represented
in upper layers as 50% of that which would have been computed for a bottom
layer of the same velocity/depth configuration. A more physically realistic
interlayer friction term would consider velocity in adjacent layers as stress
transferring momentum between layers as well as slightly dissipating mechan-
ical energy (to heat) within each layer. Eddy viscosity/diffusivity terms
have been omitted, as their use is more intimately related to numerical
schemes than to any imperative of physics. (See section on numerical smooth-
ing below.) Leendertse (10) provides a relatively detailed development and
analysis of the essential form of the differential equations presented here.
The form is analogous to that presented by Hansen (2), (3) and is substan-
tially the same as that presented by Laevastu (7), (11) and (14).
These equations for layer i of k layer system (where i=l in the surface
layer) are as follows:
Continuity;
x Momentum:
U. «J. (U.) x+V. (U.) y-fV.+r. (uX> * VV
(2)
-------
y Momentum:
, 2 2%.5-1
(3)
where:
JT —deviation of the layer surface from an arbitrary level,
U,V —horizontal components in the x,y direction, respectively, of
the vertically integrated velocity,
H —thickness of the layer,
x,y —horizontal space coordinate variables,
X,Y —horizontal components in the x,y direction, respectively, of
the external forces at the air-sea interface,
f —coriolis parameter,
r —friction coefficient,
g —acceleration of gravity,
p —density of the medium,
q —effect of external sources such as river inflow and rainfall
on layer surface deviation;
the following notational conventions apply:
""" —represents the operator 3/3t (superscript dot),
( )x —represents the operator 3/3x,
( )y —represents the operator 3/3y;
and the following conditions at upper or lower fluid boundaries apply:
t.. =0, when layer i intersects bottom topography,
{?i_1)x=(Ci_1)y=Pi_1=0, when i=l.
Numerical Difference Scheme
Computational Grid—
The finite difference equations presented in the following sub-section
are based on the staggered grid shown in Figure 1. The z grid points (inter-
sections of grid net) identify the geographic positions where t. values are
computed. The u grid points are one-half a grid length to the right of the
correspondingly indexed z grid points. The v grid points are one-half a grid
length below correspondingly indexed z grid points. Computation of U and
V velocity components occur at u and v grid points, respectively. Positive
U and V velocity components are directed to the right and upward, respective-
ly, in Figure 1. The relationship between computational coordinates and geo-
graphic coordinates is established by Figure 2. All land or sea features
must contain at least one z grid point. In particular, at least one z grid
point must separate an open boundary from a parallel closed boundary. The
last v row and the last u column become superfluous in the computational
scheme but are required for proper initialization.
-------
N=l
M=l 2 3
,z u z u - - - z • - -
.z • • • u • • • z • • • u
1
ME
\— . . .z . . . u
:\ t
v«*-(2,2)
2---U--- ,
- 1
u
SEA
Z-..U...Z...U-..2...U
\— . -z . • • u
I
I
I
V
I
I
•+U
NE
z u z u z u
•"A
z - r, computational point
u - U computational point
v - V computational point
land boundary
open boundary
z grid structure
Figure 1. Scheme of the numerical grid.
-------
Angle of Rotation
Counter-Clockwise
from North
E,90
+U,0
-V.270
ST180
Figure 2. Orientation of the numerical grid.
The grid point indexing convention used in Hansen models is for the z
grid intersection, the u point to the right and the v point below, to have
the same double subscripts with the row index first. The row index is
n:l
-------
Vt+2T(I) =
i
- . 5 ( T/Jl) { [Vfc (NM) -Vfc (NP) ] V^ ( I) + [V* (MP) -V* (MM) ] O± ( I ) }
where •-
U(I) = aU( I) + .25(l-ct) [U(NM)+U(NP) +U (MM) +U(MP)] ,
V(I) = aV(I)+.25(l-a) [V(NM)+V(NP) +V(MM) +V(MP) ] ,
UD = a^(I)+.25(l-a) [C(NM)+C(NP)+C(MM)+^(MP)1; (7)
= .25[U(I)+U(NP)+U(MM)+U(NPMM)] ,
= .25[V(I)+V(NM)+V(MP)+V(NMMP)] r <8)
Hvt+T(I) = hv.+.5{[^t+T(I)+?t+T(NP)]-[(I)+(NP)]}; (9)
where:
h —initial layer thickness,
t —time variable,
u,v —indicators of U and V computational points, respectively,
«, —half grid step,
o —smoothing parameter,
T —half time step.
Numerical Smoothing—
Explicit solution of non-linear systems requires systematic suppression
of the shortest waves through smoothing or other techniques to avoid non-
linear instability. Equations 7 (the "bar" terms) introduced numerical
smoothing into the difference solution. The differential form of these
terms was not made explicit, since their introduction into the difference
form is a numerical device, rather than a physically necessary term.
Smoothing techniques have been applied to hydrodynamic models to sup-
press higher order harmonic waves generally associated with non-linear
-------
instability. Lamb (13) applied a successive approximation analysis to the
simplified non-linear differential system:
U+U(U)x = -gU)x,
£ = -[(h+C)U]x. (10)
He demonstrated the appearance of higher order harmonics in the analytic
solution of this non-linear system. Successively higher order harmonics are
introduced at each step, limited in the difference solution by the minimum
representable wavelength. Energy is thus transferred from longer to shorter
waves and is accumulated in the shortest possible waves. Leendertse (10)
showed that the transferred energy becomes fixed in space due to the zero
velocity associated with the shortest waves.
The bar terms imply an "artificial viscosity or diffusivity" term in the
differential form. The artificial term in each momentum equation has been
interpreted physically as "horizontal or lateral eddy viscosity or turbulent
friction." These terms would appear on the right-hand side of Equations 2
22 2
and 3 as +AiV Ui and +AiV Vi, respectively, where V is the two dimensional
2 222
Laplacian operator: 9 / 9x +9 /9y . The form of the analytic solution accord-
ing to Lamb (13) shows the requirement for a corresponding "lateral eddy
diffusivity" term which would appear on the right-hand side of Equation 1 as
2
AiV £i. The coefficients of these "artificial viscosity and diffusivity"
2
terms become in the difference form (l-ai)£ /2, where ai is an arbitrary
smoothing parameter, usually chosen near .99. These artificial terms have
been retained as a tuning factor in many quasi-linear HN applications, e.g.,
Hansen (3), Laevastu (4)-(9), (11), (12) and (14). The smoothing parameter
should most properly be a function of the horizontal eddy viscosity/diffusiv-
ity coefficient and the grid/time step to avoid arbitrarily heavy smoothing
with resultant misleading dynamic impacts.
An interesting sidelight to the smoothing problem is the smoothing in-
herent in the scheme through linear interpolation. Particularly with ref-
* *
erence to the star terms (U,V), Kagan (15) states that the HN system "will
lead to extremely strong smoothing in the calculated velocities and eleva-
tions." His analysis reveals that for a grid mesh length in excess of about
30 Km the "computational viscosity" coefficient would exceed that expected
for horizontal turbulence in tidal motion. This should be a topic of
further study; however, it is doubtful that the "computational viscosity"
inherent in the star term would replace the smoothing required for non-linear
model stability.
Boundary Conditions
Differential Form—
The area of interest is defined by an enclosing boundary (T), which
includes open and closed boundary segments. The boundary conditions which
are applied at the land-sea (closed) boundary (T ) and the sea-sea (open)
boundary (F ) in a differential form are:
Condition 1: v-NJ = 0, (11)
-------
Condition 2: c|r2 =?(x,y,t), (12)
Condition 3: 3(vN)/3NJr2 = 0, (13)
Condition 4: 9H/9N|F2 = 0, (14)
Condition 5: 9?/9N|r = 0, (15)
Condition 6: 9(vT)/9N|r = 0, (16)
where:
v —vector velocity,
N —a unit vector normal to the indicated boundary,
T —a unit vector tangential to the indicated boundary,
"." —vector "dot" product,
| —may be interpreted as "evaluated on",
3N —an element of the normal to the indicated boundary.
(Note that v-T=U and vN=V along a boundary in the x direction and vice versa
along a boundary in the y direction.) Conditions 1 and 2 are Dirichlet condi-
tions, while conditions 3 through 6 are Neuman conditions. The combination
of conditions applied may be viewed as a Cauchy condition, which is precisely
the required condition for hyperbolic systems of equations.
The combination of conditions 1, 3 and 6 form a Cauchy condition on the
velocity variables. The combination of conditions 2 and 5 form a Cauchy con-
dition on the c variable. Condition 4 is required only in association with
condition 2 when ?(x,y,t) is computed dynamically. TABLE 1 shows the con-
ditions and their effect.
Condition 1 effects land-sea boundary closure by imposing zero normal
velocity at the boundary. Tidal harmonics, permanent flow characteristics,
river inflow characteristics, etc., are generally imposed via condition 2.
Condition 3 establishes a uniform velocity condition at the open boundary on-
ly. This differs from the work of Kagan (15), which suggested application
only at closed boundaries. Since only one condition may be applied at any
given boundary node, the no-flow condition. Equation 11, was chosen for all
closed normal boundaries.
When condition 2 is not applied via external information, conditions 4
and 5 complement condition 3 to form in concert a uniform flow condition at
open boundaries. Condition 5 at closed boundaries was selected to effect a
slight amplitude gradient damping, while avoiding a laminar or quasi-laminar
amplitude damping. Condition 6 implies a slip condition along both open and
closed tangential boundaries. In an extremely fine mesh model a laminar
(infinite lateral friction; zero slip) condition might be indicated. The
gross scale of HN applications indicates, however, that slight velocity
gradient damping, Equation 16, is sufficient at land-sea and open boundaries.
Difference Method of Applications—
The various methods of applying in the difference scheme conditions 1-6
specified in Equations 11-16 are presented in this section with figures for
clarification. Figures 3-8 depict computational stencils, which enclose
geographic positions for values required by the indicated computations.
10
-------
TABLE 1: BOUNDARY CONDITION CONFUSION
BOUNDARY
CONDITION
#
1
2
3
4
5
6
7
Ecj.
(11)
(12)
(13)
(14)
(15)
(16)
(33)
Fig.
3,8
10
4
5
5,6,7
11
5
5
8
8
7
11
10
BOUND
TYPE
closed
T!
open
r2
open
r2
open
F2
open
r?
closed
rl
do sec
TI
open
F2
open
r2
closed
ri
DIFFERENCE
ORDER OF
APPROX.
1st
2nd
1st
1st
2nd
1st
1st
1st
2nd
COMPUTATION
DIRECTLY
AFFECTED
c,u,v,u,v
non-linear
terms
Hu.Hv,
u,v,e
u,v,u,v,c
nonlinear
terms
C
?
U,V
nonlinear
terms
BOUNDARY CONDITION TYPE;
FRICTIONAL EFFECT
Solid wall , zero normal
flow; infinite normal
friction
External forcing
Internal dynamic
Reflective condition
with phase shift due to
light (1st order) or
moderate (2nd order)
friction effecting
normal to indicated
boundary
Reflective condition
with phase shift due to
light friction effect-
ing damping of gradient
normal to indicated
boundary
Perfect transverse re-
flection, total slip,
zero friction.
DIFFERENCE METHOD OF EFFECT-
ING BOUNDARY CONDITION
Zero normal velocity compon-
ent effects land-sea closure
Tidal harmonic/permanent
flow/river inflow/ "overt ide"
specification
Values not specified comput-
ed dynamically under Condi-
tions 3,4,5
Index modification effects
value assignment such that
same computational form may
be used for both internal
and boundary computations
Direct value modification al-
lows the same computational
form to be used for both in-
ternal and boundary computa-
tions
A combination of value and
index modification allows the
same form to be used for both
internal and boundary compu-
tations
Extrapolation of internal
transverse gradient (doubl-
ing of first order half dif-
ference computed due to in-
dex modification.)
-------
Superscript numerals in these figures identify the numbered boundary con-
ditions that are applied. Arrows where shown in Figures 3-8 indicate ef-
fective value transfer used to impose the boundary conditions.
M-l
M+2
N+l
* — *
Figure 3. U, T, (other than t) , and V
computational stencils at closed boundaries.
N=l
NE=4
M=l
'—u—
U
ME=4
2 ' J_2 ! 12
— —u— f z u- —z) • u
I
V
I
J2
u -Cz • u
V
Cz • u • z)• u
[r
[
V
'2
Cz u z>—u—
u
I
V
u
Figure 4. U, V computational stencils at open boundaries.
12
-------
NE=3
z • u •
V
-
-z u
V
1 '
z
1
1
V
1
f
z
V
I
• u
5
Y u ^"J
5
Figure 5. C and t. computational stencils at open boundaries.
Figure 6. U, V computational stencils at open boundaries.
13
-------
z—
U
v)
z •
u
.V /
-z
•u}--
V
V
z
1
Y
V
1
V
• u
1
— ^
V
- u
3
Figure 7. U, V computational stencils at open boundaries.
14
-------
N-2
N-l
N
N+l
N+2
— •£. • X
v
- Z • X
1— v
r— Z • X
•
i~(z • v
— v
|- Z • X
- V
1
>s
1 •
\
\
\
s (
\
\
\
1 •
1
1 •
z
z
V
.
z
I
z
V
z
N — '
V
• X
N
• X
•1
6?
<
b :^
• X
s
• X
"
- I
a •
a •
i •
i •
i •
z
z
f
V
z
z
z
V
• X
L
s\\\\\\\
• X
v
• X
) -
-------
Conditions 3-6 at open bounds—Equations 13-16 are satisfied in the dif-
ference scheme at open boundaries such that values representing a position:
(1) outside the defined grid mesh, (2) in the last row of V, hv or Hv values
or (3) in the last column of U, hu or Hu values, are replaced with the next
like variable value lying on an inward directed normal line from the inter-
vening boundary. (Refer to positions identified by superscript 3-6 in
Figures 5-7.)
Conditions 5 and 6 at closed bounds—Equations 15 and 16 are satisfied
in the difference scheme at closed boundaries such that values representing
positions over land (or at a normal boundary not intersected by the tangen-
tial flowline) are replaced with the next like variable value lying on the
seaward directed normal line from the intervening tangential boundary. (Re-
fer to positions identified by superscripts 5 and 6 in Figure 8.)
ADVECTION BY METHOD OF MOMENTS
The method of moments is a quasi-Lagrangian method which maintains in-
formation on the zeroth, first and second order moments of the concentration
in each cell of the grid mesh. The scheme preserves concentration gradients
in the original distribution through time based on externally supplied ad-
vection velocities while avoiding the pseudo-diffusion encountered in Eu-
lerian difference methods.
The method selected for use here is that of Pedersen and Prahm (16).
The differential equation of advection 3C/3t=-WC, holds where C is any
conservative quantity, as long as V is a non-divergent velocity field. HN
integrated velocities (vertically averaged) over a region of variable depth
do not precisely satisfy the non-divergence requirement. It has been as-
sumed that for the regions of interest these divergence effects will not be
significant. A single Lagrangian step is taken based on a velocity at the
point interpolated from the Eulerian grid. Immediately the parameters of in-
terest are decomposed from the results of the Lagrangian step back into the
Eulerian framework by summation of contributions from surrounding cells. The
parameters of interest for each cell are total conservative quantity C, cen-
ter of concentration F and width of the quantity R. The equations used for
the zeroth, first and second order moments in one dimension are:
+ .5
M = C = / C(x)dx,
+ .5
M = F = / C(x)xdx/C,
1 -.5
+ .5 2
M = R /12 = J C(x)
-------
HEAT BUDGET FORMULATION
The primary emphasis of this application of the Laevastu thermal model
is on heat budget rather than response of atmospheric parameters to ocean
surface temperature under the gradient wind, which has been discussed by
Laevastu, et al., (17), (18). A more detailed treatment of the heat budget
components, including some terms not considered below, may be found in Lae-
vastu (19). Young (20) provides quantitative values of the various heat
budget terms for the southern California shelf area, which have been used
as representative input for thermal effects.
Oceanic Heat Budget
A statement of the oceanic heat budget which is sufficient for the
purposes of this study follows:
Q = Qs-Qr-Qb-Qh-Qe+Qm, (20)
where:
Q —net heat flux,
Qs --total incident (short wave) radiative heat flux,
Qr —reflected (short wave) radiative heat flux,
Qb —net atmospheric-oceanic (long wave) radiative heat flux,
Qh —sensible heat loss to (gain from) atmosphere (convective
transfer or conduction),
Qe —evaporative heat loss (latent heat uptake by atmosphere),
Qm —heat induced by human efforts.
The preceding formulation ignores the following heat exchange processes:
transfers through ocean bottom, kinetic transformation (tidal and wind fric-
tion) , natural chemical processes, precipitation, fresh water runoff and net
circulation through region boundaries. Man-induced heat exchange processes
including power plant cooling, sewage and industrial outfalls, man-induced
chemical processes, etc., may be handled as pollutant sources.
Heat Budget Components
-2 -1
Each of the natural components of Q is presented in cal cm day units
(1 calorie = 4.18605 joules). The formulae of this section are presented
essentially following Laevastu (19). Total daily insolation is given by:
Qs = Qos(l.-.0006c3), (21)
where c is the cloudiness in tenths of the sky and Qos is the clear sky
cloudless short wave radiation. Daily values of Qos may be computed em-
pirically as:
Qos = .014A t , (22)
where A is the noon altitude of the sun (degrees) and t, is the length of
n d
the day (sunrise to sunset) in minutes. Shorter term computations modeling
the day/night cycle would require a more definitive method, such as:
Qos = 2736. sincT ' (23)
17
-------
where a is the mean altitude of the sun for the time period of interest.
The daily reflected radiation may be expressed in terms of Qs in a
linear relationship by:
Qr = .087Qs. (24)
Although presented for daily values, Equation 24 could be used with Equations
21 and 23 for a first guess approximation of short term values.
Effective back radiation, the net long wave radiation of the atmospheric-
oceanic system is computed as follows:
Qb = (297.-l.86Tw-.95U )(l-.0765c) (25)
where Tw is the water temperature, U is the relative humidity and c, as be-
fore, is the cloudiness factor. Both the net atmospheric long wave radiation
and the net back radiation from the ocean follow a fourth power law of their
respective temperatures. Equation 25 is interesting since it combines both
these long wave radiation terms into a single equation basically dependent
on the water temperature alone. See Laevastu (19) for a description of
Equation 25 development.
The evaporative heat loss (or gain) is computed as follows:
Qe = 53.9(.26+.077W)(Tw-Ta) (26)
where W is the wind speed in m/sec and Ta is the air temperature. A similar
formulation is employed for the sensible heat loss (gain):
Qh = 39.0(.26+.077W)(.98Ew-Ea) (27)
where Ew is the saturated water vapor pressure of air at the temperature of
the water surface and Ea is the water vapor pressure of air over the surface.
The man-induced component is determined externally based on the situa-
tion to be modeled.
MONTE CARLO SIMULATION
Monte Carlo techniques as a class simulate stochastic processes through
the judicious application of a random number generator. The stochastic
character of diffusive processes lends itself to simulation by Monte Carlo
methods. In order to introduce the random element, the total velocity for a
particular fluid particle is assumed to be composed of a mean flow velocity
component and a turbulent flux velocity component. The HN model provides
the mean flow velocity and the Monte Carlo scheme provides the turbulent
flux velocity. Dispersion is thus modeled by simulating the diffusion
process stochastically within the background fluid in motion. For a particu-
lar case, a sample of tracer particles representing the output of a source is
tracked in time. The number of particles in the sample is determined by the
particular case and the desired degree of statistical resolution.
Meier-Reimer (21) presented the results of several alternative ap-
proaches to applying Monte Carlo simulation to the diffusion problem and
found the Monte Carlo method superior to either Lagrangian or Eulerian dif-
ference methods in terms of results as well as cost. Using a meteorological
application Thompson (22) showed the Monte Carlo process to be simple and
direct with a good match between statistically expected and numerically
18
-------
experimental values for simple Fickian cases. Thompson also demonstrated
the ease with which the method could be generalized to more complex depen-
dency regimes and more sophisticated parametric representations.
19
-------
SECTION 4
EXPERIMENTAL EXTENSIONS OF THE HN MODEL
THE NON-LINEAR ENHANCEMENT
Most of the Hansen type HN models have neglected the convective (ad-
vective) terms since analysis indicated that these non-linear terms were not
significant for most practical applications. Further, retention of these
terms to cover the general case was not economically justifiable and com-
plicated the boundary problem. These non-linear terms are first order dif-
ferential terms which should be analytically expected to show significant
effects in fine mesh applications, especially where lateral flow specifica-
tion is required, e.g., river inflow. The incorporation of the convective
terms into the momentum equations was the first task in this effort to en-
hance the HN model for the study of near-shore processes.
Non-Linear Difference Form
Fractional z notation—
To simplify the description of the derivation of the difference form of
the convective terms, a decimal fractional z notation, illustrated in Figure
9, has been used. All points in the staggered grid have been identified by
fractional index values that have a common z grid origin. The staggered
grid point u(l,l) is represented as u(l,1.5) and staggered v(l,l) becomes
v(1.5,l). Thus fractional z notation is a positionally correct notation pro-
viding the flexibility to represent variables at any desired location.
m=l 1.5 2 2.5 3 3.5
n=l z •--u • •-z-•-u- • -z. • .u.
1.5 v o v o v o
2 z- • -u-•-z-•-u-• -z- •-u-
2.5 v o v o v o
3 z-•-u- •-z---u-•-z-•-u •
3.5 v o v o v o
z-z point, u-u point, v-v point
o - values required must be interpolated
Figure 9. Fractional z index notation.
20
-------
Derivation of the Difference Form —
The convective difference forms presented in Equations 5 and 6 are de-
rived. Boundary conditions are treated separately. This derivation treats
only those points not immediately adjacent to boundary points. Using frac-
tional z notation Equation 2 is solved at U(n,m+.5), and Equation 3 is solved
at v(n+.5,m) . The convective terms from Equations 2 and 3 have been included
here for convenience:
. . . . (28)
The difference form of Equations 28 was derived by replacing the par-
tial differentials with spatially centered differences as follows:
(Ui) (n,m+.5) - [Ui (n,m+l . ) -Ui (n,m)
(Ui) r(n,m+.5) - [Ui (n-. 5,m+.5) -Ui (n+.5,m+.5) ]/2£,
J.
(Vi) (n+.5,m) - [Vi (n+.5,m+.5) -Vi (n+.5,m-. 5) ]/2£,
X
(Vi) (n+.5,m) = [Vi (n,m) -Vi(n+l. ,m) ]/2£. (29)
Equations 29 constitute difference approximations of the respective differen-
tial factors to first order accuracy; however, u and v velocity components
are required at "z" and "o" points, as identified in Figure 9. As U and V
values will not be available at "z" and "o" points, linear interpolation be-
tween available u and v points yielded:
(Ui) (n,m+.5) -> [Ui (n,m+1.5) -Ui (n,m-.5) ]/4£,
x
(Ui) (n,m+.5) -> [Ui (n-1. ,m+.5) -Ui (n+1. ,m+.5) ]/4£,
y
(Vi) (n+.5,m) -> [Vi (n+.5,m+l .) -Vi (n+.5,m-l. ) ]/4£,
x
(Vi) (n+.5,m) -»• [Vi (n-.5,m)-Vi (n+1.5,m) ]/4£. (30)
Returning to staggered grid notation resulted in the following differ-
ence approximations :
(Ui) (I) -> [Ui(MP)-Ui(MM)]/4£,
X
(Ui) (I) ->• [Ui(NM)-Ui(NP) ]/4£,
Y
(Vi) (I) -»• [Vi(MP)-Vi(MM)]/4£,
X.
(Vi) (I) -*• [Vi(NM)-Vi(NP)]/4£-. (31)
Y
Equations 31 constitute difference approximations of the respective differ-
ential factors to second order accuracy.
Substituting the Equations 31 into Equations 29, as required, transpos-
ing terms to the right-hand side of the equation and multiplying each equa-
tion by a factor of 2x, yielded the difference forms presented in Equations
21
-------
5 and 6- as follows:
. ..-.5(T/£){Ui(I) [Ui(MP)-Ui(MM)]+Vi(I) [Ui(NM)-Ui(NP)]}...,
...-.5(T/£){Ui(I)[Vi(MP)-Vi(MM)]+Vi(I)[Vi(NM)-Vi(NP)]}
(32)
* *
where the star terms, V and U, respectively, represent the linearly interpo-
lated v velocity component at the u point, and the linearly interpolated u
velocity component at the v point. (Note that Vi(n,m+.5) and Ui(n+.5,m) are
not directly available in the difference scheme.)
Non-Linear Boundary Conditions
The boundary conditions, Equations 13 and 16, are effected by index sub-
stitution and would be the default convective computational boundary condi-
tions (unless special measures were taken) due to the procedure of computing
the convective terms. Equation 11 is also a default condition, but is con-
sidered to be correctly so (Figure 10). When index substitution occurs to
effect conditions 3 and 6, the separation of the values differenced decreases
by half due to the index modification, but the form of the difference compu-
tation causes division by the second order difference width (4£). Linear ex-
trapolation of the internal gradient would then easily be accomplished by a
factor of 2. Leendertse (10) found that the two step implicit scheme showed
characteristics of local instability when a linear extrapolation condition
(zero friction) was applied at boundaries normal to the direction of flow.
He had previously shown a theoretical tendency in this regard for a simpli-
fied equation set.
M-l
M
M+l
M+2
N+l
N+2
z-u-z-u-z-u-z-u
V \^
1
Z • 1
7- 1
V\\\\\ TT
1 •
Cv
z • u •
v^J
.
z
•
V
z
V
KW^ v x\\\\\vv^^
I . r^i
• u • z •
1 -
v7)^
• Cu • z •
KWV^S
K:
z • u • z • u • z -
V
rt
u
u
u
A;
1
^
• z - u
V
•
1
1
1
• z • u)
^^v^sss^
7
1
• z - u
V ^ V V
• 1 •
Figure 10. u, V non-linear computational stencils at closed boundaries.
22
-------
N=l
2 ME=3
I I I 1
u
Figure 11. U, V non-linear computational stencils at open boundaries.
Instability with respect to the slightly stronger (friction ^ 0) condi-
tion, Equation 13, was found when the convective terms were applied to a
Prudhoe Bay model by EPA at Corvallis. Subsequent experiments with the
Prudhoe Bay model showed that when the first order form of Equation 13 was
replaced with the second order form of the same condition (precisely twice
the implied friction of the first order form) that the instability disap-
peared (see Figure 11). This is the same result that Leendertse (10) found
successful. The source of the initial instability in the Prudhoe Bay model
was apparently local, which, as Leendertse indicated, was fed back into the
equations through the non-linear terms to drive the entire area to instabil-
ity.
No apparent instability was introduced when the first order tangential
condition, Equation 16, was modified to a second order zero friction (free
slip) condition by extrapolation of the internal gradient. Equation 33 de-
defines the free slip Neuman condition used for the non-linear terms at all
boundaries tangential to flow in lieu of Equation 16.
Condition 7:
3[v(x,y)-T]/9N| = Iim8[v(x+es,y+<5s)'T]/9N.
s-K>
(33)
The parameters e and 6 are chosed as 0, ±1, irrespectively, so that a boun-
ary point (x,y) is approached from seaward along N as s approaches 0. The
sign of the non-zero parameter would reflect the appropriate direction.
Figures 10 and 11 depict computational stencils enclosing geographic
positions for values required for the non-linear term computations. Super-
script numerals in these figures identify the numbered boundary conditions
specified in Equations 11-16 and 33 that are applied. Figure 10 should be
interpreted that the condition is applied at land edges rather than across
narrow channels. The cross channel component of the non-linear term is
zero for a single cell channel. Arrows, where shown in these figures.
23
-------
indicate effective value transfer used to impose the boundary conditions.
TABLE 1 provides a comparative description of the various conditions.
An item of further study should be the use of this free slip boundary
condition for other than convective term computations at tangential open
boundaries as a means of avoiding even the slight lateral friction currently
imposed by Equation 16.
THE LAYER DISAPPEARANCE ENHANCEMENT
Many regions in the coastal zone are characterized by the cyclic emer-
gence and inundation of extensive mud flat or marshy areas due to shallow
bathymetry and tidal range. The optimized HN model assumed fixed land-sea
boundaries based on initial depth fields. The boundry conditions associated
with those fixed boundaries were violated when tidal flat conditions oc-
curred. The second task proposed to modify the HN model to handle variable
land-sea boundaries and their associated boundary conditions. To simulate
tidal flat regions the effective land-sea boundaries were made a function of
layer thickness. The variable land-sea boundary capability was designed to
be an optional, rather than an integral feature to have minimal effect on
the basic model. An hypothetical tidal flat bathymetry was designed to fa-
cilitate procedural verification testing which would thoroughly exercise the
dynamic boundary operation. Additional verification tests were performed on
the grid representing Prudhoe Bay, Alaska, which had been prepared by EPA
at Corvallis.
Layer Disappearance Procedures
Figure 12 depicts the (n,m) cell which is defined with the z(n,m)=z(I)
point in the top center of the cell with vertical boundaries passing through
the u(I) point to the east, the u(MM) point to the west, the v(NM) point to
the north and the v(I) point to the south, and horizontal boundaries at the
upper and lower surfaces of the layer. When fluid volumes are excluded or
included in the computation they comprise a single cell or group of such
cells.
M-l
M
M+l
U- • • Z- • • U-
N+l
_._._y__._..
-j V
i...
-V-
U
•U-
•U-
•U
N
• ••z grid structure
—(n,m) cell boundary
Figure 12. The (n,m) cell.
24
-------
Cell Closure Procedure—
Following the normal computation of layer surface deviation (?) and
thickness (Hu,Hv), the average cell thickness, H, is computed at the z
point and compared to a prescribed minimum thickness, Hmin. If H < Hmin,
the cell should have been considered closed during the computation of C, Hu
and Hv. A closure status is maintained for each cell and this is changed to
indicate full closure. For each cell so indicated, those associated veloci-
ty components not bordering on a land cell are reset to zero, and the status
of adjacent, non-land, open cells is modified to show partial closure.
If H > Hmin, it is possible that some thicknesses associated with this
cell may individually be less than or equal to Hmin. For this case, a par-
tial closure of the cell is effected by setting the associated velocity com-
ponent to zero and maintaining the closure status appropriately.
Re-evaluation of Layer Surface Deviation—
The velocity components as computed in the last computational step may
have been modified by cell closure procedures. The layer surface deviations
(£) logically adjacent to the modified components must now be recomputed
based on the effectively modified land-sea boundary. The most compelling
reasons to perform the recomputation are: 1) to maintain the land-sea boun-
dary conditions described for regular land-sea boundaries; and 2) to main-
tain the layer surface deviations associated with fully closed cells consis-
tent with respect to surrounding layers. Based on the cell closure status,
actions to be effected in the single layer case follow:
1) open: no action;
t+T
2) partial closure: recomputation of £ ;
3) full closure: reset £ to the value of £
In the multi-layer case the same closure status conditions are possible
independently in each layer. The nature of the continuity equation, howev-
er, is such that a lower layer may have an effect on an upper layer, notwith-
standing the upper layer's closure status. In the recomputation procedure,
if a lower layer C. requires:
1) recomputation (partial closure);
2) non-zero "adjustment" (from the i+1 layer);
3) resetting of £. to the value of £. (full closure);
4) a combination of 1) and 2); or
5) a combination of 2) and 3);
then ?. (i-1^0) requires an adjustment 1', . The adjustment of Z' is equal to
the recomputed or reset and/or adjusted t,. less the value of £. . Thus
the recomputation procedures must occur in precisely the same order that the
original £ computations took place, i.e., lower layer ^ terms are computed
first so that the Z'. may be computed for use in succeeding layer computation.
Based on cell closure status, actions to be effected in the multi-layer case
25
-------
follow:
t+T t+T
1) open: C± = e± +
2) partial closure: £,. = 5. (recomp) + Z' ;
3) full closure: C = £ +z'+1-
Boundary Condition Consideration
Velocity components which have been set to zero by the cell closure
procedures are treated as if they were a true land-sea (closed) boundary by
the recomputation of c, as it occurs. Second order computations in £ (£)
are discontinuous since £ is not included in the recomputation. (Figure 8
does not apply to recomputation) . Further, recomputed (or changed) £ values
do not effect any changes in £ values which are adjacent to the recomputed
(or changed) £ values. Tests calculating £ with a correctly recomputed £ did
not significantly alter or improve the results, ergo £ is not recomputed.
Second order computations in U and V (U, V and the convective computa-
tions) are also discontinuous, but only for the step in which a particular
cell side is first closed (set to zero.) These second order U and V compu-
tations adjacent to a cell side (for which it is determined for the first
time that the velocity component should have been zero valued) were com-
puted using a value other than zero. All such succeeding second order com-
putations adjacent to continuously closed cell sides are computed correctly
for all subsequent steps including the step in which the side is reopened.
These second order discontinuities in the various computations have not and
should not provide any problems, since the smoothing through the remainder
of the field in subsequent steps has been designed to remove such low order
perturbations .
THERMAL ADVECTION EXTENSION
The advection of a thermal discharge is a common practical problem of
interest in the coastal zone, particularly with respect to power plant
cooling water and, to a lesser degree, waste water discharges. If one ig-
nores the density effects of the heated fluid on the dynamics of a region,
the HN computed current velocities for the top layer may be reasonably used to
advect a thermal discharge through the region. The Pedersen-Prahm (16) advec-
tion scheme has been chosen since it is a conservative scheme without the
pseudo-diffusion of Eulerian difference methods. In order to model the heat
budget effects on the thermal discharge as it is transported by the currents
through the region, the Laevastu thermal techniques (19) were selected.
A version of the Pedersen-Prahm advection routine developed by J. Hard-
ing of Naval Environmental Prediction Research Facility (NEPRF) , Monterey, was
combined with the thermal techniques used by Laevastu, et al. , (17) , (18) .
Both FORTRAN computer program systems required conversion to the Lawrence
Berkeley Laboratory (LBL) CDC 7600 and required appropriate linkages to weld
26
-------
the two systems into a single practical system of thermal advection. The
current scope of the task does not admit the study of diffusive effects.
Thermal Equations
The change of temperature T of a volume C of water induced by the quan-
tity of added heat K in addition to "ambient heat" is given by:
AT = K/cpC, (34)
where c is the average specific heat of the fluid and p is the density, as-
sumed uniform through the layer. The quantity of heat introduced over an
area A during a time period 2x is provided by the heat budget formulation:
K = QA2x = (Qs-Qr-Qb-Qe-Qh+Qm)A2T; (35)
thus
AT = (Qs-Qr-Qb-Qe-Qh+Qm) 2-r/cpd (36)
where d is the depth of the heated fluid.
These equations assume that all added heat is absorbed, stored or lost
in the upper thickness d of the surface layer. This depth of the heated
fluid is a dynamic quantity which may be determined from the advected volume
being followed and the distribution information maintained by the moment
method of advection. The distribution information also allows computation
of the average temperature of the entire cell.
Duplex Thermal Advection
The advection procedure requires a conservative quantity. The volume
of the heated fluid is one obvious conservative quantity where the tempera-
ture is not conservative. Equation 34 shows that the volume is not suffi-
cient to determine AT. Both the volume and the "added heat" must be known
to determine the change in temperature. Since both quantities must be fol-
lowed for a thermal procedure, the process is being called "duplex thermal
advection." In the duplex procedure the full method of moments is applied
to the volume of the heated fluid and only the zeroth moment values are
computed for the "added heat." The first and second moment information com-
puted for the volume is assumed to hold for the "added heat" as well.
The time step required to maintain stability in the method of moments
is determined by the velocity of the fluid rather than the Courant-Friedrich-
Lewy criterion as in the HN model. To take advantage of the longer time step
permitted, the duplex thermal advection has been developed as a separate
program. This also makes it possible to perform a series of numerical ex-
periments based on a single set of HN velocity fields.
The duplex method is composed of a sequence of procedures performed at
each time step. The operations include thermal and advective processes, each
of which is computed independently. The advective process assumes the ther-
mal process is fixed in time and the thermal process assumes that the ad-
vective process is fixed in time for the computational interval.
Subscripts refer to thermal processes: superscripts refer to advective
processes. (If a sub- or superscript is missing, it is of the same time
27
-------
order as the sub- or superscript that is present.) The sequence of pro-
cedures follows:
1. t=t , continue at procedure #7; (initialization)
2. (Tw^i0_,Ta^OT,Aw^0J=thermal (Twt,Tat,Awfc) ; (thermal step)
"t+2l'
,
3.
4.
5.
t t t+2t _
(RxRy becomes Rx ,Ry
(U,V)
)
t;
(inverse interpolation)
(added heat computation)
(advection step)
(also adjust RxRy);
(full step accomplished)
(pollution volume increment)
(pollution heat increment)
6. t=t+2T;
7. Ct=Ct+AC
(temperature change computation)
(cell temperature interpolation)
11. when t=te stop, otherwise continue at procedure #2;
9. AT =Kt/(pcC );
10. Tw = T RxRy+Aw ;
t U U
where:
Rx,Ry
Ta
te
S
Tw
AT
Aw
thermal
advect
—horizontal width factors in grid step units related to second
moment of pollutant in cell (Equation 17),
—air temperature,
—end time,
—start time,
—surface water temperature,
—temperature in excess of ambient for the heat polluted volume,
—ambient water temperature,
—thermal procedure
—duplex advective procedure.
The device of starting at procedure 7 prevents advection or thermal ad-
justment from being applied to heated material added in the same step. Pro-
cedures 7 and 8 are determined externally by the case of interest. Note that
the width factors maintained by the advective procedure must be initialized
or adjusted in procedure 7. Procedure 9 computes the temperature anomaly
(or temperature in excess of ambient temperature) for the volume associated
with the excess heat. Procedure 10 interpolates from the temperature anomaly
obtaining the horizontal mean temperature for each cell. Procedure 11 ter-
minates the computation at the appropriate point or returns to procedure 2.
The ambient temperature is maintained constant for the experiments pre-
sented in the next section. Equation 23, however, provides the key to model-
ing a diurnal cycle.
28
-------
MONTE CARLO EXTENSION
Practical diffusion problems occur more frequently than anyone would
like. Controlled and uncontrolled petrochemical, chemical and sewage spills
are particularly of interest in coastal areas due to the potential ecological
effects, property damage and human environmental factors. The HN model
provides the detailed velocity information required to give a solution to the
basic vector path problem of following a particular spill. It has been shown
by Maier-Reimer (21) that when a Monte Carlo technique is appropriately ap-
plied to the basic HN velocity components associated with a particular con-
servative element of the spill, one may expect to model the diffusive process
much more reasonably than with other ostensibly more explicit procedures.
Pure Eulerian difference schemes generally produce an extremely unrealistic
level of artificial diffusion generally associated with the interpolation in-
herent in these schemes. Pure Lagrangian schemes break down due to Lagrang-
ian grid distortion after an unreasonably short time . The Monte Carlo scheme
is pseudo-Lagrangian in that specific pollutant particles are followed in
time, but the Eulerian framework is maintained to provide the velocities for
each Lagrangian step. Each marker particle is advanced one Lagrangian step
by application of a mean flow velocity interpolated to the present position
of the particle from the Eulerian grid values. The final location of the
marker particle at the end of the step is further determined by the applica-
tion of a random velocity increment consistent with the statistical view of
diffusion. The u and v components of the velocity increment are selected in-
dependently by a stochastic process. The stochastic process is constrained
to provide each selected component within a discrete range determined by the
diffusion process and the time step between available Eulerian velocities.
The Procedure
A Monte Carlo procedure may be defined for a set of Lagrangian tracer
particles in horizontal coordinates as follows :
t+2t t . t+2T t+2T . . . _
x. =x.+(U. +u. [P])-2i,
t+2t t , t+2T t+2T r . . „
y. =y.+(V. +v. [P])-2T,
for jth particle, where the x,y are position vector components, the U,V are
mean flow velocity components, and u[P],v[P] are variates stochastically
selected from P, the total population of possible turbulent flux velocity
components. In actual computations the tracer particles are followed in di-
mensionless grid coordinate (n,m) units as:
where the U[P],V[P] are selected from P independently in order to describe a
constant diffusion coefficient.
29
-------
The Stochastic Process
The particular stochastic process employed in this Monte Carlo scheme
is defined as follows:
p. = SCR.-.5),
where p. is the jth random variate "selected" from P, S is the stochastic
coefficient and R. is a pseudo-random number evenly distributed beween 0 and
1. Successive independent selection of the x and y components results in an
even distribution of the resultant turbulent flux velocity in a square of
side S centered at the tip of the mean flow velocity vector. According to
Meier-Reimer (21), the resulting distribution of tracer particles is statis-
tically isotropic after two time steps and Gaussian after about eight time
steps. Choosing S=1.0 results in a maximum effective "diffusion velocity" of
the order of .5 cm/sec. Meier-Reimer (21) indicated that a diffusion veloci-
ty of the order of .5 cm/sec should be appropriate for most seas and oceans.
Some adjustment of the stochastic coefficient may be necessary to achieve an
average "diffusion velocity" of .5 cm/sec or some other characteristic vel-
ocity. Both Meier-Reimer (21) and Thompson (22) have presented some work
with cases where the diffusion coefficient was allowed to vary in time and
space. Extension of the present work to such cases should present no serious
Difficulties.
30
-------
SECTION 5
RESULTS AND DISCUSSION
This section presents the results of the various extensions and enhance-
ments through exposition of several figures. The figures present selected
results for a hypothetical area designed to verify the procedures implemented.
Literally thousands of figures, interesting in their own right, could have
been prepared from the computational data; however, a few representative
figures were deemed sufficient to present the case for the procedures.
NON-LINEAR TERM ENHANCEMENT
A special coastal grid with depth profile depicted by Figure 13 was de-
signed to verify the procedures introduced to implement the non-linear terms.
The 20X30 grid was characterized by longshore uniformity. The "optional"
island indicated in Figure 13 was the only break in the uniform offshore
depth profile when present in column 13 of the z grid. A combination of fac-
tors including depths surrounding the island, the longshore uniformity and
the grid length selected (1.5 Km) resulted in essentially no significant dif-
ferences between linear and non-linear cases of like bathymetry over a tidal
period.
MSL
100
X
frl
I
200
300
400
500
Figure 13.
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
COASTAL GRID ROW NUMBER (n)
MSL — mean sea level
h depth profile for coastal grid.
31
-------
The grid developed to verify procedurally the tidal flat enhancement was
selected to verify the non-linear enhancement. Figures 14 and 15 present the
tidal flat grid and associated bathymetry. TABLE 2 provides detailed grid and
run specifications for the typical computer model run involving the tidal flat
grid.
1
2
3
4
S
6
7
9
10
11
12
13
14
IS
•
1
1
1
1
e
^
C-
\i
n
v
^
•s
r-
•T
•
MI
Ci
k:
K
c-
k:
\
/
) 1
^
J
i i:
t i
(i
k!
(
V)
(-
k.
3 1
>)
^
Ml
/V
>|
y
i i
•A
^
j
4.
S 1
WH
i 1
m*
7 1
K4H
1 1
HMM
» 20
MM-
(D) identify z cells referred to in Figures 15 - 20
(4) identify bathymetric features Fl - F4, respectively
(shallow subsurface knobs referred to by TABLES 3-6)
tidal input boundary (open)
_._ . _ land-sea boundary (initial condition)
iiiiiiniiiiiiii dynamic open boundary (no tidal prescription)
thermal source cell
L..J
Monte Carlo source locations (approx.;
Figure 14. Tidal flat grid.
32
-------
m
1 2 34 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
2
3
4
5
6
7
8
9
10
11
12
13
14
15
r-
Figure 15. Tidal flat grid bathymetry (depths in cm).
33
-------
TABLE 2. TYPICAL TIDAL FIAT GRID AND RUN SPECIFICATIONS
DATA
CARD
NAME
GRID
LAYS
TIME
FACT
TIDE*
COMP
PROGRAM
VARIABLE
ICl{2)-(9)
NE
ME
DL(2£)
ROTANG
ALAT
R(r)
00
HL(1)
ALPHA (o)
RHO(p)
HMIN
MODOUT
ITO
ITS(t0)
ITE
rriNC(2T)
IPO
IPMOD
TIDSPD(l)
ITIDE(l) , (2)
ITIDE{3) , (4)
ITIDE(5)
ITIDE(l)
ITIDE (5)
MODEL
VARIABLE
DESCRIPTION
Grid description
No. of Rows
No. of Columns
Grid length
Rotation angle
Average latitude
Friction coefficient
No. of layers
Maximum depth of layer
Smoothing parameter
Density
Minimum allowable depth
Data save step
Data save start
Start computation
End computation
Computational step
Print step
Print start
Constituent speed
Min/max row spec.
Min/max column spec.
Tidal offset time
Constituent phase shift
Constituent amplitude
Grid I.D. No.
VARIABLE VALUE ASSIGNMENT
EPA TIDAL FLAT TEST GRID
15
20
1/3 Kta
180° (implies (1,1) is the SE
corner of grid)
45° (used to compute f and g)
.003
1
660 cm
.99
1.022
1.0
120 sec
120 sec (for runs saving data)
0 sec
31200 sec (8.67 hrs)
30 sec
1200 sec
1200 sec
45°/hr (false "tide"-8 hr period)
1,15
20 (western boundary specifica-
tion)
0 sec
90°
130 cm
998
*The open portions of the northern and southern boundaries were
computed dynamically based on internal information and boundary
conditions as specified previously.
34
-------
TABLES 3-6 present linear and non-linear configuration development of
the sometime tidal flat "islands" associated with bathymetric features Fl
through F4 as identified in Figure 14. Configuration development in TABLES
3-6 is shown in full cell steps. The differences between the linear and non-
linear computational modes is apparent. The non-linear model tends to "lead"
the linear model with notable exceptions. The major exceptions relate to a
change in the order of appearance of individual cells when the model is
shifted into the non-linear mode. These variations in order and relative
timing of appearance of individual cells in the development of both flats F3
and F4 imply a redistribution of the velocity field in the vicinity of these
features.
In general the amplitude of the wave incident on the eastern closed
bound (the back of the "bay") is less for the non-linear model than for the
linear model. This is in agreement with an energy shift from longer to
shorter wavelengths for the non-linear model.
LAYER DISAPPEARANCE ENHANCEMENT
A tidal flat grid and bathymetry were designed to verify thoroughly the
procedures introduced in the layer disappearance enhancement of the HN model
as shown in Figures 14 and 15. TABLE 2 provides detailed grid and run spec-
ifications for the typical computer model run involving the tidal flat grid.
TABLES 3-6 present configuration development of the tidal flat segments as-
sociated with bathymetric features Fl through F4 as identified in Figure 14.
All configuration development in TABLES 3-6 is presented in full cell steps.
These tables were designed to show the differences between linear and non-
linear model versions but quite adequately show the functional aspects of
the layer disappearance modifications in those computational regimes.
For a more detailed view of a particular closure case, Figures 16
through 19 present a series of snapshots of actual configurations showing
tidal flat development in cell side closure steps. These figures can be
directly compared to the non-linear case configurations shown in TABLE 3.
Figure 16 presents an "initial" closure configuration with specific ref-
erence to cells A, B, C and D. Initially A is closed on the north and B, C
and D are open. In Figure 17 the cells are shown at a subsequent time where
A is closed to the north and east, D is closed to the south, and B and C are
open. Figure 18 presents a later time showing A fully closed while B and D
are each closed on a single side. Figure 19 presents a still later time
where D becomes closed on two sides.
Figures 20 and 21 present HN computed time series data at time step in-
tervals for £ at A, B, C and D for separate periods. Tidal prescription
from a tidal boundary point is also included for reference. Figure 20
brackets the snapshot times presented in Figures 16-19 showing the relative
"sea height" for each of these cells through time compared to the tidal pre-
scription. Note the gradual slope change in c down to the point of closure
A
of cell A. A similar character is observed in the profile of t subsequent
to the closure of the south side of cell D (Figure 17.) Neither the closure
of a side of cell A, nor the full closure of cell A, introduced any evidence
of instability locally in cell A or adjacent cells B, C and D.
35
-------
TABLE 3: FLAT F3 FORMATION; NON-LINEAR VS LINEAR
FLAT CELL
APPEARANCE
ORDINAL
1
2
3
4
5
6
7
8
9
10
11
NON-LINEAR
(n,m)
LOCATION
(6,11)
(6,12)
(7,11)
(5,11)
(6,10)
(6,13)
(7,12)
(5,10)
(5,12)
(4,11)
(8,11)
TIME
16980
17280
17730
17760
18210
18240
18390
18480
18870
19050
19320
CONFIGURATION
0+
0++
°t+
*
mil,
O7TT
®+++
®tt+
^
tfr
*s+
sjr
LINEAR
(n,m)
LOCATION
(6,11)
(6,12)
(7,11)
(5,11)
(6,10)
(6,13)
(7,12)
(5,10)
(5,12)
(8,11)
(4,11)
TIME
17040
17430
17760
17940
18030
18300
18420
18420
19140
19260*
19410
CONFIGURATION
0+
O++
C++
*
®++
®+++
nr
i|r
$n+
f*
f*
* Connects to Flat F4 to form F34 at _.
O Represents shallow knob at (6,10) at initial condition
(Bathymetric feature F3, Figure 13.)
36
-------
TABLE 4: FLAT F4 FORMATION; NON-LINEAR VS LINEAR
FLAT CELL
APPEARANCE
ORDINAL
1
2
3
4
5
6
7
8
9
10
11
NON-LINEAR
(n,m)
LOCATION
(11,10)
(11,11)
(12,10)
(10,10)
(11,12)
(11,9)
(12,9)
(12,11)
(10,11)
(13,10)
(9,10)
—
TIME
17100
17340
17640
18060
18180
18300
18510
18570
18600
18840
19050
19320*
CONFIGURATION
O+
0++
°r
°r
O+++
e+++
$l++
?U+
fS+
$tt+
u
$p+
(no +++
change) ?++
LINEAR
(n,m)
LOCATION
(11,10)
(11,11)
(12,10)
(10,10)
(11,9)
(11,12)
(12,9)
(10,11)
(12,11)
(13,10)
(9,10)
—
TIME
17130
17400
17760
18060
18120
18210
18340
18660
18720
18990
19020
19260*
CONFIGURATION
O+
O++
OJ+
o|+
e|+
©+++
jft+
&+
fS+
f|r
fir
--+
(no ffi+++
change) +
* Flat F3 connects to F4 by accretion to form Flat F34 at "—".
0 Represents shallow knob at (11,9) at initial condition
(Bathymetric feature F4, Figure 13.)
37
-------
TABLE 5: FLAT Fl FORMATION; NON-LINEAR VS LINEAR
FLAT CELL
APPEARANCE
ORDINAL
1
2
3
4
5
6
7
8
9
10
NON-LINEAR
(n,m)
LOCATION
(7,6)
(6,6)
(7,5)
(8,6)
(7,7)
(5,6)
(7,4)
(6,5)
(8,5)
—
TIME
18150
19410
19530
20700
20910
21660
21960
22110
22560
25290*
CONFIGURATION
e
I
4
*
t jh »
1 VII I
I g\ I
1 \Lf 1
++$+
+tf+
+if+
+fr
LINEAR
(n,m)
LOCATION
(7,6)
(6,6)
(7,5)
(8,6)
(7,7)
(5,6)
(7,4)
(6,5)
(8,5)
(4,6)
TIME
17970
19110
19230
20310
20520
21120
21390
21510
21990
25050*
CONFIGURATION
9
$
4
4
i A i
T W>
+
i A i
1 W i
+
+4+
++!+
+9+
+g;
* Flat Fl is connected to Flat F234 to form "U" shaped tidal
flat F1234 at "-".
O Represents shallow knob at (7,6) at initial condition
(Bathymetric feature Fl, Figure 13.)
38
-------
TABLE 6: FLAT F2 FORMATION; NON-LINEAR VS LINEAR
FLAT CELL
APPEARANCE
ORDINAL
1
2
3
4
5
6
7
8
9
NON-LINEAR
(n,m)
LOCATION
(11,6)
(12,6)
(11,5)
(10,6)
(11,7)
(11,8)
(13,6)
(12,7)
(12,8)
TIME
18180
19350
19590
20670
20910
20970
21810
21870
21990
CONFIGURATION
®
f
+f
*
i m i
T w 1
+S+-H
+S+-H
+
+|r
+|«
LINEAR
(n,m)
LOCATION
(11,6)
(12,6)
(11,5)
(10,6)
(11,7)
(11,8)
(13,6)
(12,7)
(12,8)
TIME
18000
19050
19260
20310
20460
20550*
21270
21330
21420
CONFIGURATION
®
f
+f
+i
+®+
+I+-H
i fK i j i
T'wTTi
+
+|r
*|ttl
* Flat F2 connects by accretion to Flat 34 to form lazy "L" shaped
tidal flat F234 at "'".
0 Represents shallow knob at (11,6) at initial condition
(Bathymetric feature F2, Figure 13.)
39
-------
ffi
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Figure 16. Closure configuration TFC: t = 17700 sec.
1 2 34 5 6 7 • 9 K) 11 12 13 14 15 14 17 18 19 20
2
3
4
5
6
7
•
9
10
11
12
13
14
15
I
c?
t=18090
93.7
Figure 17. Closure configuration TFC: t
40
= 18090 sec.
-------
2
3
4
S
6
7
8
9
10
U
12
13
14
15
m
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
t =18240
» X"/
-96.6
mm m mmm. m m j
Figure 18. Closure configuration TFC: t = 18240 sec.
m
34 56 78 9 10 11 12 13 14 15 16 17 18 19 20
Figure 19. Closure configuration TFC: t = 18390 sec.
41
-------
80
fa
I
E
g
-90
100
TIDAL FLAT (6,13)
CONVECTIVE CASE
,
"***' SL£S1P—-
N=
f I Closure of indicated sided) E W
L-N-1
T Tidal Input at M=2O
B
17700
17850
18000
time (sec)
18150
18300
18450
Figure 20. Closure analysis for TFC cell (6, 13)
o
-50 -
fr.
w
I
en
o:
I
.TIDAL FLAT (6,13)
T Tidal Input at M = 2O
2*550
26700
27150
27300
26850 27000
time (sec)
Figure 21. Post-reopening analysis for TFC cell (6, 13).
42
-------
Figure 21 presents a much later period of time than Figure 20. Cells A,
8, C and D had just been reopened by the initial time of Figure 21. Although
waves of a period shorter than the "tidal cycle" are evident in the data,
there is no indication of local instability at any of the reference points.
Figure 22 depicts a relatively advanced stage of tidal flat development in-
dicating maximum seaward advance of the tidal flat array. The connection of
tidal flats F2, F3 and F4 have created a semi-enclosed tidal pool which drains
to the south and east. Later a "U"' shaped configuration, draining only to the
east, develops in both the linear and non-linear modes.
1
2
3
4
5
6
7
n 8
9
10
11
12
13
14
15
1
m
2 34 56 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1"
*
t=22022
Figure 22. Maximum seaward advance of tidal flat F234.
43
-------
THERMAL ADVECTION EXTENSION
The tidal flat grid as presented in Figures 14 and 15 was selected for
procedural testing of the thermal advection extension. TABLE 2 presents
typical tidal flat grid and run specifications as before; however, the com-
putational time step was longer than the 30 sec step for the usual HN runs and
the "continuous" heat pollution source was begun at 1200 sec. The change in
time step is due to the fact that the procedure is driven from data saved at
120 sec intervals during a prior HN computation of ?, U and V. A time step
of 120 seconds was selected when it was determined that differences from 240
second step computations were significant although not substantial.
Figure 23 presents the thermal advection plume, the region affected by
the residual heated fluid, after 8 1/3 source hours (approximately one "tidal"
cycle). The source is essentially a single cell source depending on extant
quantity of heated fluid. Heat budget computations occur to derate the excess
heat over all areas where excess heat exists considering output radiation flux,
latent heat transfer and convective sensible heat loss. Some portions of the
advected volume thereby lose all of the associated excess heat. Volumes which
have become associated with zero excess heat units, i.e., have been reduced to
ambient temperature, are removed from the advection computation to become part
of the background ambient volume. This occurs most noticeably over regions
where volumes are spread over the tidal flat. The plume is thus split by the
emergence of the tidal flat. Over several cycles the resulting plume "cast-
offs" should be completely derated to ambient (based on observed numerical
rates of cooling) or pass through the region boundary (lost to the computation
in either case). These computational results seem reasonable in spite of the
assumptions of no diffusion and no density dynamic effects.
MONTE CARLO DIFFUSION EXTENSION AND COMPARISON
The tidal flat grid as presented in Figures 14 and 15 was selected for
procedural verification of the Monte Carlo extension. TABLE 2 presents
typical tidal flat grid and run specifications as before. Here the Monte
Carlo technique was employed as a special feature during computation of £,
U and V, rather than as an after the fact procedure like the thermal advec-
tion extension. Four "continuous" Monte Carlo sources were started simultan-
eously at 1200 sec during an HN run from 0 to 31200 sec.
Figure 23 presents the Monte Carlo "plume," the approximate location of
4000 tracer points, after 8 1/3 hours (approximately one "tidal" cycle).
Note the strong relationship of the Monte Carlo sources to the thermal advec-
tion source and strong correlation of the two "plumes." In spite of the dif-
ference in time step, the fact of the randomized diffusive effects in the
Monte Carlo scheme, and the slight difference in the effective center of the
two source types, the two methods gave distributions with very strong resem-
blance. The region where the separation occurs in the thermal advection
plume is associated with a region of very low density of Monte Carlo points.
The "holes" which appear in the Monte Carlo distribution are an artifact
of the separation of the four continuous point sources which were designed to
present the outlines of a single cell source similar to the thermal advection
44
-------
source. The triangular shape of the Monte Carlo distribution is the result
of the divergence caused by the two exits to the tidal pool presented in
Figure 22.
m
4 6 • 10 12
14
16
18
20
10
12
14
31200
CONSTANTS:
T imkient - 10 C
Del. Humiditj - 85%
CloriCow - 2%
T »ir - 8*C
Monte Carlo Dispersion Plume
Thermal Advection Plume
• Monte Carlo Source (30 second step)
LJ Thermal Source (120 second step)
Resultant Thermal Advection Cell Temperature Anomaly
(AT in excess of ambient temperature) is indicated in
degrees Celsius for each cell affected by the residual
heat loss.
Figure 23. TFC continuous source "plumes."
45
-------
REFERENCES
1. Bauer, R. A. Description of the Optimized EPRF Multi-layer Hydrodynam-
ical-Numerical Model. ENVPREDRSCHFAC Tech. Paper No. 15-74, Environ-
mental Prediction Research Facility, Monterey, California, 1974. 66 pp.
2. Hansen, W. Tides. In: The Sea, Vol. I, M. N. Hill, ed. Interscience
Publishers, London, U.K., 1963. pp. 764-801.
3. Hansen, W. Hydrodynamical Methods Applied to Oceanographic Problems.
In: Proceedings of the Symposium on Mathematical-Hydrodynamical Methods
of Physical Oceanography, Institut Fur Meereskunde der Universitat
Hamburg, Hamburg, West Germany, 1962. pp. 25-34.
4. Laevastu, T. and P. Stevens. Applications of Numerical-Hydrodynamical
Models in Ocean Analysis/Forecasting. FNWC Tech. Note No. 51, Fleet
Numerical Weather Central, Monterey, California, 1969. 70 pp.
5. Laevastu, T. and K. Rabe. A Description of the EPRF Hydrodynamical-
Numerical Model, NEVPREDRSCHFAC Tech. Paper No. 3-72, Environmental
Prediction Research Facility, Monterey, California, 1972. 49 pp.
6. Laevastu, T. A Vertically Integrated Hydrodynamical-Numerical Model
(W. Hansen Type), Model Description and Operating/Running Instructions.
Part I of a series of four reports. ENVPREDRSCHFAC Technical Note
No. 2-74, Environmental Prediction Research Facility, Monterey, Cali-
fornia, 1974. 63 pp.
7. Laevastu, T. A Multilayer Hydrodynamical-Numerical Model (W. Hansen
Type), Model Description and Operating/Running Instructions. Part 2
of a series of four reports. ENVPREDRSCHFAC Technical Note No. 2-74,
Environmental Prediction Research Facility, Monterey, California, 1974.
54 pp.
8. Laevastu, T. in collaboration with M. Clancy and A. Stroud. Computa-
tion of Tides, Currents and Dispersal of Pollutants in Lower Bay and
Approaches to New York with Fine and Medium Grid Size Hydrodynamical-
Numerical Models. Part 3 of a series of four reports. ENVPREDRSCHFAC
Technical Note No. 3-74, Environmental Prediction Research Facility,
Monterey, California, 1974. 51 pp.
9. Laevastu, T. and R. Callaway in collaboration with A. Stroud and M.
Clancy. Computation of Tides, Currents and Dispersal of Pollutants
in New York Bight from Block Island to Atlantic City with Large Grid
Size, Single and Two-Layer Hydrodynamical-Numerical Models. Part 4
46
-------
of a series of four reports. ENVPREDRSCHFAC Technical Note No. 4-74,
Environmental Prediction Research Facility, Monterey, California, 1974,
79 pp.
10. Leendertse, J. J. Aspects of a Computational Model for Long-Period
Water-Wave Propagation. USAF Project RAND Memorandum RM-5294-PR, the
RAND Corporation, Santa Monica, California, 1967. 165 pp.
11. Laevastu, T. and G. D. Hamilton. Computations of Real-Time Currents
Off Southern California With Multilayer Hydrodynamical-Numerical Models
with Several Open Boundaries. ENVPREDRSCHFAC Technical Paper No. 10-74.
72 pp.
12. Laevastu, T. and A. D. Stroud. Exchange between Coastal and Continental
Shelf Waters as Computed with a Two-Layer HN Model, Manuscript Report,
ENVPREDRSCHFAC Report to the Environmental Protection Agency, Corvallis
Environmental Research Laboratory, Environmental Prediction Research
Facility, Monterey, California, 1975.
13. Lamb, H. Hydrodynamics. DOVER Publications, New York, New York, 1945.
738 pp.
14. Laevastu, T. Multilayer Hydrodynamical-Numerical Models. In: Proceed-
ings of the Symposium on Modeling Techniques, Am. Soc. Civil Engineers,
San Francisco, California, 1975. pp. 1010-1025.
15. Kagan (Cohen), B. A. Properties of Certain Difference Schemes Used in
the Numerical Solution of the Equations for Tidal Motion. IZVESTIA,
Atmospheric and Oceanic Physics, 6(7):704-717, 1970. Translated by
F. Goodspeed.
16. Pedersen, L. B. and L. P. Prahm. A Method for Numerical Solution of
the Advection Equation 8C/9t=-V'VC, Meteorological Institute of Den-
mark, as submitted to TELLUS, August, 1973. 36 pp.
17. Laevastu, T. and J. M. Harding. Numerical Analysis and Forecasting of
Surface Air Temperature and Water Vapor Pressure. Journal of Geo-
physical Research, 79(30) :4478-4480, 1974.
18. Laevastu, T., K. Rabe and G. D. Hamilton. The Effects of Oceanic
Fronts on Properties of the Atmospheric Boundary Layer. ENVPREDRSCHFAC
Technical Paper No. 2-72, Environmental Prediction Research Facility,
Monterey, California, 1972.
19. Laevastu, T. Factors Affecting the Temperature of the Surface Layer
of the Sea, Soc. SCIENT, FENNICA, Comment, Physico-Mathem., Helsinki,
1960. 136 pp.
20. Young, Chen-Shyong. Thermal Discharges into the Coastal Waters of
Southern California. Southern California Coastal Water Research Pro-
ject (SCCWRP), Los Angeles, California, 1971. 30 pp.
47
-------
21. Maier-Reimer, E. Numerical Treatment of Horizontal Diffusion and
Transport Phenomena in Marine Basins of Large Size, Institiit fur
Meereskunde der UniversitSt Hamburg, Hamburg, West Germany, as pre-
sented at IAMPA/IAPSO Assembly, Melbourne, Australia, January, 1974.
22. Thompson, R. Numerical Calculation of Turbulent Diffusion. Quart. J.
R. Met. Soc. 97(411):93-98, 1971.
23. Intersea Research Corporation. Semiannual Report of San Onofre Ocean-
ographic and Biological Monitoring Program—July 1972-January 1973,
prepared for Southern California Edison Company. La Jolla, California,
August, 1973. pp. 1-7, 17-30, 190-196.
48
-------
APPENDIX A
VERIFICATION OF THE MODEL
After the development of the various enhancements on the model using
dummy data sets to test the code an attempt was made to verify the model
using a real data case. Since data for the San Onofre Nuclear Power Sta-
tion are reported to the San Diego Region Water Quality Control Board and
were readily available, these data were selected for the test. Data are
collected for San Onofre on a bimonthly schedule and are reported semi-
annually. For 2 August 1972 the report (23) provides the tide curve, in-
termittent drogue measurements from 0800 to 1200 PDT, three infrared over-
flights, periodic wind measurements and bathythermograph measurements. Dur-
ing this period the plant was reported to be operating normally pumping
350,000 gallons of cooling water per minute and raising the temperature
approximately 10°C.
MODEL SPECIFICATIONS
A single layer 20 by 40 HN grid with 500' mesh size oriented with the
major axis along the coast was used to cover the San Onofre outfall region
as shown in Figure A-l. Depths were interpolated to tenths of feet based on
a bathymetric chart showing 6 foot depth contours from mean lower low water.
An average latitude of 33°22'30"N was used to determine the mean gravitation-
al constant and the mean coriolis effect over the grid. The usual bottom
friction coefficient of .003 and a smoothing coefficient of .99 were used.
The maximum depth at approximately 50 feet and the selected grid size re-
sulted in a maximum time step for the HN scheme of 8 seconds through the CFL
stability criterion.
The tides for the region were computed using the four major tidal con-
stituents for the area, Kl, Ol, M2 and S2, adjusted for time and position.
The tidal components were specified at three separate points on the seaward
open boundary with temporal separation so that the tidal wave would propagate
in the longshore direction toward the northwest. All additional points along
the seaward boundary between the computed values were specified by linear in-
terpolation. The tidal curve is shown in Figure A-2.
The southeastern and northwestern open boundaries to the region were not
specified from external data except at the intersections of these boundaries
with the seaward boundary. Heights at these boundaries were based solely on
the dynamic boundary conditions. Observed wind values given in Figure A-2
were specified as uniform wind fields over the entire grid. The cooling
water flow rate was used to define a "source" at the outfall positions and a
"sink" at the intake position throughout the dynamic computations.
49
-------
m
1 3 S 7 9 11 13 IS 17 19 21 23 25 27 29 31 33 35 37 39
4 intake
• outfall
generating station
land-sea boundary
depths in feet below mean lower low water
................... interpolated tidal boundary
east and west boundaries dynamically computed
Figure A-l. San Onofre outfall area grid bathymetry.
50
-------
Tide at column 23
^^^^^ ^l • 1
3 ^4^ J
is
-s _
ii i
/ 12 15 jl
•
M . .-Direction - iSWi SW
Measurement . ./ .x . .V ,«: ,»
L Speed (mph) 0 J5-10C 10
Speed used 0 i 7.5: 10
SW
10-15
12.5
is
0§
•
SSW
10
10
SW SSW -
5 s ;o
5 5:0
5
1
i
i
I
i
Low Water
i i i i 1 | """"
98848 109648 120448 131248 142048 152B48
Tide at column 23 was computed from MONACO constituents
using an offset of 160345 seconds.
The tide at column 1 leads column 23 by 260 seconds.
The tide at column 40 lags column 23 by 209 seconds.
Tides were prescribed at points 2-22 and 24-39 by
linear interpolation.
Figure A-2. San Onofre outfall tidal and wind prescription August 2, 1972.
51
-------
Both the non-linear and layer disappearance enhancements were used to
simulate the dynamics for the region. No estimate was made of the impact of
the non-linear terms on the computation; however, it was noted that the layer
disappearance impact was minimal. The dynamic computations were initiated
one full tidal cycle prior to the simulation date to assure prior proper con-
vergence of the dynamic solution prior to application of the currents to drive
the thermal advection problem. Data were saved at 120 second intervals during
the dynamic computations for the thermal advection problem.
The time required to digitize the bathymetry, set up the control cards,
verify the Phase I results and run the Phase II model was less than two man
days. The Phase II model was stopped after each 25,000 second period out to
100,000 seconds to examine the results and then the model was restarted. As
shown in Figure A-2 the model was then run for the tidal cycle that matched
the 2 August conditions.
The integrated velocities calculated during the period of the drogue
measurement agreed with the velocities measured by the drogues set at a depth
of 15 feet. The calculated velocity pattern is shown in Figure A-3. The
velocities are between 35 and 40 cm/sec while the drogues had velocities of
between 35 and 43 cm/sec with both showing the characteristic movement along
shore to the northwest during the period of rising tide. Figures A-4 and A-5
show lower velocities and finally the reversal inshore as the tide reaches a
maximum and begins to ebb.
THERMAL ADVECTION
The velocity fields saved from the Phase II model at 120 second inter-
vals were used for the thermal advection model run beginning at the model
time equivalent to 2022 PDT on 1 August 1972. Figures A-3 through A-5 show
the plume 15, 17 and 19 hours after the start and compare the computed plume
to the infrared measurements. Only two contours are compared for each fig-
ure, a contour 1°C (1.8°F) above ambient temperature and a contour 6.56°C
(11.8°F) above ambient. The computed values show the highest values down-
stream in the plume rather than the "bubble" over the outfall. Since dif-
fusion was not considered horizontally or vertically, no spreading of the
heat to surrounding waters occurs. The hottest portions of the computed
plume passed back and forth over the outfall as the tidal currents reversed
receiving a double dose of heating. Even though the radiation terms in the
model were reducing the heat content of this water, the limited surface ex-
posed did not produce sufficient cooling. The temperature of the plume did
not drop to the level of the background temperature within the area modelled
as expected. Several minor fixes were attempted to improve the agreement
between the observed and modelled results by modifying the widths of the
heated portions of the cells and changing constants in the radiation equa-
tions but all attempts simply showed that the numbers obtained reflected the
basic assumption of pure advection and that a major effort beyond the scope
of the previous tasks would be required to adequately include diffusion in
the model.
52
-------
m
1 3 S 7 9 11 13 IS 17 19 21 23 25
-31-A3-.H--K-.K
i
3
5
7
n »
11
13
IS
17
19
•71
SM
Contours: 21.7°C (71.0°F) and 27.2°C (81.0°F). Ambient: 20.7°C (69.2°F).
infrared computed —* 50 cm/sec
Figure A-3. Infrared flight 0958-1107 compared to values after 15 source hours.
1357 9..J1..j3..J5..1^.J?..2J..23._25..27..29..3.1...33..35..3?..3?
3
S
7
9
11
13
IS
17
19
--77—
WWW
Contours: 22.2°C (72.0°F) and 27.7°C (82.0°F). Ambient: 21.2°C (70.2-F).
infrared computed —*50 cm/sec
Figure-A-4. Infrared flight 1226-1344 compared to values after 17 source hours.
53
-------
m
1 3 5 7 9 11 13 IS 17 19 21 23 25 27 29 31 33 35 37 39
3
5
7
9
11
13
IS
17
19
Contours: 22.5°C (72.5°F) and28.1°C (82.5°F). Ambient: 21.5°C (70.7°F).
infrared computed >50 cm/sec
Figure A-5. Infrared flight 1503-1617 compared to values after 19 source hours.
54
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1. REPORT NO.
EPA-600/3-78-073
3. RECIPIENT'S ACCESSION NO.
4. TITLE AND SUBTITLE
Enhanced Hydrodynamical-Numerical Model
For Near-Shore Processes
5. REPORT DATE
July 1978
6. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
R.A. Bauer and A.D. Stroud
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Compass Systems, Inc.
San Diego, California 92109
10. PROGRAM ELEMENT NO.
EHE 625
11. CONTRACT/GRANT NO.
68-03-2225
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Research Laboratory-Ccrvallis
Office of Research and Development
U.S. Environmental Protection Agency
Corvallis. Oregon 97330
13. TYPE OF REPORT AND PERIOD COVERED
final 1975-197 7
14. SPONSORING AGENCY CODE
EPA/600/02
15. SUPPLEMENTARY NOTES
16. ABSTRACT
An optimized version of a multilayer Hansen type Hydrodynamical-
Numerical (HN) model is presented and discussed here as the basis for the
following experimental extensions and enhancements developed to more ap-
propriately handle near-shore processes:
•Non-linear term extension to facilitate small-mesh studies of near-
shoreT including river inflow dynamics;
•Layer disappearance extension to enable appropriate procedures in
tidal flat and marshy regions, as well as some down/upwelling cases;
•Thermal advection enhancement for treatment of thermal pollution
cases by method of moments coupled with heat budget procedures for
dynamic plume development experiments;
•Monte Carlo diffusion enhancement to deal with dispersion via sta-
tistical methods and comparison to the method of moments experiments,
Extensive efforts were invested in determining reasonable and appropriate
boundary conditions for both the basic model and the extended versions
presented here.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS C. COSATI Field/Group
numerical methods
thermal pollution
monte carlo methods
mathematical methods
near-shore
river inflow
titdal flat
marshy regions
operations research
08/C
12/A,B
18. DISTRIBUTION STATEMENT
Release to Public
19. SECURITY CLASS (This Report)
Uncla ssified
21. NO. OF PAGES
70
20. SECURITY CLASS (This page)
Unclassified
22. PRICE
EPA Form 2220-1 (Rev. 4-77)
* U. S. GOVERNMENT PRINTING OFFICE.- 1978—797-713/228 REGION I 0
55
------- |