mm

Hi

AN INVESTIGATION OF LINEAR TRANSIENTS ASSO-
CIATED WITH A TIME DEPENDENT BOTTOM SPIRAL

Norman Thomas Camp

NAVAL POSTGRADUATE SCHOOL

Monterey, California

THESIS

»

AN

INVESTIGATION OF LINEAR

TRANSIENTS

ASSOCIATED

WITH A

TIME DEPENDENT

BOTTOM

SPIRAL

by

Norman Thoma

s Camp

Thesis

Advisor

.

J

.A.

Gait

March 1972

Approved faoh puhtlc h.eA.<La£>a; dti&ubuution antariLted.

An Investigation of Linear Transients Associated
With a Time Dependent Bottom Spiral

by

Norman Thomas Camp
Lieutenant Commander , United States Navy
Ed.B., Rhode Island College, 1962

Submitted in Partial fulfillment of the
requirements for the degree of

MASTER OF SCIENCE IN OCEANOGRAPHY

from the

NAVAL POSTGRADUATE SCHOOL
March 1972

ABSTRACT

Ekman's linear equations for time dependent flow (neglecting
wind stress) are solved using a time dependent Green's function and
the method suggested by Welander (1957). The solution represents
the vertical velocity profile in terms of the local time history of
the changes in sea surface elevation determined by the divergence
of the flow in the vertically integrated continuity equation.

A fully implicit finite-difference scheme is developed to re-
present a time dependent seiche oscillating across a shallow infinite
channel. The transients associated with the formation of the bottom
spiral are clearly represented by the model and the influence of
friction and Coriolis are individually and collectively introduced.
The model allows independent calculation of velocity, volume trans-
port, sea surface elevation, bottom stress, and the total energy
balance of the system. The numerical scheme provides a method for
adequately describing and investigating certain classes of time
dependent motion, and its development suggests a mechanism for
improving numerical prediction of storm surge.

TABLE OF CONTENTS

I. INTRODUCTION 8

A. BACKGROUND 8

B. PREVIOUS RELATED RESEARCH 9

C. THESIS OBJECTIVES 10

II. THE FORMAL SOLUTION 12

A. GOVERNING EQUATIONS AND ASSUMPTIONS 12

B. METHOD OF SOLUTION 14

C. DISCUSSION OF THE FORMAL SOLUTION 15

III. THE NUMERICAL SOLUTION 17

A. FORMULATION OF THE MODEL 17

1. Establishment of a Representative

Cross Section 17

2. Establishment of Governing Non-dimensional
Ratio 19

3. Adaptation of the Formal Solution to the
Model 21

4. Initial Conditions 22

5. Boundary Conditions 22

6. The Step-by-step Scheme 22

7. The Implicit Finite-Difference Scheme 23

8. Numerical Treatment of CI and C2 28

9. Velocity Calculations 30

10. Bottom Stress Calculations 30

11. Energy Calculations 30

12. Scaling of Model Output 31

3

B. TESTING THE MODEL 32

1. Establishing a Reference Solution 32

2. Accuracy and Stability of the Model 33

C. MODEL LIMITATIONS 39

IV. PRESENTATION OF RESULTS 4l

A. TEST CASES AND GEOPHYSICAL SIMULATIONS 4l

1. Description of Test Cases --4l

2. Description of Geophysical Simulations 4l

B. DISCUSSION OF TEST RUN RESULTS 4l

1. Frictional Effects on the Flow 4l

/

2. Coriolis Effects on the Flow 51

3. Combined Effects of Friction and Coriolis 56

C. DISCUSSION OF THE GEOPHYSICAL SIMULATIONS 66

V. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER WORK 69

APPENDIX A: DETAILS OF THE FINITE- DIFFERENCE SCHEME 72

APPENDIX B: DETAILS OF REPRESENTING CI , C2 , VELOCITY

AND BOTTOM STRESS 77

APPENDIX C: COMPUTER PROGRAM DESCRIPTION 8l

COMPUTER PROGRAM 84

BIBLIOGRAPHY 90

INITIAL DISTRIBUTION LIST 91

FORM DD 1473 92

LIST OF TABLES

I. Parameters of Test Cases (Runs 1 - 25) 42

II. Parameters of Geophysical Simulations 43

LIST OF FIGURES

1. Pictorial Depiction of Model Basin 18

2. Comparison of Free Surface Elevation for Model and
Analytic Solution (PCF=0) , PCD=0) 34

3. Comparison of Free Surface Elevation for Model and
Analytic Solution (PCF=1, PCD=0) 35

4. Comparison of Energy Balance For Model and Analytic
Solution (PCF=0, PCD=0) 36

5. Changes in Energy Balance as a Function of Time Step 37

6. Effects of Friction on the Free Surface (PCF=0) 45

7. Effects of Friction on Mid-Channel Transport (PCF=0) 46

8. Effects of Friction on Energy Balance (PCF=0) 47

9. Comparison of Bottom Stress and Mid-Channel Transport
(PCF=0, PCD=1) 48

10. Effects of Friction on Mid-Channel Velocity (PCF=0, PCD=2)-50

11. Effects of Inertial Rotation on Surface Elevation (PCD=0)--52

12. Effects of Inertial Rotation on Mid-Channel Transport
(PCD=0) -- - --53

13- Comparison of Bottom Stress and Mid-Channel Transport

(PCF=.5, PCD=1) 57

14. Comparison of Bottom Stress and Mid-Channel Transport
(PCF=1, PCD=1) 58

15. Effects of Friction and Inertial Rotation on the Free
Surface 59

16. Effects of Friction and Inertial Rotation on Energy
Balance 60

17. Hodograph of Mid-Channel Velocity (t=5T) 6l

18. Effects of Friction and Inertial Rotation on Mid-Channel
Velocity (PCF=1, PCD=1 ) 63

19. Effects of Friction and Inertial Rotation on Mid-Channel
Velocity (PCF=1, PCD=1, t=3T/4) 64

20. Effects of Friction and Inertial Rotation on Mid-Channel
Velocity (PCF=1, PCD=1, t=T) 65

ACKNOWLEDGEMENTS

The topic for this thesis was suggested by Assistant Professor
Jerry A. Gait of the Naval Postgraduate School, Monterey,
California. His valuable experience in the field of ocean dynamics
and his guidance in developing the conceptual model contributed
greatly to the successful completion of this research. Assistant
Professor R.H. Franke provided helpful ideas during the research
and contributed valuable time and thought reading the manuscript.
To these people the author expresses his sincere thanks. Addi-
tionally the author wishes to express his long overdue appreciation
to his wife, who provided a needed source of encouragement and a
sympathetic ear.

I. INTRODUCTION

"The purpose of computing is insight not numbers"

(Hamming, 1962)

A. BACKGROUND

Oceanographers use a number of techniques for the study of
oceanic flow including analysis of hydrographic data, use of
hydraulic laboratory models, and mathematical treatment of the
basic equations of motion. Initially oceanic circulation was
treated as a steady state problem, and various investigators have
been successful in representing and explaining many of the gross
features of oceanic circulation using each of these techniques.

Mathematical treatment of oceanic circulation began with
Ekman's analysis (1905) of wind-driven circulation and his results
have been extended by a number of authors. (Sverdrup, 1947, Munk,
1950, Bryan, 1959) For the most part the study of large scale
ocean dynamics has been done with steady state formulations.

Mathematical treatment of time dependent motion is far less
complete in the literature than that of steady state. It is also
far more difficult to compare results due to the scarcity of
reliable data. The complex nature of the solutions to time de-
pendent boundary value problems further increases this difficulty.
These solutions take the form of varying differential, integral,
and infinite series equations whose numerical evaluation was not
feasible until the advent of the high speed computer and the de-
velopment of appropriate numerical techniques.

In the last decade these methods have been applied, through
numerical models, to the problems of explaining and predicting the
various details of fluid motion in both the atmosphere and the
oceans. Although much success has been attained with these methods,
many critical questions still need explanation. Many of these
questions relate to the structure of the various frictional boundary
layers, the details of their development, and the resulting effects
on the interior flow.

Many complicated engineering problems occur in near shore and
other shallow water regions. The time dependent nature of these
areas is well documented and the resulting flow is greatly modified
due to interaction with the ocean bottom. Therefore an understanding
of the nature of time dependent flow at all levels is fundamental to
dealing with problems in these areas. The question of how readjust-
ment takes place and how the many transients can be appropriately
represented in the bottom frictional boundary layer has not been
adequately explained. Further, it is necessary to understand the
dynamic coupling between the flow and the development of the bottom
boundary layer before accurate modeling techniques can be developed.
It is towards this understanding that this work was directed.

B. PREVIOUS RELATED RESEARCH

Work by Welander (1957) dealt with obtaining time dependent
solutions using Ekman dynamics. His results were presented in
terms of a time dependent Green's function and were related to the
problems of predicting tides and storm surges. He demonstrated
that the velocity profile, volume transport, sea level elevation,
and bottom stress could be exactly represented by his

"integr o-dif ferent ial" solution. It was suggested that this
solution could be used to improve storm surge and tide predictions.
A step-by-step numerical scheme was proposed, but the details of
the scheme were not presented in the literature.

Gait (1970) covered much of the same material that Welander
presented but used a different method of solution. The solution
allowed a more explicit formulation of results which provided more
insight into the possible transients present and the relevant time
scales of each. He discussed the need for information relating to
time dependent flow in shallow water regions, and suggested that
the numerical investigation of the interaction between a seiche
and its associated transient bottom spiral might be useful in de-
termining how the various transients present are effectively
coupled.

C. THESIS OBJECTIVES

Since the major area of interest was to be the bottom boundary
layer, the first objective was to formally solve the governing
equations neglecting wind stress. The solution would therefore not
allow a surface Ekraan spiral to develop and effect the interior
flow. It will be showed that the solution can be represented by
two equations, one integral and one differential. These equations
are similar in form to the ones developed by Welander and Gait.
Although both of these authors proposed that these equations could
be solved numerically using step-by-step procedures, it was not
initially apparent whether a stable and accurate numerical scheme
was feasible. Therefore the second objective was to formulate a
numerical model capable of simultaneously solving these equations.

10

Finally, as suggested by Gait, this model would be applied to
represent a seiche to examine the relationships which developed
between the various transients.

11

II. THE FORMAL SOLUTION

A. GOVERNING EQUATIONS AND ASSUMPTIONS

The momentum equations considered are essentially those used
by Ekraan. The non-linear and horizontal friction terms are assumed
small and neglected. The Coriolis parameter, the density, and the
coefficient of vertical friction are all assumed constant. Thus:

2

bu - . b u b£ , , .

fv _ k_ _ . g _- (1)

Oz

__ + fu . k = . g ^ . (2)

Oz

Where u and v are velocities in the x and y-dir ect ions , f the
Coriolis parameter, g the acceleration of gravity, and k the
constant coefficient of vertical friction. ?(x,y,t) is the ele-
vation of the free surface and the z-axis is considered positive
up.

The third equation necessary to solve for u, v, and § is the
vertically integrated continuity equation.

|L + |u + bV =

bt bx by K3'

where

0 „ 0

U = \udz;V = \vdz (4)

-h "-h

is the volume transport and h = h(x,y) is the depth of the water.
§ is assumed small compared to h.

12

Neglecting the effects of wind stress on the motion results in
the following linearized boundary conditions on u and v:

4

<u>x=-h " ° » <">2=-h " ° ■ <6>

For this development it was assumed that initially a state

*
of rest exists:

(u)t=0 = o ; (v)t=0 = 0 . (7)

As in Welander ' s work the following complex notation was

introduced:

w = u + iv

(8)
bn bx by

Where i = V^T. Using this notation eqs . (1), (2), (6), and (7)
can be expressed as follows:

* Urn = - 9 H (9)

(10)

(ii)
<»>z=-h " ° • <12>

Equations (9) through (12) define the initial/boundary valued
problem that must be solved to specify u and v. Although (f~~) in
eq. (9) is not externally specified it can theoretically be deter-
mined by integrating eq. (3). This makes it possible to treat it
as a time dependent forcing term (i.e. «j*- (t)).

13

bw
bt "

k

*2

b w

x 2
oz

C&.N) = °

(W)t:

=0

= 0

B. METHOD OF SOLUTION

To solve the problem begin by multiplying eqs . (9) through

ift * ift

(12) by e ' and introducing the relation w = we . Then eqs.

(9) through (12) can be written as:

I * 2 *

bw . b w _ ift b? . .

_ _ k _ _ ge g- (t) (13)

QZ

<»*>t=o ■ ° <15»

(»*)z=-h - ° • (16>

A solution to eqs. (13) through (16) can be obtained by using

a time dependent Green's function approach as described below.

Assume

w* = 0 (t) Z (z) (17)

where Z (z) satisfies the boundary conditions and 0 (t) satisfies

mv mv '

the initial conditions. Z (z) is the eigen-f unction determined by

mv ' *

applying standard separation of variables to the homogeneous part
of eq. (13) using eqs. (14) and (16). It is found to be:

Z (z) = cos(pz) (18)

wher e . _ _ .

P- ^=111 (19)

0 (t) is determined by substituting eq. (17) into eq. (13)
requiring that 0 (t) satisfy eq. (15), and using the orthogonality
of the cosine function to determine the appropriate constants.

Thus

a (tJ = 23Lg£\ e"kP <t-f)e-iff 6L(t.)dt, . (20)

n»x ph J0 0nv '

14

Therefore a solution to eqs . (13) through ( 16 ) is

(x,y,,,t) -^ l=i i-'l'-IP) ^.Vt^'l.-^^tMdf (2!)

w

For the general time dependent solution of eq. (9), dividing
eq. (21) by e yields:

(x.y.z.t) =^J, (-1 >'<">■ (P£l fVl^UX*"*') &f )«■ (22)
v "' ' ' h m=l p J_ onv ' v '

w

C. DISCUSSION OF THE FORMAL SOLUTION

It is apparent at a glance that eq. (22) satisfies the boundary

and initial conditions. Substituting eq. (22) into eq. (9) yields

as a residue:

oo m

£ . -^~) — cos(pz) = - 1 . (23)

m=l ph vt- / \ ^/

It is easily verified that the left hand side of eq. (23) is
simply the cosine expansion of the right hand side, and therefore
eq. (22), in the limit, is an exact solution to the initial
problem.

Equation (22) represents the horizontal currents at all depths
resulting from the pressure gradient produced by the slope of the
sea surface at a particular point and instant in time. The integral
term carries a running summation of the currents resulting from the
slope in the past history of the flow. These currents are modi-
fied in time by a rotation and decay represented by the weighting

function

e-(kp2+if)(t-t').

This weighting function introduces two interesting time scales
into the problem. The first is the inert ial rotation period, which

15

varies as a function of latitude. For mid-latitude locations it is
on the order of 20 hours. The second time scale is introduced by
the exponential decay (T ) represented by

Td = ^ - ^2 2 ' W

d kp k(2m-l) tt

This is a function of depth and friction, but for typical values

2
(k = .01 m /sec, h = 50m) it is on the order of 28 hours for the

first term in the series to decay to 1/e of its initial value.
This is the slowest term to decay.

The motion represented by (22) depends on the surface slope
and thus new time scales are introduced into the solution. These
depend on the horizontal dimensions represented and enter through
the integration of the continuity equation (3). Time scales re-
lated to set up of the sea surface and seiche periods will govern
the time rate of change of the surface slope. This change
initiates the transients present in the geostrophic and bottom
currents, which show the inertial rotation and exponential decay
previously noted.

16

III. THE NUMERICAL SOLUTION

A. FORMULATION OF THE MODEL

1 . Establishment of a Representative Cross Section

In formulating the numerical scheme it was decided to
simplify the model by neglecting the derivatives along the y-axis.
To accomplish this without the loss of any essential physics, it
was decided to allow the seiche to oscillate across a narrow,
shallow channel of infinite length. This prevents sea surface
slopes from developing along the length of the channel but is
otherwise non-restrictive. The simplification further reduces the
problem to one of looking at a cross-sectional distribution of
variables in the channel.

It was further decided to restrict the model by considering
a cross-section of constant rectangular dimensions. The width and
depth of the basin considered were fixed at 5000 meters and 30
meters respectively, and Figure 1 is a pictorial depiction of this
cross section.

These simplifications allow eqs . (3) and (22) to be re-
presented in the model as

M * g - 0 (35)

,(x,z,t) - 22 E , (-1)"c°«(p») rV(kp2+")(t-f ) |i(f)df (26)

h m=l p J 0xv '

17

30 m

BASIN CROSS-SECTION

5000m

z

->x

Z:o

z:-h

X=o

MODEL CROSS -SECT ION

X=L

Figure 1. Pictorial Depiction of Model Basin

18

2 . Establishment of Governing Non-dimensional Ratios

To represent the various time scales and frictional co-
efficients possible with more general time dependent motion the
model was scaled to the period (T) of a free oscillating seiche in
the model basin, and the coefficient of vertical friction determined
by Ekraan's "depth of frictional resistance (D)". Equation (26) was
scaled as follows:

T = Jk . (27)

fgh

D = TT jli ; (28)

Let,

t

=

*
Tt

t*

=

*
Tt1

dt'

=

T dt'

*

*
Lx
x = — (eleven grid points across the channel)

.. 10A vP*
b£ _ o b£

bx L " v *

Ox

* * * *

where t , tf ,5 » x , are non-dimensional quantities and A is

o

the free surface elevation at x = 0, t = 0. Substituting eqs . (27)
through (29) into eq. (26) yields:

w = 2£S (-1) CQ5(P2) o_ re-(-^- + if)(t-t') T

h m=l p L JQ 2tt

be* *

X -% dt' . (30)

bx

Recalling that p = ' m~ ' ~, and a = -—-, eq. (30) may be written

19

i„ „ rt/T r0 ,f* (2m-l) ,D, ^ .

40g AT00 . 1Nm , f -12tt(-:) -^ *— (-) + i

^ o _ (-1) cos(pz)\ v<7 ' L o vh'

e_n / ->_ t \~ Ja e

L m=l (2iq-1)tt ^0

2 _ 2

(t-t»)*

*

x ^ dt' . (31)

bx

It was determined that the two non dimensional ratios,
(f/CT) and (D/h), appearing in eq. (30) govern the time/spacial
scale and frictional dependence of the motion. For discussion and
computation the first was designated (PCF) and the second (PCD) .

The first of these ratios represents the relationship
between the earth's inertial frequency and the frequency of the
basin in the following manner .

f = 2fisin 0

where £2 is the frequency of the earth's rotation and 0 is the

latitude.

a = 2tt F ,
m

where F is the resonance frequency of the basin. Since no external
m

forcing terms were used in the solution, F is also the approximate

m rr

frequency of the motion. Thus:

f _ Q sin0

OF
m

The model was scaled for (0 £ PCF ^ 1). A value of
(PCF = 0) represents motion of a time scale which occurs too
quickly to be affected by Coriolis accelerations, while a value
of (PCF = 1) represents motion whose time scale is equal to that of
inertial rotation and subsequently greatly influenced.

20

The second term, PCD, is the ratio of the depth of
frictional influence (Sverdrup, Johnson, and Fleming, 1942) to the
actual water depth. The coefficient of vertical friction was de-
termined in the model from eq. (28) as

k = (PCDX?h)2 Xf . (3la)

2n

The model was scaled for (0 ^ PCD ^2). A value of (PCD = 0)
represents a very deep basin where frictional effects on the motion
are minimal, while (PCD = 2) represents a basin where friction domi-
nates the motion.

If the time scale of the motion being represented by the
model is controlled strictly by the resonance of the basin, then
the spacial scale of the basin is determined by a choice of PCD
and PCF. PCD will govern the approximate depth of the basin while
PCF determines the basin width through eq. (27). Certain con-
clusions about motion caused by external forces (i.e. tides,
atmospheric pressure) producing time scales independent of the
basin size were also inferred from the model and are discussed in
the conclusions.

3« Adaptation of the Formal Solution to the Model

To establish two equations that could be solved simulta-
neously for sea surface elevation and transport, eq. (26) was
integrated in the vertical using eqs . (4) and (8). Thus:

p °

and W = U + iV. Equation (32) represents the total volume trans-
port associated with the flow at a particular point in the channel.

21

The real part of eq. (32) represents the cross-channel transport,
while the imaginary part the along-channel transport.

The numerical scheme was established to solve eqs . (25)
and (32) for § and W. The velocity profile, bottom stress, and
energy balance associated with the flow were also recoverable from
the scheme, and the details for this recovery will be discussed
later .

4. Initial Conditions

The model assumes an initial state of rest in the motion

and an initial set-up of the free surface. This was represented

as

W(x,0) = 0, (33)

and

?(x,0) = Aocos &Zl } (34)

where W = W(x,t) and § = ?(x,t).

5. Boundary Conditions

Since the governing equations were solved without reference
to horizontal boundaries, these conditions had to be built into the
model. The model therefore assumes no flow across the horizontal
boundaries, or

W(0,t) = W(L,t) = 0. (35)

These conditions were satisfied by requiring the surface slope to
remain zero at these boundaries at all times. Therefore eq. (32)
is forced to meet the conditions of eq. (35).

6 . The Step-by-step Scheme

Both Welander and Gait suggested that a step-by-step
numerical scheme could be used to solve eqs. (25) and (32). As a

22

first solution the staggered central difference scheme suggested
by Gait was used and is summarized below.

t = - At Given initial conditions on W.

t = 0 Given initial conditions on §.

t = At Calculate W from eq. (32) using r* at t = 0,

t = 2At Calculate £ using an explicit finite-difference

scheme on eq. (25) and U at t = At. Calculate a new
surface slope.
t = 3At . . . etc.

The term At is the time step of the explicit scheme.

This scheme was developed neglecting Coriolis (f = 0) in
eq. (32). The scheme proved to be stable and accurate subject to
the following stability criterion:

YIh^< 1 • (36)

The expression Ax is the spacial step across the channel. It was
felt however that a more precise solution, not subject to such a
restricted stability criterion, could be formulated using a totally
implicit finite-difference scheme. The development of this method
is now presented.

7. The Implicit Finite-difference Scheme

The implicit scheme used is a modification of a method
described by Richtmyer and Morton (1967).

Start by approximating the time rate of change of the
transport as follows:

bt At

(37)

23

From eq. (32)

bw
bt

P L^O

= I. ("2fl s JL

At / h m=l :

L P

^ e-(kp2+if)(t-tM |L(t.)dt,

t+At e-(kp2+if)(t+At-f ) 6S(t.)dt

(38)

By expanding the first integral as \ + \ and regrouping terms

J0 Jt
eq. (38) may be written as

»S = JL f-22 e 1 C* e-(Kp2*«)(t-f)(e-(kp2-if)At. M.(t,)dt,\

ot At I h m=l 2 J. v ' 0xv ' J

At ( h m=l 2 ^t

At e-(kp2+if)(t+At-t') b£

bxx '

(39)

Eqs . (25) and (32) can then be represented in the following form:

^bt; At 2

bu\ /bu
bx;t l^bx

(40)

t + At J

where

i$3L = #S_% + /C2_

vbt; At vAt ; v2At'
t 2

2 .

fill +/ss

bxjt ^6x,

2 .

(41)

a».-afi. 1 f .-civ +if)(t-t»> -(kp+if)At s(t.)dt. (42)

h m-1 2 J,. v ' 0xv '

P °

00 x pt+At

= m=l ~2 \
P t

C2* = "42s , i V" """ e-(kp"+if)(t+At-t')dt,

24

Equations (40) and ( 4l ) can be represented in Matrix form as
follows :

<t?'t+At=A + f^ [(6)t + <eW] (*3)

where,

0 =

B =

0 0 C2/At
0 0 D2/At
■10 0

and,

CI = Real {ci } ; C2 = Real { 02 }
Dl = Imag {ci } ; D2 = Imag { C2 }

(Real { } and Imag { } represent the real and imaginary part of
the argument in braces).

The following implicit finite-difference equation was used
to represent eq. (43):

0n+1- 0n

D_ _

At

= A +

B

4ax

9n + * - 9n+^ + 9n , - 9n ,
j+1 j-1 3+l j-1

(44)

where 9 .
3

QnAt
JAx

25

and n = (1,2,3, •••, maximum number of time steps)

j = (1,2,3,..., J)

are the grid points in the time and x-directions respectively.

Letting

a. = -, , and

4Ax

D. = tt9n , + 9n

- <*9. _ + AtA,

eq. (44) becomes

(45)

The left hand side of eq. (45) contains the unknown values of
(9) to be calculated, whereas the right hand side (D.) are all
known values.

A recursion relation of the following form was assumed.

9 . = E. 9 + F .

D J 3+1 3

(46)

where

E. =

eJ eJ
ell e12

'21
J

'22
J

13
'23

31 32 33

and

F. =
3

26

If eq. (46) holds for all 9., then the model's left hand boundary
conditions require that

0

0

0

E =
o

0

0

0

0

0

1

and

F =
o

(47)

,n+l

Substituting the expression for 9. derived from eq. (46) into
eq. (45) yields:

D.-a F

9

n+1

a

i+a e. ,

±1

i+a e

J'1 ..

(48)

Equating the right hand side of eq. (48) to eq. (46) gives the
following recursion relationship:

E. = M 1

3 I1** Ej-i.

D.-Q! F. '
-2 111

Fj

j-l J

, j * i;

> j a i;

(49)
(50)

and I is the ( 3X3 ) identity matrix.

Using eqs . (47), (49), and (50), values of E. and F. were
calculated inductively in order of increasing j(j = 1,2,3,..., J-l).

Since the value of 9 . , is given for j = J-l by the model's right

j+1 J

hand boundary condition (U = V = 0 and %, - % ), values of

9. were then calculated inductively from eq. (46) in order of
decreasing j(j = J-l, J-2, ... , 1). This completes the calculation

y. i 1

of 9 . (j = 1,2,3, ... , J-l). The details of the mathematics for
representing E., F., D., and eq. (46) are presented in appendix A.

27

8. Numerical Treatment of CI and C2

The numerical scheme for representing CI and C2 is
extremely lengthy and is merely outlined in this section. For a
more detailed presentation the reader is referred to appendix B.

The integration of CI was performed by dividing the total
integration time into small time intervals (At) and assuming that
the surface slope and effects of inertial rotation could be re-
presented adequately by a constant throughout ^:he interval. These
values were calculated at the raid-point of each At interval and
used in the scheme as representative values throughout the entire
sub-interval. This permitted the integral to be evaluated analyt-
ically over each interval with a maximum error on the order of
these approximations. Using these assumptions CI was numerically
represented as

cl*i = -& rv-L

j_(1.e- kP2At)

M=l (^p Ikp J N=l L

i-i r ^B .. i .2

S

,^^-h, - kp At -if At ..

2

X (e"if(I"^"N)At)(e'kp (I-1-N)At

(51)

*1

for I = 2, 3, ... , and CI =0. The value of the integer I ranges

from 1 to MAXT (maximum problem time) and IAt represents the elapsed
time from initial conditions. L takes on integer values from 1 to
MAXL (maximum number of cross-channel grid points) and LAx is the
distance from the left boundary. The integer M ranges from 1 to
MAXM, and determines the maximum number of modes in the series that
are included in the approximation of CI . It was determined from
eq. (23) that MAXM = 25 would represent the solution with a maximum
error of three percent, and this value was used throughout all

28

numerical evaluations. The integer N takes on the values from 1
to (1-1) and determines, through the exponential terms, how much
rotation and decay has taken place in the recent time history of
each transient considered. Since the terms of the series decay
more rapidly with increasing friction (especially the higher modes),
it was not necessary to carry their integration as far back as
N = 1. It was decided to include only the terms of significant

magnitude, and this was accomplished by requiring that when the

2
value of kp (I-l-N)At in eq. (51) was greater than 10, the cor-
responding terms were discarded.

The integration of C2 involves only one time interval
(t to t+At) . Again the amount of inertial rotation was evaluated
at the mid-point and assumed constant throughout the interval, and
02 was evaluated as the following constant:

C2

_ MAXM . ,.,! A„_x
• = - 2g_ s e-if(?§At)

h M=l

2- (i-e-kP2At)

kp'

(52)

CI and C2 both contain the term

— — (1-e p ) , and for k = 0 this term can
. kp J

kp
not be numerically evaluated in this form. However this difficulty

2

. , . -kp At . _
was overcome by expanding e c in a Taylor series, i.e.

(i - kp2At * i^n\ us^Ati3)

(53)

and substituting this series for the term in eqs . (51) and (52)

2 -U

whenever the value of kp At was less than 10 .

The method for separating the real and imaginary parts of

* *

CI and C2 involves using the relation

e = cos(ft) - i sin(ft)

29

in eqs . (51) and (52) and gathering terras. The details are
covered in appendix B.

9. Velocity Calculations

The velocity profiles across the channel were calculated
by the model at discrete time intervals determined by an input
parameter. Noting the simularity between the integral of eq. (26)
and CI , it was numerically possible to compute and store the terms

of these integrals needed for velocity, during the calculation of

*
CI . For velocity profile computation, these stored values were

recalled and applied in eq. (26). The details for numerically re-
presenting eq. (26) is also presented in appendix B.

10. Bottom Stress Calculations

The bottom stress (T) associated with the motion was
calculated using the following relation:

' = k (If ) ...*

or fr om eq. (26) ,

r=Ms" f e-(kP2+if)(t-t-) jji(t.)dt. .

h m=l Jn bxv '

(54)

Since the integral term of eq. (54) is identical to that
of eq. (26), bottom stress is easily calculated whenever the
velocity profile is recovered. Appendix B also contains the
details for this calculation.

11. Energy Calculations ,

The total energy balance associated with the seiche is a
sum of its potential and kinetic energy. This balance was

30

calculated as follows:

PE = h p g \ %2 dx (55)

J0

KE = h P \ \ w dx dz ( 56 )

PE and KE are the potential and kinetic energies and (p) is the

is the density of the water. Equations (55) and (56) were evaluated

using Simpson's 1/3 rule (McCalla, 1967) whenever velocity was

calculated.

The energy calculations were monitored during model
testing as a calculational check since small numerical errors
were quickly apparent in increasing energy levels.
12. Scaling of Model Output

Since the governing equations for the model were scaled
to produce the ratios PCF and PCD, the output of the model must
be appropriately re-scaled to obtain real values for the variables.
The series of graphs were produced on a CALCOMP plotter with cor-
responding values from the model. The various axes were auto-
matically scaled by the computer to fit the data being represented.
These axes were then re-scaled using a non-dimensional scaling
factor derived from the continuity equation to allow real values
to be easily obtained. The exception is the figures depicting the
energy balance, which are not scaled and should only be used to
indicate relative changes in energy levels.

31

The scales for the various axes are as follows:

AXIS SCALE

Surface Elevation §/A

o

Time

Velocity

t (gh)2

2L

2w (h/g) 2

C..A
1 o

2W
Transport 1 —

C2Ao(gh)^

2T(h)3/2
Bottom Stress ' i

C3Aokg^

The values of g , h, k, A , L are determined from the
actual basin being represented and C. , C , C are scaling con-

J- £• J

stants dependent on the parameters of the prototype basin and the
scaling of the graphs by the computer. These values are calculated
for each figure and inserted in the scale. The variables, t, w, W,
§, and T are the real values associated with the basin being re-
presented and may easily be obtained from the graphs by entering
appropriate arguments for g, k, h, A , and L in the corresponding
scale.

B. TESTING THE MODEL

1. Establishing a Reference Solution

One of the major difficulties in numerical analysis is
verifying the accuracy and stability of the scheme. Since many
numerical solutions deal with problems that have no analytical

32

solutions it is common practice to test the accuracy of these
models against empirical data, or analytical solutions involving
degenerate cases. As previously noted, data pertaining to time
dependent flow is extremely scarce, and it was decided to test the
accuracy of the model against a known analytical solution for the
degenerate case where PCD = 0 (i.e. no friction). Analytical
solutions to this case are readily available in the literature and
a modification of one proposed by Prodman (1952) was used in the
verification of the accuracy of the numerical scheme.
2 . Accuracy and Stability of the Model

The accuracy of the model was determined by comparing the
numerical representation of the free surface elevation at (x = 0),
the raid-channel transport, and the total energy balance, with the
analytical values. Figures 2 through 5 highlite this comparison
and indicate that the accuracy of the model was a function of time
step (At) and the Coriolis parameter (f). It was noted that as the
time step was reduced the numerical solution appeared to converge
to the analytical solution. The phase shift in surface elevation
over five periods was reduced from 45 degrees for (At = T/20) to
24 degrees for (At = T/40) and appeared to be independent of the
value of Coriolis. This error might possibly be introduced into
the system by approximating the sea surface slope as a constant
value in the time history integration scheme. As (At) increases
this linear approximation is extended over a greater time span
and the resulting errors amplified. As the value of the Coriolis
parameter was increased, the error in representing the amplitude
of the free surface and total energy balance also increased, as

33

hX01

OljCM

cd

e

CD
N

(T5

e

u
O

2:

o

CM

O <

•H — '

>x r-l

i-H CD

c o

< 2

C

o

•rl

o

o

•H

>»

iH

c

0)

-o
o
s

w

o

Vi

c
o

A3
>

CD
iH
(X)

CD

o

Vi

u

Crt

CD ^

CD O
W

U. II

<H Q

O U

IX

c

o •>
tfl o

•H

w II

rtJ

CUfc

£ U
O (X

u —

CD

VI
•H

fa

Normalized Elevation at X - O (§/A )

34

s:

01

0>

B

•H

H

0)
N

e

O

•P

c
<

0)

o

2

h\n

Dl CM

0>

e

•H

H
*0

N

E

H

o

C

0)

XI

o

w

o
<H

c «—»

o o

•H
•H II

rt ^

> Q

0) U

H rv
(i)

0) H
O

(0 II
<H

H fc.

3 U

03 CU

0)
0>
H

Cx,

«w
O

c
o
in

•H
VI

m

c
o

iH

o

CX H

o c
u <

ro

0)

Normalized Elevation at x = 0 (§/A )

o

35

m

o

CM

H

II

•p

o

<

•H

■♦-»

>>

i-i

iH

<D

m

TD

c

0

<

2

- CO

01

o

&

u

IX

o
II
tl.
u

c
o

•H
•H

iH

o
u

•rl
4-»

fH

c

<

C
(0

0)

o

w
O

CM
QJ
U

c

CM

(1)
N

<0
E
u
O
2

PQ

vt

<u
c
w

o

c
o
</)

•H

w
rt

o
u

O
w

•H

fa.

ON

VO

CO

Model Energy Balance (joules X 10 )

36

o

CM

\
H
II

<

o

<

m

(0

c
<

c»

H^

&
w

0)

B
•H

H
<H

o

B
O
•H
■H
O
C

CM

CNCM

0)

B

•H

H

•o

(1)
N

(T3
B
w
O

(A

<T5

0)

o

c
as
l-l
(T5

VI

c

W

c

rH

CM

ON

VO

m

Model Energy Balance (joules X 10 )

c

x:
u

in

0)

w

0>

•rl

to.

37

depicted in Figures 3 and 5- The amplitude deviation was a maximum
of one percentage for (At = T/20) and the energy rise was 8.5 and
3.7 percentage for (At = T/20) and At = T/40) respectively. Errors
of this nature did not appear for solutions which neglected Coriolis,
and possibly were introduced by approximating the amount of rotation
[i.e. e * ' ] as a constant in the integration scheme. It is
not exactly clear how this error increases the energy in the model
but it was observed that it appeared linear and very systematic.

To determine whether this energy rise constituted a
numerical instability, a standard Fourier series method (Smith,
1965) of stability analysis was attempted on eq. (44). However,
due to the complicated nature of the amplification matrix, the
eigenvalues could not be readily obtained and the analysis was
abandoned. As an alternate test, the model was run with (At = T/5)
for (t = 5T) . This corresponds to fgh — > 3 , and resulted in a

l AX _t

further linear rise in the energy level. These errors were there-
fore judged to constitute a problem in numerical accuracy and not
stability, for the range in which the test cases were run.

It was further found that the magnitudes of these errors
could be predicted and controlled by varying the time step, and the
model was run maintaining a maximum allowable error of one percent
in the elevation and transport amplitude, and four percent rise in
the total energy balance. This was accomplished by using a time
step of (At = T/20) for (PCF < .5), and (At = T/40) for (PCF ^ .5).
Since the errors in phase shift were strictly a function of time
step, they were completely predictable and accounted for during
data analysis.

38

It should be noted that for cases where friction is
neglected (k = 0), the transients do not decay with time and must
be carried throughout the integration scheme. Additionally, the
bottom boundary condition degenerates to a free slip situation and
the approximation of a constant velocity profile with a truncated
cosine series becomes more inaccurate. The number of calculations
increase rapidly as friction is decreased, which may introduce
increasing errors in round-off.

Based on the comparison with the analytical solution it
was concluded that the accuracy of the model could be sufficiently
controlled such that the details of time dependent motion, as
described by the governing integral and differential equations,
could be adequately represented and examined.

C. MODEL LIMITATIONS

The most serious limitation imposed on the model were the
basic assumptions of Ekman dynamics used in the formal solution
of the governing equations. These assumptions were considered
necessary to simplify the mathematical treatment of these equations
to the extent that a numerical solution was feasible. The additioanl
simplifications imposed during the formulation of the model further
limit its capabilities. It was felt however that valuable experi-
ence in the techniques of numerically representing time dependent
motion would be gained during the formulation of this simplified
model, and insight gained from the numerical product. With this
information available it might be possible to further refine the
model and remove many of the more important restrictions. The

39

limitations inherent in the present model are listed in this
section and the more important ones will be discussed in section
(V) with recommendations for improving the scheme.

The basic assumptions made in arriving at eq. (22) have the
following limiting effects of the model's capabilities: 1) the
non-linear terms in the motion have been ignored; 2) a variation
in the coefficient of vertical friction in the layer was not
permitted; 3) the assumption of constant density neglects the
effects of stratification on the pressure gradients in the layer,
and the subsequent modification of the vertical momentum transfer;
and 4) horizontal friction was considered insignificant and
neglected.

Simplification of the model resulted in the additional
limitation: 5) the assumption of a constant basin neglects the
important effects of bottom topography on the flow. This is es-
pecially important for considering motion in shallow near shore
areas where these effects are large.

40

IV. PRESENTATION OF RESULTS

A. TEST CASES AND GEOPHYSICAL SIMULATIONS

1. Description of Test Cases

In order to gain maximum insight into the physical processes
occurring in the motion, as well as to study the effects of inertial
rotation and friction on the variables, the model was run for 25
test cases in the ranges of (0 £ PCF ^ 1) and (0 £ PCD ^ 2).
Table I summarized the various controlling parameters of each run
and includes the Central Process Unit (CPU) time used on the IBM
36O/67 system for each program.

2. Description of Geophysical Simulations

As interesting dynamic interactions developed in the test
cases, it was decided to simulate actual geophysical embayments
to determine to what extent these phenomena were present and
significant. The actual length, depth, and latitude of the basins
were measured and used to calculate (f ) , (0-), and (PCF). An Ekman
depth of 20 meters was assumed in the calculation of (PCD) . Runs
26-3O represent these simulations and the various imput parameters
are summarized in Table II.

B. DISCUSSION OF TEST RUN RESULTS

1. Frictional Effects on the Flow

To describe the conditions where friction totally dominates
the motion, the model was run with varying amounts of friction and
(PCF =0). It must be noted that by neglecting the effects of
inertial rotation, the Ekman depth is no longer specified as

hi

Table I. Parameters of Test Cases (runs 1-25)

RUN NO.

PCF

PCD

CPU(MIN-SEC)

1

.25

.5

6-2

2

.25

1.0

3-51

3

.25

1.5

3-8

4

.25

2.0

2-49

5

0

-5

5-H

6

0

1.0

3-20

7

0

1.5

2-4?

8

0

2.0

2-24

9

0

0

13-33

10

.25

0

14-58

11

.5

2.0

6-25

19

.75

2.0

6-0

13

1.0

2.0

5-48

14

.5

1.5

7-25

15

.75

1.5

6-42

16

1.0

1.5

6-25

17

.5

1.0

9-40

18

.75

1.0

8-21

1Q

1.0

1.0

7-40

20

• 5

.5

16-14

21

.75

.5

13-52

22

1.0

.5

12-29

23

.5

0

62-48

24

.75

0

62-38

25

1.0

0

62-51

NOTES: For the model basin
10 T = 583.2 Sec

2)

CT = .108 X 10_1 Sec"1

42

c
o

•p

i-H

0

•rH

co

o

•H
0)

si
&

e>

O

0)
VI

0)
+■>

0>
E
«C
w
flj

M
iH
H

u

tt

CO

1

co

rH

ON

ON

rH

z

m

CO

CO

<r

CM

t-H

i

1

I

i

1

£

m

l^

CO

r-i

CO

s„

rH

r4

r-t

rH

D

0.

U

^__

1

eg
i

CM
1

CM
1

CO

CM
1

O

m

Ui

W

w

w

(1)

vO

o

<t

CM

CM

CO

un.

<r

t-i

rH

t^

»_

m

co

r-i

m

r-i

b

•

•

•

•

•

r-\

iH

iH

rH

1

1

w

(I)

b.

m

iH

CO

o

VO

u

m

VO

oo

On

m

n.

rH

CM

in

rH

m

•

•

•

•

«

Q

00

CO

rH

co

u

CM

o

r^

ON

r^

0-

•

*

CM

•

*

•

r-t
1

i

-3-

1

i

i

i

(J

m

w

w

w

w

0)

CM

r^-

m

CO

<r

CO

VO

oo

On

t^

m

00

00

00

ON

ON

NH

•

•

«

•

•

r-l

m

•

CM

co

0

0

0

H

0

o

oo

CM

rH

<

VO

fv.

co

~d-

<r

J

ro

m

O

o

o

o

o

o

o

o

o

o

IT

o

VO

o

o

o

•*

•k

•fc

■l

w%

j

ITS

ON

m

o

o

r-i

<r

OV

co

"e

CM

o

i>

CM

r^

«— .

t^

iH

CM

CM

CM

x:

2

D

VO

r^

00

ON

o

ce

CM

CVl

CM

CM

CO

>v

>>

«J

(0

&

cq

V

CQ

a>

c

c 0

X

CD

<a

Z,

>»

<o y

«j

•H

rH

M

Q>

(0 0)

CD

W

</)

to

w

•H

a

(I)

W

<

CD

J3 U

cfl

TJ

CQ

■*->

H-> C

V)

0)

O* C

C

3 rtJ

0

.*

C 3

0

0 w

x:

(0

0 0

s

to U«

u

J

J to

I

o

I

a)

I

43

evidenced by eq. (28). Therefore, the frictional coefficients
generated by runs 1 through 4 were artificially imposed on the
model with (PCF = 0), and resulted in runs 5 through 8.

The results for these cases are depicted in Figures 6
through 8 and are simply a seiche damped by friction. As the
friction coefficient increased the damping of the motion increased
with a resulting increase in the period of oscillation. For the
maximum case of (PCD = 2), the seiche was completely damped after
4^2 periods. These results were completely predictable and are
presented to show that the model's representation of motion which
is dominated by friction conforms to previous observations.
(Pierson and Neuman, 1966)

An interesting coupling between bottom stress and transport
developed during these cases and a typical example is depicted in
Figure 9. In all cases examined a phase lag quickly developed be-
tween the transport and the bottom stress. For the case showed,
the value of this lag averaged 45 degrees through the periods of
oscillation. This phase lag can be explained by recalling that
the transport is an integration of the motion over the entire water
column, whereas the bottom stress is solely dependent on the
velocity gradient at the bottom. Since friction has its greatest
effect at the bottom, changes in magnitude and direction of the
flow nearly always initially occur next to this boundary and prop-
agate upwards in the column. The result is that changes in trans-
port always lag changes in the bottom velocity gradient and
consequently the bottom stress. This phase lag has obvious effects
on the accuracy of commonly used numerical schemes where bottom

44

O <"H CM

H II II O

Q Q Q
U U U

a, a. a.

rX*

CT CM

CD
E

N

e

w
0

2:

n

b,
U

a.

| I ft.

I i
I i

a*
o

10

HH

w

(0

0)
01

VI

fr.
0)

+■»

c
o

c
o

•rl

W

o
tf)
o

<H

(1)
w

Normalized Elevation at x = 0 (5/A )

o

45

f /

; /

\s

X

/i_ —

— : ~\

V \

V \

\ \

\ V

H^

_T>

)

/ y

t ^

/ y

f ^

j* -

-j^l-

-r \

K

v \

m

oo

CM

h\N

O rH CM C

II II II II

Q Q Q ta

V u u u

n a. cl a.

w
0
Cu
</)

c

rt3
u
H

CD

c
c
(C

Xi

o

I

T>

•r(

2

§

C
0
•H
•P
O
•H

VI

o
o

cr cm

0)
E

0)
N
'ri
r-i
(0

£
u
O
2

0)
•H

Normalized Transport » J'^i. g

A

<

46

AQ(gh)

73
<D
N
•H
iH

e

u
O

o
II

a.

V

(1)
o

c

«5

rH

(d

DQ
^

VI

0)

<§

c
o

c

o

•H

o

■H

w

(I.

O
(fl
o

0)
VI
VI

(I]

00

w

•H

Model Energy Balance (joules X 10 )

47

Normalized Bottom Stress

o

T(1.9) h

3/2

kAQg^

m

ro

rtf»

CVJ

0> (M

0)

e

H

TJ
(1)
N

<T5

e
o

52

(5 4) W
Normalized Transport J — : — t- £

.. 48

AQ(gh)

it

Q
U

O
II

u

Du

■♦->

w
O

a
<fl

c

w
H

0)

c
c
m
x:
o
i

•H

s

C
rtJ

</)

</)

QJ
w
•*->

E
O

o

CO
VI

o

c
o

U)

•H

w

o
u

ON

0)
w

3

•H

stress is parameterized in terms of the transport. This technique
places strong restrictions on possible interactions between the
flow and the bottom, and although it may be sufficient for approx-
imating steady state problems, it is unlikely that it will represent
transient motion very well.

Figure 10 represents typical velocity profiles for motion
which is totally dominated by the effects of vertical friction.
The increased period of oscillation and the mechanism causing the
phase lag between the transport and the bottom stress is clearly
evident. At times (t = T/2) and (t = 3T/4) the profiles are in
transition from being totally positive to totally negative, and the
transport and bottom stress are 180 degrees out of phase.

These velocity profiles are typical of those found in other
boundary layer flow and are similar to the wind profiles observed
at the air sea interface over smooth water. They indicate that the
stability of the flow has a strong frictional dependence. At
(t = T/4) and (t = T) the profiles are maximum and, except initially
in the upper layers, friction exerts a strong influence throughout
the column. The result is a smooth regular flow and a well de-
veloped frictional boundary layer with the region of greatest shear
at the bottom. During the transition phases (t = T/2 to t = 3T/4)
the profile changes direction and the water column appears to be
less stable. The velocities in the column are minimum just prior
to and after reversal and thus the influence of friction in the
column is reduced. This region of reduced frictional influence
moves up the column at the point where the flow changes direction.
As the flow completes its reversal and the velocities in the column

49

a

Q
N

e

u
O
2

-2

-1

Normalized Velocity w(1-9Hh/g)'

o

Figure 10. Effects of Friction on Mid-channel Velocity
(PCF = 0, PCD = 2) .

50

increase, the effect of strong frictional influences is quickly
re-established throughout the layer, the region of maximum shear
shifts to the bottom, and smooth flow once again prevails at time
(t = T).

2. Coriolis Effects on the Flow

To gain insight into how inertial rotation modifies the
transient motion, runs 10, and 23 through 25 were made neglecting
friction (PCD = 0) .

It is important to realize that the relationship between
the transport and the free surface elevation for a non rotational,
frictionless seiche, as seen in Figures 6 and 7, is similar to the
velocity/elevation relationship of the frictionless pendulum. At
(t = 0), (t = T/2) , and (t = T) , the energy of the seiche is en-
tirely potential, whereas at (t = T/4) and (t = 3T/4) it is
entirely kinetic. At intermediate times the energy balance is a
combination of potential and kinetic and is totally conserved. The
free surface elevation and transport periodically oscillate between
maximum positive and negative values of respectively equal magni-
tudes, always maintaining their constant phase shift. Figure 2
represents this oscillation for the free surface.

As the effects of Coriolis acceleration were introduced
into the model these conservative oscillations were greatly modified
as can be seen from a comparison of Figures 2, 3, 11, and 12. The
effect of a moving earth is that the total potential energy of the
surface elevation will not be available for transfer to kinetic
energy since the water particles are forced to turn due to the
Coriolis Acceleration. Although the energy balance still remains

51

m

m ii

OJ

m

r-«

•

•

•

Q
U

II

ii

ii

a.

fe

Cx.

fc<

U

u

u

c

CL

D-

a,

0

I i

N

•1-1

rH

8

w
O
2

•H

o

rtS
Vi
w

W

c
o

c

o

•H
■P

«5

O

cc

iH
AS

•rl
4->
W
0)

c

Vl
o

(J

VI

0

w

en

Normalized Elevation at x = 0 (§/A )

o

52

> D

1

! 1 1

J^—-

— ^Xr"~-— i. i^

^^^^--"'^"^ -

^^^^-"

"">>

-^j~~ "~~^"^

S

c~~~

xT~-— - ^c-~~^

*-»

^^

^^>-<r""^ "*■*"■*»

- — — ^>>

* «»^

C __ __ — - — ""

__. — — "-"""" ^^"^

^>

^-"C — — •— .

•*^"-- ""*"''"""""'

"" ./_ — "" *'"°

^<^><

_— m ~~~_

r^^^11^1^

**
\

^-*vT

-— "■"S-^"'^

—- "*" "*"^***^

_«.

<rr

^^-^^^

^•^><^"*-o

—

- .^^ -— * ^ ^"—

>«.

■^ , , — — *-*

\

. — -

* ** "^^ — "■"■■—

J

r"-"

^^*^^w^- -^^^**>

*s..^

— • **" *^<^^-^.

^*"** — — -— -e=r-""

-"""^^^^

r"

-Z_ — „» "■""

v.^

^V "*"~ « ' "**

^^** _^~

"" """"" ^^

\

^^^^ ~~ ^""""*

^^^s

■ IX^'

in

_ -d-

m

ii ii
u u

Qu Cl

— to

oj

h\m

Oil

0)

e

•H

H

TJ
0)
N
•H
i-H
<T5

e

w
O

Q

U

a.

+■>

u
O

a

CO

c

VI

H

c

c

o

•H

§

c
o

•rl
■H

v>
o

OS

VI

c

M

o

0)

■»->

o

0)

VI

M

w

O*
•H

(5.4) w
Normalized Transport J t— t-

A0(gn)

53

conservative, the relationship between the free surface elevation
and the flow changes significantly.

To understand how rotation acts to modify the flow, con-
sider first the physical processes which drive the system. These
are discussed in terms of the frictionless irrotational seiche.
The initial conditions for the model were identical for all runs;
that of a positive set-up of the free surface at the left boundary
(x - 0) and no motion in the basin (a condition of static equili-
brium). As the system is released the pressure gradients developed
by this set-up start the fluid moving in the positive x-direction
resulting in a net divergence in the left side of the channel and
convergence in the right. The free surface responds to this flow
by slumping at the left and rising at the right with equal magni-
tudes. At time (t = T/4) the surface is level (i.e. no gradients)
and the kinetic energy associated with the flow is equal to the
potential energy of the initial free surface set-up. Therefore the
surface continues to fall at the left and rise at the right until
this kinetic energy has been entirely converted into potential
energy and a state of momentary equilibrium is once again reached.
The sea surface is now set up at the right boundary equal in magni-
tude to the initial conditions, the pressure gradients equal but
opposite to the initial values, and the flow reverses direction
and the system returns to initial conditions, etc. It is important
to realize that only the divergence of the flow in the cross-channel
direction can influence the system because no gradients are allowed
to develop along the length of the channel.

54

With these processes in mind the observed effects of ro-
tation on the system can adequately be explained. With the initial
set-up previously described, the flow will start in the same manner.
As soon as the motion begins however, the Coriolis acceleration will
start to turn the flow in the negative y-direction (left rotation),
and a component of the flow will start to develop along the axis of
the channel. Therefore some of the potential energy is converted
into rotational kinetic energy. The amount of rotational kinetic
energy that the system receives is a function of the velocity of
the flow, the value of Coriolis, and the time scale of the motion.
The slumping of the surface elevation at the left and rising at the
right will continue as long as there is divergence in the positive
x-direction; however since a certain amount of energy is rotational,
the set-up at the right boundary can not achieve the magnitude of
initial conditions. When the cross-channel flow vanishes the en-
tire motion is in the along channel direction, with subsequent
rotation turning the flow back in the direction from which it
initially came. It should be noted that for a maximum value of
(PCF = 1), Figure 3 shows the free surface elevation at the left
boundary never falling below the horizontal.

An interesting relationship developed in the transport as
rotation was introduced into the system, resulting in a rectifi-
cation of the along channel component of the transport (V) . An
examination of Figure 12 showed that the cross -channel transport (U)
periodically oscillates between maximum positive and negative values
while the along channel component oscillated between zero and a
maximum in one direction. It was easily demonstrated that this

55

direction is strictly a function of initial set-up. As the flow
was initiated in the positive x-direction the initial along channel
component developed in the negative y-direction. When the cross-
channel flow reverses direction, along channel components are
generated in the positive y-direction. These components are just
sufficient to cancel the components previously established in the
opposite direction, resulting in a rectification in the along-
channel transport. For motions of large time scales, this process
implies that important changes in the flow are caused by Coriolis
accelerations and these will be discussed in the conclusions.

Another observed feature of rotation was the subsequent
reduction in oscillation period for increased values of Coriolis.
The effect of rotation in reducing the spacial scale of the motion
allows the system to oscillate with greater frequency.

It is apparent from the previous discussions that the
effects of inertial rotation can not be neglected when dealing with
large time scale motion.

3. Combined Effects of Friction and Coriolis

The remainder of the test runs were used to investigate
how friction and inertial rotation in various combinations effect
the transient motion. Some of the important consequences are de-
picted in Figures 13 through 17.

The phase shift difference between bottom stress and
transport existed throughout all values of (PCF) and (PCD), as
evidenced in Figures 13 and 14 . This phase lag occurred in varying
degrees but was a definite characteristic of time dependent motion,
vanishing only when the motion became critically damped.

56

T(1.9) h
Normalized Bottom Stress — J f

3/2

kA0g'

■y\\

// 1

if \

II It

fl f

llf

D -Yr- K *

•u

,_^

1)1 ♦> t-

// H X ^

:/ 1<_ s o

// a aw

;/ 1 (n to

// / C «

| V // H W

MV

!

\ \\

i

M\

i

V A

i

•'/ i

i

// ;

i

\

i

\ \i

N \\

* )\

//\

'/ l)

// /

' ( A

K \ //

v \ (i

> \v

\ \v

N YV

*S\

—

'/ \

.<s V

S^r f

' ( yy

f V /:

\ v //

V \//

V ^wl

>». vv.

^ ^fcS

N ^5tv.

.> V^V

y\^ V\

.*c^^ \ 1

.^z^^ \ J

.*^^ )/

'" c ^<^

■■

< \ X^V

V ^-^ /^ '

>». ^^>^l i

^•^. ^^4^

>». ^Tv^^

"»•. ^^wJN.

V. ^^-Os^

-^^^-

I J

-^. •^^-^■""'"^ <

£x^

.*"*" / ^-**"

.** c ^~<&^

<• >v^ ^*-~-*^

"rr^**— L^>-^

m

— m

CM

cnl

J

w

e

•H

H

TD
CD
N
•H
iH

m

e

w
O
2

Q
U

a.

m

a

to
a

(0

w
H

(1)

c
c

(C

O
I

n

T3

c

(0
(0

to
(V

w
•p
to

e
o
+■>
•»->
o

CQ

Sh
O

c
o

tO

•H

W

§•

o
u

ro

H

( I.U) w

Normalized Transport * ' x

AQ(gh)

57

Normalized Bottom Stress — » — '- — ' — r

3/2

kAQ(gh)

c

Normalized Transport ' ' ^

A

i

58

A0(gn)

m.

ro

h\n

eg

01 CM

T3
N

(0
E

VI

O
2

■H

I

Q
U

(x.
U
Pu

s

a

c

VI

H

O

c
c

x:
o
i

■o

•H

■O

c

CO

</>
</)

CD
VI

■«->
W

B
O
•H
•V

o

CQ

<H

O

c

o

(A
•H
VI
«J

§■
o
u

VI
•H

b.

<D

o

to

MH

M

3

(/)

OJ

(V

£

0)

.c

•H

c

0

</)

•H

iH

0

•H

W

>J

0

*"— -*

U

&

►J

01

M

C

4->

flj

0)

C

B

0

•rl

•H

H

-t->

o

TD

•H

0

H

N

fc.

•H

rH

<H

«0

0

e

w

</)

0

■P

2

O

<u

<h

«H

w

m

fa.

Normalized Elevation at x = 0 (?/A )

59

m

m

r-l eg

ii ii ii ii

Q

Q

Q

Q

U

O

U

U

Ou

CX

CX

(X

v\

#\

*»

CM

m

*■>

o

•

«

iH

II

11

ii

II

Cx,

Cx,

Cx.

Cx.

U

U

U

U

IX

CX

CX

CX

I-

CO

CM

o

C

rt

iH

rt
CQ

>*

en

VI

0)
<§

§

c
o

•r(
•P

o

•H
W
Cx.

c

«J

c
o

•H
4->
«5

OljCM -H

o

cx

b «j

•H .rl

H

hX01

•a

N

B

H

o

H

c

VI

o
o

CD
Vl
Vl

w

\0
iH

u

•H
Cx,

Model Energy Balance (joules x 10 )

60

o

H
in

il

•♦-»
•H

o

o

H
>

^ IT» ^

O • r-«

II II II

£ £ j:

v. \ \

N N N

I I I

0

c

c
ns
.c
o
I
•o

•rl

o

x:
a

05
w

o

TD
0

K

I

X <"

r-i

0)

M

•H

61

The rectification of the along channel flow was also a pre-
dominate characteristic of the flow where inert ial rotation was
considered. Figures 13 and 14 again show this relationship, as
modified by friction, and Figure 15 depicts the combined effect on
the free surface. It can be seen that for an extreme case (PCD = 1,
PCF = 1) the free surface only very slowly slumps to a zero potential
because friction retards the cross-channel flow long enough for ro-
tation to reverse its direction before the free surface reaches the
hor izontal .

Figure 16 indicates that with rotation the energy level is
damped with relatively fewer oscillations. It must be realized
however that as (PCF) increases, the model is describing motion on
an increasing time scale.

Figure 17 represents the path that would be taken by fluid
particles (or suspended matter) at three levels in the column over
five periods of oscillation. In addition to the rectification pre-
viously noted, there appears to be a net drift of the motion in the
initial cross-channel transport direction. A closer examination of
Figures 13 and 14 also show this drift and is a combined result of
rotation and friction. Since the motion is being continuously
damped by friction, the cross-channel transport can not return the
fluid to its initial position after each oscillation and a net
transport in the initial direction of motion results. The net
drift is then strictly a function of initial set-up conditions of
the free surface. The velocity profiles depicted in Figures 18
through 20 show some of the fine details of the flow with the
combined influences of friction and inert ial rotation.

62

Figure 18. Effects of Friction and Inert ial
Rotation on Mid-channel Velocity
(PCF = 1, PCD = 1) .

<5

N

£1

a

Q

V
0)
N
•H
■H

£

w
O

2

-2

-1

0

Normalized Velocity 2ILl21I2iM

T/4
T/2

63

-2

-1

■p

a

0)
Q

TD
0)
N

•i-l
rH

e

w
O

z

Norma

i- ^ „ -. •+ w(11.4)(h/g):
lized Velocity — i " *"-

Figure 19. Effects of Friction and Inertial Rotation
.. on Mid-channel Velocity (PCF = 1, PCD = 1,
t = 3T/4).

64

-2

-1

Norma

.c
+->

9*

o

•o

o

N
•H
iH

<T5

e

u
0

z.

w(2.86)(h/g):
lized Velocity ■* ;' 2-i-

Figure 20. Effects of Friction and Inert ial Rotation on

Mid-channel Velocity (PCF = 1, PCD = 1, t = T)

65

The controlling influence on the stability of the water
column continues to be the friction, however the interactions of
the fluid are now more complicated due to the introduction of
inertial rotation. Under the influence of increasing frictional
dominance towards the bottom, the left rotation induced in the flow
by Coriolis is increasingly damped, resulting in the familiar Ekman
spiral. During periods of the motion previously described as
smooth, the spiral remains fairly constant as the column rotates
to the left. However, during the transition periods (Figure 19)
when frictional influence is reduced due to decreasing flow velocity,
the spiral unwinds from the bottom and var ing amounts of twisting of
the column results. This is again evidence of the transitional
nature of the fluid during these periods. It should be noted that
the velocity scale of Figure 19 is an order of magnitude smaller
than that of Figure 18, indicating how small the velocities of the
column are during transition. As the flow once again becomes uni-
form, the spiral is quickly re-established under the stabilizing
influence of bottom friction, as depicted in Figure 20.

It is interesting to note from Figure 18 that during the
transition periods the along channel component of velocity does in
fact change direction for brief periods. Because these periods are
of such short duration, this velocity reversal along the channel does
not contribute to any significant transport, in this direction.

C. DISCUSSION OF THE GEOPHYSICAL SIMULATIONS

All of the previously discussed processes characteristic of
time dependent motion, were observed during runs 26 through 30.

66

However, since the magnitude of these processes are Coriolis and
frictional dependent, they occurred on a much smaller scale than
observed for the test cases.

The model generated the time dependent variables produced by
the seiche on a time scale equivalent to the resonance period of
the basin. For these geophysical simulations, the periods were
all too small to allow inertial rotation to have any significant
effect on the motion. However, for shallow seas of large hori-
zontal extent, values of (PCF) approaching .25 are entirely pos-
sible and test runs 1 through 4 showed that inertial effects would
be significant and must be considered. Additionally, motion whose
time scale is driven by external forces, independent of the re-
sonance basin, are likely to be significantly altered by rotational
dependence. For example the long period waves propagating up and
down the continental shelf probably experience all the rotational
characteristics associated with the time dependent motion described
in the previous section. The observed phenomenon of transport
rectification might occur in these systems and suggests a mechanism
for modifying bed load transfer across the shelf. This process is
a function of initial conditions and appears to require a period of
static equilibrium just prior to initial motion. A region where the
seasonal migration of atmospheric pressure systems are predominantly
in one direction, might satisfy these requirements, with a resulting
net rectification in one direction. Although this process is highly
speculative, it certainly appears possible and a modification to
the model to facilitate its study is suggested in Section V.

The frictional effects observed during runs 26 through 30 were
also minimal. It should be noted that Ekman depths of greater than

67

20 meters are possible and thus these test runs represent a minimum
frictional dependence.

The exception to these minimal effects was the observed phase
lag between transport and bottom stress. This lag was pronounced
in all cases and for the simulation of South San Francisco Bay,
for example, averaged 40 degrees. This indicates that even small
amounts of friction are important and re-emphasizes the need for a
mathematical representation of time dependent transients that allows
independent calculation of variables such as bottom stress and
transport .

68

V. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER WORK

The model demonstrated that eq. (22) was a suitable mathematical
representation of the time dependent motion described by the governing
equations and the initial/boundary conditions. In this formulation
the time dependent variables were uniquely determined by the local
time histories of the surface at each point and this allowed inde-
pendent representation of the various transients. The model further
demonstrated the feasibility of developing a totally implicit finite-
difference scheme for approximating eqs . (3) and (22) with accept-
able accuracy. With the insight gained from the results it should
be possible to improve the accuracy of the system by removing some
of the important restrictions imposed by the Ekman dynamics and
the simplified model.

The physical processes observed with time dependent motion
showed the relative importance of friction and inertial rotation
in determining the nature of the variables. It was apparent that
for the larger time scales (PCF on the order of .25 or larger),
Coriolis accelerations are important to the flow and must be con-
sidered. However, shorter time scales will probably be adequately
represented neglecting rotation with little accuracy loss. During
the data analysis it appeared that friction, even in small amounts,
exerted considerable influence on the characteristics of the vari-
ables and their relationships, and should therefore not be neglected
in representing time dependent motion.

69

Since frictional effects are predominant in the system, the
simplified assumption of a constant coefficient of vertical friction
throughout the layer, is probably one of the most limiting factors
on the model's accuracy. A technique for improving the repre-
sentation of the influence of friction might be the assumption of
constant stress layer (Prandl layer) underlying the Ekman layer.
The solutions in these layers would be matched at their interface;
that is, at some height above the bottom determined as a function
of the amount of friction present and the roughness of the bottom
boundary. This method has been used in air -sea boundary layer
models with considerable improvement in accuracy, and merits con-
sideration at the ocean boundary in this scheme.

The assumption of constant depth does not account for modifi-
cations of the flow caused by changes in the potential vorticity
induced by the shrinking and stretching of the water column over
varying topography. This restriction is probably most serious in
shallow regions where time dependent flow is significant and changes
in the bottom extreme, for example the flow in the vicinity of the
Monterey Submarine Canyon.

Certain large scale motions merit representation and investi-
gation using this formulation. One of these is the study of
seiches, especially the trapped modes, on the continental shelf
since their buildup can cause storm surge (Pierson and Neuman,
1966). Also, as previously mentioned, the large scale oscillations
on the continental shelf might make important contributions to
sediment transfer across the shelf. Therefore it might prove use-
ful to remove the horizontal boundary conditions from the model

70

and drive the flow by an externally specified time dependent
forcing term. This would indicate how effectively the motion was
coupled to these terms and show whether the previously noted trans-
port rectification is of sufficient magnitude to be significant.

An examination of Tables I and II indicates how costly this
type of numerical representation is in terms of computer time. It
was not determined how much loss in resolution would occur by in-
creasing the spacial grid size and reducing the length of the local
time history considered. It is clear that considerable reduction
in computing time would result and that such considerations are
worth investigating before a similar model is used for extensive
production runs.

71

APPENDIX A
DETAILS OF THE FINITE- DIFFERENCE SCHEME

Consider first the calculation of the components of the vector

D., recalling that
3

D. = GiQ*) , + o" - ao" , + AtA

(Al)

where

a =

At
4Ax

0 0 C2/At
0 0 D2/At
-10 0

(A2)

and

on =

A =

D. =
D

/d?

(A3)

(A4)

(A5)

Expanding eq. (Al) using (A2) through (A4) yields

72

At

3 4Ax

0 O C2/At
0 0 D2/At
-10 0

3

vn

3

At

4ax

0

0

C2/At

0

0

D2/At

-1

0

0

«"-!

n

vj-i

(A6)

or

D. =
3

(C2/4Ax)(?^ + 1 - 5^) + U% CI
(D2/4Ax)(§^ + 1,- Sj.x) + v" +

(-At/4Ax)(U^+1- U^) * I1]

(A7)

Now consider the calculation of the matrix

E. =
3

eD eD eD
11 12 13

eD e3 eD
21 22 23

eJ e^ eD
31 32 33

, J ^ 1

(A8)

Recall eq. (49) and express it in inverse form as

Ej = (I + aEj_1)*1 0£, j * 1

(A9)

where I is the (3X3) identity matrix. Using the appropriate

values for I, a, and E. ,, it is easily showed that

3-1 *

73

(I + aEj._1) =

Ate

3-1

11

4ax

Ate

12

4ax

Ate

l -

j-ll

13

4ax

= y

j-i

(A10)

Therefore,

(I + «Ej_1)

-1

DETCy-5- )

AOJCy3"1)

(All)

Where DET(y ) is the determinant of the matrix (y^ ) and
ADJ(y ) is its adjoint. The computation of eq. (All) are ex-
tremely lengthy and are not included. The resulting product of
applying eq. (All) to eq. (A9) is as follows:

Ej " DET

j-l

C2e

33

4Ax

D2e

33

4ax

A

+C2e

j-l'
31

4ax

0

0

C2e

3-1

C2

-

"13

At

4

Ax

D2

-

D2

•S1

At

4

Ax

cz4]

■i \

\

4ax

j * i

(A12)

wher e ,

DET. .
3-1

= 1 -

Z^Ate^-rae^V

C2At j-1 j-1 j-1 j-1

(4AX)2(G33 Gl1 "G31 Gl3 '

(A13)

This completes the calculations of the components of E . .
Finally consider the calculations of the components of the vector

74

F. =
D

(A14)

Recall eq. (50) and express in its inverse form as
F. = (I + «E. , )_1(D.- Q!F. , ) , j > l

(A15)

Again the computations of eq. (All) are involved and are not
included. The results of applying eq. (All) to (A15) are as
follows :

F . = — —
j DET ,

j-l

1 - Ate

u1] Q1" (c2dillQ3

"5 Ax / V~4"ax

lei"1 + D2At(e^'1e^1-e^;1e?71) ^

31

4ax

'31 13
(4Ax)

'e^1 ■ C2At(e^1e^1-e^1e^1)

D2e

"31

ZfAx

j-l
33

4ax

Q3

(4Ax)'

(A16)

, j*l

Ate

j-l
11

k Ax

Ql

1 + C2e

j-l
31

4ax

Q3

where,

C2f

j-l

Ql = d^ -

Q2 = d^ -

Q3 = &>

4ax

d24'

■l

4 Ax
Atf^"

• l

4ax

(A17)

75

Note that the middle terra of the vector F. is the entire quantity
in the brackets.

The implicit scheme for computing U, V, and § can now be
constructed. Writing eq. (46) as

0n+! = E. ,9n+1 + P. _ (A18)

j-1 j-1 j j-1

and substituting the quantities previously derived for E. . and
F . . yields :

«£ - eJ-V+1 ♦ eJ"1 5^ ♦ i{* (A19)

«T-l - 4>T + e^ s-+1 + *2_1 (A20)

-n+l eJ-l„n+l + J-1 f"+l + fJ"l fA21x

5j-l e31 J e33 &j f3 (A X)

Equations (A7) , (A12), (A13), (Al6), (A17) , and (A19) through
(A21) represent the expressions used to implicitly calculate values

for 9. at each grid point in the manner described in the text.

3

76

APPENDIX B

■X- ¥r

DETAILS OP REPRESENTING Cl , C2 , VELOCITY AND BOTTOM STRESS

Consider the representation of Cl by first noting that the
limits of integration have the following meaning:
t = 0 represents the initial condition time.
t = t represents the present time.

t = t + At represents the time at which the value of the
various variables are to be calculated.

The integration of eq. (48) was performed numerically by
dividing the total time into small intervals (At) and evaluating

p*(tM as a constant at the mid-point of each interval. This

*
allowed Cl to be numerically represented as follows:

l h_, _2 h_, [ bxj

2
-kp At -if At
e e - 1

M=l p N=l \ /L

r(I-l-N)At .^2 -ifft-tM

X \ - e *P (t X }e lt(t X }dt' (Bl)

J(I-N)At

With the further assumption that the effect of inertial ro-
tation was constant in each interval, e was approximated

* *v. -,* • -if(I-^-N)At

at the mid-point as e x ' and treated as a constant.

Analytical integration of the remaining integral resulted in

eq. (51). To solve for Cl and Dl the following relation was used,

and the real and imaginary terms collected.

elft = cos (ft) - i sin(ft) (B2)

77

This allowed CI and Dl to be represented as

2

MAXM 1-1

CI1 = Al S -f S A3

M=l p N=l

- cos (f(i-^-N)At}

MAXM 1-1

Dl1 = Al S ^§ £ A3

M=l p N=l

;"kp Atcos{f(I+%-N)At}

a4

-kp At .

(B3)

sin{f(I+^-N)At}

- sin {f(I-^-N)At} A4

(B4)

where

Al

-2o.
h

A2 =

kp'

(1 "

-kp At)

A3 =

)x

N-^

<£>

a4 =

-kp (I-l-N)At

Using the same assumption regarding constant inertial rotation
in each time interval, C2 readily reduced to eq. (52). Applying
eq. (B2) to eq. (52) yields:

MAXM

C2 = Al cos (^f At) E

M=l p

A2
2

(B5)

MAXM

A2

D2 = - Al sin(^fAt) E

M=l p

(B6)

78

Consider now the scheme formulated to recover the velocity
profile. Since velocity was calculated at time t + At, eq. (26)
was represented by:

w

_ 2g ^^ (-l)mcos(pz)(ct e-kp2(t + At-t')e-if-(t + At-f) 6g
h M=l P Oo ' bx

dt

X't+At e-kp2(t + At-f)e-if(t + At-t») b£

s.

bx

( t ' ) dt '

(B7)

Applying the same assumptions and relations to eq. (B7) as

were applied to eq. (Bl) yields:

MAXM / 1-1 r .2 -,

u1 = - Al S ^ A2 S A3 e"1^ At cos{ f (1+%-N) At} A4
L M=l P I N=l

+ cos(^fAt)(A2)(A6)

MAXM - / 1-1
v.1 = - Al S ^2 ■ A2 S A3
L M=l P C N=l

-e"kp At sin{f(I+^-N)At}

A4

- sin(^fAt)(A2)(A6)

)

where

M

A5 = (-1) cos(pz)

A6 = (I1)
v Ox'

i-k

(b8)

(B9)

The term (A5) represents the depth dependence and (A6) the
slope in the most recent time interval. (A6) was calculated and
stored by the model after U, V, and § were computed.

The final quantity needing representation is bottom stress (T)
Recalling eq. (54) and noting that the integral is identical to

79

that of eq. (26), bottom stress is easily represented as:

<Tx»L =

MAXM
kAl 2
M

AXM ( 1-1 [ .2

\k2 S A3 U'^5 cos{f(I+3g-N)At}
= 1 ^ N=l L

a4

■}

cos(%fAt)(A2)(A6)

MAXM / 1-1 w,2A+

(T ) = - kAl S JA2 E A3 -e'^5 sin{f ( I+^-N) At}

y L M=l L

a4

(BIO)

- sin(?§fAt)(A2)(A6)

(Bll)

Where (T ) and (T ) are cross-channel and along-channel components
x y

of stress respectively.

During the calculation of CI, Dl , C2 , D2 , the model stores
those common quantities needed for velocity and stress calculations
as defined by eqs . (B3) through ( B6 ) , and (B8) through (Bll ) . The
program recalls these quantities and applies them to eqs. (B8)
through (Bll) whenever it is required to calculate velocity and
bottom stress.

80

APPENDIX C
COMPUTER PROGRAM DESCRIPTION
The computer program was written in FORTRAN IV and was used
with the IBM 36O/67 computer system at the W.R. Church Computer
Center, Naval Postgraduate School. The program has been generalized
such that basin size, run time, initial free surface elevation, and
number of modes considered in the solution can be varied in the data
statements. Additionally, the model allows variable grid spacing,
requiring only changes in the imput parameters and DIMENSION state-
ments. The influence of inertial rotation and vertical friction are
introduced into the model by the two non dimensional ratios dis-
cussed in the text.

FORTRAN IV symbols used for the major model parameters are
defined below.
A Model time (real time normalized to the free oscillating

seiche period)
DH Spacial step in depth (fixed at one meter)
DX Spacial step in across -channel direction
DT Time step
ETA Free surface elevation
ETAX Sea surface slope
ETAXS Array for storing the time history of the sea surface

slopes
F Coriolis parameter
FREQ Frequency of the free oscillating seiche

81

G Acceleration of gravity

H Depth of water

HI Free surface elevation at left boundary

MAXH Integer representation of water depth (must be even)

MAXL Maximum number of cross-channel grid points (must be odd)

MAXP Interval in time step when velocity, energy, and bottom

stress is calculated

MAXT Maximum number of time steps

PCD Ratio of Ekman depth to depth of the water

PCF Ratio of Coriolis parameter to depth of the water

PCT Number of time steps per free wave period

PE Potential energy of the seiche

PERD Period of a free oscillating seiche

R Coefficient of vertical friction

RHO Density of the water

SE Kinetic energy associated with the y-component of velocity

SIGMA Angular frequency of a free oscillating seiche

STRESC y-component of bottom stress

STRESR x-component of bottom stress

j Real time

TE Kinetic energy associated with the x-component of velocity

TOT Total energy balance of the seiche

UT x-component of volume transport

VELC y-component of velocity

VELR x-component of velocity

VT y-component of volume transport

82

The procedure for initialization of the program is as follows:

1) Modify the initial DATA statements to correspond to the
parameters of the motion being represented.

2) Change DIMENSION statements accordingly.

3) Insure that appropriate storage and computer time require-
ments are satisfied. The program required a maximum of 100K BYTES
storage on the IBM 360/67 system and variable CPU usage time as
depicted in Tables I and II.

4) A minor modification to the program is necessary when
frictionless motions are being represented. If (PCD = 0) or

(PCF = 0), remove the following three expressions from the program

K = (1-1) - 10.0/CC

IF(K.LE.0)K = 1

DO 40 J = K, IM1
Replace these expressions by

DO 40 J = 1, IM1

5) Because of the formulation, the program can not compute
velocity at I = 1; therefore a value of MAXP = 1 should not be
used. To recover the velocity profile at each time step (I j- 1),
initialize the DATA with MAXP = 2 and insert the statement
(MAXP = 1) immediately before statement number 300.

The program, as appearing in the next section, is initialized
in the MKS system and changes in the imput should be made ac-
cordingly. The major computational sections of the program are
prefaced by appropriate comment statements, and a reader with
FORTRAN IV experience should have little difficulty following the
program.

83

COMPUTER PROGRAM

DIMENSION ETAK11) .STA2( 11) .ETAX( 11 ) .ETAXKlDi
lETAXSd 1,201 )t UT1(11),UT2( 11 ) ,VTH 11 ), VELR ( 11, 31) t
2VELC(11 ,31) ,FUNR(11 ,2 5 ) ,FUNC(11 ,25 ) , P ( 25 ) , P2 ( 25 ) ,
3P3(25).P4(25)tPK(25).3(31).C(31),TIN(ll),PIN(ll),
4P0T(11) ,E11 (10 ),i21 (10 ) ,E31(1C ),cl3(10), E23( 10),
5E33( 10) ,D1( 10) ,D2( 10) ,D3(10) ,F1 ( 10) ,F2 (10 ) .F3 (10 ) .
6DET(10),CR1( 11),CC1(11),STRESR(11) , STRE SC ( 11)

C INPUT INITIAL DATA

C

DATA MAXT,MAXL , MAXH , MAXP, MAXM/ 200, 11 , 30, 5 , 25/

DATA WIDTH, HI , PC T, DH , RH 0/5000. .1 . #40. tl • •1027./

CAT A PIE, G,PCF, PCD/ 3. 1415927, 9. 8, .5,0.0/
C
C ECHO CHECK DATA

WRITE(6T310)

WRITE (6, 320) MA XT ,MAXL, VAXH, MAXP ,M AXM

WRITE (6, 330 ) WI DTH, HI ,DH •RHO. PI E ,G .PCT.PCF.FCD
C

DX=WIDTH/FL0AT(MAXL-1 )

H=FLOAT(MAXH)

PERD=(2.0*WIDTH)/SQRT(G*H)

FREQ=1.0/P£RD

SIGMA=2.0*PIE*FREQ

F=PCF*SIGMA

R=(PCD*H)**2*F/( 2.0*PIE**2)

CT = PERD/PCT

Q1=DT/(4.0*DX)

Q2=4.0*DX

C3=-2.C*G/H

WRITE(6,350)H,WIDTH,DT,DX,PERD,SIGMA,F,R
C

C COMPUTE AND STORE INITIAL CONDITIONS ON ETA ,E TAX, UT, VT

C

DO 10 L=1,MAXL

X=(L-1)*DX

ETA1(L)=HI*C0S(PIE*X/WIDTH)

ETAX(L) = (-1.C*PI5*HI/WIDTH)*SIN(PIE*X/WIDTH)

UT1(L)=0.0

VTKL) =0.0
10 CONTINUE

A = 0.0

ETAX(MAXL)=0.0
C

C PRINT INITIAL CONDITIONS

C

WRITE(6,360)A

WRIT £(6,370) (ETAKL ),L = 1,MAXL)

WRITE ( 6,3 80) (ETA X(L) , L = l , MAX L)

WRIT£(6,400)(UT1(L ),L=1,MAXL)

WRITE (6, 40 5) (VT1 (L) ,L=1,MAXL)
C
C CALCULATE AND STORE MODAL CONSTANTS

C

DO 20 M=1,MAXM

P(M)=(2.0*FL0AT(M)-1.0)*PIE/(2.0*H)

P2(M)=P(M)**2

PK(M)=P2(M)*R

P4(M)=EXP(-PK(M)*DT)

C5=PK(M)

C4=C5*DT

IF(C4.LT.1.E-4)G0 TO 15

P3(M)=(1.0-EXP(-C4) )/C5

GO TO 20
15 P3(M)=(DT*(1.0-(DT*C5/2.0)*( 1 .0-C5*DT/3.0 ) ))
20 CONTINUE

84

c

C CALCULATE CR2 ANC CC2

C

FDT=p*DT

C1=C0S( . 5*FDT)

C2=SIN( .5*FDT1

CR2=0.C

DO 30 M=1,MAXM

CR2=CR2+P3(M)/P2(M )

30 CONTINUE
CC2=-C3*C2*CR2
CR2=C3*C1*CR2

C

C BEGIN TIME LOOP

C

DO 300 I=1,MAXT

T=I*DT

A=T/PERD

IFd.NE.DGO TO 35

DO 31 L=1,MAXL

CR1(L)=0.0

CC1(L)=0.0

31 CONTINUE
C

C SET LEFT BCUNDRY VALUES FOR MATRICES

C

E11(1)=0.0

E21(l)=0.0

E31(l)=0.0

E13(l)=0.0

E23( 1)=0.0

E33(l)=1.0/.951

Fl(l)=C.O

F2( 1)=0.0

F3(l)=0.0

DET(l) =1.0

GO TO 6 5
C

C (I.NE.l) CALCULATE CR1 AND CC1
C

35 CO 60 L=1,MAXL

C6=0. 0

C7=0.0

DO 50 M=1,MAXM

FACR1=0.0

FACR2=0.0

FACC1=0.0

FACC2=0.0

CC=PK(M)*DT

IM1=I-1

K=(I-1)-10.0/CC

IP(K.LE.0)K=1

DO 40 J=K,IM1

C4=ETAXS(L,J)

C5=EXP(-CC*(I-1-J) )

FACRl=FACRl+C4*C0S(FDT*U + .5-J) )*C5

FACR2=FACR2-C4*C0S( FDT*< I-. 5-J) )*C5

FACC1=FACC1-C4*SIN(FDT*( I+.5-J ) )*C5

FACC 2 = FACC 2+C4* SI N < FDT* ( I-. 5-J ) )#C5
40 CONTINUE

FACR1=FACR1*P4(M)

FACC1=FACC1*P4(M)

C5=P3(M)

FUNR(L,M)=FACR1*C5

FUNC(L»M)=FACC 1*C5

C8=1.0/P2(M)

C6=C6+(FACR1+FACR2 )*C5*C8

C7=C7+(FACC1+FACC2)*C5*C8
50 CONTINUE

CR1(L)=C3*C6

CCKL)=C3*C7
60 CONTINUE

85

c

C CALCULATE VALUES OF MATRICIES

C

65 MM1=MAXL-1

DO 70 L=2,MM1

DHL) = (CR2/Q2)*(ETA1(L+1)-ETA1 (L-l) )+UTl (LI+CRKL)

D2(L) = (CC2/Q2)*(ETAHL + 1HETA1(L-1 ) )+VTl ( L ) +CC1 ( L )

D3(L)=-Q1*(UT1 (L+l )-UTl (L-l ) )+ETAl (L )

Q3=Q1/DET(L-1)

Ell (L)=Q3*CP 2*533 (L-D/Q2

E13(L)=Q3*(CR2/DT-CR2*E13(L-1 )/Q2)

E2KL )=Q3*CC2*£33(L-1)/Q2

E23(L)=Q3*(CC2/DT-CC2*E13(L-1)/Q2)

E31(L)=-03M 1.0+CR2*E31(L-1)/Q2)

E33(L )=Q3*CR2*cll(L-l) /Q2

DET(L)=1.0-(1.0/Q2 )* ( DT*E1 3 ( L )-CR2*E31 (L ) )+(CR2*DT/Q2
1**2)*(E33(L)',E11(L )-c31(L)*E13(D)

Q4=D1(L)-CR2*F3(L-1 )/0 2

05=D2(L)-CC2*F3(L-1 )/Q2

06=D3(L )+DT*FKL-l)/Q2

Fl (L)=( (l.C-DT*E13 (L-1)/Q2)*Q4-(CR2*E33(L-1)/Q2)*Q6)/
lDET(L-l)

C4=(-CC2*E31(L-1 )/02 + CC2*DT*(E31(L-l)*E13(L-D-
1E33(L-1)*E11(L-1 ) ) /Q2**2)*Q4

C5=(1.0+(CR2*E31(L-1)-0T*E13(L-1) ) /02-CR2*DT*( E31 ( L-l )
1*E13CL-1)-E33(L-1)*E11(L-1) )/Q2**2)*Q5

C6=(CC2*E33(L-1) /Q2)*Q6

F2(L)=(C4+C5-C6)/D2T(L-1)

F3(L)=( (DT*cll (L-l )/Q2)*Q4+(1.0+CR2*E31(L-l)/Q2)*Q6)/
lDET(L-l)
70 CONTINUE
C

C COMPUTE UT2, VT1, ETA2

C

ETA2(MM1)=F3(MM1 )/(1.0-E33 (MM1 J/.951)

ETA2(MAXL)=ETA2(MM1)/.951

UT2(MM1 >=E13 (MM1 )*ETA2 (MAXL ) +FKMM1 )

VTKMM1) =E23(MM1 )*ETA2 (MAXD+F2 (MM1)

UT2(MAXL )=0.0

VT1(MAXL)=0.0

MM3=MAXL-3

DO 80 L=1.MM3

K=^AXL-L

C4=UT2(M)

C5=ETA2(M)

ETA2(M-1 )=E31(M-1)*C4+E33(M-1 )*C5+F3(M-1 )

UT2(M-1 )=£11 (M-l )*C4+E13(M-1)*C5+F1(M-1)

VT1(M-1)=E21(M-1 )*C4+E23(M-1)*C5+F2(M-1 )
80 CONTINUE

ETA2U )=£TA2(2)/.951

UT2(1)=0.0

VT1(1)=0.0
C
C COMPUTE ETAX1

C

DO 85 L=2,MM1
ETAX1(L)=(ETA2(L+1)-ETA2(L-1) )/(2.0*DX)

85 CONTINUE

ETAX1(1)=0.0

ETAX1(MAXL)=0.0
C

C RESET ETAXS ARRAY AND RESET ETA1 , ETAX, AND UT1

C

DO 90 L=1,MAXL

ETAXS(L,I)=(ETAX1(L)+ETAX(L) )/2.0

UT1(L)=UT2(L )

ETA1(L)=£TA2(L )

ETAX(L)=ETAX1(L)
90 CONTINUE
C

C PRINT OUT UT1, VT1, ETA1, AND ETAX

C

86

WRITE(6,410 )A

WRITE (6,420) (UTKL) tL=l,MAXL)

WRITE(6,421XVT1(L )tL-l,MAXL)

WRITE (6, 42 5) ( ET Al ( L ) , L = l , MAXL )

WRITE( 6,429) ( E TA X( L ) . L=l . MA XL )
C

IF(MOD(I ,MAXP) .EQ.O)GO TO 100

GO TO 300
C
C I IS A MULTIPLE OF MAXP, CALCULATE VELOCITY

C

100 CO 115 L=1»MAXL
C6=ETAXS(L,I )
DO 110 M=1,MAXM

FUNR(L,M)=FUNR(L»M)+C1*P3(M)*C6
FUNCC L,M)=FUNC(L,M)-C2*P3(M)*C6

110

CONTINUE

115

CONTINUE

DO 140 L=1,MAXL

MH1=MAXH+1

DO 130 J=1,MH1

Z=(J-1)*DH

SI=-1.0

C4=0.0

C5=0.0

DO 120 M=1TMAXM

C6=P(M)

C7=C0S(C6*Z) /C6

C4=C4+SI*FUNR(L,M)*C7

C5=C5+SI*FUNC< L,M)*C7

SI=-1.0*SI

120

CONTINUE

V£LR(L,J)=-C3*C4
VELC(L, J)=-C3*C5

130

CONTINUE

c
c
c

140

CONTINUE

CALCULATE BCTTOM STRESS

C6=C3*R

DC 142 L=l ,MAXL

C4=0.0

C5=0.0

DO 141 M=1,MAXM

C4=C4+FUNR(L,M)
C5=C5+FUNC(L,M )

141 CONTINUE
STRESR(L)=C4*C6
STRESC(L) = C5*C6

142 CONTINUE
C

C CALCULATE ENERGY BALANCE USING SIMPSON1 S RULE

C

145 DO 165 L=ltMAXL
DO 150 J=1,MH1
B( J)=VELR(L,J)**2
C( J)=VELC(L,J)**2

150 CONTINUE

EVSUM1=0.0

EVSUM 2=0.0

0DSUM1 = 0.0

0DSUM2=0.0

ENSUM1=B( 1)+B( MH1)

ENSUM2=C(1 )+C(MHl )

MH2=MH1-1

DO 155 J=2,MH2,2

EVSUM1=EVSUM1+B( J)

EVSUM2=E VSUM2+C( J)

155 CONTINUE
^H2=MHl-2
DO 160 J=3,MH2,2
0DSUM1=0DSUM1+B(J )

87

c
c
c

c
c
c

0DSUM2=0DSUM2+C( J)
160 CONTINUE

TIN(L)=0H/3.0*(ENSUM1+2.0*0DSUM1+4.0*EVSUM1)

PIN(L)=DH/3.0* (ENSUM2+2.0*0DSUM2+4.0*£VSUM2)

P0T(L)=ETA1(L)**2
165 CONTINUE

EVSUM1=0.0

EVSUM 2=0.0

EVSUM3=0.0

0DSUM1=0.0

0DSUM2=0.0

0CSUM3=0.0

ENSUM1=P0T(1)+P0T( MAXL)

ENSUM2=TIN( 1 )+TIN(MAXL )

ENSUM3 = PIN(1 )+PIN(MAXL)

DO 170 L=2,MM1,2

EVSUM1=EVSUM1+P0T(L)

EVSUM2=EVSUM2 + TIN( L)

EVSUM3=EVSUM3+PIN( L)
170 CONTINUE

MM2=MAXL-2

DO 175 L=3,MM2,2

0DSUM1=0DSUM1+P0T( L)

0DSUM2=0DSUM2+TI N( L)

0DSUM3=0DSUM3+PIN(L)
175 CONTINUE

PE=.5*RH0*G*(DX/3. 0)*(ENSUM1+ 2 . 0*0DSUM1 +4 .0*EVSUM1 )

TE=.5*RH0*(DX/3.0)*(£NSUM2+2.0*ODSUM2+4.O*EVSUM21

SE=. 5*RH0* (DX/3.0)*(ENSUM3+2.0#0DSUM3+4.0*£VSUM3)

TOT=P£+TE+SE

PRINT VELOCITY

WRITE(6,430 )

DO 180 J=1,MH1

JM1=J-1

WRITE (6, 440) JM1, ( V5LR ( L, J ) , L=l , MAXL )
180 CONTINUE

WRITE(6,435)

DO 190 J=1,MH1

JM1=J-1

WRITE(6,440)JM1, (VELC(L,J) ,L = 1,MAXL)
190 CONTINUE

PRINT PE,TE,SE,TOT,STRESR, AND STRESC

WRI TE< 6,45 0) PE ,TE, SE,TOT
WRITSC6. 45 5) ( STRESR(L ) ,L=1, MAXL)
WRITE (6 ,460) (STRESC (L) ,L=1,MAXL )

300 CONTINUE

STOP

310 FORM AT ('l1

320 FORMAT(»0»

330 FORMAT( '0'

350 FORMAT CI1

l'FGLLOWSS/aX,' D
2'WIDTH=',F7. 1, IX,

•ECHO CHECK OF INPUT DATA',/)

» s
.9

360
]
370

3 80
400
405
410

4 20
4 21
425
429

31X,,DX=« ,F6.1, IX
4/,lX,« SIGMA=« ,E11

11. 3)

•SOLUTION TO

• nsPTH=i tF5

•METERS' ,

SEICHE

FORMAT (
F4.1 ,1X
FORMAT(
FORMAT (
FORMAT (
FORM£T(
FORMAT (
FORMAT(
FORMAT (
FORMAT(
FORMAT(

PROBLEM
,., "METERS
1X,'DT = ' »F5
•METERS •,/, IX, «PERD=»,
.3,/,lX,lF=« , Ell. 3,/, IX, •
CONDITIONS', /,56X,

WITH CORIOLI S*

,/, IX,
1.1X.1

SEC ,/

55X, • INITIAL
•SECONDS' )

•ETA'/( 11E11.3) )

•ETAX'/( 11E11.3) )
55X,'TRANSPCRT(UT)'/ (11 El 1.3)

,55X
t55X,

55X
55X

• TRANSPORT* VT) •

NORMALIZED T IME= ' , F8 .2, / )

/(11E11
i

3) )

R=',E11.3,//
•TIME=I ,

,55X,' TRANSPORT (UT) »/ (11511.3) )
,55X, • TRANSPORT! VT) •/( 11E11.3) )
,55X v'ETAa/(llE11.3) )
,55X, •ETAX,/(llcll.3) )

88

430 FORMAT CO
435 FORMATCO
440 FORMAT* *
450 FORMAT ( '0
1E15.7,8X,
455 FORMAT( «0
460 FORMATCO
END

•VELOCITY! VELR) »,/, 3X, 'Z'

t «Z«
2.( 11=11.3) )
■t 2Xt E15«7i 8X,
,E15.7)
•STRSSS(R) ■

5 5X
f55X,« VELO

» 2X.I
, »PE=

TOT='
55X

»/)
,/)

•TE=,,E15.7,8X, «SE=

55X, 'STRESS(C) •

11F11.3) )
11E11.3))

89

BIBLIOGRAPHY

Bryan, K. , "A Numerical Investigation of Certain Features of the
General Circulation," Tellus . v. 11, no. 2, p. I63-I74,
March, 1959.

Gait, J. A., "Linear Transients Associated with a Time Dependent
Wind Over Shallow Water," Unpublished Paper, Naval Postgraduate
School, 14 pp. , 1970.

Hamming, R.G., Numerical Methods for Scientists and Engineers,
McGraw-Hill, 1962 .

McCalla, T.R., Introduction to Numerical Methods and FORTRAN
Programming , p. 270-271, Wiley and Sons, 1967.

Morton, K.W. and Richtmyer , R.D., Difference Methods for Initial-
Value Problems . p. 198-201, Inter science, 1957.

Munk, W.H., "On the Wind-Driven Ocean Circulation," Jour . of
Meter ology, v. 7, no. 2, p. 79-93, 1950.

Neuman , G. and Pierson, W.J., Principles of Physical Oceanography,
p. 374-375, Prentice-Hall, 1966.

Proudman, J., Dynamical Oceanography, p. 243-245, Dover, 1952.

Smith, G.D., Numerical Solution of Partial Differential Equations,
p. 70-71, Oxford University, 1965.

Sverdrup, H.U., "Wind-Driven Currents in a Baraclinic Ocean; with
Application to the Equatorial Currents of the Eastern Pacific,"
Proc. Nat. Acad. Sci., v. 33, no. 11, p. 318-326 , 1947.

Sverdrup, H.U., Johnson, M.W., and Fleming, R.H., The Oceans ,
pp. 1087, Prentice-Hall, 1942.

Welander , P., "Wind Action over Shallow Sea: Some Generalizations
of Ekman's Theory," Tellus , v. 9, no. 1, p. 45-52, 1957.

Welander, P., "Numerical Prediction of Storm Surges," Adv. in
Geogpy. , v. 8, p. 315-379, 1961.

90

INITIAL DISTRIBUTION LIST

No . Cop i es

1. Defense Documentation Center 2
Cameron Station

Alexandria, Virginia 22314

2. Library, Code 0212 2
Naval Postgraduate School

Monterey, California 93940

3. Department of Oceanography 3
Naval Postgraduate School

Monterey, California 93940

4. Oceanographer of the Navy 1
The Madison Building

732 N. Washington Street
Alexandria, Virginia 22314

5. Dr. Ned A. Ostenso 1
Code 480D

Office of Naval Research
Arlington, Virginia 22217

6. Asst. Prof. J. A. Gait, Code 58G1 4
Department of Oceanography

Naval Postgraduate School
Monterey, California 93940

7. Asst. Prof. J.J. vonSchwind, Code 58Vs 1
Department of Oceanography

Naval Postgraduate School
Monterey, California 93940

8. Asst. Prof. K.L. Davidson, Code 5lDs 1
Department of Meter ology

Naval Postgraduate School
Monterey, California 93940

9. Asst. Prof. R.H. Franke, Code 53 Fe 1
Department of Mathematics

Naval Postgraduate School
Monterey, California 93940

10. Asst. Prof. R.S. Andrews, Code 58Ad 1
Department of Oceanography

Naval Postgraduate School
Monterey, California 93940

11. LCDR Norman T. Camp, USN 1
USS OPPORTUNE (ARS 4l)

FPO, New York, N.Y. 09501

91

Unclassified

Security Classification

DOCUMENT CONTROL DATA -R&D

(Security clas si lication ol title, body of abstract and indexing annotation must be entered when the overall report Is classified)

I . ORIGINATING ACTIVITY (Corporate author)

Naval Postgraduate School
Monterey, California 93940

2l. REPORT SECURITY CLASSIFICATION

Unclassified

2b. GROUP

3 REPOR T TITLE

An Investigation of Linear Transients Associated With a Time Dependent
Bottom Spiral

4. DESCRIPTIVE NOTES (Type ol report and, inclusive dates)

Master's Thesis; March 1972

5. authORiSi (Fifsi nams, middle initial, last name)

Norman Thomas Camp

6- REPOR T D A TE

March 1972

7«. TOTAL NO. OF PAGES

93

76. NO. OF REFS

13

Sa. CONTRACT OR GRANT NO.

6. PROJEC T NO.

9a. ORIGINATOR'S REPORT NUMBERIS)

Ob. OTHER REPORT NO(S) (Any other numbers that may be assigned
this report)

10. DISTRIBUTION STATEMENT

Approved for public release; distribution unlimited,

II. SUPPLEMENTARY NOTES

12. SPONSORING MILI TAR Y ACTIVITY

Naval Postgraduate School
Monterey, California 93940

13. ABSTRAC T

Ekman's linear equations for time dependent flow (neglecting wind stress)
are solved using a time dependent Green's function and the method suggested
by Welander (1957). The solution represents the vertical velocity profile
in terms of the local time history of the changes in sea surface elevation
determined by the divergence of the flow in the vertically integrated
continuity equation.

A fully implicit finite-difference scheme is developed to represent a time
dependent seiche oscillating across a shallow infinite channel. The transients
associated with the formation of the bottom spiral are clearly represented by
the model and the influence of friction and Coriolis are individually and
collectively introduced. The model allows independent calculation of velocity,
volume transport, sea surface elevation, bottom stress, and the total energy
balance of the system. The numerical scheme provides a method for adequately
describing and investigating certain classes of time dependent motion, and
its development suggests a mechanism for improving numerical prediction of
storm surge.

DD

FOR- I47O

1 NOV 68 ■ "T / Sj

S/N 0101 -807-681 1

(PAGE 1 )

92

Unclassified

Security Classification

A-31408

Unclassified

Security Classification

key wo R OS

Ekraan dynamics

Numerical model

Seiche

Shallow water flow

Storm surge

Time dependent motion

DD,F.T..1473 <back

S/N 01 Ot -807-682 I

93

Unclassified

Security Classification

A- 3 t 409

134471

Thes

C1932

c.l

Camp .

An investigation of

linear transients asso-
ciated with a time de-
pendent bottom spiral.

Thesi s

C1932

c.l

Camp 13 £ i 71

An investigation of
linear transients asso-
ciated with a time de-
pendent bottom spiral

thesC 1932

An investigation of linear transients as

3 2768 002 08469 1

DUDLEY KNOX LIBRARY