UFL/COELTR/080
NUMERICAL MODELING OF WAVEINDUCED
CURRENTS USING A PARABOLIC WAVE EQUATION
by
Harley Stanford Winer
Dissertation
1988
NUMERICAL MODELING OF WAVEINDUCED 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 soulmate 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 Setup.............................
5.2 Longshore Currents .......................
TIONS
.o...
o. ..
. .
. 68
. . . . . . ii
............................
..............
..............
..............
..............
..............
..............
..............
..............
..............
5.3 Waves On a RipChannel ...................
5.4 ShorePerpendicular 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 setup and setdown. 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 Startup 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
setup 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 setup 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 Nondimensional longshore velocities vresus the nondimensional dis
tance offshore from the analytical solution of LonguetHiggins. Reprinted
from LonguetHiggins (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 nonbreaking 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 onoffshore
flow. Grid spacing equal to 20 meters. . . . .... 85
5.18 Experimental setup for the tests conducted by Hales (1980). Source:
Hales (1980). ................... .............. 86
5.19 Experimental measurements of the middepth averaged velocities obtain
by Hales. Measurements at one foot intervals. Source: Hales (1980). 87
5.20 Ripcurrent 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 semiinfinite 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 Setup 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 Setup 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 WAVEINDUCED 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
refractiondiffraction 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 setup and wave generated currents in the lee of a breakwater or headland,"
Proc. 14th Coastal Eng. Conf., Copenhagen, 1974, pp. 19761995).
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. LonguetHiggins and Stewart (1963) give an analytic solution for
wave setup 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. LonguetHiggins (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 setup 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 twodimensional models of fluid flows when the reality is threedimensional. 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 threedimensional flow field using
the NavierStokes 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 prespecified 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 timemarching 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 semiinfinitely 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 nonlinear
advective acceleration terms and the lateral diffusion of current momentum. Although Liu
and Mei produced results for the wave field, currents, and setup resulting from the diffrac
tion of w. ve Yv hrckwaters, their study is not a wavecurrent 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 waveinduced 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=N1
J=N2
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 startup.
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 ydirection 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 CrankNicolson scheme. This will be further discussed in
section 4.2 which describes the finitedifferencing 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 anychanges
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 startup 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 setup.
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 doubleelimination 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,M1
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 gy 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 stillwater 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 refractiondiffraction 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.12.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 twodimensional flow has
the distinct advantage of reducing an essentially intractable problem to a simpler tractable
form.
; Some of the terms in Eqs. (2.12.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 mudwater 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 LonguetHiggins
(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 LonguetHiggins 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 setup 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 LonguetHiggins 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 LonguetHiggins (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. LonguetHiggins wrote
( p av) (2.23)
= % Bz /
in which e. and ey are mixing length coefficients.
LonguetHiggins 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 LonguetHiggins 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
fn1.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 LonguetHiggins 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 NewtonRaphson 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 = ESZ = 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
= 44 + #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.282.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 shortcrested 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.462.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
ig cosh k(ho + z) (2.52)
= B ei (2.52)
a cosh kh
r B = BeiCw (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 LonguetHiggins 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.282.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 LonguetHiggins'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 setup and setdown. 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 setup starting at the break point while the solid line in the lower graph shows
that setup 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 threedimensional 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)
Jhe Jho ,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 Ih. 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 Ih. 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 1h, (3.28)
Using Eq. (3.20) this becomes
f# lho.= f 7 Vhho Ih +tf2 h Vh. Vhho + ef Vhf Vho Ih, (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 1h, 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 lh= 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(pUk U2) 1+ 2) (3.61)
Booij points out that the splitting matrix approach employed by Radder (1979) is equivalent
to the choice
yR = k 1pU
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(PU') (1+ k2(p U2))} = ikk (p U')i + k2(p U2) ) (3.62)
By expanding the pseudooperators 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}= pM + 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 nonlinear term which makes the equation
applicable for weakly nonlinear 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 offshore 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 IcoTIdzfj 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 whib 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.12.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" D6,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 nonlinearity
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 pFUnU"n+ and the friction term as suggested by LonguetHiggins 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 + VIj+1) U.j+ Ui1. (4.24)
2AX 4 2AY
and in the y direction equation they are
(Uij + Uil,, + Uij+ + Ui,j+) V+lV + 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 leapfrog
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
leapfrog schemes can exhibit instabilities. In the present study a two time level scheme
is used and the wave forcing terms as well as the nonlinear 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 CrankNicholson 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(1cos2 )
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 nonreflecting noninterfering
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 nonreflective 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=j1
y
i I=i1 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 = 1N (4.38)
This condition necessitates that the offshore boundary be sufficiently far from any setup or
setdown. Other possible offshore conditions are a noflow 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 setdown
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 startup transients. This condition was applied at one time in the
present model but it did not preserve the offshore zero setup 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 startup 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: Startup 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 startup 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+ imY
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,N1) and A(I+1,N)
at the upper boundary and the governing wave equation is written for J=2,N1, 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 setup 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 setup, comparing the results
with the analytical solutions of LonguetHiggins 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 LonguetHiggins
(1970a,b) and experimental results of Visser (1982).
The model was run for the case of a ripchannel 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) whodid 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 Setup
LonguetHiggins and Stewart (1964) gave a simple but elegant derivation for the wave
induced setup 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 setup 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
S0.5
400 300 200 100 0
DISTANCE FROM STILLWATER LINE
IN BEACH, x CMS
Figure 5.1: Comparison of the numerical solution of Vemulakonda et al. (1985) for setup
with experimental data of Bowen (1968). Source: Vemulakonda et al. (1985)
a constant. For the numerical results the setup 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 setup 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 setup 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 setup on a plane beach with normally incident
waves differ somewhat slightly from the analytical solution as given by LonguetHiggins 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 setup 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 setup. Although the problem of where the setup 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: Nondimensional longshore velocities vresus the nondimensional distance off
shore from the analytical solution of LonguetHiggins. Reprinted from LonguetHiggins
(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
LonguetHiggins (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 LonguetHiggins. In the figure
the ordinate V is the nondimensional 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 nondimensional 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 SetUp 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 LonguetHiggins. 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 RipChannel
Both Birkemeier and Dalrymple (1976) and Ebersole and Dalrymple (1979) examined
the case of waves incident upon a shore with the presence of a ripchannel. 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 ripchannel 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+ AeS20 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 I25m
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 ripchannel 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 setup was produced. However with a gap in
the bar lateral flow is allowed and setup 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 ShorePerpendicular Breakwater
mIYh ',n.is t.q tfd hpe f tc ,'* of waves interacting with a shoreperpendicular
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
nonbreaking 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 currentinduced 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 setup 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 onoffshore 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 setup for the tests conducted by Hales (1980). Source: Hales
(1980).
87
1.0 FT/SEC NGLE OF INCIDENCE = 30 DEGREE
WVE 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 IDDEPTH 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 middepth 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
middepth 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' nierll '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: Ripcurrent 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
middepth 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
