Numerical modeling of wave-induced currents using a parabolic wave equation

UFL/COEL-TR/080

NUMERICAL MODELING OF WAVE-INDUCED
CURRENTS USING A PARABOLIC WAVE EQUATION

by

Harley Stanford Winer

Dissertation

1988

NUMERICAL MODELING OF WAVE-INDUCED CURRENTS USING A PARABOLIC
WAVE EQUATION

By

HARLEY STANFORD WINER

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

1988

ACKNOWLEDGEMENTS

The writing of a dissertation is a long process involving interaction with and support

from many people. The author would like to thank all those who in one way or another

offered encouragement, support, advice and helped in the problem solving process. The

contributions of those who are not mentioned by name are in no way diminished by the lack

of individual acknowledgement.

First, the author wishes to express deep appreciation to the members of his committee:

the chairman, Dr. Hsiang Wang, who enticed him to come to Florida and enroll in the

doctoral program, and the cochairman, Dr. James Kirby, who shared his knowledge of wave

equations and numerical modeling, along with Dr. Robert Dean, Dr. Bent Christensen and

Dr. Joseph Hammack. Each of these men has been both a teacher and a friend. The author

wishes to also thank Dr. Peter Sheng who at one time served on the supervisory committee

but due to a prior commitment was unable to be at the defense and officially be on the

committee.

The author would like to thank the staff of the Coastal and Oceanographic Engineering

Laboratory at the University of Florida and the staff of the department of Coastal and

Oceanographic Engineering who helped and aided the author in innumerable ways. Thanks

are also extended to Ms. Lillean Pieter who helped prepare the figures and Dr. Claudio

Neves, a fellow student, who manipulated the IATEX software so as to produce a document

in thesis mode.

A special acknowledgement is given to the author's parents who encouraged and pushed

him to complete this doctorate. The author would also like to thank his teachers, friends

and fellow students who offered encouragement and interesting discussions.

Lastly the author thanks his soul-mate and partner, Esther DeJong, who managed to

be patient with each extension of the final date. Baby cakes, it's Phinally Done. I love you.

**4

TABLE OF CONTENTS

AtKNOWLEDGEMENTS

LIST OF FIGURES ...................................

ABSTRACT ............

CHAPTERS

1 INTRODUCTION .......

1.1 Literature .........

1.2 General Description of the

2 CIRCULATION MODEL .

2.1 Governing Equations ...

2.2 Radiation Stresses ....

3 WAVE MODEL ........

3.1 Governing Wave Equation

Model

. .. *

. .*

....

3.2 A Parabolic Approximation .......

3.3 Energy Dissipation Due to Wave Breaking

4 FINITE DIFFERENCE SCHEMES AND BOUNDARY CONDI

4.1 Numerical Solution of the Governing Equations .......

4.2 Finite Differencing of the Parabolic Wave Equation .....

4.3 Boundary Conditions in the Circulation Model .......

4.4 Boundary Conditions in the Wave Model ..........

5 RESULTS .........................

5.1 Wave Set-up.............................

5.2 Longshore Currents .......................

TIONS

.o...

o. ..

. .

. 68

. . . . . . ii

............................

..............

..............

..............

..............

..............
..............
..............

..............

..............

5.3 Waves On a Rip-Channel ...................

5.4 Shore-Perpendicular Breakwater . . . .

5.5 Shore Parallel Emerged Breakwater . . . .

5.6 Comparison with Experimental Results of Gourlay . .

6 CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER

6.1 Improved Wave Formulations .................

6.2 Including Wave Reflection ...................

6.3 Extending the Model to Include Sediment Transport ....

APPENDICES

STUDY

... 70

. .. 76

. .. 88

. 92

. 103

... 104

. 105

. 105

A DERIVATION OF THE DEPTH AVERAGED EQUATIONS ..........

B DERIVATION OF THE RADIATION STRESS TERMS .............

C LAGRANGIAN DERIVATION OF THE GOVERNING WAVE EQUATION .

D OBTAINING EQUATION (3.79) .........................

BIBLIOGRAPHY .. .................................

BIOGRAPHICAL SKETCH ...............................

LIST OF FIGURES

1.1 Computational domain of the model. ....................... 8

1.2 Flow chart of computer program. .................. ... 11

1.3 Flow chart of circulation model . ... ... ....... .. 13

1.4 Flow chart of wave model. .......................... 14

2.1 Longshore current profiles for different friction formulations for an 11
second, 1.028 meter wave at 17.28 degrees to a 1:25 plane beach. . 22

2.2 Experimental results for wave set-up and set-down. Reprinted from
Bowen, Inman, and Simmons (1968), Journal of Geophysical Research,
vol. 73(8) page 2573, copyright by the American Geophysical Union. .. 30

4.1 Definition sketch of grid nomenclature and coordinate system. . 56

4.2 Start-up function according to Eq. 4.43 with C1=30 and C2=3. . 60

4.3 Waves passing through a lateral boundary unaffected by the boundary. 61

5.1 Comparison of the numerical solution of Vemulakonda et al. (1985) for
set-up with experimental data of Bowen (1968). Source: Vemulakonda
et al. (1985) ................ .............. 65

5.2 Comparison of numerical results from the model with experimental re-
sults of Hansen and Svendsen (1979). Wave set-up and wave heights for
a 2 second wave on a beach slope of 1:34 and incident wave height of 7
centimeters ................... .............. 67

5.3 Non-dimensional longshore velocities vresus the non-dimensional dis-
tance offshore from the analytical solution of Longuet-Higgins. Reprinted
from Longuet-Higgins (1970b), Journal of Geophysical Research, vol.
75(33) page 6793, copyright by the American Geophysical Union. . 68

5.4 Comparison of experimental results of Visser (1982) and numerical result
for longshore currents due to a 1.1 second wave on a slope of 1:20 for
S'O ,-*r, ponditir- -- oA and go = 20 degrees. ...... ..69

5.5 Top: Depth contours, in meters, using Eq. 5.7 with the coordinates as
shown. Bottom: Depth contours as used in the present model. Grid
spacing is 5 meters. ............... ........ ... 71

5.6 Velocity vectors from the results of Ebersole and Dalrymple for a wave
angle of 30 degrees. Scale for velocity vectors is 1 inch = 2.87 m/sec.
Source: Ebersole and Dalrymple (1979). . ... .. ...... 72

5.7 Velocity vectors obtained from the present model using the same input
conditions as for Ebersole and Dalrymple. ........... ....... 73

5.8 Mean water surface elevation contours showing that the water surface is
lower above the trough. Values indicated are in millimeters. ...... 74

; 5.9 Depth contours in meters for the case of a bar with a gap. ...... 75

5.10 Velocity vectors for the case of waves over a bar with a gap. Dashed
lines indicate depth contours. Grid spacing equal to 10 meters. . 75

5.11 Depth contours for the case of normally incident waves upon a shoal.
Grid spacing is .05 meters. Contour interval is .02 meters. . ... 77

5.12 Velocity vectors for the case of waves breaking upon a shoal. . 78

5.13 Maximum velocity across the top of a shoal versus incident wave height
for waves over a shoal on a plane beach. Incident wave height less than
9.1 centimeters are non-breaking waves. Incident wave heights greater
than 9.1 centimeters produced waves breaking over the shoal. . 79

5.14 Position of the groin in relation to the grid rows and columns. . 80

5.15 Results for the numerical model of Liu and Mei (1975). Deepwater
wave height equal to 1 meter, deepwater wave angle equal to 45 degrees.
Source: Liu and Mei (1975). ..................... . 83

5.16 Results from the present numerical model using the same conditions as
were used by Liu and Mei. Offshore boundary condition of constant wa-
ter surface level which precludes tangential flow at the offshore boundary. 84

5.17 Results from the present numerical model using the same conditions as
were used by Liu and Mei. Offshore boundary condition of no on-offshore
flow. Grid spacing equal to 20 meters. . . . .... 85

5.18 Experimental set-up for the tests conducted by Hales (1980). Source:
Hales (1980). ................... .............. 86

5.19 Experimental measurements of the mid-depth averaged velocities obtain
by Hales. Measurements at one foot intervals. Source: Hales (1980). 87

5.20 Rip-current velocities adjacent to the down wave side of the jetty. Com-
parison of experimental and numerical results. . . ... 88

5.21 Vel.,:.y *. c4: -- .- .: .- I-,... j for the same con-
ditios as tor the Hales experiment. Gid ,pacing equal to one foot. .89

5.22 Definition sketch of grids for the case of a shore parallel breakwater. .. 90

5.23 Wave crests in the lee of a semi-infinite breakwater in a region of constant
depth . . .. . .. . . . . . 91

5.24 Velocity vectors obtained from the model for a 100 meter long breakwater
100 meters offshore for a 1 meter, 6 second wave for increasing angles of
incidence. ................... ... .. .. ....... 93

5.25 Stream function contour lines obtained from the Liu model for a 10
second 1 meter wave normally approaching a 700 meter long breakwater
350 meters offshore on a 1 on 50 slope. Source: Liu and Mei (1975). 94

4 5.26 Flow lines obtained from the present model for a 10 second 1 meter wave
normally approaching a 700 meter long breakwater 350 meters offshore
on a 1 on 50 slope ................... ........... 95

5.27 Velocity vectors obtained from the present model for a 10 second 1 meter
wave normally approaching a 700 meter long breakwater 350 meters
offshore on a 1 on 50 slope. ......................... 95

5.28 Set-up contour lines obtained from the present model for a 10 second
1 meter wave normally approaching a 700 meter long breakwater 350
meters offshore on a 1 on 50 slope. ...... . . ....... 96

5.29 Physical layout for the experiment of Gourlay. Source: Gourlay (1974),
Proceedings 14th Coastal Engineering Conference, copyright by the Amer-
ican Society of Civil Engineers, reprinted with permission. . ... 97

5.30 Computational domain for the experiment of Gourlay. Grid spacing
equal to .1 meter.................. .......... 98

5.31 Comparison of experimental and numerical results for the experiment of
Gourlay: Wave Set-up Contours. Top figure is from Gourlay (1974) Pro-
ceedings 14th Coastal Engineering Conference, copyright by the Amer-
ican Society of Civil Engineers, reprinted with permission. The bottom
figure is the computer results . . . . . 100

5.32 Results from the experiment of Gourlay: Contours of current veloci-
ties and streamlines. Source: Gourlay (1974), Proceedings 14th Coastal
Engineering Conference, copyright by the American Society of Civil En-
gineers, reprinted with permission. . . . ... .. .101

5.33 Results from the numerical model: Streamlines of the flow. ...... 101

5.34 Results from the numerical model: Contour lines of the current veloci-
ties. ...................... .... .. ........ 102

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

NUMERICAL MODELING OF WAVE-INDUCED CURRENTS USING A PARABOLIC
T. WAVE EQUATION

By

HARLEY STANFORD WINER

August 1988
Chairman: Dr. Hsiang Wang
Major Department: Aerospace Engineering, Mechanics and Engineering Science

Prediction of waves and currents in the nearshore region in the presence of proposed

breakwaters is needed in the design process in order to optimize the placement of the

structures. A numerical model is developed which solves for the wave field, the depth

averaged mean currents and the mean water surface elevation within a given computa-

tional domain with the presence of breakwaters. The model iterates between a wave model

and a circulation model. The wave model uses a parabolic approximation to a combined

refraction-diffraction mild slope equation which includes currents. The circulation model

uses an alternating direction implicit method solution to the depth averaged equations of

momentum using the radiation stress terms obtained through solution of the wave model as

the forcing terms. Each of the terms in the depth averaged equations of momentum involve

assumptions that are discussed in the report. The model is constructed in modular form so

that improvements can easily be incorporated into the model.

The model was employed for several test cases such as rip channels, bar gaps, and

both shore perpendicular and shore parallel breakwaters. Good agreement was obtained

when the model was used to simulate a i: ^ '

("Wave set-up and wave generated currents in the lee of a breakwater or headland,"

Proc. 14th Coastal Eng. Conf., Copenhagen, 1974, pp. 1976-1995).

Suggestions are made for improving the model and for including sediment transport in

the model.

*1

CHAPTER 1
INTRODUCTION

As waves propagate they are transformed by variations in currents and depths, along

with interaction with obstructions. The major processes of wave transformation include

shoaling, reflection, refraction, diffraction and dissipation. Shoaling is the change of wave

energy, propagation speed and wave height due to changes in the water depth and/or current

field. Refraction is the change in the direction of wave propagation due to the spatial

gradient of the bottom contours and/or current field. Diffraction is the phenomenon of

wave energy being transmitted laterally along a wave crest. Dissipation is the process of

energy loss due to a multitude of possible causes such as bottom percolation, damping

due to viscous effects, turbulence, and wave breaking. As waves approach a shoreline and

break, momentum exchanges can produce currents and mean water level variations. These

currents can be in the form of longshore currents, rip currents, or circulation cells in the lee

of obstructions. They in turn modulate the incident wave field. The process can be quite

complicated depending upon the complexity of the bathymetry, the ambient current field,

and the presence of obstructions.

Coastal engineers and planners have many reasons for wanting to obtain accurate pre-

dictions of waves and currents and water elevations in the coastal zone. For instance, pre-

dictions of extreme conditions are needed for site planning and erosion mitigation studies.

Harbor authorities need wave predictions for structure design and for operational decisions.

This dissertation develops a numerical model for obtaining predictions of waves, currents

and mean surface elevations in a given domain with the offshore wav -' '

bathymetry as the input. The model uses a linear wave equation and depth integrated

equations of momentum and mass conservation. Comparison with experimental results

2

indicates that the model produces useful and valid predictions of the waves and currents in

the nearshore region in the presence of obstructions such as groins and breakwaters. The

predictions obtained through use of the computer model herein reported are not intended

to be the sole resource for the engineer. A good engineer or planner will seek answers from

multiple sources and cross check them, recognizing the limitations of each method used to

obtain the answer.

Predictions of waves and currents can be obtained for different cases analytically,

through experiments or by numerical solutions. Analytical solutions are obtained only

for simple cases involving idealized beaches using simplified assumptions which are not al-

ways satisfied in nature. Longuet-Higgins and Stewart (1963) give an analytic solution for

wave set-up for steady state conditions on a plane beach with straight and parallel con-

tours, assuming linear waves approaching the beach normally, and spilling breakers where

the wave height shoreward of the breaker line is a constant times the total depth of the

water column. Longuet-Higgins (1970a) gives a cross shore profile of the longshore current

for a plane beach with straight and parallel contours, for linear waves approaching the beach

at an angle, spilling breakers, and linearized bottom friction, while neglecting the influence

wave set-up has upon the total depth of the water column. In a companion paper, Longuet-

Higgins (1970b) extended the solution to include lateral diffusion of the current momentum.

There are other analytical solutions for waves and currents but all are limited to very simple

cases.

The literature is replete with laboratory investigations of waves and currents. Physical

models have advantages and disadvantages. The advantages of experimentation are that

results can be obtained for complex situations where there is uncertainty about the exact

form of the governing equations or about the precise form of any of the component terms

of the equation. Any reasonable amount of detail can be included. Visualization of the

physical processes provides an opportunity to identify problems which otherwise miu,-t no.

be evident.

3

However there are limitations to the use of laboratory experiments. There is the problem

of scale effects which means that all scaling requirements cannot be satisfied simultaneously;

thus, all features may not be interpreted properly to the prototype. For example, the

time scale dictated by Froude modelling laws in which gravitational forces dominate is

different from that obtained from Reynolds modelling laws in which viscous forces are most

important. Another disadvantage of physical models is the cost of labor to construct and

run large scale experiments. Also, because of the cost of construction, it is not always very

easy to change details in the model. Thus, it is not always possible to isolate one parameter

for study. For example, an investigation of the influence of the size of the gap between

segmented breakwaters might require many reconstructions of the rubble mound structure.

The large size of a model means that it is not easily stored once the initial experiment

is completed. The options are to either dismantle the model upon completion of testing

or to incur some maintenance costs. The significance of this is that the model cannot

always be easily retested if, for example, additional input information becomes available or

man or nature makes a substantial alteration of the prototype after the model has been

dismantled. Also the results of a model for one prototype are not always transferable to

other prototype situations. Often the available space for the physical model will require

the use of a distorted model in which the vertical scale is different from the horizontal

scale. This presents the additional problem that for a distorted model, when modelling

a situation where both refraction and diffraction are equally important, the laws of scale

modelling dictate different time scales for refraction and diffraction. For refraction the

time ratio is Tr = -/v, where the subscript r is for the ratio of the value in the model

to the value in the prototype and T and Lv are time and vertical length respectively. For

diffraction T, = VA in shallow water and Tr= \/Vlj in deep water, where Ly, represents
LV?
the horizontal length scale.

Numerical models do not have many of the disadvantages found with physical models.

Any mechanism for which the governing equations are known can be modeled numerically.

4
There are no limits in terms of scale effects. Numerical models, once developed, can be

adapted easily to many different prototype situations. There is no problem in the storage or

duplication of the model. Today computer facilities are far more available than laboratory

facilities. Computer models easily allow for sensitivity testing of the importance of an

individual variable and any specific feature of the model can be monitored.

!. Numerical models also have their disadvantages. In fact some of their advantages

are also their disadvantages. For example, being able to monitor any feature also means

that the computer will deliver only the information that is explicitly specified. Thus an

investigation of one aspect will not indicate problems with other aspects. Though any

mechanism for which the governing equation is known can be modeled, the model is only as

good as the assumptions involved in the equations. For example, the forces exerted on the

water column by the complex interaction between the fluid motion and the bottom are not

totally understood and are often greatly simplified. Another example is the use of depth

integrated two-dimensional models of fluid flows when the reality is three-dimensional. Both

of the above are incorporated in the present study. Knowing the correct equation does not

always make for successful computer modelling since the complexity of the equations may

make for difficulties (for example solving for the complete three-dimensional flow field using

the Navier-Stokes equations). Computer time and cost for large programs can often be

significant.

1.1 Literature

This section is a brief history and literature review of the attempts of the coastal

engineering profession to obtain predictions of waves and currents in the nearshore region

given the offshore wave conditions.

The extension of Snell's law which governs optical wave refraction to the analogous

problem of water waves over straight and parallel contours has been known for some time.

Johnson, O'Brien, and Isaacs (1948) reported graphical methods for constructing wave

refraction diagrams. Two methods were reported, one producing wave fronts and the second

5
producing wave rays. Pierson (1951) showed that the method using wave fronts introduced

more error than the wave ray method. These methods proceeded upon the assumption that

the contours are locally straight. The wave rays obtained were important in that the wave

height was determined to vary inversely as the square root of the spacing between the wave

rays. Using these ray construction techniques O'Brien (1950) was able to demonstrate that

thp April 1930 destruction of the Long Beach Harbor breakwater on a seemingly calm day

most probably resulted from the focusing of long wave energy from a South Pacific storm.

Munk and Arthur (1952) showed a method to obtain the wave height directly from the

curvature of the ray and the rate of change of the depth contours. The wave ray tracing

techniques did not include the effects of currents and had the disadvantage (from a numerical

modeling point of view) that wave heights and directions were obtained at positions along

a ray and not at pre-specified grid positions. Noda et al. (1974) solved this problem by

devising a numerical scheme that solved a wave energy equation and a statement for the

irrotationality of the wave number to obtain wave height and direction at grid points. This

scheme also included the effects of current refraction but neglected diffraction.

Numerical models of wave induced currents were developed by Noda et al. (1974) in

which time dependent and advective acceleration terms were ignored. Birkemeier and Dal-

rymple (1976) wrote a time-marching explicit model which neglected lateral mixing and

advective acceleration. Ebersole and Dalrymple (1979) added lateral mixing and the non-

linear advective acceleration terms. Vemulakonda (1984) wrote an implicit model which

included lateral mixing and advective acceleration terms. All of the above used the numer-

ical scheme of Noda et al. (1974) to obtain wave heights and wave angles. Each of these

models neglected diffraction.

Diffraction of water waves represents some very difficult mathematical problems. His-

torically most solutions were based upon the assumption of a flat bottom, since to quote

.ewian (1972),
The restriction to constant depth h is important.... Problems involving wave
propagation in a fluid of variable depth are of great interest, and the resulting

6
diffraction effects are among the most important of all ... but these problems
are less tractable than those involving a fluid of constant depth. (page 2)

Penney and Price (1944) showed that Sommerfeld'L solution of the optical diffraction prob-

lem as described by Bateman (1944) is also a solution of the constant depth water wave

diffraction problem. Penney and Price (1952) solved the constant depth diffraction problem

for breakwaters assuming a semi-infinitely long, infinitesimally thin barrier. The solution

is m terms of Fresnel integrals which need to be evaluated numerically. The Shore Protec-

tion Manual (SPM,1977) produced many graphical solutions based upon Penney and Price

(1952) that give the wave height and direction in the lee of breakwaters.

The problem of combined refraction and diffraction according to the SPM could be

handled by assuming a constant depth for a few wavelengths so as to obtain the diffraction

effects and then continue with the refraction effects. Perlin and Dean (1983) used essentially

the same procedure in a numerical model which calculated the longshore thrust in terms of

the breaking wave height and angle.

Liu and Mei (1975) solved for the wave field around shore parallel and shore perpendicu-

lar breakwaters by using a curvilinear coordinate system where one coordinate corresponded

to the wave ray and the other to a curve of constant phase. Their solution was for "an am-

plitude modulation factor D((, C) which accounts for diffraction" (page 72), where C and C

are the two coordinates. An equation for D(f, C) was obtained and a boundary value prob-

lem was established by matching conditions in the incident zone and the shadow zones. For

constant depth the solution for D(C, C) is in terms of Fresnel integrals. Upon obtaining the

wave field the radiations stresses were obtained and used to force a circulation model that

solved for the stream function and mean water surface level while ignoring the non-linear

advective acceleration terms and the lateral diffusion of current momentum. Although Liu

and Mei produced results for the wave field, currents, and set-up resulting from the diffrac-

tion of w. ve Yv hrc-kwaters, their study is not a wave-current interaction model since

there are no current terms in the wave equations and no feedback of current terms from the

circulation model.

7

Berkhoff (1972) derived an equation governing combined refraction and diffraction with

the assumption of a mild slope. Radder (1979) applied the splitting method of Tappert

(1977) to obtain a parabolic equation governing the propagation of the forward propagat-

ing component of the wave field. Booij (1981) extended the work of Berkhoff to include

the effects of currents and also obtained a parabolic approximation. Kirby (1984) made

corrections to the equation of Booij.

Parabolic wave equations have recently been used in wave-induced circulation models.

Liu and Mei (1975) solved for the wave field by matching conditions between the shadow

zones and the illuminated zone in the lee of a structure. The wave information was then

used to force a circulation model which formulated the governing depth integrated time

averaged momentum equations in terms of a stream function. However this model did not

have current feedback in the wave equation. Yan (1987) and the present study both include

a current feedback for the parabolic wave equation and solve directly for the mean currents.

Prior to the use of parabolic wave equations to obtain the wave forcing in circulation

models, refraction wave models using a numerical scheme derived by Noda et al. (1974)

were used. Birkemeier and Dalrymple (1976) produced an explicit numerical solution of

the depth averaged equations of momentum and continuity using the wave model of Noda

et al. (1974). Their model did not include the nonlinear advective acceleration terms

or any mixing. Ebersole and Dalrymple (1980) included the advective acceleration terms

and lateral mixing in an explicit model which also used the refraction scheme of Noda

et al. Vemulakonda (1984) Produced an implicit model which also was based upon the

wave refraction scheme of Noda et al. Yan (1987) used a parabolic approximation to the

mild slope equation which accounts for combined refraction and diffraction. The present

model Improves upon Yan's model in numerical efficiency, by introducing structures and

most importantly by obtaining the radiation stress terms directly from the complex wave

.111tuufr

J=N
J=N-1
J=N-2

LL)
cccc

Lo

Y
LJ=3
J=2
J=l

8

LPTERAL BOUNDARY

LATERAL -
BOUNDARY : X
II II II

Figure 1.1: Computational domain of the model.

1.2 General Description of the Model

This section is a description of the numerical model that has been constructed to calcu-

late the wave and current field within a given domain. This domain usually is bordered by

a shoreline along one boundary and an offshore region along the opposite boundary. The

other boundaries are sufficiently far from the region of interest (i.e. a structure or shoal,

etc.) so as to not affect the region of interest. This will be discussed further in the section

on boundary conditions. The computational domain for the purposes of this investigation

consists of a rectilinear coordinate system. Figure 1.1 shows the computational domain.

The model consists of two parts, a circulation model and a wave model. The circulation

model computes the mean water surface level and the depth averaged mean currents using

th- radit.,.n serass* s from the wave the f, -ng terms and includes bottom friction,

advective acceleration terms and the ILeral diffusion of current momentum. The wave

model determines the radiation stress terms by first solving for the complex amplitude as a

9
function of the total depth and the depth averaged mean currents, and then computing the

radiation stresses from the complex amplitudes. The model starts with initial conditions

of zero mean currents and zero mean water surface elevation. The offshore incident wave

height is initially zero and then is smoothly increased up to its given value using a hyperbolic

tangent curve. This is done so as to minimize any transients due to a shock start-up.

The model iterates between the wave model and the circulation model until steady state

convergence is obtained. Each of the two parts of the model will be described separately.

The circulation model uses an alternating direction implicit (ADI) method of solution.

Three equations are solved for three unknowns. The unknowns are fj the mean water surface

level, and U and V, the two components of the depth averaged mean current. The equations

are the continuity equation and the x- and y-direction momentum equations. Each iteration

of the circulation model solves for the three unknowns on the next time level. In other words

each iteration advances the values of q, U and V by a time increment, At. For the purposes

of this investigation the model solves for steady state conditions using the time marching

character of the circulation model. The model can be used to solve for time dependent

wave fields, but this would require the substitution of a time dependent wave equation in

the wave model. The ADI method first solves, grid row by grid row, for U at the next time

level and an intermediate value of q, by solving the x direction equation of motion and the

continuity equation using an implicit numerical scheme. This solution is in terms of the

known values of U,V,and 4 at the present time level. Once all the grid rows are solved in

the x direction, the y equation of motion and the continuity equation are solved, again by

an implicit scheme, grid column by grid column, for V and f at the next time level in terms

of U at the new time level, q at the intermediate level and V at the present level. For each

of the implicit schemes in the two directions, the governing equations yield a tridiagonal

matrix. The tridiagonrl matrix and the method of solution will be discussed later.

The wave i.udel "iso yields a tridiagonal mati---: quatio,. The governing equation is

a parabolic equation for the complex amplitude A. The complex amplitude contains phase

10

information meaning among other things that y variations of the wave field are contained

within the complex amplitude. The term "parabolic equation" indicates that there are

first order derivatives in the x direction and second order derivatives in the y direction.

The method of solution is to solve for the amplitude on the I+1 grid row in terms of the

known amplitude values on the I grid row. The solution scheme starts on the first grid row

using the prescribed offshore incident wave and then advances grid row by grid row to the

last grid row. To advance to the next grid row the amplitude values on the I+1 grid row

are implicitly formulated using a Crank-Nicolson scheme. This will be further discussed in

section 4.2 which describes the finite-differencing scheme for the parabolic wave equation.

Once the complex amplitudes are known the radiation stress terms are computed in terms

of the absolute value of the amplitude and the gradients of the amplitude. This is described

in section 2.2.

The circulation model and the wave model are the heart of the program; however there

is much more. In all there are some 27 subroutines in the program. These subroutines cover

everything from input and output to convergence checks and flooding of the shoreline. The

main program is compartmentalized into many subroutines so as to facilitate any-changes

or future upgrading. For example, the bottom shear stress, a term of importance in the

equations of motion, currently uses a crude approximation. If a better approximation is

developed that is also computationally efficient, it will be an easy task to replace the present

subroutine that computes the bottom shear with the newer subroutine. There are many

places where the model can be improved with better assumptions and better approxima-

tions. These will be discussed in detail in the final chapter outlining recommendations for

further study. The modular construction of the model also allows for ease in substituting

one wave model for another.

Figure 1.2 shows a flow chart of the computer program. The subroutine INPUT is called

to initialize the values of thi anknuwns and to input the wate .*n, ahw offshore wave

height and direction and the position of any structures. The three subroutines WAVE,

Figure 1.2: Flow chart of computer program.

12

SHEAR, and LATMIX compute the radiation stresses, the bottom shear stress, and the

lateral mixing coefficients. These terms are the forcing terms in the equations of motion

and are used in the subroutine CIRC, which is the circulation model, to solve for q, U and V.

There is a decision tree which directs the program to the three subroutines WAVE, SHEAR,

and LATMIX at selected intervals. This is purely arbitrary and is done in the interest of

computer speed. It is assumed that the forcing terms do not change at a great rate and can

thus be updated periodically. This updating is more frequent during the start-up time while

the offshore incident wave height is being increased. The CHECK subroutine is to check the

model for convergence. Since steady state conditions are assumed, this subroutine indicates

convergence when the absolute value of the percentage difference between the updated values

of U, V, and q and the previous values is less than an arbitrarily small specified value. The

UPDATE subroutine simply updates the values of the unknowns for the next iteration. The

FLOOD subroutine allows for the addition or elimination of grid rows at the shoreline to

allow for beach flooding due to set-up.

Figure 1.3 is a flow chart of the circulation model. The subroutines XCOEF and

YCOEF determine the coefficients of the tridiagonal matrix equation for the general cases

in the x and y directions. The subroutine TRIDA solves the tridiagonal matrix equation

using the double-elimination scheme described by Carnahan, Luther, and Wilkes (1969).

Subroutines XCOJ1 and XCOJN are adaptations of XCOEF to determine the coefficients

of the tridiagonal matrix for the special case along the lateral boundary where J=1 or

J=N. Likewise YCOIM is an adaptation of YCOEF for the I=M grid row which borders

the shoreline. These adaptations are necessary because boundary conditions necessitate

different numerical formulation of some of the derivative terms. Similarly with the presence

of an emerged impermeable barrier bordering the shore side of the JJJ grid row, special

subroutines YCJJJ and YCJJJ1 are used on the grid rows I= JJ and I=JJJ+1 so as to

properly include the boundary conditions for these grid ,ows, and for shore perpendicular

groins subroutines XCJGR and XCJGR1 incorporate the boundary conditions for a groin

Figure 1.3: Flow chart of circulation model.

Do Loop I=M

I=
1,M-1

CALL AMPCO
CALL CTRIDA
Amp(I+l,J) = Vector(J)

CALL RAD L

(Return

Figure 1.4: Flow chart of wave model.

located between the J = JGR and the J = JGR+1 grid columns.

Figure 1.4 is a flow chart of the wave model. Subroutine WAVNUM determines the wave

number k, the group velocity C, and the intrinsic frequency a at each of the grid centers.

These terms are used in the governing parabolic equation. Subroutine AMPCO determines

the coefficients of the tridiagonal matrix and CTRIDA solves the matrix equation for the

complex amplitude on the next grid row. CTRIDA is essentially the same as TRIDA except

that it has been modified to handle complex numbers. The subroutine WW is called by the

AMPCO subroutine in order to determine the dissipation coefficient following the work of

Dally (1980,1987).

Chapter 2 describes the circulation model in greater detail. Individual sections w;

cover the governing equations and the formulation of each of the component terms. Chap-

ter 3 describes the wave model, the derivation of the linear wave equation, the parabolic

15
approximation, the manner of including wave energy dissipation due to wave breaking and

the method of obtaining the radiation stress terms directly from the complex amplitude.

The fourth chapter gives details of the numerical methods used to obtain solution of the

governing equations and the boundary conditions employed. The fifth chapter presents the

results of applications of the model and comparisons with other studies both numerical and

experimental, as well as some analytical solutions to some idealized cases. The final chapter

summarizes the study and makes recommendations for future research.

CHAPTER 2
CIRCULATION MODEL

2.1 Governing Equations

The governing equations for the circulation model are the depth integrated time aver-

aged horizontal equations of momentum

aU vU au aq 1 1
-~+ + + g + ; rTz
Bt Qz ay ox pD pD
1 as,, as,, 1a8
+- aI= + + = (2.1)

av 9av + av aqi 1 1
a ax g-y y+ +pDTbv- *p
1 as ,, as, 1lat
+ + + + = =O (2.2)
pD azx ay p 8a

and the continuity equation

a. a a
-a+ (UD)+ -(VD) = 0 (2.3)
Tt az ay

where

U = x component of the mean current

V = y component of the mean current

q = mean water surface elevation

p = water density

D = the total water depth = ho + *j

ho = local still-water depth

r, = the lateral shear stress due to turbulent mixing

ri. = x component of the bottom shear stress

17
-- = y component of the bottom shear stress

r,, = x component of the surface shear stress

r,y = y component of the surface shear stress

and S,=, S., and Sy are the radiation stress components which arise from the excess

momentum flux due to waves. It should be noted that U and V include the wave induced

miss transport.

These equations are obtained by integrating the local x and y momentum equations and

the continuity equation over the depth of the water column and then time averaging the

results. This is demonstrated in Appendix A, which follows closely the work of Ebersole and

Dalrymple (1979). For greater details of the derivation the reader is referred to Ebersole

and Dalrymple (1979) or Phillips (1969).

These are the same equations that have been used by other modelers such as Birke-

meier and Dalrymple (1976), Ebersole and Dalrymple (1979), Vemulakonda (1984), and

Yan (1987). Birkemeier and Dalrymple used an explicit numerical scheme and neglected

the advective acceleration terms and the lateral mixing. Ebersole and Dalrymple added the

advective acceleration and lateral mixing terms but still used an explicit numericalscheme.

Vemulakonda introduced an implicit numerical formulation for these equations. Each of

these investigators used a refraction wave model based upon the work of Noda et al. (1974)

which solves for the wave height and wave angle at given grid points. Yan (1987) used a

combined refraction-diffraction wave equation to determine the complex wave amplitude

(magnitude and phase) and then solved for the wave angle in terms of the slope of the wa-

ter surface. Each of these studies obtained the radiation stress components in terms of the

wave height and wave angle. The present study differs from the previous ones in that the

radiation stress terms are obtained directly from the complex amplitude and its gradients.

This w;U be discussed in section 2.2.

The use of Eqs. (2.1-2.3) involves many assumptions. The primary assumption is

that the flow field can be represented in two dimensions using depth and time averaged

18
equations. This assumption implies that vertical fluctuations of the flow which are zero in

the time average can be ignored. Many details of the flow may be missed by this assumption

of purely horizontal flow. Of course reducing the problem to a two-dimensional flow has

the distinct advantage of reducing an essentially intractable problem to a simpler tractable

form.

; Some of the terms in Eqs. (2.1-2.3) also involve certain assumptions. These terms

are the bottom friction, the surface wind stress and the lateral mixing. These are all

phenomenon that are not completely understood so that assumptions are made concerning

the complex forces that are represented. These assumptions usually involve the use of

empirical coefficients. In this study the wind forces are ignored so the assumptions involved

with the surface stress will not be discussed in detail here. It would be a relatively simple

task, though, to add a wind surface stress to the model, using a quadratic shear stress as

formulated by Birkemeier and Dalrymple (1976) and also used by Ebersole and Dalrymple

(1979) following the work of Van Dorn (1953) in determining the value of the empirical

coefficient.

The interaction of an oscillatory fluid flow with a real bottom is very complex. Many

physical phenomena may be at work. With sand bottoms there can be percolation into the

sand and frictional effects resulting in wave energy dissipation. With mud bottoms where

the mud is essentially a dense fluid layer there can be damping and internal waves along

the mud-water interface. Wave energy losses due to bottom effects are not included in the

model but can easily be incorporated into the wave dissipation subroutine. Resuspension

and settling of particles also constitutes a complex phenomenon with density changes and

momentum transfers. All of these effects are glossed over in the model by assuming a rigid

impermeable bottom. As with all of the models discussed above the bottom friction term is

represented in the quadratic form. Noda et al. (1974), Birkemeier and Dalrymple (1976) and

Vemulakc.'a (1984) all used the weak current formulation as proposed by Longuet-Higgins

(1970a) in which

T- = pF juIb U (2.4)

S= pF u",bl V (2.5)

where F is a drag coefficient (of the order of .01) and juorbl is the time average over a wave

period of the absolute value of the wave orbital velocity at the bottom. From linear wave

theory it is found that
2H
|Uorbl (2.6)
1 T sinh kh(2

where H is the wave height and T is the wave period, k is the wavenumber, and h = D is

the total water depth. Ebersole and Dalrymple used a more complete quadratic formulation

rb = pFtlI(j (2.7)

where the overbar represents the time average over a wave period and tr is comprised of both

the mean current and the wave orbital velocities. The total velocity vector ti< is expressed

as

ut= (U + cos O)7+ (V + fsin 0)j (2.8)

so that its magnitude is given by

|(l = %VU2 +V2 + U2 + 2Ui cosO + 2Vsin 0 (2.9)

The wave orbital velocity C is expressed as

S= a cos at (2.10)

where . is the maximum wave orbital velocity at the bottom which is found to be
jrH
T sinh kh (2.11)

The x and y components of the time averaged bottom friction are thus

of /(U + Q, cos coat) U (t d(at) (2.12)

T p = f (V + sin 0co. ot) lu d(at) (2.13)
2T Jo

where Itru is written as

Itj = /U2 + V2 + cos2 t + 2Um cos at cos + 2Vm cos at sin 0 (2.14)

The integrals in Eqs. (2.12) and (2.13) are then evaluated by a 16 point summation using

Simpson's Rule.

In the present study several different formulations are used and compared. The first is

the formulation given by Longuet- Higgins (1970a) and used by several previous investiga-

tors. The second is a modification of the Longuet-Higgins formulation with partial inclusion

of current strength in the friction formula.

rb = pF jub U + pFIUIU (2.15)

r, = pF juoM V + pFIUIV (2.16)

It is thought that this formulation has some value since in the lee of structures where there

is relatively weak wave activity, there may still be relatively strong currents generated by

the gradients of the set-up outside of the sheltered region. The third formulation is similar

to that of Ebersole and Dalrymple in that the full quadratic expression is used. However the

wave orbital velocities at the bottom are expressed as gradients of the complex amplitude

and the resulting integral is determined using an eight point Gaussian integration. It is

found that an eight point Gaussian integration of a cosine squared integration over an

interval of 27r is accurate to within 1 percent.

The quadratic formulation of the bottom friction is used

rb = pFuiit (2.17)

where tr is composed of both the mean current and the wave orbital velocities. The wave

orbital velocities are expressed as the real part of the gradient of the complex velocity

; ,enti.,i P which is expressed as

S= A(,y) cosh k(ho + z) (2.18)
o cosh kh

21
where A is the complex amplitude, w is the absolute frequency, o is the intrinsic or relative

frequency following the current, and ho, h and k are as previously explained. Equations

governing and a parabolic approximation solution for A are presented in the next chapter.

The total velocity vector ti is expressed as

[^ ( SI h (2.19)
C = [U + R ( ) + [V + ( )] (2.19)

where R indicates the real part, so that its magnitude is given by

(2 +V2+)2U +2VR ( + + 2 +O(2.20)

The x and y components of the time averaged bottom friction are thus

r = PF f U +R -)l .'tdt (2.21)

rby = PF [V +R( R )]-jdt (2.22)

As in Ebersole and Dalrymple, the assumption is made that only the gradients of A are

retained and the gradients of the depth and spatially varying quantities such as k and a

are ignored. This is essentially the assumption of a flat bottom and very slowly varying

currents.

Figure 2.1 plots the magnitude of the longshore current as a function of the distance

offshore for the various friction formulations discussed above. These velocity profiles were

obtained numerically for an eleven second wave with a deep water wave height of 1.028

meters (1 meter at offshore grid of the computational domain, using a AX of 10 meters

and 30 grid rows) and a deep water wave angle of 17.28 degrees (10 degrees at the model

boundary) approaching a 1:25 plane beach. A friction factor, F, of .01 was used for each case.

Using a friction formulation that ignores the contribution of the waves results in a very weak

friction and thus a large longshore current. Using the Longuet- Higgins formulation gives a

friction slightly ... ...-. r-,, on, while a combination of

the Longuet-Higgins formulation and the strong current model (no wave contribution) yields

a slightly stronger friction. For practical numerical considerations the latter formulation is

tN
0 ..... =pFIUI

3. '.......... r = pFlui| C

o

S- TUor\+
s" .." u.U
C .y.

J / T .."j....'....-U
00
I-

6C D BREAK
,/ "POINT ..-.....

o 50. 100. 150. 200.
DISTANCE FROM SHORE (M)

Figure 2.1: Longshore current profiles for different friction formulations for an 11 second,
1.028 meter wave at 17.28 degrees to a 1:25 plane beach.

most often used since it requires less computation time and more importantly the Gaussian

integration involving the real part of the gradient of complex amplitudes at times produces

stability problems.

The lateral mixing term in Eqs. (2.1) and (2.2) also involve several simplifying assump-

tions. The first assumption is found in the derivation of the two depth integrated equations

of motion where it is assumed that the turbulent shear stress term -pu'v' is independent

of depth. The other assumption is that the process of momentum transfer due to turbulent

fluctuations can be represented by a product of a mixing length parameter and derivatives

of the mean current. This is what is done in the present study and also in previous studies.

Ebersole and Dalrymple (1979), Vemulakonda (1984), and Yan (1987) all used the formu-

lation given by Longuet-Higgins (1970b) in which he developed an analytical expression for

the cross shore distribution .. .--- ---- .. .. rized friction

formulation as described above, linear waves at an angle to the normal of a plane beach,

23
and a spilling breaker assumption. Longuet-Higgins wrote

( p av) (2.23)
= % Bz /
in which e. and ey are mixing length coefficients.

Longuet-Higgins reasoned that the mixing length parameter c, must tend to zero as

the shoreline was approached. He assumed that e, is proportional to the distance from the

shoreline, X, multiplied by a velocity scale which he chose as the shallow water celerity

V/7. Thus e, is written as

e, = NX Vg (2.24)

where N is a dimensionless constant for which Longuet-Higgins chose the limits as

0 < N < 0.016 (2.25)

Obviously, for physical reasons there needs to be some limit on the scale of these eddies.

Ebersole and Dalrymple arbitrarily chose the value of es near the breaker line to be the

maximum value and for e to be constant offshore from the breaker line. The same limit

upon c, is used in this study. The coefficient e, is assumed to be a constant and it is

arbitrarily assigned the maximum value of e,.

In the present study the further assumption is made in computing rT that the cross

derivatives are negligible (a common assumption in fluid mechanics). Therefore, the lateral

mixing terms become
1 r( a ( ) (2.26)
p ay ay 'ay
in the x direction and

Sar, -= a (2.27)
p ax ax ax
in the y direction. This means that the lateral mixing term is represented by a purely lateral

diffusion of momentum. That the lateral mixing term is a lateral diffusion term has some

f--n1.cr.:- -;i iat the aitiusion 'e.-n 13 expressed exoDicitly r -i thL,-:. -esults a time

step limitation to an implicit model. This will be discussed in section 4.1 dealing with the

finite differencing of the equations governing the circulation model.

24
2.2 Radiation Stresses

In a series of papers in the early 1960s Longuet-Higgins and Stewart (1962, 1963, 1964)

introduced and developed the mathematics of the radiation stresses. These stresses which

can be conceptually thought of as the "flux of excess momentum" can be obtained through

integration of the momentum equations over the water column as is shown in Appendix A.

Using linear wave theory the radiation stress components are found to be

SC, = E [n(cos' + 1) (2.28)

S, = En(in' +1) (2.29)

Sy = Encosesine (2.30)

where E is the wave energy equal to spgH2, n is the ratio of the group velocity to the wave

celerity and 0 is the angle the wave rays make with the x axis.

If the wave height and wave angle are known at each grid point, the radiation stresses

and their gradients are readily determined for use in the governing equations. To determine

the wave height and angle at each grid point Birkemeier and Dalrymple (1976), Ebersole

and Dalrymple (1979) and Vemulakonda (1984) each used the following scheme which is

based upon the irrotationality of the wave number and an equation for the conservation of

wave energy. The irrotationality of the wave number determines changes in the wave angle

due to refraction while the energy equation determines changes in the wave height due to

shoaling and refraction. This scheme was first introduced by Noda et al. (1974).

The irrotationality of the wave number

V x k=O (2.31)

where

S= Ikl cos 07+ kl sin OY (2.32)

is expanded in Cartesian coordinates to obtain

80 a8 1 8 kl 1 alk =
cos +sin T + sin T cos1 =0 (2.33)
zx By rk\ oz \k\ 9y

25

For waves on a current the dispersion relation is

w = a + VU = gk tanh jklh + Ik cos + Ikl in P (2.34)

Equation (2.34) is differentiated to eliminate and "I from Eq. (2.33). Using a

forward difference derivative in x and a backward difference derivative in y, Oij is obtained

as a function of cos 0, sin 9, U,V, k,h and their derivatives as well as 0 at the adjoining grids.

The sin 0 and cos 9 terms are expressed in terms of the surrounding four grids using a 2nd

order Taylor expansion. The solution is an iterative process alternating solution of the 0

equation with solution of the dispersion relationship (2.34) for k = [Ik which is accomplished

by a Newton-Raphson method until acceptable convergence is obtained for k and 0.

Next an energy equation is solved for the wave height. A conservation of energy equation

is given by Phillips (1969 ed. Eq. 3.6.21) as

aE a aU-
-T+ a {E(Ui + C,,)} + S,-i= 0 (2.35)
at axi axi

Assuming steady state conditions and expanding in Cartesian coordinates gives

1 8E 1 8E 8
(U + C, cos 8) + (V + C, sin O) + (U + C cos 9)
E E Uay ax V V
a au au av av
+-(V + C, sin )+ + +,, z + f + zy e 0 (2.36)

where

1 1
a = E-SZ = n(cos' +1) -
1 1
= EI
uyy = SVV = n(sin' 0+1)-2
1 1
= ncossin

Defining

Qi~j= -( + Ccos ) + a (V+c, sin9)+ (2.3"'

where

= 4-4- + #y + r tj (2.38)
at ayay

26
and noting that U,V,C,, and a are known quantities; taking forward difference derivatives

in x and backward difference derivatives in y and solving for the wave height, Hij, at the

i,j grid yields
(V+Cr in )H,i,_t (U+Ccos )HRi+,j
HiJ = A U (2.39)

The computation is a row by row iteration proceeding from large i (offshore in the

coordinate system of the above cited authors) to small i. During the computation of each

Hij the breaking height is also calculated and if Hij exceeds this value, then the breaking

wave height is used. Several breaking height indices are possible. Miche (1944) gave as the

theoretical limiting wave steepness the condition

Hb
= .142 tanh(kahh) (2.40)

where the b subscript indicates breaking conditions. Le MWhautW and Koh (1967) indicated

from experimental data that the limiting steepness is

= .12 tanh(kshb) (2.41)

or

Hb = .12tanh(kbhb) (2.42)

This is the breaking height criterion used by Noda et al. (1974). Perhaps the simplest

breaking height criterion is the spilling breaker assumption

H1 = .78h1 (2.43)

which is used by Vemulakonda (1984) and Vemulakonda et al. (1985).

The scheme described above of determining the wave angle 0 from the irrotationality

of the wave number (2.33) and the dispersion relation (2.34) and the wave height from

the energy equation (2.35) has its limitations. Prinreoally ;s scheme does not take i- '

account the effects of diffraction. This limitation is surmounted by use of a parabolic wavy

equation. However this presents another problem in that the wave angle is not determined

27
in a parabolic equation solution. Only the complex amplitude, which contains the mag-

nitude and the phase of the wave, is obtained. Yan's (1987) solution to this problem is

to obtain 0 directly from the complex amplitude by assuming that the slope of the water

surface resulting from the waves indicates the direction of wave propagation. This value

for 0 and twice the absolute value of the complex amplitude for the wave height are then

used in Eqs. (2.28-2.30) to obtain the radiation stresses. The value thus obtained for 0 is

then also used in the dispersion relationship. For plane waves this works fine; however, in

diffraction zones where short-crested wave are encountered problems may arise, since the

instantaneous slope of the wave at a particular location can be other than in the direction

of wave propagation.

In this report the radiation stresses are obtained directly from the complex amplitude

and its gradients using equations developed by Mei (1972) and repeated in his book (1982).

The basis of the derivation is to use the definition of the velocity potential to substitute for

the wave induced velocities in the definition for the radiation stresses.

For an inviscid incompressible fluid the instantaneous equations of motion are

p L+ V') = -Vp pg9 (2.44)

V.f = 0 (2.45)

where represents the velocity components, p is the pressure and e, is the z or vertical unit

vector. Denoting the horizontal components of q as Ua, a = 1, 2 and the vertical component

of q'as w, the boundary conditions are

P= 0 (2.46)

S+ u = w (2.47)
,t aBy
on z = q, where r) is the instantaneous free surface, and the bottom boundary condition

ah

q'is assumed to be composed of a time averaged part q and a fluctuating part q'; i.e.

4= q(z, ,) + (x, z, t) (2.49)

28
Integrating a horizontal component of (2.44) vertically over the depth of the water column

from z = -h to z = q), using the boundary conditions (2.46-2.48) and then taking time

averages (the same procedure as in appendix A) the following definition of the radiation

stresses is obtained:

S0p = P j audz + 60 pdz 2(q + h)2 + ()] (2.50)
I Jfh [ -h 2 2
where the overbar is used to indicate the time averaged values and Sap is the Kronicker

delta function

1 for a (2.51)
= 0 for a

Mei (1972) was able to express Eq. (2.50) in terms of the general complex amplitude

function B by relating it to the velocity potential and the free surface height by

i-g cosh k(ho + z) (2.52)
= B-------- ei (2.52)
a cosh kh
r B = Be-iCw (2.53)

An expression for p is obtained by integrating the vertical component of Eq. (2.44)

aw aw 8w ap
-t + PU +pw- = T P9 (2.54)

Integrating from any depth z to the surface, using Leibnitz's rule for interchanging the

order of differentiation and integration, employing the two surface boundary conditions,

Eqs. (2.46) and (2.47) and the continuity Eq. (2.45), Eq. (2.54) yields:

p(z, Y, z,t) = pg( )+pj wdz+ pa up"wdz pw (2.55)

Taking time averages gives

y, ) = pg(I z) + p dz p(w)2 (2.56)

Substitution of Eqs. (2.52),(2.53) and (2.56) into (2.50) and making the assun .,on .-a.

horizontal derivatives of k, h and a can be neglected (this is essentially a mild slope as-

sumption), the following is obtained after straightforward but tedious calculations which

are shown in Appendix B:

S, L9 B IaBB( 1 1 2kh )
4 tZ= ~Zap k2 +sinh 2kh
Srlsinh2kh W+
S[IBji2nh2jkh l+ (vBi k2 B2) (2khl coth2khL 1)J} (2.57)

where B* denotes the complex conjugate of B. In obtaining Eq. (2.57) use was made of the

flat bottom Helmholtz equation

V'B + k'B = 0 (2.58)

The general complex amplitude function B is related to the complex amplitude A which

is solved for in the wave model by

B = Aei(f coI d) (2.59)

There are some problems with Eq. (2.57). One may question the assumption of a

flat bottom; however the inclusion of the slope and gradients of the wave number and the

intrinsic frequency leads to intractable mathematics. The use of the localized flat bottom

is also implicit in the derivation of Longuet-Higgins and Stewart to obtain Eqs. (2.28- 2.30)

since in carrying out the integration involved it is assumed that horizontal derivatives of

the velocity potential are limited to derivatives of the phase function. For plane waves

using a spilling breaker assumption, results using Eq. (2.57) differ little from those using

Eqs. (2.28-2.30); however in diffraction zones Eq. (2.57) is superior for reasons mentioned

previously.

The major problem with the use of Eq. (2.57) is also a major problem associated with

the use of Longuet-Higgins's formulation. This problem is that in both formulations the

initiation of wave breaking implies the dissipation of energy and thus the forcing of set-

up and longshore currents starting at the breaker line. The physical reality is somewhat

different in that in the initial stages of breaking, energy is not immediately dissipated but

rather is transformed from organized wave energy to turbulent energy and to a surface roller

(Svendsen 1984). There is a space lag between the initiation of breaking and the dissipation

30

-30

2a

sACH, z.0

40
S- I.S

MEAN WATER LEVEL, 9
S10
0

POINT

PLUNGEY
ENVELOPE OF E HET *
BEAEAK
POINT
| PLUNGE
POINT
ENVELOPE OF WAE HEIGHT .t| S

WAVE CREST, .' *

-----------===^ ^^.---- 0
8 /
.. ..* ,.* "too.

WAVE TROUGM
I I 1 I 1 -.4.
400 300 200 100 0
DISTANCE FROM STILL WATER LINE ON BEACH, x (cms)
Figure 2.2: Experimental results for wave set-up and set-down. Reprinted from Bowen, In-
man, and Simmons (1968), Journal of Geophysical Research, vol. 73(8) page 2573, copyright
by the American Geophysical Union.

of energy. This can be seen in figure 2.2 which gives the results for the experiments of Bowen

et al. (1968). Analytical and numerical solutions for the mean water surface elevation will

give a set-up starting at the break point while the solid line in the lower graph shows

that set-up measured in the experiment begins shoreward of the break point. This will be

discussed further in Chapter 5.

CHAPTER 3
WAVE MODEL

3.1 Governing Wave Equation

Several methods have been used to derive the governing equation for the propagation

of waves over varying topography and varying currents. Among the methods used are per-

turbation techniques, multiple scales, a Lagrangian and the use of Green's second identity.

With the assumption of a mild slope and weak current conditions each of these methods

should yield the same equation. For illustrative purposes the method of using Green's sec-

ond identity will be used here to derive a linear equation. A Lagrangian derivation of the

same equation is presented in Appendix C. The linear equation is used since second order

waves exhibit singularities as the water depth goes to zero.

A velocity potential 0 is assumed such that the water particle velocities are given by

V4, where V is the three-dimensional gradient operator

V = + -3 + (3.1)
TV By Tz

and s,3, and k are the unit vectors in the x,y, and z directions, respectively. The velocities

are a superposition of steady currents and wave induced motion. The velocity potential is

given by

0= + ef4 (3.2)

where o is defined such that its gradient yields U, V, and W which are the steady current

terms, and f is the wave component of the velocity potential. An undefined scaling e

so as to separate the wave and steady current part of the velocity potential. The

wave part of the velocity potential is composed of two parts: f, which contains the depth

dependency and 4, which is a function of the horizontal coordinates and time. The depth

dependency f is given by linear theory as

Scosh k(ho + )(3
cosh kh

where ho is the depth of the bottom referenced to the still water line, and h is the total

water depth and k is the wave number vector. It should be noted that f is also a function

of the horizontal coordinates since k and ho vary with horizontal position.

For incompressible, irrotational flow the governing equation for the velocity potential

is the Laplace equation

V20 = 0 ho << r1 (3.4)

where q is the instantaneous water surface elevation.

Using a horizontal gradient operator Vh defined by

Vh= -~ + T- (3.5)

the Laplace equation is written as

V2 + ,, = o ho < z < 7 (3.6)

The governing equation is complemented by the following boundary conditions. On the free

surface the kinematic free surface boundary condition is given by

l+ Vh* Va =, z = 17(3.7)
Ft

The dynamic free surface boundary condition with zero atmospheric pressure

O, + -(V,)o + g = 0 z = (3.8)
4'+i 4)2+z=O z=w, (3.8)

and the bottom boundary condition

VA Vhho = ,. z = -ho (3.9)

The derivation technique is to write Green's second identity for 4 and f as defined by

,-4. 4 )ui ) respectively.

f f4,dz Of,- dz = [f -h (3.10)
h.ft- h~o

33
Then Eq. (3.6) is used to substitute for 4,, in the first term of Eq. (3.10) and the boundary
conditions are used to substitute for 4, at the mean free surface and at the bottom in the
third term of Green's identity. Equation (3.2) is used and terms of order e only are retained
in order to obtain the governing equation for 4.

In order to obtain ~, at f it is assumed that

S= i + eq. (3.11)

The kinematic free surface boundary condition, Eq. (3.7), is then expanded in a Taylor
series about z = q7,

a[ + V h -, + e + V 1,0V V J = 0 (3.12)
at a, Qz at ],=17
solving for 4, f#j using Eqs. 3.2 and 3.11 while retaining only terms of order e yields

-, l= L + U. VP + V,. .h fW,] (3.13)

A total horizontal derivative
D a
= + U. Vh (3.14)
Dt at
is introduced so that the first two terms on the right hand side of Eq. (3.13) can be written
as the total horizontal derivative of i. The final manipulation of Eq. (3.13) is to recognize
that at the surface is the horizontal divergence of the horizontal current (= Vh U).
Thus

S\= + V h" Vhq + qVh" C (3.15)
Dt
Sis eliminated from Eq. (3.15) by expanding the dynamic free surface boundary condition,
Eq. (3.8), in a Taylor series about z = f.

[t + (V4)' + gz = t+ + + (vo)

+Ec [ + 4.t + -(V4) = 0 (3.16)

Using Eq. 3.2 to sua O ute for 4 the order c terms yield a solution for i

1=- + p v ] = (3.17)
g g Dt

Substituting into Eq. (3.15) gives

1Dt2 1 D(31
O In= g DO + V h~*V (3.18)

The first term of the Green's identity using Eq. (3.6) becomes

Sf..dz = fV dz = -J fVL. (Vh)dz (3.19)
J-he J-ho ,-ho

Now
Vh = + ef VA# + eVAf (3.20)

and
Vh* (VhA) = Vh. 1 + eVh (f Vh)) + eVh Vhf + e hVf (3.21)

Retaining only terms of order c and dropping the last term since it is second order in the
slope the integral is equal to

L f Vh (f V )dz f V Vhf dz (3.22)
-ho -ho

These two terms are handled using Leibnitz's rule of integration.

Vh f(fVhV)dz = Vhf (fVh)dz+ f Vh(fVh)dz
-ho -ho -ho
+Vhlf2 I. Vh, + Vhhof2 I-h. Vh* (3.23)

Using this relationship and recognizing that f I| is equal to one and the integral over the
depth of f2 is CCg/g, where C is the wave celerity, C, is the wave group velocity and g is
the gravitational constant, the first term is equal to

Vh (-C Vh) + Vhq Vh) + Vhhof2 I-h. Vh$ (3.24)

where the e has been dropped.
The second term of Eq. (3.10) is evaluated quite easily.

fzdz= #of.dz- f#fdz (3.25)
-ho -ho ho

35
The first integral above is dropped since there are only steady current terms. The remaining
term, dropping the v, is

=- f 4 f,,dz = -ki f2'dz = -k2 (3.26)

The third term in Eq. (3.10) is

f #, l"-.= 1, I~ -/1, I 1-. (3.27)

The first term on the right hand side is equal to 4#, j~ since f I,7= 1 and has already been
determined. The second term on the right hand side of Eq. (3.27) is easily determined using
the bottom boundary condition Eq. (3.9)

fh, -A.,= fVh V*ho 1-h, (3.28)

Using Eq. (3.20) this becomes

f# l-ho.= f 7 Vhho I-h +tf2 -h Vh. Vhho + ef Vhf Vho I-h, (3.29)

The first term is dropped since it not order c and the third term is clearly higher order in
the slope. Thus, again dropping the e

f, -h.= f2 1-h, V,' Vhho (3.30)

This term cancels the bottom boundary term resulting from the Leibnitz integration above.
The fourth term of Eq. (3.10) is also easily determined.
sinh k(ho + z) 2
4'f l-h= -k c'sin h k lz = -4k tanh kh LO= (3.31)

Evaluating # at the mean water surface and retaining only terms of order e yields

f I= (3.32)

Collecting terms from (3.24),(3.26),(3.18),(3.30) and(3.32) yields the following linear
governing equation for 4.

D (V ) -t Vh(CC,V) + (or k2CC,) = 0 (3.33)
This equation was first derived by Kirby (1984).

36
3.2 A Parabolic Approximation

Equation (3.33) is a hyperbolic wave equation. It can be altered to an elliptic equation

by assuming that any time dependency is solely in the wave phase. Before doing this an

energy dissipation term DD is introduced as follows. It is assumed that Eq. (3.33) holds

for the case of no energy dissipation; then for the case with dissipation the left hand side is

balanced by any energy that is dissipated. The governing equation is then written as

Dt- + (V.- ) D Vh(CC,Vh) + (0' k CC,) = DD (3.34)

In the derivation of the governing hyperbolic equation using a Lagrangian as presented in

Appendix C the DD term can be viewed as a virtual work term.

Following Kirby (1983) the DD term is assumed to be of the form

DD = -w- (3.35)
Dt

where w is an undefined coefficient indicating the strength of the dissipation. Eventually

the w term will be related to the energy dissipation due to wave breaking following the work

of Dally, Dean and Dalrymple (1984).

Using Eq. (3.35) in Eq. (3.34) and making the assumption that the only time dependency

is in the phase, i.e.

t = -iwo (3.36)

reduces the hyperbolic equation (3.34) to the elliptic form

-2iwU Vh + U Vh(y VI) + (Vh. )(U.7 V ) Vh (CCVhj)

+{U2 w k2CC, iw(VA U7)}} = iowq (3.37)

where only the phase contribution to the horizontal derivative of is retained in obtaining

the term on the right hand side of (3.37).

Before deriving a parabolic approximation to Eq. (3.37), the need for such an approxi-

mation should be discussed. There are two major computational drawbacks to numerically

37

solving an elliptic equation. The first is that solution is required simultaneously for each

grid. This means that for a large computational domain, say m by n grids, it is necessary

to invert an m x n by m x n matrix. If m and n are sufficiently large, problems may be

encountered in terms of the internal memory storage, and even if the available storage is

not a problem, the number of computations required to invert the matrix may make for

anzexceedingly slow solution. In some situations this may be acceptable; however for the

purposes of a wave induced circulation model in which repeated solution for the wave field

is necessary the program would take a very long time to run. The other disadvantage in-

volved in the numerical solution of an elliptic equation is that boundary conditions must be

specified at all of the boundaries. In many applications the only known boundary condition

is the offshore, or incident wave amplitude. For example, the domain of interest may be the

wave field in the lee of an obstruction where the only down wave condition is that the waves

freely pass through the down wave boundary. Or in the present study in which the wave

amplitude is known to go to zero at the shoreline a problem ensues in that the shoreline is

(for purposes of the circulation model) at the grid edge and not at the grid point.

A parabolic equation eliminates both of these problems. Since the solution proceeds grid

row by grid row where each solution uses the results of the previous grid row, no down wave

information is required. Since only one grid row is solved at a time, the solution requires

only that a tridiagonal matrix equation be solved to obtain values for the grid row. The only

required boundary conditions are the conditions on the first grid row (usually the offshore

boundary) and lateral boundary conditions which are usually a derivative condition that

allows for the wave to pass through the boundary without any reflection at the boundary.

There are several ways to transform the elliptic Eq. (3.37) to a parabolic equation. The

most direct is to make certain assumptions concerning the length scale of variations in the x

direction. At the heart of the parabolic approximation is the assumption that the direction

of wave propagation is essentially along the x axis. For waves propagating at an angle to

the x axis the proper form of is

S= -ig z e)(f kcod+f k.ndy) (3.38)

and the proper form of the dispersion relationship is

w= o + kcoe s U + ksin V (3.39)

where 0 is the angle of the wave propagation relative to the x axis. The essence of a parabolic

approximation is to assume that the angle that the wave rays make with the x axis is small

so that the ksin 0 term can be neglected in Eqs. (3.38) and (3.39) and cos 0 is assumed to

be unity. The manner of doing this varies with different authors. Yan (1987) retains the

complete dispersion relation (Eq. (3.39)) but deletes the y dependency in Eq. (3.38) in the

derivation of his governing wave equation. He solves for 8 by determining the slope of the

instantaneous water surface. For long crested waves this method is satisfactory although

requiring additional computations. However, this method may produce erroneous results in

diffraction zones where the presence of short crested waves will produce a horizontal com-

ponent of the slope normal vector that is different from the direction of wave propagation.

Kirby's (1983) method is to make the assumption that

g= _A(-z, ei(f kdz) (3.40)

so that the e' f sin*dy part of the phase function is absorbed into the amplitude function,

A.

Using Eq. (3.40) directly in Eq. (3.37) and further assuming that derivatives of As are

small compared to derivatives of the phase function (i.e. that << ikA, ) will yield a

parabolic equation. For illustration the case of linear diffraction in the absence of currents

and in a region of constant depth is examined. Substitution of Eq. (3.2) into the Laplace

-:'V: ,is; the Helmholtz equation

Vi + k= 0

(3.41)

39
This equation is also obtained directly from Eq. (3.37) for the stated assumptions of no

currents, no dissipation, and constant depth. Substitution of (3.40) into the Helmholtz

equation gives

A,, + 2ikA, + Ayv = 0 (3.42)

and employing the above mentioned assumption concerning derivatives of A, yields the

lihnar parabolic equation

2ikA, + Avy = 0 (3.43)

Another method of obtaining a parabolic equation is to use a splitting method in which

it is assumed that < is composed of forward moving and backward moving components that

are uncoupled.

S= + + 4- (3.44)

Uncoupled equations for 4+ and 4- are developed and since the equations are uncoupled

only the equation for + is retained. In the following the work of Kirby (1983) and Booij

(1981) is closely followed.

The second order elliptic equation

a ( ))+- = 0 (3.45)

can be split exactly into equations governing forward and backward disturbances Q+ an I-

which are given by

a5+
-- = +{ (3.46)
ox
'- = -i'- (3.47)

(3.48)

Booij chose the relation

(3.49)

40
where C is an as yet unknown operator on S. Substituting into Eq. (3.45) and neglecting
derivatives of C alone leads to the result

&Z+ + + 7'=0 (3.50)

The elliptic equation (3.37) is now put into the form of Eq. (3.50). Letting p = CC,
for notational convenience and expanding some of the terms of Eq. (3.37) in their spatial
coordinates gives

(Ps.)t + (p,)y + k2P + (w' 02 + iW(Vh U. )) + iuow

-(U2&), (V2'y), (UV&),v (UV ,), + 2iwO Vt, = 0 (3.51)

Booij chose to neglect terms involving the square of the mean current components assuming
that they are smaller than terms involving p and then arranged the above equation as

& + P(ip + 2iwU), + k2 I + M = 0 (3.52)

where
M4 = (w2 o2 + iw(Vh C) + iow); + 2iwV9y + (p5v,)v (3.53)

Kirby points out that this choice by Booij leads to his inability to derive the correct wave
action equation. Kirby instead does not drop any of the squares of the mean current
components until after the splitting is complete. By assuming that the waves are oriented
in the x direction the dispersion relation is written as

o = w kU (3.54)

leading to the result
W a = 2wkU (kU)2 (3.55)

Usim; Ibis relationship Kirby there arranged Eq. (3.51) as

,, + (p U)-(p U,),, + k2 ( + k2(p U2) = 0 (3.56)

where

M, = {2wkU + iw(V,. 0) + i ,w) (UV,), (UV) ,)v

+ (p- V'), + 2iwO V, (3.57)

Comparison with (3.50) gives

I 72 = k2 1+ k+(p U2) (3.58)

= (p- U) (3.59)
-y

or solving for and C

S= (+ k2(p U2) (3.60)

(= k(p-Uk U2) 1+ 2) (3.61)

Booij points out that the splitting matrix approach employed by Radder (1979) is equivalent
to the choice

yR = k 1p-U
2k2(p U2)

Using the results (3.60) and (3.61) in Eq. (3.46) along with the definition (3.49) gives the
parabolic equation

I {k-(P-U') (1+ k2(p- U2))} = ikk (p- U')i + k2(p- U2) ) (3.62)

By expanding the pseudo-operators in the form
( M ) + 1+ M $* (3.63)
+ k(p- U) 4k(p U2) (3.63)

( M U + + ( 3 M )+ (3.64)
k2(p Us) 4 k2(p U))
the result c. otL expansions may be s;n.-.-. I :' .. -;.nc generall equation

a_ { -_ U) (1+ P}= p-M + i (+ P2M ( (3.65)
z U (+ k2(p U2)) ikk U) + k2(p U2) (3.65)

where
Io Radder (3.66)
pil I B6oij (366)
and where P2 = PI + I for both approximations. Kirby chose the approximation with
PI = 0 since for his derivation which was for weakly nonlinear waves it should correspond
to the information contained in the nonlinear Schr6dinger equation. With Pi = 0 Eq. (3.65)
becomes

k(p U') + {k(p U')} = ik2(p U2)4 + 2M^ (3.67)

where Mr+ is defined by Eq. (3.57). An equation for the complex amplitude A is then
obtained by making the substitution

= -ig (A) eif d (3.68)

which after dropping squares of the components of the mean current results in

a c,+ V
(C, + U)A + A+ )VA.+ (V2 )A
2 0 0 2 o

i a, + A = 0 (3.69)

Kirby then further modified the equation by a shift to a reference phase function using the
relation
A = A'ei(I- kd) (3.70)

where k can be taken as an averaged wave number, resulting in

(C, + U)A, + i(k k)(C, + U)A' + ( A'

+ VA+ -+ A' = 0 (3.71)

The shift to a reference phase function is intended to incorporate the phase information
in the complex amplitude so th' t the instantaneous water surface can be easily recovered
in the form

7(r, Y) = R {A'eEd} (3.72)

43
where S designates the real part. This is an important consideration especially when ob-

taining the radiation stresses.

It is important that in obtaining Eq. (3.69) that squares of the products of the current

components be retained until the last step. The reason for this is seen in that terms such
as kUV and kU2 are combined with terms of wV and wU using the dispersion relationship

to give terms of aU and oV.

Equation (3.71) has been successfully used by Kirby in the form given and with modifi-

cations. The modifications are the addition of a non-linear term which makes the equation

applicable for weakly non-linear waves. However for present purposes Eq. (3.71) presents a

certain problem which will shortly be corrected. The problem is that for waves propagating

at an angle to the x axis it is not possible to obtain a correct form of the conservation

of wave action. The case of no current on a beach with straight and parallel contours is

given as an illustration. The equation for the conservation of wave action reduces to the

conservation of wave energy flux. In mathematical terms this is

[C,coso0a] =0 (3.73)

where a is the wave amplitude. For waves approaching the beach in a direction normal to

the beach Eq. (3.71) reduces to

C,A' + j(C,),A'= 0 (3.74)

Multiplying Eq. (3.74) by the complex conjugate of A', and add to it the complex conjugate

equation of (3.74) multiplied by A', yields

[C,IA1] =0 (3.75)

which is Eq. (3.73) for the special case of cos = 1. However it is not possible to obtain

Eq. (3.73) from Eq. (3.71) for the general case of waves at angle 6 on a beach with straight

and parallel contours.

The importance of this inability to give the correct form of Eq. (3.73) is that if energy

flux outside the breaker line is not conserved then gradients of the S,y radiation stress

44
will generate currents off-shore of the breaker line in the circulation model. Although
such currents in general will be of a small magnitude, they can cause some numerical
computational problems. The problem though is easily corrected by using as the dispersion
relationship
w = + kcos U (3.76)

and the substitution
= _iA(,y)ei(f kot. ) (3.77)

Rather than using a splitting method, the procedure will be to go directly from Eq. (3.51)
using Eq. (3.77) and writing

w2 a2 = 2wkcos U (kcos U)2 (3.78)

in place of Eq. (3.55). Dropping the ( ), term and all terms containing products of the
current components that cannot be reduced as explained above results in the following
parabolic equation.

SCcos + U A+V v aV
(Co cos0 + U)A, + Cos + U A+VA+ ) A

-k,(1A- CC.(),] + A = 0 (3.79)

Now adding a phase shift
A = A'e (fl IcoTIdz-fj ck d) (3.80)

so that the free surface can be recovered by

r (z, y) = A'ei(f Ioidz) (3.81)

where i is the wave angle obtained by Snell's law in the absence of any obstructions or y
variation of the depth or current. The following governing equation is thus obtained.

(Cy cos O + U)A', + i( cos # kcoe )(C, cos + V}A' + f C'O ] A

+VA + () A'- kC,(1 cos O)A C,(),+ A' =0 (3.82)
i -01 i i a Y T '111

45
Details of the algebra and assumptions involved in obtaining Eq. (3.79) are given in Ap-

pendix D. For computational purposes it will be assumed that cos = cos 0.

From Eq. (3.82) the correct form of the conservation of wave action can be obtained.

For the case of no y variation of the wave number Eq. (3.82) reduces to

a C, co0 + U +V A V+
(C, co 0 + U)A', + [0CA+VAI+ A'

SkC,(1- cos2 0)A' C,(-)v + tA = 0 (3.83)

Multiplying Eq. (3.83) by the complex conjugate of A and adding to it A times the complex

conjugate equation of (3.83) yields

[(C, coso+U)iA .+ [V AI =0 (3.84)

which is a closer approximation to the exact expression which should have a (C, sin 0 + V)

term instead of V in the y derivative term.

3.3 Energy Dissipation Due to Wave Breaking

Kirby (1983) found a method of relating the dissipation coefficient w to the work of

Dally (1980) who proposed that the decay of energy flux in the surf zone was proportional

to the amount of excess energy flux. The excess energy flux is the difference between the

actual flux and a stable energy flux. In mathematical terms this relationship is expressed

as

(EC,), [EC, (EC,),] (3.85)

where

EC, = the energy flux

E = wave energy= -pgH9
ao
C, = the rate of energy propagation whi-b is -

K = constant

h = the water depth

46

(EC,), = the 'stable' energy flux.

H = the local wave height

The stable energy flux is obtained by assuming it to be a function of the height of a wave

propagating over a flat bottom after breaking ceases, and that this height is a constant
times the water depth, (h) = 7. The constants K and were determined by analysis of

experimental data obtained from Horikawa and Kuo (1966).

Kirby related Daily's work to the dissipation term w by deriving a conservation equation

for the wave action Garrett (1967) and Bretherton and Garrett (1968) showed that for

waves propagating in a moving and variable medium, the quantity that is conserved is the

wave action (), which is the wave energy divided by the intrinsic frequency. Kirby derived

a wave action conservation equation from the linear hyperbolic wave Eq. (3.34) using the

definition of W as given by Eq. (3.35):

S+ (Va. *) Vh(CC,Vh) + (2 k2CCO) = ioaw (3.86)

Kirby follows the approach of Jonsson (1981) assuming that 4 can be represented as

i = Re'e (3.87)

where R and E are real. From linear theory it is obvious that

R = PAI. (3.88)

To lowest order the derivatives of E are approximated as

Ot = -w (3.89)

Vhe = c (3.90)

Substituting into (3.86) gives the complex equation

D2R DR
2+ (V D) VA(CCVhR) i(oaR + 2aR1)-
Dt- Dt
i {V. ( U)R + 2(U) VhR + Vn *(kCC,)R+ 2kCC,. VhR}
-iw$4 = 0 (3.91)

47
Since R and e are real, the real and imaginary parts of (3.91) must individually equal zero.
Taking the imaginary part and multiplying by R gives

(oaR2) Of VhR' Vh (oOf)R2 CC,k VhR'
-VAh (kCC,)R2 owR = 0 (3.92)

which using the substitution kCC, = aC, is written as

(R2'), + V (1R'(0 + ,)) = -owR' (3.93)

Using (3.88) R2 may be replaced by

R2 = g' A 2g((3.94)

which yields the wave action conservation equation

(f) ( + ,) = ) (3.95)

Assuming a time steady wave field and neglecting currents, (3.95) becomes

(EC,), = -wE (3.96)

Comparing (3.96) with (3.85) and noting that (C,), = C,, w is written as

C, E KC,( (H.) (397)
h ( E h H2
where H is the local wave height given by 21AI.
In obtaining Eq. (3.97) Kirby neglected the effects of the currents. However, the same
expression (3.97) is obtained when currents are included in the derivation. Dally (1987)
extended his work to investigate wave breaking on currents and proposed the following
equation:

(C, + U) = [EC, (EC,).] (3.98)

For the case of a steady state wave field propagating in the x direction the conservation of
wave action Eq. (3.95) reduces to

+ U) W (3.99)

48
Setting the right hand sides of Eqs. (3.98) and (3.99) equal to each other results in Eq. (3.97).

At first glance it may appear unusual that the same expression for the dissipation

coefficient results for the case of currents as for the case of no currents. This can be

explained in part by the fact that Daily's work was with waves where breaking was depth

limited. For the case of current limited breaking (for example a stopping current), another

expression is necessary. In this study Eq. (3.97) will be used, since in general the current

field will be moderately weak. Such an assumption of weak currents was invoked in the

derivation of the governing parabolic equation.

CHAPTER 4
FINITE DIFFERENCE SCHEMES AND BOUNDARY CONDITIONS

1 4.1 Numerical Solution of the Governing Equations

The finite differencing of the governing equations in the circulation model is derived

by a matrix analysis as suggested by Sheng and Butler (1982). A rational analysis of the

complete Eqs. (2.1-2.3) including all the nonlinear terms is not possible. However guidance

can be obtained through an analysis of the linearized equations. For this purpose the linear

equations for the special case of a flat bottom and without any of the forcing terms will be

examined.
a + au av
a0 D9U DBV
-+ + D- = 0 (4.1)
at ax By
aU a9
+ gz -= 0 (4.2)
aV al9
-7 + g = 0 (4.3)
Tt ay
Let W be a vector of j,U and V

W = U (4.4)
V
Using Eq. (4.4), the governing equations are written as

Wt+AW,+BWV =0 (4.5)

where A and B are the following matrices
O D 0 0D
A= g 0 0 and B= 0 0 0 (4.6)
0 0 0 g 0 0
To find an implicit finite difference scheme for the governing equations, Eq. (4.5) is

finite difference as
wn+l W"
W+- + AW,+' + BW,+" = 0 (4.7)
AT

50
This equation is solved for W"+1 in terms of Wn

(1 + 2A, + 2A,)W"+1 = W" (4.8)

where
1 AT 1 AT
A = -A- 6 and A, = -B-- (4.9)
2 AX 2 AY
4AAYW"+' is added to the left hand side of the equation and 2AW" is added and sub-

tracted to the right hand side of Eq. (4.8) so that the equation can be factored to get

(1 + 2A,)(1 + 2AV)W"+1 = (1 2A,)W" + 2A,W" (4.10)

An intermediate value W* is introduced so that the above equation can be written as

two equations

(1 + 2A,)W = (1 2Ay)W" (4.11)

(1+ 2Ay)W"+l = W* + 2AyW" (4.12)

Note that by eliminating W* from the two Eqs. (4.11) and (4.12), Eq. (4.8) is recovered
with the addition of an error of 4AA,X,W"n+ which has been introduced.

Writing out the component equations of Eq. (4.11) gives

AT AT
1* + D- 6,U' = Q" D-6,V" (4.13)
AX AY
ST
U'+g- = U" (4.14)
AX
AT
V*= V"- g y" (4.15)

(4.16)

While writing out the component equations of Eq. (4.12) gives

AT AT
f"+l + D ,V"n = + D 6,V" (4.17)

Un+1 = U' (4.18)

Vn+l + "n+1= V +g ,T (4.19)

51
Using Eq. (4.15) to eliminate V* in Eq. (4.19) and (4.18) to eliminate U* in Eq. (4.14) the

following are obtained. For the x direction

AT AT
+ D 6, U"1 = q" D 6,V" (4.20)
AX AY
AT
U"+1 +g O 6S* = U" (4.21)

and for the y direction

AT AT
"q+l + D A-,Vn+ = q* + D n (4.22)

V" + gT 6. = V" (4.23)

The above equations are numerically stable for any combination of AT, AX, and AY.

However as a circulation model they are incomplete in that they do not contain any of the

forcing terms that drive the circulation. In particular in this model wave forcing represented

by gradients in the radiation stresses and bottom friction are of concern. Bottom friction is

especially important because at steady state conditions gradients in the radiation stresses

must be balanced by either gradients in the water surface elevation or by bottom friction

forces. In the absence of friction forces the model could yield continuously accelerating

currents. The method of adding the forcing terms is somewhat arbitrary in that the rational

derivation of the finite differencing scheme of Sheng and Butler cannot be expanded to

include these essentially non- linear terms. The addition of these terms to the model also

limits the latitude of the ratio of in a manner that is not analyzable. That both the

friction and the radiation stresses are nonlinear should be recognized. The non-linearity

of the velocity friction term is quite evident since it is represented by the quadratic term

pFlUiU. The radiation stress terms, which represent excess momentum flux induced by

waves, are functions of the wave height squared, based on linear wave theory. In shallow

water the wave height H is a function of the total water depth, ho + i so that H2 is a

All the forcing terms and the nonlinear acceleration terms are added at the present or

nth time level with the exception of the bottom friction term. The quadratic friction term

52
is finite difference as pF|UnU"|n+ and the friction term as suggested by Longuet-Higgins is
difference as pF luor, U"+1, where luo.,r is in terms of values at the n time level. However,

when using the integral quadratic expression, it is evaluated entirely in terms of the present

nth or known values.

The nonlinear advective acceleration terms are added at the nth time level using central

difference derivatives. In the x direction momentum equation these terms are

Urj U+l,y- U- + U(ij + -L. + Vij+1 + V-Ij+1) U.j+ Ui-1. (4.24)
2AX 4 2AY
and in the y direction equation they are

(Uij + Ui-l,, + Uij+ + Ui-,j+) V+l-V + V + -1 (4.25)

In the model of Vemulakonda (1984) which is based upon the WIFM (Waterways Ex-

periment Station Implicit flooding Model) (Butler 1978) model, a three time level leap-frog

scheme is used. This gives the advantage that all of the forcing terms (and the nonlinear
acceleration terms) are added at the time centered level. However this has some disad-

vantages in that additional computer memory and computational time are required and

leap-frog schemes can exhibit instabilities. In the present study a two time level- scheme

is used and the wave forcing terms as well as the non-linear acceleration terms are added

at the present time level. Since, for all of the applications in this report, the model is run
until steady state conditions are obtained, it is reasoned that there is no great necessity to

properly time center all of the terms. If the use of the model is extended to time dependent

applications, this may need to be corrected.

4.2 Finite Differencing of the Parabolic Wave Equation

Dividing Eq. (3.82) by (C, cos0 + U) yields

A! + i(k cos kcos )A' + A'+ V A'
2 (C, cos + U) (C, cos + U)
+ V) A' i kC,(1 cos )A
2(Ccos9+U) \. 2 (C, cos9+U)

-2(CcosO+) [C( + W( U A' = (4.26)
2(Cpcos 0 + U) [CC,( IV 2(Cgcos 8 + U)

53
The following notation is used to describe the finite differencing of Eq. (4.26): Fj
represents the value of F(i,j) where F is any term and i and j designate the grid row
and grid column. Because of the staggered grid system used for the circulation model as
previously explained, in which the velocity components are specified at the grid edges, the
velocity terms U and V are averaged across the grid so as to produce a grid centered term
for use in the wave model. Thus Uj in the wave model is actually (UV + UJ+) and similarly
for V and other grid designations.
A Crank-Nicholson finite differencing scheme is used for Eq. (4.26). All terms not
involving an x derivative are an average of the value of the term on the i grid row and on
the i +1 grid row. Terms having y derivatives are centered at i, j on the i row and i + 1,
(A+'-A')
on the i+ 1 grid row. The one term with an x derivative is written as A' Thus all
terms are centered at i + ,j. The finite differencing of Eq. (4.26) is written as
A+1 A'
AX
+i [(k~ cos, + cos +F)A+ + (I cos' k; cos#')A]
1 [(C, cos0+ U)/ol' [(C, cosO + U)/l~ .+ + A
+2r X +( (C, cos + + U)/c] J
1 V Itj+ Ai Awt1 ( Va V A\ A -1)
2 (C,cose + U). 2AY + (Ccos+ U) ) 2AY

+ V+1 I(VyU)+) (V j-| A+1
( 2(Ccos +U) I 2AY J

( V 1/)+-(/o)i .
+ 2(C, cos 0 + U) 2AY
i kC,(1 cos2 0) W+ 1 \kC,(- cos2 0) A 1
(C, cos + U) (C, cos a + U)

+ (, )o + 1 U
4 (C, cos + U) +1 (C cos+U)'

1 f *1 A+1 w A ,
+ (C, co j + ( C = 0 (4.27)
S(C, cos a + U){+1 (C, coe 0 + U)

54
where p = CC, for notational convenience and the y derivative of p (A) is written as

( +. [(A -(A)'
[ I ( )2(AY)2
In order to evaluate w.+1 by the method to be explained in section 3.3 an intermediate
value of A+'1 is obtained by the explicit formulation
A.+1 At
SA -- L + i(' cos i k} cos )A
+1 \(c, cos 0 + U)/o]j+ [(C, cos 0 + U)/oj )
AX ((C, cos e + U)/a]f1 + [(C, cos 9 + U)/orJ
+ V V [Al AX_,]
S((C, cos + U)/o). 2AY
0 (V c/of )( (v/) 1 kC(1-cos2 )
2(C cosOe+U ) 2AY 2 (C cosO +U) j .

--- } +] y = 0 (4.29)
2(C, cos + U) ) 2(C,cos + U) 0 (
When Eq. (4.27) is written out using Eq. (4.28) and wi+1 obtained through use of
(4.29) there results for any particular values of i and j an equation for the three unknown
values A!41, Ay+1, and A t. The equation is put into the form

ai A+i + b, A'1+ +- ci Aj;t = dc (4.30)

where the coefficients a,, bi, and c, are known quantities involving the properties of the
environment (depth, current, celerity, etc.) and the dj involves the environmental properties
and the known values of the amplitude on the i grid row. When Eq. (4.30) is written for
all the j values, j equal to 1 through N, for a specified i value, the ensemble of N equations
takes the form of a tridiagonal matrix equation. For illustrative purposes let N be equal to
four. The resulting matrix equation is

a[ b, cl A1+ di
6s 63 cZ A!+ d (4.31)
a3 b3 c3 A'+1 dJ
a4 64 C4 A 4 d4
6A'+

55

Equation (4.31) contains four simultaneous equations with six unknowns. Two bound-

ary conditions should be specified. At each boundary a functional relationship is established

between the grid inside the domain and the grid outside of the domain. For the wave prop-

agation model there are two types of boundary conditions that are readily available for use

!in a computer model. These are a totally reflective boundary condition which would be

used to model a laboratory experiment in a wave basin, and non-reflecting non-interfering

boundary which will allow waves to freely pass through the boundary. For the reflective

boundary, the mathematical statement of the boundary condition is = 0. At the lower

boundary this is finite difference as a = 0 which gives the functional relationship

Al = Ao. (4.32)

At the upper boundary the functional relationship for a reflective boundary is

AN+1 = AN. (4.33)

For the non-reflective boundary the mathematical statement is 9A = -imA where m is

the longshore component of the wavenumber. Since the lateral boundaries are sufficiently

far from any diffraction effects, only refraction and shoaling effects are present. Thus

using Snell's Law m = k sin 0 = k, sin 8o. The non- reflective boundary condition is finite

difference at the lower boundary as

A1 Ao im
A [A + Ao) (4.34)
EY 2

which gives

Ao = Al (4.35)
(1+ imAY)
Similarly at the upper boundary the functional relationship used to eliminate AN+I is

The use of either pair of boundary conditions, (4.32) and (4.33) or (4.35) and (4.36) to

eliminate Ao and AN+I will then yield an N by N tridiagonal matrix equation comprising N

-4 + I

J=j+1

J=J UiJ- .Aj
-4--

J=j-1
y

i I=i-1 I=i I=i+1

Figure 4.1: Definition sketch of grid nomenclature and coordinate system.

equations for N unknowns. Returning to the illustrative example where N is equal to four,

the matrix Eq. (4.31) becomes

bF + 7al cl A,+1 di

as b3 C3 A da
a4 4 + TNC4 A+1 d4
where 71 and 7N represent the functional relationships as described above on the j=1 and

j=N boundaries.

Equation (4.37) is solved easily for the unknown A+'1 using a double sweep method.

In the computer program the subroutine CTRIDA uses an algorithm of Carnahan, Luther,

and Wilkes (1969) to efficiently accomplish the double sweep solution of eq. (4.37).

4.3 Boundary Conditions in the Circulation Model
1
As shown in figure 4.3 a staggered grid system is used so that velocities are expressed

at the grid edge and 4 and the wave properties are defined at the grid center. So that the

wave model will have full grids it is desirable to place the boundaries of the model at the

grid edges on all sides, with the exception of the off -

are an input or boundary condition for the wave model. Figure 1.1 can also be referred to.

The offshore boundary condition that is used in the circulation model is that the mean

57
water surface level is not changed by any of the circulation within the model, mathematically

ETA(1,J) = 0 J = 1-N (4.38)

This condition necessitates that the offshore boundary be sufficiently far from any set-up or

set-down. Other possible offshore conditions are a no-flow condition or a radiation condition.

Thr no flow condition would conserve water within the computational domain so that the

volume of water involved in set- up would be compensated by an equal volume of set-down

offshore. In an open ocean environment the open ocean can be conceived of as an infinite

reservoir. However in modeling laboratory experiments a no fow condition as the offshore

boundary condition would be the most appropriate. A radiation boundary condition as

proposed by Butler and Sheng (1984) imposes the condition

t + C = 0 (4.39)
Ot (z

This condition allows long waves to freely propagate out of the computational domain

without any reflection. This tends to eliminate any seiche within the computational domain

which is induced by the start-up transients. This condition was applied at one time in the

present model but it did not preserve the offshore zero set-up condition and was changed in

favor of the more rigid zero water surface elevation condition. If the model is adapted for

use with time varying wave trains (an adaptation that the model can easily accommodate)

then the radiation condition will become imperative.

At the shoreline the boundary condition is that there is no flow into or out of the beach.

Since the model allows for flooding of the beach the index of the shoreline changes. The

model allows for differential flooding along the shoreline, for example, due to differences of

wave intensity in the lee of structures. For each grid column in the y direction the last grid

is the IWET(J) grid, where IWET(J) is allowed to change with changes -f the mean water

level. The numerical statement of the no flow boundary conditic at tie shoreline is

U(IWET(J)+1,J) = 0 J=1,N

(4.40)

58
This boundary condition is used in other models, though not all models allow for differential

flooding of the shoreline. Yan (1987), whose model was used to investigate waves and

currents around a tidal or fluvial discharge, used the no flow condition away from the

discharge and defined U(IWET(J)+1,J) as equal to the discharge per grid divided by the

water depth at the grids where discharge occurred.

Other conditions are also imposed at the shoreline boundary. These are that the wave

height and the longshore velocity are both zero. These conditions are imposed at the

shoreward grid edge of the IWET(J) grids even though the wave height and longshore

velocities are specified a half grid space form the shoreline. This is done by writing x

derivatives of the radiation stress terms and of the V velocities that need to be positioned

at the center of the IWET(J) grid in terms of three points. The three points are at the

IWET(J)-1 and IWET(J) grids and the zero value at the shoreline at the edge of the

IWET(J) grid.

For the lateral boundary conditions several alternatives are available. Birkemeier and

Dalrymple (1976) and Ebersole and Dalrymple (1979) both used a periodic or repeating

beach. This meant that computations were performed between the second and Nth grid,

and values outside of the computational domain are determined according to

Q(I, 1) = Q(I, N) (4.41)

Q(I, N + 1)= Q(I,2) for all I (4.42)

for any quantity Q.

At one point this was used in the present model until the problem of what should have

been transient edge waves generated by wave obstructions during the start-up were found

to be trapped in the model. The periodic boundary condition allowed for disturbances

propagating out of the domain at one lateral boundary to enter the domain at the other

boundary. For straight and parallel contours, the periodic beach condition .

the present model but with the presence of breakwaters another condition had to be used.

For the lateral boundary conditions there are two choices, either closed or open bound-

59

aries. Closed boundaries are easily applied by the condition of no flow across the boundary.

Such a choice of boundary conditions wouid accurately duplicate the physical conditions ex-

isting in a laboratory basin (or for prototype conditions within a harbor). Open boundaries

are obtained by specifying derivative conditions. In the present study an open boundary

.is obtained by specifying that there are no gradients of the depth and the wave proper-

ties across the lateral boundaries. This necessitates that the lateral boundaries are placed

sufficiently far from any obstructions or any y gradients caused by the obstructions. An

alternative open boundary condition would be that there is no gradient of the flow. This

would, however, put the position of the boundary condition at the center of the bordering

grid rather than at the edge of the grid and could possibly introduce error at the boundary.

4.4 Boundary Conditions in the Wave Model

In the wave model only lateral boundary conditions are needed for computational pur-

poses. The offshore condition is merely the specified wave amplitude and wave angle. No

condition is needed at the shoreline since the parabolic wave model marches grid row by

grid grow and needs only the lateral conditions.

The offshore condition is that each time the wave model is called the wave amplitude

on the first grid row is specified as

A(1,J) = C(t) Ao(J) J=1,N (4.43)

where Ao(J) is defined in the input subroutine according to

Ao(J) = eiksin"oJa (4.44)

where ko is the wave number along the first grid row and 0o is the angle the wave ray

makes with the x axis. This offshore condition necessitates that the water depth be constant

along the offshore grid row. This is another reason for choosing the offshore condition in

iLir :;'':;uunion moaei to De a cc' --,a.: ~ -it tn e presence ,iI wave obstructions and

current deflectors such as jetties it is possible that with other offshore conditions the mean

water surface elevation could vary along the offshore grid row. The use of one condition in

L)5

50. 100. 150. 200. 250.
TIME(SEC.)

Figure 4.2: Start-up function according to Eq. 4.43 with C1=30 and C2=3.

the wave model which is contrary to the condition in the circulation model can introduce

instabilities in the overall model. The interaction of the two models is more sensitive than

either of the models individually.

The C(t) term in equation 4.43 is a start-up function which smoothly goes from zero

to unity so as to bring the incident wave height from zero to H. and minimize any shock

loading of the model. Mathematically it is given by

C(t) = .5 [ + tanh(Z C2)] (4.45)

where C1 and C2 are arbitrary constants to adjust the slope and offset of the curve. Fig-

ure 4.2 shows equation 4.45 for C1 and C2 equal to 30 and 3 respectively.

The lateral boundary conditions are either open or reflective. The reflective condition

is expressed as
aA
S= 0 (4.46)
ay
while the open condition is expressed as
aA
= imA (4.47)
9y

where m is the longshore component of the wavenumber which in the absence of diffraction

effects remains constant

m = kosin 0o = ksin (

(4.48)

- I

Direction of
W Wave Propagation
-- --x

Figure 4.3: Waves passing through a lateral boundary unaffected by the boundary.

according to Snell's law. Equation 4.46 or 4.47 are used to eliminate the amplitude outside

of the domain from the tridiagonal matrix equation. To illustrate equation 4.47 is finite

difference at the lower lateral boundary as

A(I,1) A(1,0) imA(I,1) + A(I,0)
=im (4.49)
AY 2

Solving for A(I,0)

A(I,0) = 1 A(I,1) (4.50)
1+ im-Y

Thus A(I,0) is eliminated in terms of A(I,1). similarly at the upper boundary A(I,N+1) is

eliminated in terms of A(I,N). The same procedure is used with equation 4.46 which yields

the simple relationships

A(I,0) = A(I,1) (4.51)

and

A(I,N+I) = A(I,N)

(4.52)

62
It is important that the boundary conditions 4.46 and 4.47 be used as prescribed

above. The parabolic equation gives N equations containing the N+2 unknowns A(I+1,J)

for J=O,N+1. The two boundary conditions are thus used to reduce the number of un-

knowns to a number equal to the number of equations. If equation 4.47 is written in terms

of A(I+1,1) and A(I+1,2) at the lower boundary and in terms of A(I+I,N-1) and A(I+1,N)

at the upper boundary and the governing wave equation is written for J=2,N-1, then errors

may be introduced at the boundaries. This is especially true for the case of a reflective

boundary. Figure 4.3 shows waves freely passing through the lateral boundaries without

any noticeable reflection for the case of a plane wave with a period of 10 seconds propagating

over constant depth at an angle of 10 degrees to the x axis.

CHAPTER 5
RESULTS

The numerical model combining the iterative interaction of the circulation model and

the wave model as described in the previous chapters was tested for several cases in which

either analytical or experimental results are available. The model was also run for cases

where previous numerical results are available.

The two cases involving analytical solutions are for set-up induced by normally incident

waves on a plane beach and longshore currents generated by incident plane waves at an angle

on a plane beach. The first section of this chapter covers set-up, comparing the results

with the analytical solutions of Longuet-Higgins and Stewart (1964) and the experimental

results of Bowen et al. (1968) and Hansen and Svendsen (1979). The second section covers

longshore currents and compares results with the analytical solution of Longuet-Higgins

(1970a,b) and experimental results of Visser (1982).

The model was run for the case of a rip-channel on a plane beach. Comparisons with the

numerical results of Birkemeier and Dalrymple (1976) and Ebersole and Dalrymple (1979)

are given in the third section. A rip current was also generated for the case of normal waves

breaking on a shore parallel bar with a gap. No comparisons are available for this case.

The model was run for the cases of shore perpendicular and shore parallel breakwaters

represented in the model as infinitesimally thin impermeable barriers. For the shore per-

pendicular barrier, the results are compared with the experimental results obtained at the

Waterways Experiment Station as reported by Hales (1980). For the shore parallel barrier
the numerical results of Liu and Mei (1975) who-did not

include a wave current feedback mechanism.

In the final section the model was used to simulate the laboratory experiments of

64
Gourlay (1974) in which the shoreline in the lee of a breakwater was made to curve and

meet the breakwater as occurs with tombolo formation. The results of the model were in
close agreement with Gourlay's experimental values.

5.1 Wave Set-up

Longuet-Higgins and Stewart (1964) gave a simple but elegant derivation for the wave

induced set-up on a plane beach assuming linear shallow water theory. For the steady state

case of normally incident waves on a plane beach without any y variations and without any

currents the depth integrated equations of motions reduce to

a0i 1 aS,
a + = (5.1)
am pD az

In shallow water for normally incident waves the radiation stress term from linear wave

theory, Eq. (2.28), becomes

33
S.- = E(n + 1)= 2E = -pgH (5.2)

A spilling breaker assumption is made so that the wave height H is a constant ratio to

the water depth D = ho + qi,

H = x~(h, + f7) (5.3)

Substituting equations (5.2) and (5.3) into Eq. (5.1) gives

aq 3 aBPC(ho + )2 3 (aho + ao
at = 1-(ho + O) a = --IC - + (5.4)
ax 16(h0 + q) Bx 8 Qx Bx)

Solving for | gives the result
a4= WK2 aOh (5.5)
am 1+ x2 am
that the slope of the water surface inside the breaker line is a constant times the slope of

the beach.

S" ,.-..?' -. c'" "a* e results of the analytical expression but to
date none, including the present, have matched the experimental results. In figure 5.1 it

is seen that the slope of the mean water surface produced by the wave set-up is nearly

S3.0

S2.5

NUMERICAL SOLUTION 2.0
o EXPERIMENTAL DATA
OF BOWEN, ET AL /> A 1.5s

Break c 10
Point o
Io
S0.5
SWL 0

S-0.5
400 300 200 100 0
DISTANCE FROM STILL-WATER LINE
IN BEACH, x CMS

Figure 5.1: Comparison of the numerical solution of Vemulakonda et al. (1985) for set-up
with experimental data of Bowen (1968). Source: Vemulakonda et al. (1985)

a constant. For the numerical results the set-up starts at the breaking point, whereas in

the experimental results the water surface elevation is nearly constant for some distance

shoreward of the breaking point and then there is a change in the slope of the water surface

elevation. The numerical results follow the analytical solution. This is an area that needs

further investigation. The reason for the difference between the analytical and experimental

results is that the analytical formulation is predicated upon the radiation stress term being

a function of the wave height, so that as the wave height is reduced after the breaking point

so is the radiation stress term. As pointed out by Svendsen (1984) the radiation stress term

remains nearly constant after breaking for a distance roughly equal to five or six times the

breaking depth. This distance is termed by Svendsen as the outer or transition region in

which there is a rapid change of wave shape. In the inner region the wave shape according

to Svendsen remains nearly constant and is characterized by a surface roller or bore. By

incorporai .g th. surface roller Svendsen is ablh wo obtain close agreement for wave height

and set-up to the experimental results of Hansen and Svendsen (1979). However this close

66

agreement is just for the inner region. For the outer region there is no viable solution other

than that the radiation stress term and the water surface elevation remain nearly constant.

Further there is no analytical separation of the outer and inner regions other than the a

posteri definition of where there is a sharp change in the slope of the water surface elevation.

This definition is not useful in a model that seeks to determine the water surface elevations.

i Svendsen's work shows promise and should be pursued; it is not however at this point

viable for inclusion in a numerical circulation model. For one it needs to be expanded for the

case of waves at a angle (this need will be discussed further in the next section). Another

limitation is that Svendsen's model is for monotonically decreasing depths, which is a severe

restriction (for example excluding barred profiles).

Svendsen conceptually addresses the question concerning the spacial difference between

the initiation of the set-up in the analytical and experimental results. The analytical model

assumes that once breaking commences energy dissipation also begins. Svendsen states that

the organized potential energy of the wave prior to breaking is transformed into forward

momentum in the outer region. In keeping with the definition of the radiation stresses as

the flux of excess momentum, i.e. momentum that would not be present in the absence of

the wave, it is easy to conceptualize that the radiation stress term should remain nearly

constant just after breaking. The problem is in assuming that the wave form in the breaking

zone is characterized the same as outside the breaking zone. The author, at present, offers

no solution other than to recognize the existence of the problem.

In the present study results for the set-up on a plane beach with normally incident

waves differ somewhat slightly from the analytical solution as given by Longuet-Higgins and

Stewart. This is because the present model does not use the spilling breaker assumption.

Instead wave energy is dissipated according to a model suggested by Dally (1980) which is

explained in section 3 3. Th'~ia he slope of the mean watr,-w sawrfe :v not a constant times

the bottom slope. A result from the experiments of Hansen mnd Svendsen (1979) show this

clearly in figure 5.2. The figure shows a slight curvature in the wave height after breaking

mm
140

120

100
0 Model
80

60

40 Experiment

20

0 5 10 15 20 /25

.20.

Figure 5.2: Comparison of numerical results from the model with experimental results of
Hansen and Svendsen (1979). Wave set-up and wave heights for a 2 second wave on a beach
slope of 1:34 and incident wave height of 7 centimeters.

and in the slope of the set-up. Although the problem of where the set-up begins is still

present, the curvature of the mean water surface is the same for both cases. It is also quite

clear that the diminution of the wave height is not constant as is assumed by the spilling

breaker assumption.

In figure 5.2 the breaking wave height in the experiment is much greater than in the

model. This is because linear wave theory cannot adequately predict wave shoaling. Ven-

gayil and Kirby (1986) were able to duplicate the shoaling recorded in the experiments of

Hansen and Svendsen (1979) by numerically solving coupled equations for the many har-

monic components. Their results sl i. :- -, _'- -... making,

wave energy was transmitted to the higher frequency components. Their work however is

not adaptable for use in the current model since it is not known how to dissipate energy

P=0

0.0 0.5 1.0 1.5 X 2.0

Figure 5.3: Non-dimensional longshore velocities vresus the non-dimensional distance off-
shore from the analytical solution of Longuet-Higgins. Reprinted from Longuet-Higgins
(1970b), Journal of Geophysical Research, vol. 75(33) page 6793, copyright by the Ameri-
can Geophysical Union.

from the various frequency components once breaking is initiated.

5.2 Longshore Currents

Longuet-Higgins (1970b) gives an analytical solution for the longshore current generated

by obliquely incident waves on a plane beach. The solution involves a spilling breaker

assumption and assumptions concerning the form of the bottom friction and the lateral

mixing. For steady state conditions the y equation of motion reduces to a balance between

the gradient of the Sxy term and the bottom friction with the lateral mixing serving to

distribute the current. Details of the solution can be found in the cited paper.

Figure 5.3 gives the results of the analytical solution of Longuet-Higgins. In the figure

the ordinate V is the non-dimensional velocity v/vo, where v is the actual velocity and vo

is the velocity at the breaking point obtained for '- .'

5x, sin 0
Vo = -(ghs) (5.6)

The abscissa X is the non-dimensional distance z/Zb, where xb is the distance from the

V(mrsec)

* r = 20 1/sec
S Qr = 25 1/sec
S. r = 45 1/sec

S 1 2 *3 4 5 6
Wave Set-Up Line Breaker Line

Figure 5.4: Comparison of experimental results of Visser (1982) and numerical result for
longshore currents due to a 1.1 second wave on a slope of 1:20 for the deep water conditions
Ho = .065m and 00 = 20 degrees.

shore to the breaker line.

Previous models have reproduced the results of the Longuet- Higgins solution by incor-

porating all of the same assumptions (i.e. spilling breakers and bottom friction). Since the

present model does not contain the same assumptions comparison will not be made with the

analytical results of Longuet-Higgins. Instead comparison is made with the experimental

results of Visser (1982). For both the analytical and the numerical results it should be noted

that the same inadequacy is present as was discussed in the previous section. In both cases

it is implied that wave energy dissipation commences as soon as the wave begins breaking.

Thus both analytical and numerical solutions initiate a longshore current starting at the

breaker line. It is reasonable to expect that the same holds for the y component of the

onshore flux of excess momentum as for the x component. That is, that the Sx, term also

remains nearly constant over a certain unspecified distance as the organized momentum flux

of the non- breaking wave becomes turbulent momentum upon breaking. This is illustrated

in figure 5.4 which shows a comparison of the experimental results of Visser (1982) and the

model results for the same test conditions. The maximum '.

both the model and the experiment but the maximum is further offshore in the model. This

can be attributed to the initiation of a longshore current in the model at the breaker line.

70
5.3 Waves On a Rip-Channel

Both Birkemeier and Dalrymple (1976) and Ebersole and Dalrymple (1979) examined

the case of waves incident upon a shore with the presence of a rip-channel. The bottom

topography was given according to a formula first used by Noda et al. (1974)

Sh(x, y)=m zl+ Ae- 20i' sino (y y- xtan l) (5.7)

where 6 is the angle the rip-channel makes with the x axis which is, in their coordinate

system, directed offshore from the still water line, m is the beach slope, A is an amplitude

of the bottom variation and A is the length of the repeating beach. The origin of the

coordinates is at the intersection of the still water line and the trough as is shown in

figure 5.5 which gives the depth contours for the case where the slope is .025 and P is equal

to 30 degrees.

Due to the restriction of the present model that there be no y gradients at the lateral

boundaries the above formulation is modified to be

h()= mzl+ Ae-S20 sinl (y- tan#) (y- tan)< (5.8)
mz J Jx( (y tan #)> (.
Also shown in figure 5.5 are the contours obtained through use of the modified Eq. 5.8

for 6 equal to 30 degrees and a .025 slope. The two formulas differ only in that the repeating

rip channels in the Noda formulation are removed. For application in the present model the

presence of the trough in the farthest offshore grid row is removed so that the depth in the

offshore grid row remains constant in compliance with the requirements of the boundary

conditions as discussed in sections 4.3 and 4.4.

Figure 5.6 shows the results obtained by Ebersole and Dalrymple for the case of a wave

angle of 30 degrees a wave height of 1.0 meters and a period of 4 seconds. Figure 5.7 gives

the results for the present model. Even though the scales of the two figures are diff '-ent it

is evident that there is agreement between the two. The length scale of the eddy is about

30 meters for both cases. In the caption to the figure of his results, Ebersole indicated that

the lateral mixing coefficient may have been large. The present author does not believe

Figure 5.5: Top: Depth contours, in meters, using Eq. 5.7 with the coordinates as shown.
Bottom: Depth contours as used in the present model. Grid spacing is 5 meters.

72

0 6 9 9 9 9 9 9 9 9 9 9 6 9 0 9 9 9

0 9 0 0 0 0 9 0 0 0 0 0 0 0 a 0

0 0 0 ** 0 0 ** **
0 0* 0* 9 ** 0 00 d 0 6 9 0 S 0- 0.** **

*^ -*****'00 # .0* 4 -

24m
Figure 5.6: Velocity vectors from the results of Ebersole and Dalrymple for a wave angle
of 30 degrees. Scale for velocity vectors is 1 inch = 2.87 m/sec. Source: Ebersole and
Dalrymple (1979).

73

. . . . . . . . .-

...... .. .. .*..... .*. .....*

--- ---.-. -.- -------*- ------ -----*--
- - - -- -, -. - --.

0.635m/s -I-25m---

Maximum Vector

Figure 5.7: Velocity vectors obtained from the present model using the same input condi-
tions as for Ebersole and Dalrymple.
- - - - - - -- - I -
- - - ----- - - - - -

- ------------- -- --- -- -------~- -----
----------------------------------~4-----
------------------ -------- --\~-~~----

- - - - - -- - -

- - - -- - -

Figur 6.7 Veloity ector o~aned fom t- - - - ondi

74
I i I I I I 1 1 1 1 1 1 1 I 1111111 1 1 i IiI

..............................-...............

".................................................,. .,t -o **** .......... ..............

.......... .......... .... .....:::::::::::::::::...

----100m ----
Figure 5.8: Mean water surface elevation contours showing that the water surface is lower
above the trough. Values indicated are in millimeters.

this to be the case. In the present study the lateral mixing coefficient was reduced to as

low a value as 2 meter2/second without substantially changing the magnitude of the eddy

adjacent to the rip trough. However Ebersole did find a significant difference between the

case of no lateral mixing and with lateral mixing. Results for the present model for the

case of no lateral mixing are not available since the model sometimes exhibits an instability

without the stabilizing influence of the lateral mixing terms. Figure 5.8 shows the mean

water surface elevation contours obtained from the model which show that the water mean

surface is lower above the trough.

The agreement between the results of the Ebersole model and the present model for

t -' -- if -,- q incident upon a rip-channel should not be surprising since diffraction in

this case is small compared to current and depth refraction.

Another test result is the case of normally incident waves breaking on a bar in which

Figure 5.9: Depth contours in meters for the case of a bar with a gap.

. . . 0.622m /s
-'- "-" -- Maximum Vector

S .. .
. . .. Maximum Vector .

-* a -. -

Figure 5.10: Velocity vectors for the case of waves over a bar with a gap. Dashed lines
indicate depth contours. Grid spacing equal to 10 meters.

76
there is a gap in the bar. Figure 5.9 gives the contour lines for this test case. The topography

was created by adding an inverted parabola on to a plane beach with a slope of 1:40. The gap

was created by suppressing the bar with a symmetric smooth tangent function. Figure 5.10

shows the velocity vectors resulting from a 6 second wave with an incident wave height of 1

meter. The contours of the bottom topography are represented by dashed lines. The driving

force of the circulation is the return flow through the gap in the bar of the water that is

carried over the bar by the breaking process. These results may not be totally correct given

the arguments of section 5.1. If there were no gap in the bar the waves breaking on the bar

would push water shoreward until sufficient set-up was produced. However with a gap in

the bar lateral flow is allowed and set-up cannot be maintained.

Kirby (conversation with the author) has remarked that he has observed the initiation

of a very strong current in the direction of the waves produced by waves breaking on a sand

mound in a small wave basin. This was tested in the model using normally incident waves

on a slope of .02 with a shoal described by the following equation which is the shoal used

in the Delft experiments of Berkhoff et al. (1982).

.45 z < -5.82m
h(x, y)= .45 .02(5.82 + z) ( s)2 + ()2 > 1. (5.9)
.45 .02(5.82 + ) .51 (__)' ()2 + .3 (.)2 + (V)' < 1.0

Figure 5.11 shows the depth contours for this test case. Figure 5.12 shows the velocity

vectors obtained for the case of a barely breaking wave. For the case of a non- breaking

wave the velocity vectors plotted to the same scale are not discernable. Figure 5.13 plots

the max velocity obtained in the shoal region versus the incident wave height. At a wave

height of about 10 centimeters there is a sharp change in the slope of the curve indicating

increased velocities once breaking commences.

5.4 Shore-Perpendicular Breakwater

mIYh ',n.-is t.q tfd hpe f- tc ,'* of waves interacting with a shore-perpendicular

breakwater. The breakwater is represented as an infinitesimally thin impermeable barrier

located in between the J = JBRK and J = JBRK+I grid columns and extending from

Figure 5.11: Deoth contours for the case of normally incident waves upon a shoal. Grid
spacing 05 .. ---- --v, 1 r .' 'e'g.

78

.- .( . .

a. . . .

. . .
.Maxim....m Vector... .

Figure 5.12: Velocity vectors for the case of waves breaking upon a shoal.
9 * * * *
* * * * *

* * * * *

* * * * *

* * * * *
* * * * *

* * . . . .
* * * * *

* * * * * 4

* * , , ,

* * * * * * a
* * * * * * *
* * * * * * *
* * * *

79

-o

-Jo
LLJo

X

o
0
.i

-0.00 5.00 10.00 15.00
WN VE HEIGHT (CM)

Figure 5.13: Maximum velocity across the top of a shoal versus incident wave height for
waves over a shoal on a plane beach. Incident wave height less than 9.1 centimeters are
non-breaking waves. Incident wave heights greater than 9.1 centimeters produced waves
breaking over the shoal.

J= JeBX*
JJ JORK

r

Figure 5.14: Position of the groin in relation to the grid rows and columns.

the shore boundary to the I = IBRK grid row. Figure 5.4 illustrates the position of the

breakwater. Mathematically, this thin impermeable barrier behaves as an internal solid

boundary that prevents any activity on one side of the barrier from being transmitted to

the other side. This requires modifications to both wave and current models.

In the wave model an internal boundary condition is imposed so that at the breakwa-

ter face wave energy is totally reflected. This is accomplished by imposing the reflection

condition

aA

(5.10)

on both sides of the breakwater from grid row I = IBRK to the shoreline, or

A(I,JBRK) = A(I,JBRK+1)

for all I > IBRK

When the wave equation is expressed in finite diff -n, ~, ,* ...

is on the other side of the barrier, is eliminated by virtue Eq. 5.11. Likewise when the wave

equation is applied to the other side of the barrier A(I, JBRK) is also eliminated. Further

- I -- -t- 7 q q

(5.11)

. .

I I A I I I I I I

1

81
reflective conditions are imposed at the barrier by assuming that the derivatives of any

property perpendicular to the breakwater has a zero value at the breakwater.

In the circulation model the internal boundary condition of impermeability at the bar-

rier is imposed by setting the velocity component normal to the barrier equal to zero.

V(I, JBRK+) = 0 for all I > IBRK (5.12)

To insure that current-induced momentum is not diffused across the barrier, subroutines

were introduced in the computation for the x sweep for the two grid columns adjacent to

the barrier.

Two sample cases were run for comparison purposes. The first case is a breakwater

extending 400 meters from the shoreline on a beach with a composite slope of 1 on 10 for

the first one hundred meters and a slope of 1 on 100 beyond that point. This is a geometry

selected by Liu and Mei (1975). A deep water wave height of 1 meter and a deep water

wave angle of 45 degrees was used by Liu and Mei. Using linear shoaling and refraction

this gives a wave height of .84 meters and a wave angle of 31.15 degrees as the wave input

conditions at the offshore grid row 800 meters from shore in a depth of 16.95 meters.

The solution method of Liu and Mei was to solve for the wave field in the absence of

currents and mean water level changes and then use the wave field as a constant forcing in the

solution of the circulation. By ignoring the lateral mixing and advective acceleration terms

and assuming steady state conditions, and by expressing the two horizontal components of

the depth integrated flow as derivatives of a stream function b, the two horizontal equations

of motion were reduced to one Poisson equation

V = /(, y) (5.13)

where the f term contains the forcing due to water level gradients, friction and wave forcing.

This equation represents a boundary value problem with tw .an a, s y ana q. upul.

specifying boundary values an initial guess for the water surface q is made and then Eq. 5.13

is solved for ik. With this new solution for 0 the value of t7 is updated by solving steady

82
state equations of motion. This process is repeated until convergence.

Figures 5.15 and 5.16 show the results for Liu and Mei and the present study respec-

tively. Figure 5.15 shows the contours of the stream function whereas figure 5.16 shows

flow lines which give the direction of the flow but not the intensity of the flow. The present

model shows a separation eddy at the tip of the breakwater which is absent in Liu's results.

Both models show a series of rip cells in the up wave side of the jetty which result from

the interaction of the incident and reflected waves. There are not as many rip cells in the

present model as the weaker ones farther upstream are buried by the stronger longshore

current. Neither the Liu model nor the present one show a rip current at the downwave side

of the jetty. In figure 5.16 the flow lines in the offshore region are directed in the onshore

- offshore direction which is different from the stream function contour lines in figure 5.15.

This is because of a difference in the offshore boundary condition in the two models. In

the Liu model the offshore condition is a constant 4 which means that the flow must be

tangential to the offshore boundary. In the example used in the present model the offshore

condition is a constant water surface (q = 0) which forces the tangential velocity to be zero

at the offshore boundary. An attempt to change the offshore condition to a zero normal

flow condition did not accurately conserve water within the computational domain, even

though convergence was achieved. Figure 5.17 shows the flow lines resulting from using no

onshore flow at the offshore boundary.

The second case is to simulate the experiments conducted at the Waterways Experiment

Station and reported by Hales (1980). A jetty perpendicular to the shoreline extended 15

feet out from the still water line on a plane beach with a slope of 1 on 20. The plane beach

extended beyond the jetty tip until a depth of one foot was reached and further offshore

the laboratory basin was a constant depth of one foot. Waves were sent at an angle by

use of a moveable wave paddle and wave guides which were set up to follow the w e rays

at the ends of the wave paddle. Figure 5.18 shows the experimental set-up tor these tests.

These wave guides in essence made the experiments to be in a closed basin and greatly

I 0 0 i yt A 4 U U ^ o N -
o00 o 00 o o "o o 00 o 0
oo o c o o o o c o o 0

ro 0I

o
0 'S -o
0 0

1 ,

0
l'

Figure 5.15: Results for the numerical model of Liu and Mei (1975). Deepwater wave height
equal to 1 meter, deepwater wave angle equal to 45 degrees. Source: Liu and Mei (1975).
1 -i
-8

1 1 1 1 1 1 1 1 1 1 1 1

equal to 1 meter, deepwater wave angle equal to 45 degrees. Source: Liu and Mei (1975).

i45

0
C
ar

0-

Figure 5.16: Results from the present numerical model using the same conditions as were
used by Liu and Mei. Offshore boundary condition of constant water surface level which
precludes tangential flow at the offshore boundary.

0
ST

figure 5.17: Results from the present numerical model using the same conditions as were
used by Liu and Mei. Offshore boundary condition of no on-offshore flow. Grid spacing
equal to 20 meters.

4 4?5

EL 0.0 MWL

WI

S j / / WAVE GAGE ARRAY

a;/ // /

0 50 IIT
^ ZGD c -I L ^ D EPTH = 0.75
VANES

EL -1.0

11_
F(1L8AT (E)-- .0
PLUNGER FOR M;E
WAVE GENERATOR 0:

I ^^N&4 /VANES

O ^FLAT (EL 1.0)

WAVE ABSORBER
50 FT

Figure 5.18: Experimental set-up for the tests conducted by Hales (1980). Source: Hales
(1980).

87

1.0 FT/SEC NGLE OF INCIDENCE = 30 DEGREE
W-VE PERIOD -.0 SECONDS
Hales. Measurements At one foot intervals. Source: Hales (1980).

S t o i o o
a

mapping techniques that are not yet included in the capability of the model.
the magnitude of the rip current adjacent to the down wave side of the jetty for the cast of
S
A A a

RVERPGE ID-DEPTH VELOCITIES
----- CO4SINED REFRACTION QND DIFFRQCTION'
1.0 FT/SEC ANGLE OF INCIDENCE = 30 DEGREES
WOVE PERIOD a 1.0 SECONDS
WAVE HEIGHT NEAR GENERATOR a 0.139 FEFT
Figure 5.19: Experimental measurements of the mid-depth averaged velocities obtain by
Hales. Measurements at one foot intervals. Source: Hales (1980).

affected some of the circulation patterns away from the jetty. Figure 5.19 shows the average

mid-depth velocities recorded on the down wave side of the jetty. Note the strong offshore

current next to the boundary. Since the boundaries are less than two jetty lengths from the

jetty, one must ask whether the circulation near the jetty was also affected. The numerical

simulation of the Hales experiments was performed using open lateral boundaries, since

to fully represent the experimental conditions in which the lateral boundaries are neither

straight (i.e. curving wave rays) nor perpendicular to the shoreline boundary would require

mapping techniques that are not yet included in the capability of the model.

Figure 5.20 shows a comparison of the experir^ntaJ i' ni-erl-l '-rwlts obtain

the magnitude of the rip current adjacent to the down wave side of the jetty for the cast of

a 1 second wave with an angle of incidence of 30 degrees and a one foot depth wave height

S--- Middepth
S, --- Bottom
0.3
S----- Model Results

0 0.2
O
-J
wu Shoreline -- ,
Wm 0.1 /

I I .
% 0 2 4 6 8 10 12 14
DISTANCE FROM SHORELINE (ft)

End of Breakwater (15 ft from shoreline) -

Figure 5.20: Rip-current velocities adjacent to the down wave side of the jetty. Comparison
of experimental and numerical results.

of .139 ft. One would expect the depth averaged velocity to be somewhere between the

mid-depth velocity and the bottom velocity. The model appears to slightly underpredict

the rip current. In figure 5.19 a circulation cell in the lee of the groin is barely discernable.

In figure 5.21 the circulation cell is very obvious.

5.5 Shore Parallel Emerged Breakwater

In order to solve for the wave and current field resulting from the interaction of waves

and an emerged breakwater several boundary conditions at the breakwater are incorporated

into the model in addition to the boundary conditions already discussed in sections 4.3

and 4.4. The breakwater is represented as being infinitesimally thin and impermeable.

Figure 5.22 shows the position of the breakwater on the shoreward side of the I = IBRK grid

row. The impermeability condition implies that neither waves nor currents pass through

the breakwater. The no current condition is easily imposed in the 0lation mnndel by

making U(IBRK+I,J) equal to zero for the entire length of the breakwater. In addition the

circulation model is modified so that momentum diffusion and derivatives do not take place

89

S.. . . . . . . 11

. . . . .

. . 1111

. .I.. .llll S

..... .. .. .ll.. 51
S .. . . . l

........... . . . . .' i i
. .
.o .* V s l l l lJ 1

. ........ ~.r 5. . .. "

as fo th Ha...e.s. .e.x.p.e .G.i.d spac. to on foot.
............ ............. ............* '

... ....... ..... .. ....... S %

.r~~~. a he1. Hes exeie1 Gr pn eu 5 o *o.
Figure 5.21:. V elocity vectors o th currents o i n rica y for the sa m e cond tions

s fr te H e e. Grid spacing equal to one foot.
: '* * * * : : ; : : : : : 1 '. -
. . . * * * * ,

. . . * * * *
. . * * * '*' "

. * * s s
' ' ^ . * * ^ "
' ' ^ . .. r ~ ,L 1
~~~~ r *r rr l l ************* I
' ' . * * *l
* * . . . * il l
' ' . . * * ll l

MILLISECOND	CLASS.METHOD	MESSAGE
0	sobekcm_page_globals.constructor
0	sobekcm_page_globals.constructor	Application State validated or built
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.constructor	Navigation Object created from URI query string
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.display_item	Retrieving item or group information
0	sobekcm_page_globals.get_entire_collection_hierarchy	Retrieving hierarchy information
0	sobekcm_assistant.get_entire_collection_hierarchy
0	cached_data_manager.retrieve_item_aggregation
0	cached_data_manager.retrieve_item_aggregation	Found item aggregation on local cache
0	item_aggregation_builder.get_item_aggregation	Found 'all' item aggregation in cache
0	system.web.ui.page.page_load (ufdc.page_load)
0	sobekcm_page_globals.constructor.on_page_load
0	html_echo_mainwriter.add_style_references	Adding style references to HTML
0	html_echo_mainwriter.add_text_to_page	Reading the text from the file and echoing back to the output stream
62	html_echo_mainwriter.add_text_to_page	Finished reading and writing the file