i

f 1
zZ—~ —u—
n=l [| :
v v
| : :
| : .
j : u 1Zeeeu +Z- u- Zz u
2 Lh NT:
Sv ve(2,2) v SSN LANDS y SSSSS
| : :
Ze. vou
3 |
Y
|
|
Z— —u— — u
« [ :
v v v v
zZ- & computational point SSK land boundary
u - U computational point ——— open boundary
v- V computational point 3 ..... z grid structure

Figure 1. Scheme of the numerical grid.


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/coc, (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 21 is provided by the heat budget formulation:

K = QA2t = (Qs-Qr-Qb-Qe-Qh+Qm) A2t; (35)
thus
AT =  (Qs-Qr-Qb-Qe-Qh+Qm) 2t/cpa (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


Cc = sup -z--u}-z-u

1
Vv
t :
|
z-lyjoz-u-
»L i
v @

1
de -u- ~z- -u- -u
=f

v v

i
if
_
4

<
if —N— =<} —

Vv
UY

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: a[v(x,y)*T]/aN|, = Lima [vGctes,y+6s) +71 /2N. (33)
s>0

The parameters € 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


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.

123 4 5 6 7 8 9 011 12:13:14 15 1 17 B19 20

a a a 5
af J
srt N
s 1 4
A
6
7H » i
s 4
“Hi
ert
) 4) iz
af
3
“ tl slams e e
s i «| a + i

®® Ss identify z cells referred to in Figures 15 - 20
®® @ identify bathymetric features Fl - F4, respectively
(shallow subsurface knobs referred to by TABLES 3-6)
sateeteeteee eee tidal input boundary (open)
== seesem jand-sea boundary (initial condition)
+H dynamic open boundary (no tidal prescription)

C} thermal source cell

TM Monte Carlo source locations (approx.)

Figure 14. Tidal flat grid.

32


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. Haines for preparing the
illustrations and tables and to D. Bauer for typing the report.


LAYER SURFACE FROM MSL (cm)

LAYER SURFACE FROM MSL (cm)

TIDAL FLAT (6,13) . , CONVECTIVE CASE
“© WS open, |
. ~~ — N
‘Sle SE
LNSey, ‘ J
~ siosto_ J,
-90 4
mee
.
~s. TM“B
—
N ane,
=100 B rs c
Closure of indicated side(s) w t D
6lAi Cc " tnd E
7\/D T Tidal Input at M=20
L L 1 1
1700 17850 18000 1150 18300 18450
time (sec)
Figure 20. Closure analysis for TFC cell (6, 13).

T ~—
| TIDAL FLAT (6,13) CONVECTIVE CASE |

T Tidal Input at M=20
1

“
o

26550 26700

Figure 21.

26850

27000

time (sec)

42

27150 27300

Post~reopening analysis for TFC cell (6, 13).


12 3 4 5 6 7 8 9 WH 2 3 4 1 6 7 WB WY 2

_ 7 t=17700

Figure 16. Closure configuration TFC: t = 17700 sec.

Fa A

—s om “= et ed—)- t=18090

Figure 17. Closure configuration TFC: t = 18090 sec.

40


REFERENCES

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.

Hansen, W. Tides. In: The Sea, Vol. I, M. N. Hill, ed. Interscience
Publishers, London, U.K., 1963. pp. 764-801.

Hansen, W. Hydrodynamical Methods Applied to Oceanographic Problems.
In: Proceedings of the Symposium on Mathematical-Hydrodynamical Methods
of Physical Oceanography, Institut Fiir Meereskunde der Universitat
Hamburg, Hamburg, West Germany, 1962. pp. 25-34.

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.

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.

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.

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.

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.

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


13) 5 7 9 11 13 15 17 «19-2123 25-27 29 31 33 35 3789

0 woo 3000 FEET
300 METERS,

+ intake ® outfall @ generating station

ad

Wes land-sea boundary
depths in feet below mean lower low water

senanecaesensennes interpolated tidal boundary

east and west boundaries dynamically computed

Figure A-l. San Onofre outfall area grid bathymetry.

50


TABLE 5: FLAT Fl FORMATION; NON-LINEAR VS LINEAR

FLAT CELL NON-LINEAR LINEAR

APPEARANCE (n,m) (n,m)

oRDINAL | LocaTzon | Trme | ConFIGuRATION | LocaTIon | TIME | CONFIGURATION
1 (7,6) | 18150 e (7,6) | 17970 ®
2 (6,6) | 19410 $ (6,6) | 19110 $
3 (7,5) | 19530 +6 (7,5) | 19230 +b
4 (8,6) | 20700 +§ (8,6) | 20310 +8
5 (7,7) | 20910 ++ (7,7) | 20520 +o

+
6 (5,6) | 21660 he (5,6) | 21120 ++
+
7 (7,4) | 21960 as (7,4) | 21390 +g
+
8 (6,5) | 22110 tb (6,5) | 21510 +g
+f +f
9 (8,5) | 22560 +g (8,5) | 21990 +g
+
ai
10 -- 252904 ie (4,6) | 250504 +

* 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



Condition 2: clr, = 0(x,y,t), (12)

condition 3: aivem/an|T, = 0, (13)
Condition 4: aH/an|T, = 0, (14)
Condition 5: ac/en|T = 0, (15)
Condition 6: a(ver)/am|r = 0, (6)
where:
v --vector velocity,
N --a unit vector normal to the indicated boundary,
? --a unit vector tangential to the indicated boundary,
wen --vector "dot" product,
| --may be interpreted as “evaluated on",
ON --an element of the normal to the indicated boundary.

(Note that v°T=U and v*N=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 ¢ 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


order as the sub- or superscript that is present.) The sequence of pro-
cedures follows:

1. tty. continue at procedure #7; (initialization)
t t _ .
2. (Deore Tegan ME ear) “thermal (Tw, Ta, Aw.) ; (thermal step)
t . : *
3. Teor Mega MM eaar) / (RERY) ; (inverse interpolation)
t t,t .
4. Rear Pee AT yor? (added heat computation)
t+2T t ot +
5. (Cc 1K 97) =advect (C Regar) (u,v) * (advection step)
+2 +:
crkre becomes RE T Ry ary
6. t=t+2T; (full step accomplished)
7. ctectsack, (also adjust rid); (pollution volume increment)
8. KHK TAK (pollution heat increment)
9. AT =K,/ (pec®) ; (temperature change computation)
t
10. Tw. 7 RERY AM, (cell temperature interpolation)

11. when t=te stop, otherwise continue at procedure #2;

where:

Rx, Ry ~-horizontal width factors in grid step units related to second
moment of pollutant in cell (Equation 17),

Ta --air temperature,

te --end time,

ty --start time,

Tw --surface water temperature,

aT temperature in excess of ambient for the heat polluted volume,

Aw ambient water temperature,

thermal thermal procedure

advect --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


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

iii


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:

++ -tUi (Ui) xtVi(Ui)y...,
«+ -tUi (Vi) xtVi(Vi)y.... (28)

The difference form of Equations 28 was derived by replacing the par-
tial differentials with spatially centered differences as follows:

(Ui) (nymt.5) = [Ui(n,m+l.)-Ui(n,m)]/22,

Rv

(Ui) (n,m+.5) (Ui (n-.5,m+.5)-Ui (nt+.5,m+.5)]/22,

(vi) (nt.5,m) = [Vi (n+.5,m+.5)-Vi(n+.5,m-.5)]/22,

(vi) (n+. 5,m) = [Vi(n,m)-Vi(ntl.,m)]/22. (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) (nsm+.5) > [Ui(n,m+1.5)-Ui(n,m-.5)]/42,

(Ui) | (nym+.5) + [Ui(n-1.,m+.5)-Ui (n+l. ,m+.5) 1/42,

(vi) (nt.5,m) > [Vi(n+.5,m+1.)-Vi(n+.5,m-1.)]/42,

(vi) (n+ .5.m) > [Vi(n-.5,m)-Vi(n+1.5,m)]/42. (30)

Returning to staggered grid notation resulted in the following differ-
ence approximations:

(Ui) CD) > {Ui (MP) -Ui (MM) 1/42,
(vi) > [Ui (NM) -Ui (NP)]/42,
(vi) (D> (vi (MP) -Vi (MM) 1/42,
i) > [Vi (NM) -Vi(NP)]/42. (31)
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 21, yielded the difference forms presented in Equations

21


Angle of Rotation
Counter-Clockwise
from North

E,90

-U,180

W,270 | ~~

V,270 7
W270 S180

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<n<NE and the column index is m:l<m<ME. In the optimized model the con-
vention is retained; however, a single index (I)=(n,m) is defined and is
computed using FORTRAN array subscript convention, as I=(m-l1)*NEt+n. (Col-
umns are stored by consecutive elements; arrays are stored by consecutive
columns.) The following additional indexing abbreviations are used:

MM --index (n,m-1)=I-NE=m minus

MP (n,m+1)=I+NE=m plus
NM (n-1,m)=I-1,n minus
NP (n+1,m)=I+tl=n plus
NPMM (n+1,m-1)=I+1-NE=n plus, m minus
NMMP (n-1,m+1)=I-1+NE=n minus, m plus

Finite Difference Form of the Equations--

The finite difference form of Equations 1-3 that follow are essentially
the same as in W. Hansen's single layer model except the terms dictated by
the presence of k layers are added and the convective terms usually omitted
in the basic model have been included. Symbols are as previously defined
except as noted. The k layer model is:

ttt _ ttt _ t-T t —pyt7t t
oT) = Bo (1) (1/2) Fak 5 (2)0, (2) -Hws* UF Cm)
t-T t t-T t
- } ;
+HVS (MM) V5 (NM) Hv, My (I) 3427) 5 (4)


A304
pL, 2, ffs
ryt = 2> Cus uy!
Nel UIE : if
“Tv v v4
' . wv
ld]. ou Zz uf). u
2 |v ; 1
La v iY
3,4 : H 5
ref 09 Ga mu eh a =)
NE=3 t : ie
WB 4 Vv Vv
Hs
U

Figure 5. ¢ and Tt computational stencils at open boundaries.

x .
Figure 6. U, V computational stencils at open boundaries.

13


instability. Lamb (13) applied a successive approximation analysis to the
simplified non-linear differential system:

G+v(U)x = -9(0)x,
B= -[(hte)Ulx. (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

and 3 as +AiV-UL and +aiV-vi, respectively, where v is the two dimensional

Laplacian operator: 3°/ ax-43°/ay?. 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

Aiv-ti. The coefficients of these "artificial viscosity and diffusivity"

terms become in the difference form (1-ai) 2772, 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 (0,¥) , 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 (I), 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 () in a differential form are:

Condition 1: ven | Ty = 0, (11)


5 and 6 as follows:
*
oe 7-5 (1/2) {Ui (1) [U4 (MP) -U4 (MM) ]+Vi (I) [Ui (NM) -Ui(NP)]}...,
*
2 eee 5 (1/2) {UL (T) [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 (42). 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-1
C7
Zu
maf
v
aS
Z-u-
N [ z N
@

SAE @ 2 ful za)
vw f Pf Se
v a ee FS¥S

7
Z-u-.z-u-z--z-u
we [ ; ;
Vv Vv Sov Vv

Figure 10. U, V non-linear computational stencils at closed boundaries.

22


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:
eS a 4 t HG U) MARYS) Hay, (1)
x Momentum:

6.40, (U,) V7, (0,). -£v, 4, (u24v) “Sut,
iviei’x ivivy  iiviii i

1

°.
VC) ra (5 EX, (2)


21.

22.

23.

Maier-Reimer, E. Numerical Treatment of Horizontal Diffusion and
Transport Phenomena in Marine Basins of Large Size, Instittit fur
Meereskunde der Universit&t Hamburg, Hamburg, West Germany, as pre-
sented at IAMPA/IAPSO Assembly, Melbourne, Australia, January, 1974.

Thompson, R. Numerical Calculation of Turbulent Diffusion. Quart. J.
R. Met. Soc. 97(411) :93-98, 1971.

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


Number
21
22

23

Number

Post-reopening analysis for TFC cell (6, 13). ....
Maximum seaward advance of tidal flat F234. .....
TFC continuous source "plumes"... . . 2. ee eee
San Onofre outfall area grid bathymetry .......
San Onofre outfall tidal and wind prescription August

Infrared flight 0958-1107 compared to values after 15
hours. 2 ee ee ee ee ee ee ee

Infrared flight 1226-1344 compared to values after 17
hours. 6 6 6 ee ee ee ee ee ee

Infrared flight 1503-1617 compared to values after 19

hours. 2 ee ee ee ee ee
TABLES
Boundary Condition Confusion. .......-....

Typical Tidal Flat Grid and Run Specifications. .. .
Flat F3 Formation; Non-linear vs Linear. ......
Flat F4 Formation; Non-linear vs Linear .......-
Flat Fl Formation; Non-linear vs Linear .......

Flat F2 Formation; Non-linear vs Linear. ......

2, 1972

source

source

source

Page
42
43
45
50

51

53

53

54

36
37
38

39


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


TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)

11, REPORT NO. 2. 3. RECIPIENT'S ACCESSION NO.
EPA-600/3-78-073
4. TITLE AND SUBTITLE 5. REPORT DATE
July 1978

Enhanced Hydrodynamical-Numerical Model

6. PERFORMING ORGANIZATION CODE
For Near-Shore Processes

7. AUTHOR(S) IB. PERFORMING ORGANIZATION REPORT NO.
R.A. Bauer and A.D. Stroud

I9. PERFORMING ORGANIZATION NAME AND ADDRESS 10. PROGRAM ELEMENT NO.
EHE 625
Compass Systems, Inc. 5
San Diego, California 92109 11. CONTRACT/GRANT NO.
68-03-2225
12, SPONSORING AGENCY NAME AND ADDRESS 13. TYPE OF REPORT AND PERIOD COVERED
Environmental Research Laboratory-Ccrvallis final __1975~-1977

* 14, SPONSORING AGENCY CODE
Office of Research and Development

U.S. Environmental Protection Agency

EPA/600/02
Corvallis, Oregon 97330

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-

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.

7. KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS b.IDENTIFIERS/OPEN ENDED TERMS |c. CosaTi Field/Group
numerical methods marshy regions 08/c
thermal pollution operations research 12/A,B
monte carlo methods
mathematical methods
near-shore
river inflow
titdal flat
118. DISTRIBUTION STATEMENT. 19. SECURITY CLASS (This Report) 21.NO. OF PAGES
Unclassified 70
20. SECURITY CLASS (This page) 22. PRICE
Release to Public Unclassified

EPA Form 2220-1 (Rev. 4-77)
# U.S. GOVERNMENT PRINTING OFFICE: 1978-797-713/228 REGION 10

55


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 aCc/dt=-v*VC, holds where C is any
conservative quantity, as long as V isa 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
mM =c=f  C(x)dax, (17)
°
~.5
+.5
mM =F =f C(x)xdx/c, (18)
=.5
2 +5 2
M, = R°/12 = J Cx) (x-F) “dx/c. (19)
=.5

A more detailed description of the Pedersen-Prahm advection technique may be
found in (16).

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-

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.


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.


y Momentum:

¥.40, (V,)_4V, (V,) +60, +r, (02402) “PHT V,
1 1 1x 1 ly 1 a 1 1 1 1
Pan. Py-1
+g(i- — (oi) ta 7 (iD yr¥y (3)
i i
where:
c --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,
XY --horizontal space coordinate variables,
X,Y --horizontal components in the x,y direction, respectively, of
the external forces at the air-sea interface,
£ coriolis parameter,
x friction coefficient,
g acceleration of gravity,
0 density of the medium,
a, effect of external sources such as river inflow and rainfall

on layer surface deviation;
the following notational conventions apply:

we represents the operator 3/a3t (superscript dot),
( )x represents the operator 3/2x,
Cy --represents the operator 9/dy;
and the following conditions at upper or lower fluid boundaries apply:

B47 when layer i intersects bottom topography,

(051) x= (G51) YFP 5 =O” 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 ¢ 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 2 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.


Number

10
11
12
13
14
15
16
17
18
19

20

FIGURES

Scheme of the numerical grid. .....--

Orientation of the numerical grid .. .-

* = *
U, G (other than ~), and V computational stencils at closed

boundaries... . 6 eee ee eee

U, V computational stencils at open boundaries. ......

ct and zc computational stencils at open boundaries .... .
* *

U, V computational stencils at open boundaries. .....- -
U, Vv computational stencils at open boundaries. ..... .

Fractional z index notation. .....

U, V non-linear computational stencils at
U, V non-linear computational stencils at
The (n,m) cell. . 2. - - - ee ee ee ee
hy depth profile for coastal grid... .
Tidal flat grid... .. 2s. ee ee ee
tidal flat grid bathymetry (depths in cm)
Closure configuration TFC: t=17700 sec. .
Closure configuration TFC: t=18090 sec. .
Closure configuration TFC: t=18240 sec. .
Closure configuration TFC: t=18390 sec. .

Closure analysis for TFC cell (6, 13)

viii

Ge U, v computational stencils at closed boundaries .

closed boundaries

open boundaries .

12

12

13

13

14

15

20

22

23

24

31

32

33

40

40

41

41

42


TT

TABLE 1: BOUNDARY CONDITION CONFUSION

BOUNDARY BOUND [DIFFERENCE COMPUTATION | BOUNDARY CONDITION TYPE; | DIFFERENCE METHOD OF EFFECT-
CONDITION | TYPE ORDER OF DIRECTLY FRICTIONAL EFFECT ING BOUNDARY CONDITION
Eg. |Fig. APPROX. AFFECTED
3,8 closed ist ¢,0,7,0,¥ Solid wall, zero normal Zero normal velocity compon-
(11) non-linear | flow; infinite normal ent effects land-sea closure
10 Ty 2nd terms friction
4 open External forcing Tidal harmonic/permanent
(2) Tp ist Hu,lv, flow/river ; inflow/"overtide"
UW, . specification —
5 ur Internal dynamic Values not specified comput-
ed dynamically under Condi-
tions 3,4,5
5,6,7| open ist OV,0,0,6 Reflective condition Index modification effects
(13) T2 nonlinear with phase shift due to | value assignment such that
il 2nd terms light (1st order) or same computational form may
open moderate (2nd order) be used for both internal
(14)] 5 T2 lst id friction effecting and boundary computations
damping of gradient
open normal to indicated Direct value modification al-
5 rT _ boundary lows the same computational
(1s) sizseq lst 5 form to be used for both in-
8 T eae ternal and boundary computa-
1 Reflective condition tions
losedj with phase shift due to A combination of value and
8 ctose light friction effect- . es .
(16) Ty ist ay ing damping of gradient index modification allows the
7 open ' normal to indicated same form to be used for both
T internal and boundary compu-
2 oundary tations
11 open Perfect transverse re- Extrapolation of internal
(33) qT nonlinear flection, total slip, transverse gradient (doubl-
10 | closed) 2nd terms zero friction. ing of first order half dif-
Ty ference computed due to in-
dex modification.)



United States Environmental Research EPA-600/3-78-073
Environmental Protection Laboratory July 1978
Agency Corvallis OR 97330

Research and Development

SEPA Enhanced Hydrodynamical-
Numerical Model for
Near-Shore Processes

yy SS


TABLE 2.

TYPICAL TIDAL FLAT GRID AND RUN SPECIFICATIONS

DATA PROGRAM VARIABLE
CARD VARIABLE DESCRIPTION VARIABLE VALUE ASSIGNMENT
NAME _|1C1(2)-(9) Grid description EPA TIDAL FLAT TEST GRID
GRID NE No. of Rows 15
ME No. of Columns 20
DL (2k) Grid length 1/3 Km
ROTANG Rotation angle 180° (implies (1,1) is the SE
corner of grid)
ALAT Average latitude 45° (used to compute f and g)
R(r) Friction coefficient +003
LAYS (k) No. of layers 1
HL(1) Maximum depth of layer 660 cm
ALPHA (a) Smoothing parameter 99
RHO (p) Density 1.022
HMIN Minimum allowable depth | 1.0
TIME MODOUT Data save step 120 sec
ITO Data save start 120 sec (for runs saving data)
ITS (to) Start computation 0 sec
ITE End computation 31200 sec (8.67 hrs)
ITINC (21) Computational step 30 sec
IPO Print step 1200 sec
IPMOD Print start 1200 sec
FACT | TIDSPD(1) Constituent speed 45°/hr (false "tide"-8 hr period)
TIDE* |ITIDE(1),(2) | Min/max row spec. 1,15
ITIDE (3), (4) | Min/max column spec. 20 (western boundary specifica-
tion)
ITIDE (5) Tidal offset time 0 sec
ITIDE (1) Constituent phase shift | 90°
ITIDE (5) Constituent amplitude 130 _cm
COMP |MODEL Grid I.D. No. 998

*The open portions of the northern and southern boundaries were
computed dynamically based on internal information and boundary

conditions as specified previously.

34



=
for

= 2
w

N

OQ Voss

tre g
z

DEPP Sa

uv
ulP],v[P]
U,V

ea al
<+ <|

--zeroth order moment of conservative quantity
--first order moment of conservative quantity

--second order moment of conservative quantity

--a unit vector normal to the indicated boundary

--stochastic variate

--population of stochastic variates

--effect of external sources on layer surface deviation

--net heat flux

--net atmospheric-oceanic (long wave) radiative heat flux

~--evaporative/condensative heat flux

--net atmospheric-oceanic sensible heat flux (convective trans-
fer or conduction )

--heat flux resulting from human efforts

--reflected (short wave) radiative heat flux

--total incident

--total clear sky incident (short wave) radiative heat flux

--friction coefficient

--width of a conservative quantity or a random value between zero
and one

--x,y widths of a conservative quantity

--distance variable

--stochastic coefficient

--time variable

--length of a day (sunrise to sunset in minutes)

--initial time (start of a computation )
--temperature variable

--unit vector in the tangential direction

--surface air temperature

--mean surface water temperature for a cell

--indicators of points of computation for U and V, respectively

--u,v turbulence component assigned the value of p

--vertically integrated velocity components for the x,y direc-
tion, respectively

--the "smoothed" U and V components, respectively (result of
eddy viscosity)

--linearly interpolated u,v velocity components at the u,v
points, respectively

--relative humidity (%)

--wind speed (m/sec)

--horizontal space coordinate variables

--horizontal components in the x,y direction, respectively, of
external forces including the air-sea interface

--indicator of points of computation for ¢

xi


LIST OF ABBREVIATIONS AND SYMBOLS

ABBREVIATIONS

F1-F4 --bathymetric features identified by notes 1 - 4 on Figure 14;
also tidal flat segments associated with bathymetric features

F1234 --tidal flat associated with combined features F1-F4

F234 --tidal flat associated with combined features F2-F4

F34 --tidal flat associated with combined features F3 and F4

HN --Hydrodynamical-Numerical

I --the array index (n,m)=(m-1)NEtn

m,M ~-the column index (m=1,2,...,ME-1,ME)

MM --the array index (n,m-1)=I-NE(m minus)

MP --the array index (n,m+1)=I+NE(m plus)

n,N --the row index (n=1,2,...,NE-1,NE)

NE --number of rows for array variables

NM --the array index (n-1,m)=I-1(n minus)

NMMP --the array index (n-1,m+1)=I+NE-1(n minus, m plus)

NP --the array index (n+1,m)=I+1(n plus)

NPMM the array index (n+l1,m-1)=I-NE+1(n minus, m plus)

TFC --Tidal Flat Convective (nonlinear tidal flat case)

A --area variable

An --noon altitude of the sun in degrees

Aw --ambient water temperature

Cc --cloudiness factor in tenths of the sky
c
Cc
a

--mean specific heat
--conservative quantity (concentration or volume)
--depth (thickness) of the heated fluid

Ea --moist vapor pressure of air

Ew --Saturation vapor pressure of air

£ --coriolis parameter

F --center of concentration (one dimensional case)
Fx,Fy --x,y coordinates of center of conservative quantity
g --acceleration of gravity

h --initial layer thickness

H --dynamic layer thickness

H --mean thickness at z point

i --layer indicator (i=1,k)

j -~general index

k --number of layers in the HN model

K --new heat quantity (added heat)


10.

ll.

12.

13.

14,

15.

16.

17.

18.

19.

20.

of a series of four reports. ENVPREDRSCHFAC Technical Note No. 4-74,
Environmental Prediction Research Facility, Monterey, California, 1974.
79 pp.

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.

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.

Laevastu, T. and A. D. Stroud. Exchange between Coastal and Continental
Shelf Waters as Computed with a Two-Layer HN Model, Manuscript Revort,
ENVPREDRSCHFAC Report to the Environmental Protection Agency, Corvallis
Environmental Research Laboratory, Environmental Prediction Research
Facility, Monterey, California, 1975.

Lamb, H. Hydrodynamics. DOVER Publications, New York, New York, 1945.
738 pp.

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.

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.

Pedersen, L. B. and L. P. Prahm. A Method for Numerical Solution of
the Advection Equation 3C/3t=-V+VC, Meteorological Institute of Den-

mark, as submitted to TELLUS, At gust, 1973. 36 pp.

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.

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.

Laevastu, T. Factors Affecting the Temperature of the Surface Layer
of the Sea, Soc. SCIENT, FENNICA, Comment, Physico-Mathem., Helsinki,
1960. 136 pp.

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


The Stochastic Process

The particular stochastic process employed in this Monte Carlo scheme
is defined as follows:
. = S(R,-.5)
P5 5 ,
where Ps 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


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 ina
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.-1.86Tw-. 95) (1-.0765c) (25)
where Tw is the water temperature, U5 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-Fa) (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


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 oy down to the point of closure

of cell A. A similar character is observed in the profile of BD 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


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+Om, (20)
where:
Q --net heat flux,
Qs --total incident (short wave) radiative heat flux,
Or --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),
fen --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

Each of the natural components of Q is presented in cal om? 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 = Q0s(1.-.0006¢°), (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:
=. 22
Qos O14A tyr (22)
where AL is the noon altitude of the sun (degrees) and ta is the length of

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.sina : (23)

17


m
13 5 7 9 NBS 7 W NM 2 2 27 Ww 3H 33 35 37 39

3 > > > > > » > > > »

5

7 > > > > > > > > > >
n?

Contours: 22.5°C (72.5°F) and 28.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


Foreword . - - - - e+ ee ee et ee :
Preface. . 6 6 6 ee ee ee ee ee .
Abstract 2. 6 6 ee ee ee ee .
Figures. . . 2 6 ee ee ee ee ee .
Tables . . ee ee ee ee ee ee ee -
Abbreviations and Symbols. . ......- .
Acknowledgment ©... 6 6 ee ee ee ee .
1. Introduction... . +. eee .
2. Conclusions and Recommendations. :
3. Numerical Methods and Techniques .
4. Experimental Extensions of the HN Model.
5. Results and Discussion... . - .
References - - - - - e+e ee ee ee ee .
Appendix
A. Verification of the Model. .. - :

CONTENTS

vii

49


ott py = TEC fare, smut (2p) pot (1) 2H" (2) 79 doe (a) #20 (TD
1 1 1 1 1 1 1
t t t t t et

= 25 (0/8) { (0, (MP) U5 (MM) JU S(T) + (05 (xm) -U (NP) ]V5 (1) }

og (0/2) ((1-(0,_4/04)N [op MP) -25 (I

t+T

ttt
gap) l}4etxy” (Ds (5)

ttt
+(0,_/0;) (55) (MP)-%

_ ot ttt t 1.2 et). 2).5,.t 4 et
(1) = vy-{ (tr ,/ivs (1) IV(2) ws (1) Ww, (1) -21£0, (7)

t t t t t at
~-5(t/2) CV, (NM) -V (NP) Iv, (1) +1, (MP) -V, (MM) JU; (1) }

-a(t/2{[1-(0, 4/0, ) eg (05 Ne)

t+T t+T t+T
+6; 1/05) (05 4 (D)-o 5 (NP) } }42T¥ (I); (6)
where:
, U(qI) = aU(I)+.25(1-a) [U (NM) +U (NP) +U (MM) +U (MP) 1],
V(I) = aV(I)+.25(1-a) [V (NM) +V (NP) +V (MM) +V (MP) 1,
U(I) = af (1)+.25 (1a) [ (NM) +5 (NP) +0 (MM) +0 (MP) }; (7)
*
U(I) = .25[U (I) +U (NP) +U (MM) +U (NPMM) ] ,
vn) = .25[V(I)+V (NM) +V (MP) +V (NMMP) 1]; (8)
ttt t+T tt+T t+T t+T
muy (1) = hu,+.5{(tp 0 (1) +e; (mP)]-(05 7, (1) +0, 5 (MP) 1},
ttt = ttt ttt eprttt t+T .
Hv, (1) = bvy+.5€06p0 (1) 40,0 (NP) tory (+554, (NP): (9)
where:
h --initial layer thickness,
t --time variable,
u,v --indicators of U and V computational points, respectively,
& --half grid step,
a --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


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.

2 4 6 8 10 2 “ 16 18 20
J T= 31200

CONSTANTS:

12

~ Monte Carlo Dispersion Plume
oO Thermal Advection Plume
@ Monte Carlo Source (30 setond step)
aa) 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


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


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
(S) 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;
2) partial closure: recomputation of ot,
3) full closure: reset ttt to the value of 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 a requires:

1) xrecomputation (partial closure);
2) non-zero "adjustment" (from the i+l layer);

T T

3) resetting of tf*" to the value of cf" (full closure);

4) a combination of 1) and 2); or
5) a combination of 2) and 3);
then oy 1 (-140) requires an adjustment aye The adjustment of 25 is equal to

T

5 t+ ttt
the recomputed or reset and/or adjusted oF less the value of oy . Thus

the recomputation procedures must occur in precisely the same order that the
original ~ computations took place, i.e., lower layer oy terms are computed

first so that the Zi 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:
ttt t+T

20. = G. Sart
1) open ae a + Zin
. ttt _ tt+T ay
2) partial closure: i = oy (recomp) + Ziyi
ttt t-T -
2 C. = . +Z7 4.
3) full closure uF aS 41

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 ¢ as it occurs. Second order computations in ¢ (2)

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 T aia
not significantly alter or improve the results, ergo Tt is not recomputed.

Second order computations in U and V Ww, 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


m
12°43 nuuu swe view»

Tavares rs
FACT ~7 18240

96.6

ee ee

Figure 18. Closure configuration TFC: t = 18240 sec.

je A uo2 3 4 1 6 7 Bw 19 20

Tv x) te1a90

Leen LAL.

Figure 19. Closure configuration TFC: t = 18390 sec.
41



3
1, 2
Paes tis —s
. hf
ul. z+ fal. {2}. ou
. Uf
@& a =)
: Vv
to | -u
u -
\ Vv ly|3
6
U

Figure 7. U, v computational stencils at open boundaries.

14


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 t, 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


m
Vi dece See Jene Senet 19 S19 2h 2301250270273) 9335, 37D

7 — > >
n? 71.

n cS" > — > —»

3 NS Ya SIN

x} >» >» 3» 3 @- res

Contours: 21.7°C (71.0°F) and 27.2°C (81.0°F). Ambient: 20.7°C (69.2°F).

infrared 2 000 -----"- computed — > 50 cm/sec

Figure A-3. Infrared flight 0958-1107 compared to values after 15 source hours.

m
L305 7 9 1138171920 232827 2931333537,
3 > —> > > >
5
7 _—-— orc o—TMVm" >
9

yo = aan OO? eee

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


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.

ii


6 7 8 9 0 1 2 «#13 4 «1 1% 17 8 OW DO

Figure 15.

Tidal flat grid bathymetry (depths in cm).

33


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

Mf

DEPTH (cm)
8

Y

Yy,
WM Uf

123 4 5 6 7 8 9 10 11 12:13 4 15 1 17 18 19 20

COASTAL GRID ROW NUMBER (n)
MSL -- mean sea level

Figure 13. hy depth profile for coastal grid.

31


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.


Figure 21 presents a much later period of time than Figure 20. Cells A,
B, 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.

m
12 3 4 5 6 7 8 9 © 1) 122 13:4 «15 % 17 18 19 2

a es we ie ee te22022

3

J
ee ee ed

Elta aath esata teal head toedteandttes 7

Figure 22. Maximum seaward advance of tidal flat F234.

43


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


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.


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


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-1 M M+1
ee ee

+2 grid structure

+---(n,m) cell boundary

Figure 12. The (n,m) cell.

24


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-1 M M+2

* = *
Figure 3. U, ¢ (other than ¢), and V
computational stencils at closed boundaries.

ME=4
pam
~o~ fe -u- GSE:

Nel
-u

2
-u

3
_u

Figure 4. U, V computational stencils at open boundaries.

12


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
oianned to foster technology transfer and a maximum interface in related fields
The nine series are

Environmental Health Effects Research

Environmental Protection Technology

Ecological Research

Environmental Monitoring

Socioeconomic Environmental Studies

Scientific and Technical Assessment Reports (STAR)
Interagency Energy-Environment Research and Development
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

mom

DNONERW

This document is available to the public through the National Technical Informa-
tion Service. Springfield, Virginia 22161.


TABLE 3:

FLAT F3 FORMATION; NON-LINEAR VS LINEAR

FIAT CELL NON-LINEAR LINEAR
APPEARANCE | (nym) (am)
orprnat | Location | rmme | conricuratron | tocarton | TIME | CONFIGURATION
1 (6,11) | 16980 oF (6,11) | 17040 ot
2 (6,12) | 17280 OH (6,12) | 17430 ott
3 (7,11) | 17730 oft (7,11) | 17760 ott
+ +
4 (5,11) | 17760 of+ (5,11) | 17940 of+
+ +
5 (6,10) | 18210 ot+ (6,10) | 18030 of
+ +
6 (6,13) | 18240 of++ (6,13) | 18300 ote
7 (7,12) | 18390 ett (7,12) 18420 ttt
8 (5,10) | 18480 dy (5,10) | 18420 oh
9 (5,12) | 18870 othe (5,12) | 19140 oH
10 (4,11) | 19050 rseay (8,11) | 192604 oe
£
11 (8,11) | 193204 (4,11) | 19410 at.
4 Bs

* Connects to Flat F4 to form F34 at __.

© Represents shallow knob at (6,10) at initial condition

(Bathymetric feature F3, Figure 13.)

36



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--

fo 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(1,1) is represented as u(1,1.5) and staggered v(1,1) becomes
v(1.5,1). Thus fractional z notation is a positionally correct notation pro-
viding the flexibility to represent variables at any desired location.

m11.5 2 2.5 3 3.5

n=1 Zz se Zee

1.5 v ° v o vio Z-z point, u-u point, v-v point

2 : : i . o - values required must be interpolated
2.5

3 .

3.5

Figure 9. Fractional z index notation.

20


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 damagé 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:

5 Taxi coy as [P])-2t,

t+2t_ ot t+2T +27

. =y.t(V_ tv,

J 3 3 3
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:

[P])*2t,

mj 7 =my+ cos ust [P])*/2,
ny Pens wy ays [P])*/2,

where the U[P],VI[P] are selected from P independently in order to describe a
constant diffusion coefficient.

29


M-2 M-1 M M+ M+2
rm 7

! 1 f 1 f 1
— Zz u z u Zz us Zz u Z-u
N-2 4
Lov a v ADS N Vv v
rz u Zz u Zz u z u
N-1 64
— Vv 1c v v
N lv
op 6
- z u Zz a z u z ‘ul - 2 u
N : . -
— Vv yy Vv v a
U
-@- u-z u Zz-@- Zz wf z at
NeL : AS
Lv Vv v.
rz: ou zjo-u- 6.
W Ue zs
Nt2 : :
Lv v v

Figure 8. Tu, U, v computational stencils at closed boundaries.

Condition 1 at closed bounds--By definition of the grid structure, land-
sea boundaries are associated with positions of zero depth. Velocity compon-
ents associated with zero depth values are initialized and remain equal to
zero. As a result, Equation 11 is satisfied in the difference scheme without
modification of the derived difference forms adjacent to normal land-sea
boundaries: the zero velocity components on land-sea boundaries are simply

—* ot
used for computations of ¢,U,V,V,U as required. (Refer to positions identi-
fied by superscript 1 in Figures 3 and 8.)

Condition 2 at open bounds--Equation 12 is satisfied in the difference
scheme by providing ¢ values from external information, where available, or
by computing a dynamic ¢ value constrained by conditions 3-5. (Refer to
positions identified by superscript 2 in Figure 4.)

15


NDA PM ORI R

ent

area

--smoothing parameter

--the mean altitude of the sun for a short period
~-directional parameter

--directional parameter

--half grid length

--half time step

--density of fluid (assumed uniform through layer)
--deviation of the top of the layer from an arbitrary level

surface
--smoothed value

of C(result of eddy diffusivity)

--vector velocity (resultant of U and V)

--vector "dot" product

--superscript dot represents the operator 3/9t

--interpreted as

“evaluated on"

--an element of the normal to the indicated boundary

--represents the
--represents the

--represents the

--the difference
“approximately

--represents the
represent unit
tively)

operator 3/3x
operator 3/dy

operator 97 /ax743°/ay-
analog to 2
equal to"

vector operator id/dx+jd/dy (where i and j

normal vectors in the x,y direction,

xii

respec-


Tide at column 23—s,
. T POT Hours.
12 iB i
i
g
7 a tie
Direction - sw sow -
Measurement [ sad (mph) 0 : 10-15 is io
Speed used 0 i 125 is io
f T T < ss Mode! Time Sec.
98848 0964s = 20848131248 142088 1548

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


TABLE 4:

FLAT F4 FORMATION; NON~LINEAR VS LINEAR

FLAT CELL NON-LINEAR LINEAR
APPEARANCE | (n,m) (n,m)
ORDINAL LOCATION | TIME | ConFIcuRATION | LocaTIon | TIME | CONFIGURATION
rT (11,10) | 17100 oF (11,10) | 17130 oF
2 (11,11) | 17340 Ott (11,11) | 17400 Ott
3 (12,10) | 17640 of+ (12,10) | 17760 oft
4 (10,10) | 18060 ott (10,10) | 18060 of+
5 (11,12) | 18180 oft (11,9) | 18120 eft
6 (11,9) | 18300 of++ (11,12) | 18210 ote
7 (12,9) | 18510 gi (12,9) | 18340 gin
8 (12,11) | 18570 ett (10,11) | 18660 ett
+4 7
9 (10,11) | 18600 ett+ (12,11) | 18720 ett+
4
lo (13,10) | 18840 tts (13,10) | 18990 @
+ ¢
ty +
ll (9,10) | 19050 gtt+ (9,10) | 19020 t+
+
(no q rs
-- 193204 sarge) Sti -- 192604 (no @ft+
hange)

* Plat F3 connects to F4 by accretion to form Flat F34 at "—",
© Represents shallow knob at (11,9) at initial condition

(Bathymetric feature F4, Figure 13.)

37



TABLE 6:

FLAT F2 FORMATION; NON-LINEAR VS LINEAR

PLAT CELL NON-LINEAR LINEAR
APPEARANCE (n,m) (n,m)
orprnaz | Locatton| Time | conricuration | LocaTron | TIME | CONFIGURATION
1 (11,6) | 18180 ® (11,6) | 18000 ®
2 (12,6) | 19350 g (12,6) | 19050 g
3 (11,5) | 19590 +9 (11,5) | 19260 +9
4 (10,6) | 20670 +$ (10,6) | 20310 +§
5 (11,7) | 20910 ++ (11,7) | 20460 +b+
6 (11,8) | 209704 the (11,8) | 20550* ten
7 (13,6) | 21810 agen (13,6) | 21270 then
+ +
8 (12,7) | 21870 +44 (12,7) | 21330 +40
+ +
9 (12,8) | 21990 “ben (2,8) | 21420 Bey

* Flat F2 connects by accretion to Flat 34 to form lazy "L" shaped
tidal flat F234 at "I".
© Represents shallow knob at (11,6) at initial condition

(Bathymetric feature F2, Figure 13.)

39



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