WAVECURRENT INTERACTION AND
QUASITHREEDIMENSIONAL MODELING IN NEARSHORE ZONE
By
JUNG LYUL LEE
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
1993
I
ACKNOWLEDGMENTS
The author would like to express the sincere appreciation and gratitude to my
adviser, Professor Hsiang Wang whose encouragement and guidance with his sage
advise and infinite patience over the years have been a source of inspiration without
which this work would not have been possible.
Gratitude is extended to the members of my dissertation committee, Professor
Michael K. Ochi, Professor Y. Peter Sheng and Professor Ulrich H. Kurzweg for their
guidance and suggestions. Professor Robert G. Dean, Professor Donald M. Sheppard
and Professor Ashish J. Mehta are also acknowledged for the review of the manuscript.
Special thanks go to Dr. LiHwa Lin and Becky Hudson for providing timely
advice and comfortable circumstances. The author owes a lot of debt to Mr. Chul
Hee Yoo, Mrs. Yoo and my fellow students for their hospitality and companionship
during my years in Gainesville.
Finally, the author wishes to truly thank my parents for their love and support
which sustained me throughout studying abroad.
TABLE OF CONTENTS
ACKNOWLEDGMENTS ................................ ii
LIST OF FIGURES ................................ vi
ABSTRACT ........ .. ........................ ix
CHAPTERS
1 INTRODUCTION ............... ..................... 1
1.1 Literature Review .............. ........... .... 3
1.1.1 WaveCurrent Interaction Model .................... 3
1.1.2 Nearshore Circulation Model ................... 5
1.2 Summary of Contents ............................... 7
2 DESCRIPTION OF NEARSHORE HYDRODYNAMIC MODEL ..... 10
2.1 W ave M odels .... ....... .................... 12
2.2 Circulation Models .... ..... .................. 14
2.2.1 Vertical Current Profile ........................ 14
2.2.2 Quasi3D Model .................... ........ 15
3 ENERGY PROPAGATION THROUGH VARYING CURRENT FIELD 17
3.1 Governing Equations and Boundary Conditions ............ 21
3.2 Wave Energy Equation ........... ................. .24
3.2.1 TimeAveraged Wave Energy .................. 26
3.2.2 TimeAveraged Wave Energy Flux ............... 27
3.2.3 TimeAveraged Wave Energy Equation ............. 28
3.3 Wave Action Equation .......................... 28
3.4 Excitation Due to WaveCurrent Interaction .......... .... 29
4 MATHEMATICAL WAVE MODELS ................. .... 32
4.1 Mild Slope Equation ........................... 33
4.2 Derivations of Governing Equation of Each Model ........... 36
4.2.1 HyperbolicType Model I (HM I) ............. ... 36
4.2.2 HyperbolicType Model II (HM II) ............... 38
4.2.3 EllipticType Model I (EM I) .................. 39
4.2.4 EllipticType Model II (EM II) ................. 40
4.2.5 ParabolicType Model (PM) ............ ..... 41
5 NUMERICAL SCHEMES AND COMPARISONS OF WAVE MODELS 44
5.1 Numerical Scheme of Wave Models ..................... 44
5.1.1 Hyperbolic Model I .................. ....... 44
5.1.2 Hyperbolic Model II ....................... 47
5.1.3 Elliptic ModelI ................. ........... 48
5.1.4 Elliptic M odel II ......................... 51
5.1.5 Parabolic Model. .............. ........... 52
5.2 Comparison of Wave Models ..... ......... ........... 53
5.2.1 Computational Difficulty ....................... 53
5.2.2 Wave Shoaling and Refraction .................... 54
5.2.3 W ave Diffraction ......................... 56
5.2.4 Wave Reflection .......................... 56
5.2.5 WaveCurrent Interaction . . . . ... 62
5.2.6 Summary ............................. 62
6 SURF ZONE MODEL ........................... .. 70
6.1 TimeAveraged Wave Energy Equation in Surf Zone . ... 71
6.2 Wave Action Equation in Surf Zone . . . . ... 73
6.3 Wave Height Transformation in Surf Zone . . . ... 74
6.4 Surface Currents in Surf Zone . . . . ... 80
6.5 SetUp and SetDown ........................... 85
7 MATHEMATICAL MODEL FOR WAVEINDUCED CURRENTS .... 88
7.1 TurbulenceAveraged Governing Equations . . . .. 89
7.2 Horizontal Circulation Model .. ............... .. 90
7.2.1 DepthIntegrated and TimeAveraged Equation of Mass .. 90
7.2.2 DepthIntegrated and TimeAveraged Equations of Momentum 92
7.2.3 Roles of New Surface Advection Terms in Radiation Stress 95
7.3 Vertical Circulation Model ........................ 98
7.3.1 Theoretical Undertow Model . . . . ... 100
7.3.2 Longshore Current Model . . . ...... 104
7.3.3 Estimation of Eddy Viscosity . . . ... 105
7.3.4 Theoretical Solutions . . .......... 107
7.3.5 Model Adoption for General ThreeDimensional Topography 112
7.4 Quasi3D M odel .............................. 116
7.4.1 Modification of 2D DepthIntegrated Equations . ... 116
7.4.2 Estimation of Integral Coefficients for Vertical Velocity Profile 120
8 NUMERICAL SCHEME AND RESULT OF CIRCULATION MODEL 122
8.1 Numerical Scheme of Quasi3D Circulation Model . . ... 122
8.1.1 Fractional Step Method . . . . ... 123
8.1.2 FiniteVolume Discretization . .... . .. 124
8.1.3 Approximate Factorizations . ..... . .. 126
8.1.4 Formation to Tridiagonal Matrix . . . ... 130
8.2 Comparison with Gourlay's Experiments .... . .. 133
9 CONCLUSION AND FURTHER STUDY ............... 144
9.1 Conclusions ................................ 144
9.2 Further Study .. ... ........................ 145
APPENDICES
A WAVEAVERAGED ENERGY DENSITY ................. 147
B WAVEAVERAGED ENERGY FLUX ........................ 149
C GREEN'S IDENTITIES CONSISTENT WITH ENERGY EQUATION .151
D NUMERICAL DIFFERENCE OPERATORS ................ 152
BIBLIOGRAPHY ............... ..................... 153
BIOGRAPHICAL SKETCH .................. ........... 158
LIST OF FIGURES
2.1 Structural overview of nearshore hydrodynamic model. . 11
3.1 Time histories measured at the offshore (upper) and inlet points
(lower) during flood tide and ebb tide. . . ... 18
3.2 Change in absolute frequencies and wave heights measured in
physical model........ ................ . 19
5.1 Staggered grid system for hyperbolic model I. . . ... 45
5.2 Subgrid system for elliptic model I. . . . .... 49
5.3 Shoal configuration for comparison of CPU time (concentric cir
cular contours of h/Li). ..................... .. .. 55
5.4 Comparison with the laboratory data of Ito and Tanimoto (1972). 55
5.5 Comparison of wave shoaling. . . . ...... 57
5.6 Comparison of wave refraction (wave height). . . ... 58
5.7 Comparison of wave refraction (wave angle). . . ... 59
5.8 Comparison of wave diffraction for semiinfinite breakwater (00)
between analytic solutions (dotted line) and numerical solutions
(solid contour line of 0.8, 0.6, 0.4 and 0.2 diffraction coeff. from
left)..... . . . . . . . 60
5.9 Comparison of wave diffraction for semiinfinite breakwater (300). 61
5.10 Wave reflection tests against wall. . . . .... 63
5.11 Wave reflection tests against bottom slope. . . ... 64
5.12 Condition of collinear wavecurrent interaction. . ... 65
5.13 Condition of shearing wavecurrent interaction. . ... 65
5.14 Comparison of collinear wavecurrent interaction. . ... 66
5.15 Comparison of shearing wavecurrent interaction (height). . 67
5.16 Comparison of shearing wavecurrent interaction (angle) ...... 68
6.1 Effect of a parameter / on the dimensionless wave height in the
surf zone . . . . . .. . . 77
6.2 Comparison with laboratory experiments presented by Horikawa
and Kuo (1966). ......... ................. 78
6.3 Comparison between theoretical and numerical results. . 79
6.4 Effect of a parameter / on the dimensionless surface onshore cur
rent in the surf zone........................ .. .. ..82
6.5 Effect of a parameter P on the dimensionless surface longshore
current in the surf zone ...................... 83
6.6 Comparison of longshore current with laboratory experiments pre
sented by Visser (1991). ................... .... .. 84
6.7 Effect of a parameter / on the dimensionless setup in the surf zone. 87
7.1 Vertical profile of crossshore currents measured by Hwung and
Lin (1990) ............ ....... .. . . 99
7.2 Vertical profiles of crossshore current. . . ... 109
7.3 Comparison with experiments presented by Hansen and Svensen
(1984). . . . . . . . .. 110
7.4 Effect of the advection term in the undertow model. ...... 111
7.5 Vertical profiles of longshore current. . . . ... 113
7.6 Comparison with experiments presented by Visser (1991). . 114
7.7 Combined threedimensional profiles. . . . ... 115
7.8 Vertical profiles of crossshore current by using bottom shear stress.117
7.9 Comparison with experiments presented by Hansen and Svensen
(1984) . . . . . . . .. 118
8.1 Computational grid of finite volumetric scheme. . ... 125
8.2 Grid system of circulation model. . . . ... 131
8.3 Physical layout of Gourlay's experiment (1974). . ... 135
8.4 Wave height and setup contours represented by Gourlay (1974). 136
8.5 Measured current pattern. . . .... . . 137
8.6 Comparisons of wave height contours. . . . ... 138
8.7 Comparisons of setup contours. . . . ... 139
8.8 Comparisons of depthmean current pattern resulted from quasi
3D m odel ............................ 140
8.9 Comparisons of depthmean current pattern resulted from depth
averaged model. .......... ................ 141
8.10 Comparisons of mean surface current pattern. . . ... 142
8.11 Comparisons of bottom current pattern. . . ... 143
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
WAVECURRENT INTERACTION AND
QUASITHREEDIMENSIONAL MODELING IN NEARSHORE ZONE
By
JUNG LYUL LEE
August 1993
Chairman: Dr. Hsiang Wang
Major Department: Coastal and Oceanographic Engineering
This dissertation represents a study of wavecurrent interaction, surfzone hydro
dynamics, and waveinduced three dimensional current pattern. Based on this study,
a quasithree dimensional circulation model has been developed.
When waves propagate through a region with varying current, an instability in
wave heights has been observed in the tidal inlets. This fact is clarified by both
the surface action and energy transformation which arises due to the change in the
absolute frequency in the varying current field. Consequently, the varying current can
excite the wave field to lead the time variation of wave heights. Two new equations
governing the excitation of wavecurrent field are derived, and form the cornerstones
of this dissertation. These equations are termed here 'wave action equation' and 'the
kinematic conservation of intrinsic angular frequency.'
Five representative wave models are used here; two hyperbolic types, two elliptic
types, and one parabolic type. At this moment, no single model has been proven
perfect or has clearly outperformed the others. The governing equation of each model
is derived from the mild slope equation of hyperbolic type, with varying degrees of
approximation.
I
A simple surf zone model is presented. The model is based on the consideration
of wave energy balance and wave action conservation. It is noteworthy that wave
height, surface onshore current and surface longshore current across the surf zone are
presented in analytical forms for the twodimensional graduallysloped bottoms. This
surf zone model provides the surface current pattern of the vertical circulation model,
and consequently, significantly contributes to solving the threedimensional current
pattern for the general topography.
A hydrodynamic model for the nearshore zone is developed by combining the cir
culation models and the five wave models; one wave model is selected from them ac
cording to the problem of interest. The model is based on the mathematical models of
wavecurrent interaction, and solved by using a fractional step method in conjunction
with the approximate factorization techniques leading to the implicit finite difference
schemes.
I _
CHAPTER 1
INTRODUCTION
The purpose of this dissertation is to present 1) a connected account of the mathe
matical theory of gravity waves and current motions in the nearshore region extending
into the surf zone where the effects of turbulence due to wave breaking become im
portant, and 2) a numerical model applicable to the nearshore zone.
Wave motion has been one of the more challenging problems among the hydrody
namic phenomena since it is basically unsteady and nonlinear. When waves interact
with currents, their interaction induces wave instability as often was observed near
entrances of tidal inlets and harbors. The effect is particularly pronounced when
incoming waves encounter strong opposing currents, producing unstable but nearly
group wave environment. This poses extremely hazardous navigation condition as can
be tested by the mariners. The very nature of the dynamic fluid motion also induces
active sediment transport creating unpredictable shoals and causing large nearshore
morphological changes. Thus, the problem of wavecurrent interaction in nearshore
zone has great engineering significance.
The common assumption that wave energy on a spactially nonuniform current
propagate with a constant velocity, which is equal to the vectoral sum of the current
and wave group velocity, is shown here to be incorrect. In fact, the excitation of
wavecurrent field has been observed near the region where waves meet a nonuniform
current, which might be initiated by shifting of wave frequency. Thus, the wave
motion slowly modulates in the presence of a nonuniform current which influences
the rate of wave energy transport. This slowing time varying wave energy transport
in turn affects the current fields. This interaction eventually could lead to instability
I
2
as will be shown in this dissertation. Under such circumstance, large radiation stress
gradient appears. Under normal circumstance, the primary wave field is unable to
compensate for the effect; the excess wave energy radiates toward far field in the form
of progressive long waves. The radiation of the forced waves is made possible through
surface action, satisfying the surface boundary conditions. In surf zone, the excess
wave energy could be disspated due to turbulence, therefore, the flow field could be
steady as often observed. The excess radiation stress is balanced here by the wave
setup and the waveinduced current such as the longshore, undertow, etc.
To solve the wavecurrent interaction problem two equations governing the dy
namics are to be satisfied, the wave action equation and the wave energy equation.
They are formally derived first. Another basic equation governing the wave kine
matics named here the conservation of intrinsic wave frequency is proposed. The
relationships between the newly derived equation and that of the traditional forms
are examined.
One of the main objectives of this dissertation is to construct a numerical model
that is able to reproduce the threedimensional features of nearshore zone flow such
as surface onshore flow, undertow, longshore current, etc. In the past two decades, we
have witnessed remarkable progress in modeling nearshore hydrodynamics by numer
ical techniques. Commonly, the vertically averaged differential equations have been
employed to predict dynamics in coastal waters and shallow seas. This approach ap
pears to yield reasonable results in both circulation patterns and wave field beyond
the surf zone but fails to produce useful information inside the surf zone. This is be
cause of the highly threedimensional flow structure within it. At present, the vertical
structure of the flow field inside the surf zone has been investigated only for simple
case of twodimensional profiles. Most of the proposed surf zone models whether the
oritical or numerical were also restricted to twodimensional. The model developed
here allows for applications to more general threedimensional topographies.
I
3
The surfzone problem is extremely difficult to formulate although it is a very
narrow zone on oceanic scale and may appear to be uncomplicated. In particular,
the wave breaking and the consequent turbulent flow are exceedingly complex pro
cesses to deal with theoretically. The current pattern in the this zone is certainly
three dimensional. The crossshore component in the surf zone contains a shore
ward mass transport and an offshoredirected return flow. This circulation system
is considered to be the main cause of shore erosion and profile changes under storm
waves and of formation and maintenance. The longshore component is known to be
responsible for long term shoreline changes. Therefore, if one is interested in predict
ing nearshore morphological changes the hydrodynamic model must appropriately
describe this threedimensional current features. Detailed description of the fluid
motion including the wave breaking phenomenon and the subsequent turbulent flow
field is clearly too complicated. The approach used here is to decompose the motion
into different time scales. The large scale motion such as the mean flow can then be
obtained by time averaging over smaller scale motions. In this manner, the effects of
small scale motions are retained in an integrated manner even though the small scale
motions can not be explicitly described. For example, the effects of turbulence on
wave motion are manifested by energy dissipation and the influence of wave motion
on mean flow is the residue flow, excess radiation stress, etc. The surf zone model
is then constructed by the applications of the three fundamental governing equations
derived in this dissertation.
1.1 Literature Review
1.1.1 WaveCurrent Interaction Model
Wavecurrent interaction problems have interested a considerable number of math
ematicians beginning apparently with LonguetHiggins and Stewart (19601961), and
Whitham (1962). LonguetHiggins and Stewart (1960) introduced the concept of ra
diation stress against which the variation in wave energy corresponds to work done by
I
4
the current. This has been a singificant contribution to the advancement of coastal
science from both the theoretical and the practical point of view. The concept of
radiation stress has become fundamental to modern the analysis of water waves in
nearshore region. Further notable contribution was made by Bretherton and Gar
rett (1969) who showed that wave energy flux is no longer conserved owing to the
nonlinear interaction of the radiation stress with the current and introduced the con
servative quantity called wave action. So far, the concept of wave energy flux balance
and the conservation of wave action are the primary tools for analyzing wavecurrent
interactions. The problems of wavecurrent interaction can generally be classified into
two categories; the wave field as affected by the presence of current and the current
field induced by the waves.
The literature concerning wave models is very extensive. The Boussinesq equation
which is based on the timedependent, depthintegrated equations of conservation of
mass and momentum, for instance, has been applied to a number of practical engi
neering problems. The prediction of nearshore waves, however, took a new dimension
with the introduction of the mild slope equation by Berkhoff (1972) which is capa
ble of handling the combined effects of refraction and diffraction. The mild slope
equation is formulated by linear monochromatic waves in areas of moderate bottom
slope. Since then significant progress has been made in computational techniques as
well as model capabilities, notably by Radder (1979), Copeland (1985), Ebersole et
al. (1986), Yoo and O'Connor (1986a), and Dalrymple et al. (1989). The approach
based on the mild slope equation is presently preferred over the more comprehen
sive approach based on the Boussinesq equation since the former can be numerically
solved with faster and more efficient algorithms.
The combined refractiondiffraction effects of wavecurrent interaction were in
cluded by Booij (1981) and Kirby (1984), who derived wave equations by using vari
ational principle. Kirby (1984) made corrections to the equation of Booij (1981) so
I
5
as to yield the wave energy equation suggested by Bretherton and Garrett (1969).
Ohnaka et al. (1988) extended the equation set of hyperbolic type suggested by
Copeland (1985) to wavecurrent interaction; Yoo and O'Connor (1986b) extended
their equation set of hyperbolic type; Kirby (1983) derived the elliptic and parabolic
wave equations; and more recently, Jeong (1990) applied the wavecurrent interac
tion to Ebersole's model (1986). The main differences of these wave models are their
governing wave equations and the associated numerical methods. However, all wave
equations can be derived from the mild slope equation given by Kirby (1984), with
varying degrees of approximations.
1.1.2 Nearshore Circulation Model
The effort in nearshore circulation modeling can generally be broken down into
three related areas, namely, the crossshore circulation and the long shore current
generation, both inside the surf zone and the circulation outside the surf zone. Efforts
of combining them, such as the present dissertation, are also appearing.
In the crossshore surf zone modeling, the usual elements include the fixation of
breaking criterion, the energy decay rate across the surf zone, the wave setup and
finally the crossshore current patterns. Of the various breaking criteria proposed,
Miche's criterion (1951) is still the most durable and widely used for its simplicity.
His model can also extended into the surf zone for energy decay simply assuming
the local wave height to be proportional to the local water depth. This assumption
has also been used extensively for the analytic approach to surf zone models. Le
Mehaute (1962) first suggested the physically appealing approach assuming that the
rate of energy dissipation in a breaking wave is the same as a propagating bore, and
similarly Battjes and Janssen (1978) developed a bore model inspired in the energetic
dissipation produced in a tidal bore. Mizuguchi (1980) introduced a twoparameter
eddy viscosity model which satisfies the two requirements for a stable wave height
in constant water depth and for an approximately constant ratio of wave height to
I
6
water depth on plane beaches. The turbulence model was developed by Izumiya and
Horikawa (1984), based on the similarity law for the structure of turbulence induced
in wave breaking. Dally, Dean, and Dalrymple (1984) presented a wave decay model
that takes the wave reformation into consideration. They assumed that the dissipation
rate is simply proportional to the difference between the local energy flux and that
in the reformation zone, divided by the local water depth. The model also has two
empirical parameters. Their model has been widely used.
The vertical profile of crossshore current pattern has been studied during the
past ten years for the simple twodimensional case on the plane beach. Since the
postulation of the driving mechanism by Nielsen and Sorensen (1970), the analytical
treatment was first appeared by Dally (1980). Svensen (1984) proposed a theoreti
cal undertow model using the first order approximation technique in describing the
breaking waves, and Hansen and Svensen (1984) further considered the effect of the
bottom boundary layer in the undertow. More recently, as the study progresses, sev
eral ideas have been added. Okayasu et al. (1988) estimated undertow profiles based
on the assumed mean shear stress and eddy viscosity, and Yamashita and Tsuchiya
(1990) developed the numerical model which consists of surface and inner layers.
Waveinduced longshore currents within the nearshore zone may be generated
by a number of mechanisms including an oblique wave approaching the shoreline, a
longshore variations in wave breaking height, or the combination of the above. The
need for the development of an adequate theory for longshore currents has attracted
a large number of contributors since the pioneer paper of Putnam and Arthur (1945).
The theory for longshore currents by LonguetHiggins (1970) was a rather elegant
analysis based on the concept of radiation stress. His model and the concept behind
it are widely accepted. Most of the later attempts were more or less modifications
on the original model. The longshore current profile has been proposed by a number
of investigators (de Vriend and Stive, 1987; Svensen and Lorenz, 1989); most of
I
7
them assumed the profile to be equivalent to the logarithmic velocity profile found in
uniform steady stream flows.
Prediction of nearshore waveinduced currents beyond the surf zone has also ad
vanved considerably since some of the earlier developments by Noda et al. (1974) and
Ebersole and Dalrymple (1979). Both of these earlier models were driven by a wave
refraction model but with no current feedback. More recently, Yoo and O'Connor
(1986b) developed a waveinduced circulation model based upon what could be clas
sified as a hyperbolic type wave equation; Yan (1987) and Winer (1988) developed
their interaction models based upon parabolic approximation of the wave equation.
The nearshore circulation, which was predicted by all models listed above, deals pri
marily with the verticallyaveraged longshore current. Recently, however, de Vriend
and Stive (1987) improved the nearshore circulation model by a quasithree dimen
sional approach, which employed a combined depthintegrated current model and a
vertical profile model. They introduced a primary and secondary profile for velocity
variation over depth and assumed the primary velocity profile to be the same in the
crossshore and longshore direction. Svensen and Lorenz (1989) assumed that the
equations in the crossshore and longshore directions could be decoupled and solved
the crossshore and longshore motion independently of long coast with straight bot
tom contours.
1.2 Summary of Contents
It has already been stated that this dissertation presents a set of governing equa
tions for wave and current motions for incompressible fluids with a free surface. The
main force is the gravity but forces induced by turbulence due to wave breaking and
by bottom friction are also included where they become important such as in surf
zone and shallow water, respectively. A numerical model is then developed.
Chapter 2 gives an overview on the main structure of this nearshore hydrody
namic model and briefly introduces the governing equations and numerical schemes
8
involved. In Chapter 3 the governing equations are derived with physical interpreta
tions on the phenomenon of wavecurrent interaction. It begins with a brief descrip
tion of basic wave characteristics and continues with the fundamental equations for
estimating energy and energy flux of flows of invicid and incompressible fluids. Two
dynamic equations governing wave energy transport and wave action conservation are
presented.
Chapter 4 derives the hyperbolic wave equation for the treatment of combined
wave refractiondiffraction. Based on this equation, five different versions of wave
equation suitable for numerical modeling are given as they have their own merit
depending upon the problems to be solved. In Chapter 5, the numerical methods
of the five types of mathematical models are described. These models are evaluated
through mutual comparisons using a number of bench mark cases.
The theory for wave transformation in the surf zone is introduced in Chapter 6
an analytical wave decay model is presented for the simple plane beach cases. The
solution is compared with experiments conducted by Horikawa and Kuo (1966). Based
on the proposed wave action equation and the wave decay model, analytic solutions for
surface currents in onshore and longshore directions are derived for the case of plane
beach. In Chapter 7, the crossshore and longshore circulation model is established
by solving the time averaged mass and momentum equations. Finally, a quasi3D
model is constructed by combining the horizontal circulation model and the vertical
circulation model.
In Chapter 8, The numerical scheme solving the quasi3D model is briefly pre
sented. The governing equations are solved implicitly using a fractional step method
in conjunction with the approximate factorization techniques. The equation of each
step is discretized by the finite volume scheme which yields more accurate and con
servative approximations than schemes based on finite differences. Examples of com
puted nearshore current patterns are presented to demonstrate the applicability of the
I
9
model for typical situations through comparison with laboratory experimental data.
Chapter 9 concludes with the discussion of a few restricted problems overlooked by
the assumptions employed in the model.
I
CHAPTER 2
DESCRIPTION OF NEARSHORE HYDRODYNAMIC MODEL
In the present model, the wave model and circulation model can be combined or
separated through the control of the main program. Therefore, the model is basically
applicable to problems of shallow water wave propagation and nearshore cirlulations
driven by tides, wind and/or wavesinduced radiation stress.
This chapter describes the nearshore hydrodynamic model for wave and current
field in the nearshore zone. The main model is the circulation model for comput
ing mean water surface and mean currents. The wave model is a submodel used
to determine the radiation stress which is required in the circulation model as the
forcing mechanism. The wave model takes offshore wave conditions as input offshore
and propagates into nearshore zone while accounts for various nearshore processes
including shoaling, refraction, diffraction, reflection, and surfzone wave breaking.
The basic wave model is developed for an irrotational flow field with steady mean
currents. Various modifications are incorporated to accommodate for breaking and
bottom friction effects. Figure 1 illustrates the main structure of the waveinduced
nearshore circulation model.
The governing equations for the nearshore circulation model are divided into three
major computational steps; 1) advection, 2) diffusion, and 3) propagation, according
to the fractional step method. Each time step is composed of x and ydirectional
computations by the alternating direction implicit method which speeds up the algo
rithm.
The wave field is affected by mean surface elevation and currents which are the
outputs of the circulation model. As mentioned earlier, five different wave models
Figure 2.1: Structural overview of nearshore hydrodynamic model.
CIRCULATION MODEL
1) vertical circulation model
2) horizontal circulation model
2.1) advection step
2.2) diffusion step
2.3) propagation step
* MAIN PROGRAM *
* read input
* call wave model selected
* call circulation model
* write output
WAVE MODELS INCLUDED
1) hyperbolic model I
2) hyperbolic model II
3) elliptic model I
4) elliptic model II
5) parabolic model
12
are tested. They are categorized into 1) hyperbolic model I, 2) hyperbolic model II,
3) elliptic model I, 4) elliptic model II, and 5) parabolic model. At this moment,
no single model clearly outperforms the others. Therefore, the wave model should be
selected according to the problem of interest. Each of them is programmed separately
and can be chosen to couple with the circulation model.
The computational boundaries and the boundary conditions are automatically
specified by input depth data. Therefore, the model requires no adjustment from
run to run. Detailed description will be given in a separate operation manual. The
flooding and drying of beach face are also taken into account in a twodimensioal
sense as will be explained later.
The following sections provide a brief description of each models and the com
putational methods. All mathematical symbols are to be defined later so that the
structural overview is more clearly illustrated.
2.1 Wave Models
Five representative types of wave models are coupled with the quasithree dimen
sional circulation model. The main differences of the five models are their governing
wave equations and the associated numerical methods. Of the five wave models four
of them were selected from existing literature and one is developed by the author.
Hyperbolic Model I (HM I):
[1+ ( 1)1 + V(U)+V. )CCgV) 0
w[ C at g
+ wgV = 0
at a
Qoy 71
S+ wgV = 0
Hyperbolic Model II (HM II):
aK K k__ COg _
S+ (CgK + U) VK + K VU + s Vh V[V2( )] = 0
at k sinh 2kh 2a a
Sa2 K 2a
( )+ .[(g + U) = 0
at 0 k 01
Elliptic Model I (EM I):
iw{2U V + (V U)} + (U V)(U V)) + (V U)(U Vr)
V (CCgV4) + ( w2 w2 k2CCg) = 0
Elliptic Model II (EM II):
Kc a
V[(CgK + U)] = 0
CCga K2 V. (CCgVa) k2CCga = 0
VxK=0
Parabolic Model (PM):
BA' 1
(Cg, + u)  + i(ko k,)(Cg. + u)A' + a x[(Cg, + u)]A' =
i 9 OA' wOAv A A'
S(CCg ) A' wv
2 .y y 2 0y y
Model Unknowns Numerical scheme
HM I 77 and Vq FDM on a staggered grid system
HM II a and K FDM on a staggered grid system
EM I complex < Combined Gragg's methodFDM
EM II a, K and 0 Generalized LaxFriedrich FDM
PM complex A' CrankNicholson FDM
The wave breaking model nested in to each wave model can be applied by replacing
the CgPCgb to the group velocity in a breaking zone, where f can be estimated by
the dimensionless wave height arriving to the shoreline on a uniform slope, and a
value between 1.1 to 1.4 is suggested depending on the slope, bed condition, etc.
14
2.2 Circulation Models
2.2.1 Vertical Current Profile
A crossshore circulation model is developed to account for the effect of vertically
nonuniform currents. Based upon theoretical solution of simple cases as well as lab
oratory measurements, the velocity profile may be approximated by a second order
parabola.
u = Cz'2 + C,2z' + C,3
v = CylZ'2 + Cy2z' + Cy3
where C1, C2 and C3 are determined in terms of discharge, Q, wave height, H, total
depth, h + c7, turbulentinduced bottom shear stress, TB,tb, and wind stress, Tw as
follows:
77c + h i g H2 w IBxtb
C = h H2ez 16 ax p P
C.2 r= c + h g QH2 rw}
Cz2 +
E, 16 Ox p
3 Q. Cl C.2
1 + h 3 2
C77 + h g OH2 Wy By,tb
C +
yl 2e, 16 9y p p
2 r + h g OH2 rwy
E, 16 Oy p
C3 Qy Cyl C,2
y3 + h 3 2
where,
ez = N Iuorb
77, + h
TB,tb = FIUorb I^
P V (KH/k)
S24(1 7)/L (P Cg/Cgb)
P V. (KH/k)
F C
4rPL (P Cg/Cgb)
where Iorbl is defined as gH/2C.
2.2.2 Quasi3D Model
The quasi3D depthvarying circulation model is governed by the continuity and
momentum equations integrated over depth by parameterizing the vertical structure
of the currents. The model described herein may be considered as a approximation
to the full threedimensional model.
The continuity equation
c Q gH(2k 8L gH2k
7+ + (Q + )+a (Q,+ 8 )=
Wt TX 8o Ty 8o
The xdirectional modified momentum equation
8Q. 0 Q QQy
S+ [ + (h + ce)T.X] + yQ + (h + rc)TYy]
at Ox h + z7, ay h +17,
1 aS,, 1 ias, 8iO 7Wz OTx
+ + + g(h + 77) + = 0
p ax p oy ox p p
The ydirectional modified momentum equation
aQ, o QQ __ o
+ I + (h + 1 c)Tyv1+ [ + (h + "e)T,,y
at ax h + 77 ay h + 77c
1 as 1 as,, 877 TWy TB Y
+ + + g(h + ic) + = 0
p a p ay ay p p
where,
QX = udz, Q = vdz
1 Us
S = E'[n(cos2 + 1) + 2cos 0
2 v
S.y = Sy, = E'[sinO(cosOn + ) + cos 0]
1 Vs
Sy = E'[n(sin2 0 + 1) +2sin 0B
2 C
4C1+ C22 C+C1C.2
S 45 12 6
S= y 4Cxj1Cy C.2Cy2 C+ l Cy2 Cxz2Cyl
y ~ 45 12 12 12
T 4C1 + CylC2
y 45 12 6
16
As noted below, the bottom friction consists of turbulent shear stress at the
bottom, 'B,tb, and bottom friction due to viscous and streaming flows, rB,bf:
TB = TB,tb + B,bf
= F ,wuorblU + FclJorbl(UB + Ustrm)
where F, is the wave friction factor varying over the wave breaking zone as given
earlier, while F, is the current friction factor assumed constant here.
The lateral shear stress is added to the momentum equations as
au 9v
The mixing length coefficient, e, is assumed to be proportional to the distance
from the shoreline, lxz, multiplied by the shallow water phase speed as suggested
by LonguetHiggins (1970).
E, = NlxI\Vgd
It was suggested that a dimensionless constant N, should be less than 0.016. The
ydirectional mixing length coefficient, sy, is assumed to be a constant everywhere.
The main numerical technique follows that proposed by Rosenfeld et al. (1991)
for solving timedependent, threedimensional incompressible NavierStokes equations
in generalized coordinate systems. The governing equations are implicitly solved by
using a fractional step method in conjunction with the approximate factorization
techniques. The equation of each step is discretized by the finite volume scheme which
yields more accurate and conservative approximations than schemes based on finite
differences. The use of a finitevolume method allows one to handle arbitrary grids
and helps to avoid problems with metric singularities and nonconservative nature
that are usually associated with finitedifference methods.
CHAPTER 3
ENERGY PROPAGATION THROUGH VARYING CURRENT FIELD
A common property of waves is their ability to transport energy without the
need of any net material transport. In gravity waves, energy is propagated through
the fluid media via the oscillations of the potential and the kinetic energies. When
waves propagate through a region with currents, their energy is also transported by
the moving fluid. The general appearance of the waves including wave height, length
and period will also be altered. It is commonly observed that when the currents and
waves are in the same direction, waves are lengthened but with reduced wave heights.
Opposing currents, on the other hand, shorten the waves but with increased wave
heights. This latter situation is particularly hazardous for navigation. Moreover,
our recent field and laboratory wave measurements near an inlet entrance seemed to
indicate that waves could become unsteady, or modulated, in a nonuniform current
field. This unsteadiness is more pronounced if waves counter the current. Figure 3.1
shows the time histories of waves measured at offshore and near inlet in the laboratory
(from Wang et al., 1992). It can be seen that under ebb tidal conditions, the wave
field is heavily modulated. Figure 3.2 shows the changes in wave periods (estimated
from crest points) and wave heights under flood and ebb conditions. The reference
wave period and wave height are 1.03 sec and 3.36 cm, respectively, in the far field.
The phenomena described above have not been previously investigated.
In this chapter, two fundamental dynamic equations governing the behavior of
wavecurrent interactions are derived. They are known as the wave energy equation
Flood Tide
Ebb Tide
I' : i j iii I I ''ti Iriii *':'
1.7
Ii t H !
)il 1:I 1 1! !11; i! :!! Ill. : lili!'i I l ii, i
7_ IF.
'7 ,. ,l I
,I I I k I ii 1 !
I ''j
ISLE, "IT
77
1:I 1111~; '''~1~~1al
Figure 3.1: Time histories measured at the offshore (upper) and inlet points (lower)
during flood tide and ebb tide.
Wave Period
0.0 10.0
0.0 10.0
Variation
20.0 30.0 40.0 50.0 60.0
Wave Height Variation
20.0 30.0 40.0 50.0 60.0 70.0
Time (sec)
Figure 3.2: Change in absolute frequencies and wave heights measured in physical
model.
70.0
20
and the wave action equation. The wave energy equation has the final form of:
a ,+ V [(Cg+O ]k)E =0
and the wave action equation that of
OB
+ Vh [V B] = 0
at
where E and B represent the wave energy and wave action, respectively, defined as
S H2 B Pg H2
E=P H, B=
8a 8 0
These two equations differ from the commonly accepted forms such as presented in
Phillips' monograph (1977). In there, the wave energy equation is given in terms of
E' = pgH2/8 as
E + Vh [(Cg + f)E'] + Radiation Stress Term = 0
and in terms of B as
t + V [(Cg + )B] = 0
In Section 3.1, basic problem formulations and appropriate boundary conditions
are given. In Section 3.2, the depthintegrated wave energy equation is derived and the
exact forms of wave energy and wave energy flux are presented in waveaveraged quan
tities. Section 3.3 introduces the wave action equation and the conservation equation
of intrinsic frequency. It will be shown that the wave energy equation originates from
the conservation of energy whereas the wave action equation is derived from the free
surface boundary conditions. The relation between these two new equations and the
traditional forms, along with the differences, are also discussed in Section 3.4. These
two equations form the cornerstones of this dissertation upon which various specific
wavecurrent interaction models both inside and outside surfzones are constructed
in later chapters.
21
3.1 Governing Equations and Boundary Conditions
A velocity potential 0 is assumed to exist such that the water particle velocities
are given by VO where V is the 3dimensional differential operator
V + i j+ k
ax 9y az
The kinematic and dynamic boundary conditions to be satisfied at the free surface,
z = 77, are, respectively, z = r7,
t + Vh Vh77 = 0 (3.1)
Ot + 2(V)2 + gz = C(t) (3.2)
where C(t) may depend on t, but not on the space variables. We may take C(t) 0
without any essential loss of generality. The subscripts t and z indicate the differen
tiations with respect to time and zaxis, respectively. Vh is the horizontal differential
operators defined as
a a
Vi = + 
ax ay
The cartesian coordinate system is used with origin at the still water level, x(x, y)
in the horizontal plane and z directed vertically upwards. The velocity vector,
U(u, v, w), is related to by
a'9 ay, aq
U= x, v =, and w 
9x' 9y Qz
The velocity potential and the free surface displacement are assumed to be com
posed of current and wave components.
(x,z,t) = c(x; t, z) + ,(x, ,t) (3.3)
7(x,t) = (x; t) +6e,(x,t) (3.4)
where e is an undefined factor used to separate the current (such as tidal current,
waveinduced current, etc.) from the wave parts of the velocity potential. The '; t, z'
22
in the current part is used to recognize that q, and r7, may vary slowly over time much
longer than the wave period and it could also accommodate small vertical variations
in currents. Equations (3.1) and (3.2) are then expanded in Taylor series to relate
the boundary conditions at the mean water level z = rc,
[,t + Vh VhA ? (+]Z=,?c + ,7w (77t + Vh VAh7) + v] +... = 0
z= 7C
[,t + 1(V)2 + gz]z=_, + E77w [t + 1(Vr)2 + gz]=,. + .. = 0
Equations (3.3) and (3.4) are substituted into the above equations to give
[(,7c)t + Vhq Vh7.l],=,ze +
e[(77w)t + Vhc Vh,7 + Vhw VhTec ( w)z + 7wV c]z=0, + = 0 (3.5)
[(4c)t + 2(Vc)2 + g7c]Z=,o + E[()w)t + Vc. Vw + 97w]z=,n + = 0(3.6)
2
The above equations are separated for current and wave parts and then truncated to
retain only the significant terms:
(OC) + Vhc Vhr? (3.7)
1 (Vohc)2
77C = K[()t + (3.8)
g 2
(W) Dt + (Vh&)rw + eVh4 Vh77, (3.9)
1 Dw (3.10)
g Dt
where D/Dt /O9t + Vhc Vh.
It should be noted here that the terms retained in the above set of equations are
not necessarily of the same order of magnitude for all general conditions. For instance,
the last term in Eq. (3.9) is, in general, a higher order term than the first two and
only becomes significant when wave diffraction occurs. The (0c)t term in Eq. (3.8),
as will be shown later, becomes significant only under counter current condition. For
slowly varying water depth, the wave part of the velocity potential may be written as
O,(x, z,t) = f(z)q,(x,t) + Enonpropagating modes
(3.11)
23
where f(z) = cosh k(h + z)/ cosh k(h + r77) is a slowly varying function of x, k is
a real value wave number and q, denotes the velocity potential at the mean water
level, termed as 'surface potential.' For progressive waves, the velocity potential and
the free surface displacement can be written in terms of the waveaveraged, slowly
varying quantities as
w,(x,z,t) = f(z)A(x;t)ie'i (3.12)
r7,(x,t) = a(x;t)e'i (3.13)
where a is commonly defined as wave amplitude. The phase function is defined as
0=(K x wt), where K is a wave number vector including the diffraction effects
owing to the retention of the 3rd term in Eq. (3.9), and w is an absolute frequency.
All slowly varying quantities are given here as real numbers. The relation between a
and A can be established by the dynamic free surface boundary condition specified
in Eq. (3.10), which, after substituting Eqs. (3.12) and (3.13) into it, yields
Dt
gae' = { +U.0V}{Aie'O}
!A
= dAe' + ie'{ + 0 VA} (3.14)
at
where ad is the intrinsic frequency including the diffraction effects, defined as ad =
w U K, and U defined as VhOc.
This equation states that a and A should have a phase difference unless we impose
the condition
BA
+ V A = 0 (3.15)
Then, the relation between A and a can be given by the following familiar form
A= (3.16)
Ord
24
Similarly, substituting Eqs. (3.12) and (3.13) into the kinematic free surface boundary
condition given by Eq. (3.9) yields
qi = gk tanh k(h + ,) VA Vc7 (3.17)
Ba
+ V (Ua) + AK V7r = 0 (3.18)
Again, the last term in the above equations reflects the wave diffraction effect and,
under normal circumstances, is of a higher order.
3.2 Wave Energy Equation
In this section, the Eulerian expression of energy equation will be presented. The
Eulerian expression governs the local balance between the rate of change of energy
and the divergence of energy flux at a point. The Euler equation for incompressible
and inviscid fluid flow is
DU
P = V(p + pgz) (3.19)
Dt
Taking the scalar product of U(u, v, w) with the respective terms in Eq. (3.19) and
then summing the three components yields
D [q= U. V(p +pgz) (3.20)
Dt 2
where q2 = (u2 + v2 + w2)/2. With the use of continuity equation, the mechanical
energy conservation equation becomes
S[2] + V. [U( 2 + p + pgz)] = 0 (3.21)
at 2 2
Integrate over water depth,
S + V PU( + p + pgz) dz = 0 (3.22)
The first term in the integrand represents the volumetric rate of energy change and the
second term gives the energy flux through the enclosing surface. The depthintegrated
25
energy equation has been derived from Eq. (3.22) by many including LonguetHiggins
and Stewart (1961) Witham (1962). A brief account is given here.
Using Leibnitz's rule (f2h D fdz = D f h fdz fl,Dr + flhDh), Eq. (3.22) can
be written as
7 +[ 2 j]dz + V. 7 [UQ + p + pgz)]dz
h 2, Jh[ 2
2 2
[U(2 + p + pgz)], V [U(2 + p + pgz)]_h Vh
+[( +p + pgz)l'" = 0 (3.23)
Substituting the kinematic boundary conditions at the free surface and at z = h,
w U Vh77 = 0 (3.24)
0t
WIh + UVhh = 0 (3.25)
and letting p be zero at the water surface, the above equation becomes
Sh + + Vh / [ U + p + pgz)]dz = 0 (3.26)
&J 4t 9t Jh[ ^Pq
Letting the total energy, and the total energy flux, F, as
E= [ dz + 9'2 (3.27)
UF = U + p + pgz)dz (3.28)
Jh 2
Equation (3.26) yields the wellknown energy equation given by LonguetHiggins and
Stewart (1961) and Witham (1962),
at
+ Vh" = 0 (3.29)
Here, is the energy density in the water column per unit surface area and F is the
energy flux through the vertical surface enclosing the water column.
3.2.1 TimeAveraged Wave Energy
The energy per unit surface area given in Eq. (3.27) is expanded in a Taylor series
with respect to the mean water level z = 77c,
S= P 2_ [(V 1 )2 + (z)2] dz + P9[(1c + e,)2 h2] + p Ej [(Vh)2 + (00)2]^
= p [(Vidc + eVhq5)2 + (#cz + ew2] dz + 1pg[(7c + e1w)2 h2]
+P 1E [(Vh, + CVh+t)2 + ( + eq4.)2] (3.30)
Taking time average over the wave period and collecting terms associated with cur
rent( 0(1) terms) and wave ( 0(e2) terms) separately, we obtain the following pair
of equations,
Ec = p/ l[(Vhc)2 + (c)2]dz + pg[ h2 (3.31)
h 2 2
E = p [(Vh )2 + (W) dz + pg9? + pr/[Vh Vh~ ]e (3.32)
h 2 2 32
where Ec and E are, respectively, the mean values of current and wave energy. The
mean wave energy density, which is of primary interest here, can be expressed in
terms of the slowly varying quantities by substituting Eqs. (3.12), (3.13) and (3.17);
E = PwH2 (3.33)
8
where H is the wave height defined as twice the wave amplitude a. Detailed deriva
tions of Eq. (3.33) from Eq. (3.32) are given in Appendix A.
A few general remarks regarding the definition of wave energy density are made
here:
1. E is the energy density directly associated with the 0(e) flcutuation motions
only. It does not take into account contributions associated with mean water
level change. The contribution due to mean water level changes simply moves
the reference from the mean water level to that of no wave condition.
27
2. The last term in Eq. (3.31) represents the contribution due to wavecurrent
interaction. LonguetHiggins and Stewart (1961) and Phillips (1977) all in
cluded this term in the mean flow energy rather than in the mean wave energy.
However, this is truly a O(e2) term from the fluctuation motion.
3. In the absence of current, E reduces to the conventional definition of energy
density in a wave field.
3.2.2 TimeAveraged Wave Energy Flux
By virtue of Bernoulli equation for unsteady flow the energy flux term given in
Eq. (3.28) can be written as
= dz (3.34)
Introducing the current and wave components defined in Eqs. (3.3) and (3.4) and
expanding in Taylor series with respect to the mean water level z = jc, we have
a a
= p Vh(4( + )0.)t(Oc + E6f.)dz pe77[Vh(4O + e)N(4& + e)],
(3.35)
Collecting the terms of O(e2) and taking time averaging over the wave period, we
have the mean wave energy flux, F,
F = p ,7hw dz p.[VOrc + VhVw o] (3.36)
fFh at at at ,,
This mean wave energy flux can be expressed in terms of the slowly varying quantities,
utilizing Eqs. (3.12) and (3.13);
F= [Cg + k] E (3.37)
Detailed derivations of Eq. (3.37) are given in Appendix B.
Again, the quantity of mean wave energy flux given here differs from the conven
tional one in that the mean wave energy density is different and that an additional
term, ctk/w, is included. It will be shown that this additional term is responsible
for the unsteadiness of energy transport under a nonuniform current field.
28
3.2.3 TimeAveraged Wave Energy Equation
By substituting Eq. (3.37) into Eq. (3.29) the timeaveraged wave energy equa
tion is obtained
+ Vh(Cg + k)E = 0 (3.38)
Here the mean wave energy density, E, is as defined in Eq. (3.33).
3.3 Wave Action Equation
In this section, the wave action equation and the conservation equation of intrinsic
frequency are derived from the surface boundary conditions using the same definitions
given in Eq. (3.12), (3.13) and (3.17). Substracting Eq. (3.10)xpgrij from Eq. (3.9)
xpf,, and ignoring the higher order term due to diffraction effect, we obtain
aJ(Pw ) + V (,JP ) P wz g9, = 0 (3.39)
where U, is current velocity of the mean flow at the water surface level. Substituting
Eqs. (3.12) and (3.13) into Eq. (3.40), the following equation is obtained:
9 2k iV ,Pgk tanh k(h + r7:)
(Bie2) + V (U,Bie2') + aBei2 { anh + 1 = 0
where B is defined as
B = pg H (3.40)
8 a
Expanding and separating the harmonic motions,
ie2i' + V. (OB) + .Be2' 2 2 + k tanh (h + l)1 = 0
which yields the dispersion relation, o2 = gk tanh k(h + rc7), and the following wave
action equation:
B+ Vh [,B] = 0 (3.41)
The above equation can also be derived directly from Eqs. (3.16) and (3.18).
29
Now if we substitute Eq. (3.15) into Eq. (3.18), the following equation is obtained,
DaA
aA+ V [UoA] = 0 (3.42)
at
Eliminating A from Eqs. (3.42) and (3.16), we arrive at the final equation that governs
the change of the intrinsic wave frequency in a current field:
+ Vh [,a] = 0 (3.43)
which is termed here as 'the kinematic conservation equation, or simply the conser
vation equation, of intrinsic frequency.'
3.4 Excitation Due to WaveCurrent Interaction
Three basic equations governing wavecurrent interactions have been derived in
this chapter. They consist of two dynamic conservation equations and one kinematic
conservation equation of the following forms:
Energy Conservation Equation
[+ Vh (Cg+ = ck) =0
T
Wave Action Conservation Equation
aB
+ Vh. [i,B] = 0
Kinematic Conservation of Intrinsic Frequency
+ Vh.[UI] = 0
where
E = WH2 B P H2
8o 8 "
The wave energy conservation equation is different from the conventional forms pre
sented by previous investigators. Under the nocurrent condition, this equation reverts
to the familiar form of
aE'
+ Vh [CgE'] = 0
30
where E' = pgH2/8. In the presence of current, it is known that the wave energy
as defined by pgH2/8 is not conserved (LonguetHiggins and Stewart, 1961) but,
instead, the wave action, given by pgH2/8a, is conserved (Bretherton and Garrett,
1969). From Eq. (3.37), we can see that this is only true if the absolute wave frequency
remains constant. Garrett (1969). Although the wave action conservation equation
presented in this chapter deals with a quantity identical to that defined by pgH2/8a
as wave action, the meaning of the equation is very different. It is a surface property
governed only by the surface current condition.
We now like to demonstrate that the absolute wave frequency does not, in general,
remain constant in the presence of nonuniform current. Eliminating B from the
wave energy equation and wave action equation, the following equation governing the
change of absolute frequency results:
S+(Cg + OU k). Vw +Vh. (Cgck)B =0 (3.44)
at + 5 j+
If the absolute frequency is constant over all domains, the first and second terms in
the above equation become zero; then we have
Vh. (Cg ctk)B] =0 (3.45)
For the waves propagating in xdirection,
(Cg t)B = CgoB (3.46)
Ca
Differentiating with x,
9u 0 C
S= (CgB CgoBo)]
where Cgo and ao are values at the area far from the zone of nonuniform currents.
The right hand side in the above equation may not be equal to zero in the nonuniform
current zone. Therefore, the assumption of the constant absolute frequency may yield
the longterm fluctuation of currents. Consequently, when we consider the Doppler
31
equation, in which the intrinsic frequency and wave number are determined by the
kinematic conservation of intrinsic frequency and dispersion relationship, respectively;
w= +Uk
we realize the assumption of the constant absolute frequency may yield the inconsis
tency with the Doppler equation.
Therefore, if the absolute frequency is a variable, we have three unknowns, namely,
the wave height, the absolute wave frequency, and the intrinsic wave frequency in
the above set of equations. The wave number is determined by the intrinsic wave
frequency by use of the dispersion relationship; the vector can be obtained by irro
tationality of wave number vector. Therefore, in theory, the wave properties in the
current field can be solved with given wave conditions at the boundaries. In practice,
the problem is more complicated and simplifying assumptions are to be made as will
be presented in later chapters.
CHAPTER 4
MATHEMATICAL WAVE MODELS
In the past two decades, prediction of nearshore waves took a new dimension
with the introduction of the mild slope equation by Berkhoff (1972) which is capable
of handling the combined effects of refraction and diffraction. Since then significant
progress has been made in computational techniques as well as model capabilities, no
tably by Radder (1979), Copeland (1985), Ebersole et al. (1986), Yoo and O'Connor
(1986a), and Dalrymple et al. (1989). The five types of wave equations are concerned
here since no single model has been proven to be perfect or has clearly outperformed
the others at present; hyperbolic type I by Copeland (1985), hyperbolic type II by
Yoo and O'Connor (1986a), elliptic type I by Berkhoff (1972), elliptic type II by
Ebersole et al. (1986), and the parabolic type by Radder (1979). Type 'I' implies
that the unknown to be solved is expressed in terms of wave mode, while type 'II' in
terms of all waveaveraged quantities, and the parabolic type is expressed in terms of
waveaveraged quantities for the main propagating axis, and wave mode for the the
other axis.
The extension of the mild slope equation to the wavecurrent interaction has been
successfully done by Kirby (1984), and the corresponding equation eventually yielded
the energy equation suggested by Bretherton and Garrett (1969):
D + (V U)i V (CCgV^) + (a2 k2CCg)< = 0
Dt2 Dt
where ^ is the complex velocity potential at the water surface. The equation set
of hyperbolic type I suggested by Copeland (1985) was extended to wavecurrent
interaction by Ohnaka et al. (1988), the equation set type of hyperbolic type II by
33
Yoo and O'Connor (1986b), the elliptic I and parabolic types by Kirby (1984), the
elliptic type II by Cheong (1990). Each can be derived from the mild slope equation
revised by Kirby (1984), with varying degrees of approximations.
In Section 4.1, the mild slope equation revised by Kirby (1984) will also be directly
obtained from the energy equation, assuming that the steady state condition of wave
and current field could be retained for the waveinduced currents which are usually
accompanied by the waves. Discussion will be given to the equation revised by Kirby
(1984) with the alternative mild slope equation of elliptic type. In Section 4.2, the
governing wave equation of each model will be derived from the mild slope equation.
4.1 Mild Slope Equation
The mild slope equation has been derived from Green's second identity. Green's
first and second identities corresponding to the mechanical energy conservation equa
tion and the mild slope equation, respectively, are presented in Appendix C. In this
section, however, the mild slope equation will be derived directly from the energy
equation as described below. Introducing the velocity potentials and free surface dis
placements composed of current and wave components as given in Eqs. (3.33.4), Eq.
(3.26) becomes
9 /'7(V + eVo)2 9
/ h e )2]dz + g(rq, + eq,) ( + Ew)
Vh, [(Voc + eV)(c + es),ldz = 0 (4.1)
Taking Taylor expansion about z = 7c,
S (VOf + eV,)2 9, (Vo + eV.)2 a
JH[ ]dz + I[E1 ]g(77 + g erl ) (qC + e()
at t 2 at
Vh, [[(f + eV+ )(o c + e0)]dz + eq7[(Vc + eVq.)(S0 + e,)t] = 0
Collecting the terms of O(e2) and ignoring the longterm fluctuation of current field:
at [h( 2]d + [7,voc V0.1 + g,7w
Vh [ f[Vm,,]dz + Vcqwswt] = 0 (4.2)
L/h
Substituting Eq. (3.11) into Eq. (4.15),
at fh 2 + f ed f~ [ + [7V a 7 C VO ] + 9n
81 A 2 a t 2azt t
Vh .[ f2 dzV~b^,^ + V = 0 (4.3)
h h
where
If2dz = CCg (4.4)
._h g
1C7 22 k CCg
ffdz dz = k g (4.5)
h g
Therefore, Eq. (4.3) becomes
CCg a (Vh^$)2 a2 k2CCg O (^)2 .
[ I+ 1 + 1[77V
2g at 2 2g Ot 2 t] t
Vh CCgsV. w., + Vc ww, = 0
9
Hereafter, we omit the subscript denoting the wave mode, w.
CCgVhq (Vh))t + qt(02 k2CCg)S + g(7VqC V),t + g279 t
(tVh (CCgV ) CCgVq Vht gVh (Vq,,t) = 0
The first and sixth terms offset each other and expanding the third and last terms,
<^$t(2 k2CCg)} + g77tVjc Vq + g97V7 Vt + 2777

and also offsetting the third and seventh terms,
tVh (CCgV) gAt Vh(V,) = 0
Substituting the dynamic free surface boundary condition of O(e) given in Eq. (3.10)
to q of the third term,
^t(22 k2CCg) g + g77tV c V g(4t + V c. Vq),t
VA (CCgV^) gt Vh(V(OO) = 0
35
^t(,2 k'CCg) gr7t qtVh (CCgV9 ) gt Vh(V 7e1) = 0
Combining the second and fourth terms by use of the kinematic free surface boundary
condition of (e) given in Eq. (3.9),
(0r2 k'CCg)4 Vh (CCgV7 ) g, = 0 (4.6)
Now we consider two kinds of expressions of qz. The first is obtained by substituting
the dynamic free surface boundary condition, Eq. (3.10), into the kinematic boundary
condition, Eq. (3.9),
+ (V ) (4.7)
g Dt2 Dt]
and the second expression is simply obtained from Eq. (3.11) as
z = 0 (4.8)
which eliminates the slowly varying motion through the mean surface. With the first
relation, we finally get the mild slope equation as same as derived by Kirby (1984);
D + (V. ) V. (Cg ) + (CC +(2 k2CCg) = 0 (4.9)
Dt2 Dt
and with the second expression, we get the mild slope equation of elliptic type:
Vh (CCgVq) + k2CCg^ = 0 (4.10)
which can be regarded as the wave equation of elliptic type having the same form as
the case of no current. This may imply that the wave energy is eventually conserved
by the relative group velocity, Cg, in the field of semisteady state although they
may be transported by the absolute group velocity, Cga = Cg + U, causing the
unsteady motion of waves and currents. This hypothesis might be examined through
the laboratory experiments. If the wave and current field were retained in the steady
state, both equations would be consistent with each other by the help of the wave
action equation which describes the slowly varying motion at the mean water level.
However, we select the mild slope equation of first type as a governing equation of
each model, which has been believed to be the one describing the wave motion.
36
4.2 Derivations of Governing Equation of Each Model
In this section, the governing wave equations of the five numerical models selected
are derived from the linearized mild slope equation Eq. (4.9) for waves interacting
with currents. It is rewritten below with the symbol description; the five governing
wave equations include two hyperbolic types, two elliptic types, and one parabolic
type.
D + (V I) V (CCgV^) + (02 k2CCg)^ = 0
Dt2 Dt
where, q is the two dimensional complex potential
C is the relative phase velocity (r/k)
Cg is the relative group velocity (acr/k)
a is the intrinsic frequency (02 = gk tanh kh)
w is the absolute frequency
k is the wave number
h is the water depth
V is the steady current velocity vector (u,v)
V= 9/8x + 9/Oy, omitting the subscript h
where c and k are determined by the Doppler relation, w = c + U k.
4.2.1 HyperbolicType Model I (HM I)
The first governing equation of the hyperbolictype is given as a pair of firstorder
equations which constitute a hyperbolic system similar to those used for the solution
of the shallow water equations. The basic idea in deriving the equation came from
an approach of Ito and Tanimoto (1972) and completed by Copeland (1985) through
the mild slop equation of hyperbolic type given by Booij (1984). Ohnaka, Watanabe
and Isobe (1988) derived a set of timedependent mild slope equations extended to a
wave and current coexisting field. In the presence of strong currents, however, this
governing equation seems to be invalid. Thier equation will be written in the slightly
altered form.
The complex velocity potential at the water surface, q, can be related to the
surface elevation ,ri, by the dynamic boundary condition at the water surface as
follows.
DS
g= D (4.11)
Dt
Assuming (x, t) = q(x)e'O where V(x, t) is a phase function of the wave defined as
= k x wt, we have obtained
g
io.
Differentiating with t gives
= g r (4.12)
Substituting Eqs. (4.11) and (4.12) into Eq. (4.9), we obtain
a Cg dy GCgVq
[1 + ( 1)] + V. (O+) + V (C ) = 0 (4.13)
w C at g
where Vq yields another equation given by taking V on Eq. (4.12) and differentiating
with t,
S+wgV~ =0 (4.14)
at a
0+wgV, =0 (4.15)
Eqs. (4.1315) are a set of time dependent wave equations for linear water waves in
teracting with currents which are derived from the mild slope equation and expressed
in terms of r), V,, and V,$.
Substituting Eqs. (3.12) and (3.13) to the first order PDE set yields the following
wave energy equation, which isn't consistent with the conservation of wave action
given by Brethorton and Garett (1968) (see Equation 4.29).
a [ +V(Ua) + V[Cg] =0
0at aJ c
This inconsistency is caused by Eq. (3.15) which satisfies a slowly varying property
of the dynamic motion.
4.2.2 HyperbolicType Model II (HM II)
The second governing equation of the hyperbolic type is based on the kinematic
and dynamic conservation equations which are defined in terms of the waveperiod
and wavelength average properties of wave motion. While the governing equation
of HM1 reflects the effects of wave reflection over the mild slope, this equation isn't
able to involve any reflection effect due to the assumptions induced by waveaveraged
quantities of progressive waves.
The kinematic conservation of progressive waves can be expressed as
8K
+ Vw = 0 (4.16)
where w is the apparent angular frequency, and K is the wave number vector modified
by diffraction effects as approximated by (details are shown in 'Elliptictype model
II'),
K2 = k2 + ,V2(k) (4.17)
a 01
ignoring the derivative of CCg. a is the wave amplitude and k is a wave number
defined by the Doppler relation which can be approximated by
w = d+ U.K
ScT+UK= gktanhkh+ .K
Substituting this into Eq. (4.16),
8K
K + Var + K VU + U VK + K x (V x U) + 0 x (V x K) = 0 (4.18)
where
V V + Vh = gVk + Vh (4.19)
k h sinh 2kh
Equation (4.19) is now differentiated with respect to x to obtain,
Vk = K VK V2a)] (4.20)
k 2a 0a
39
ignoring the derivative of CCg again. Equation (4.20) is now combined with irrota
tionality of the modified wave number vector to yield the conservation equation of
the wave crest,
4K
+ (Cg + U) VK =F (4.21)
where
ka CCg
F = K VU + K x (V x U) + Vk Vh V[V2()]
sinh kh 2a a
Equation (4.21) states that the rate of change of wave number following a point
moving with the group velocity riding on the convective velocity is equal to F; if F
is zero, namely, U and h are constant, then K is constant following such points. The
timedependent dynamic conservation equation is generally given by
Oa o a2
+ 2V [(Cg + ,)a] =0 (4.22)
0t 2a a
4.2.3 EllipticType Model I (EM I)
The governing equation for the first elliptictype model can be obtained directly
from Eq. (4.9), expressing the two dimensional velocity potential in the linear sta
tionary wave field as follows.
O(x, t) = b(x)e' (4.23)
where q(x) is the surface potential in steady state. Substituting the above equation
into Eq. (4.9) gives
iw{2U V + (V .U)} + (o V)(J V<) + (V )(o ) V)
V (CCgV)) + (o2 w k' CCg)o = 0 (4.24)
The above equation is always valid as long as the wave motion is timeharmonic.
4.2.4 EllipticType Model II (EM II)
If the twodimensional velocity potential is assumed through use of the wave
period and wavelength average property, A, as
O(x, t) = A(x; t)ie'O
the linearized dynamic free surface boundary condition yields, as given in Eq. (3.15),
a
A= g
Urd
where a(x; t) is an amplitude function expressed as 7 = a(x; t)e*v, and od is the
intrinsic frequency including the diffraction effects, defined as ad = w U K. The
linearized kinematic free surface boundary condition also yields, as given in Eq. (3.17),
aO = gk tanh kh VA V.
A
Substituting Eq. (3.12) into Eq. (4.9) yields the real and imaginary parts shown as
follows:
Real part
a aa aa
a )++ OOV.() + +V.(Oa)+
N ad at
CCgK V() + K V(CCg) + CCgV K = 0 (4.25)
U9dd Od
Imaginary part
CCgaK' V (CCgV ) oda + (a2 k2CCg)a = 0 (4.26)
Od od Od
Multiply a/ld to Eq. (4.25) to obtain
2 a a2
( ) + V. [(Cg + U)] =0 (4.27)
at ad Ud
This equation can be expressed in terms of energy. Neglecting the diffraction effect
on the intrinsic frequency,
)E E
() + V [(Cg + 0)] = 0 (4.28)
5 t a 0
41
It states that the wave action (E/a) riding on the convective velocity is conserved by
the group velocity. From Eqs. (4.25) and (4.26), the so called wave energy equation
(real part) and Eikonal equation (imaginary part) become finally
8 (a a2
S) + V [(Cg + )] = 0 (4.29)
CCga K2 V (CCgVa) k CCga = 0 (4.30)
1a r a
The term V (CCgV!) can be expanded as
V (CCgVa) = CCgV2(a + VCCg. V(a) (4.31)
01 a a0
The second term in Eq. (4.30) may be negligible in the mild slope approximation.
Since we have three unknowns, a, K, and K,, we introduce one more equation which
expresses the irrotationality of the gradient of the wave phase function.
Vx K = 0 (4.32)
4.2.5 ParabolicType Model (PM)
The parabolic approximation to the elliptictype equation of harmonic wave mo
tion (Eq. 4.24) is derived by splitting the velocity potential into two components,
S= + + < (4.33)
composed of a transmitted field q+ and a reflected field q. Eq. (4.24) can be written,
through the transformation carried out by Booij (1981), in the following form,
( ) + y = 0 (4.34)
ax Y ax
which is split exactly into two equations,
4 = +i"/+, = ir (4.35)
ox ox
where ( can be chosen to be ( = ( .
Under the assumption that the waves are oriented in xdirection so that ky, 0,
expanding Eq. (4.34),
9o 1 9 
+ (g(CCg U")
aX2 (CCg u2) ax ax
+k 1( +
M )
Mk( = 0
k(CCg u2))
where
Me = [2wk.u +iw(V U)] x(uy) 
Ox oy
+[(CCgv2) + 2iwUV
and matching with Eq. (4.34), we get the two relations
and matching with Eq. (4.34), we get the two relations
(UVL)
9, 0 o,
Nty,
7 = + 2 1/2
 = k k(CCCg u2)/ +
V= kcgu) lk(CCg u2))
(4.37)
(4.38)
Substituting the above relations into Eq. (4.35) yields, for the transmitted field q,
S{k1/2(CCg u2)1/2
ik1k1/2(CCg u2)1/2 ( +
ik~vk (l
By expanding the pseudooperators in the form
( 1
M 1/4
(CCg u2))
M \1/2
+ k(CCg u2)
M2(s u)
+2k2(CCg U2)
results in the simple form
kI {k/2(CCg u2)1/2+} = ikk/2(CCg u2)1/2
TXI ik.
+' M
(l+2k(CCg u2))
(4.40)
where, k,(CCg u2) kCCg u(o + w) = a(Cgk./k + u) wu.
Therefore, expanding Eq. (4.40) and multiplying kl/2(CCg u2)1/2 yields
a(Cg, +u) 8x
S10
28[,l,(c, +u,)1+ =
(4.36)
M 4
k (CCgu2)
(4.39)
S1+
M 1/4
k2(CCg u2)
43
ika(Cgs. +u) u) (+ { UV ) (u )
+ 8(+ v +
+ [(CCg v) + iw + 2iwv } (4.41)
dy By By ay
where Cg, = Cgk/,k.
The above equation gives the modification to the equation of Kirby (1983) to obtain
the more correct energy conservation. Now the equation is modified by introducing
an reference wave number which can usually taken as an average wave number, k, as
=+ = A'eik
where the ydirectional part of the phase function, e f kydy, and the xdirectional
phase modulation due to the difference between the local and averaged wave num
bers, ei(kfk'd), are absorbed into the amplitude function of surface potential, A'.
Substituting this into Eq. (4.41) and after dropping squares of the components of
mean current, the above equation is simplified to
.A' 1
(Cgx + u) + i(k kx)(Cgx + u)A' + x[a(Cg + u)A' =
i 0 OA' wOv 9A'
(Cg A' w (4.42)
2,9y y 2 ay 'y
The above equation admits the instantaneous wave surface to be recovered in the the
following form:
7 = Im{A'eik} (4.43)
9
CHAPTER 5
NUMERICAL SCHEMES AND COMPARISONS OF WAVE MODELS
5.1 Numerical Scheme of Wave Models
We have derived the mild slope equation of hyperbolic type in the last chapter,
however, in terms of practical applications, the equation is not only accurately solv
able, but also attractive since it can only be used for the small domains due to the
present computer capacities. In this chapter, four of the five wave models were se
lected from existing literature and one was developed by the author. These will be
described later. The five models have their own numerical methods depending on
the type of governing equations. All numerical methods applied here fall under the
category of the finite difference method. All spatial finite difference operators shown
in this chapter are given in Appendix D.
5.1.1 Hyperbolic Model I
On the staggered grid system as shown in Figure 5.1, the wave displacement is
specified at the center of each grid, while the gradient of surface velocity potential is
situated at the center of the side of the grid. Eqs. (4.134.15) are solved explicitly for
the xy Cartesian coordinate system as shown in Figure 5.1 as follows.
+ +wgD)() = 0 (5.1)
At
+1 V, +0wgD(!) =0 (5.2)
At 0
Equation (4.13) will be solved by the following explicit numerical equation which
was suggested by Ohnaka et al. (1988) to avoid the numerical damping on the con
45
> > > J>
A A A A
> 0> > 0 >
A A A A
A A A A
> > > >
A A A A
w
a
> V
e
s
.4
0: 7 > : V3, A: Vy<
Figure 5.1: Staggered grid system for hyperbolic model I.
vective term.
0_ Cg 7"+ 77n
{1+ ( 1)} + Dx(u7) + 9,(vM)
wC At
+[Dv(CCgV) + VY(CCgV,)] = 0 (5.3)
g
9
There are three kinds of boundaries. The first is the upwave boundary which
generates the wave field, the second is the nonreflective boundary which passes the
waves without reflection, and the last is the reflective boundary such as breakwater,
jetty, seawall, etc.
The upwave boundary requires some refinement in order that reflected waves
traveling back towards the upwave boundary pass out of the model. This boundary
condition is given by the incident surface elevation and the reflected surface elevation
passing out of the model as follows:
(5.4)
77t(X, yo) = ri (Xo, yo) + 77(xo + Ax, yo)
where,
t7(xo, yo) = a, sin(k cos Ox0 + k sin Oyo wt)
7rt(Zo + Az, yo) = tr(xo + AX, yo) + a, sin[k cos O(xo + Ax) + k sin 0yo wt]
S k cos 0
where a, is an incident wave amplitude and 0 is an angle measured from clockwise
from xaxis.
The nonreflective boundary is satisfied either by characteristic method for the
downwave boundary or by Neumann condition based on Snell's law for the lateral
boundary
lt(xye)M = 77''(x, AXi,)
Dy77 = iky77
The three unknowns, 7, V~$ and Vy,, are given complex to accommodate the appli
cation of the lateral boundary condition and calculation of wave amplitude and wave
angle. Similarly, the reflective boundary is also satisfied by
,t(xe,ye) = Rt7`(xA x,Y,)
Dy77 = iky7]
where R is a weighting factor expressed in terms of reflection coefficient, and ky can
be specified as zero for the perfect reflection. The subscript e implies the end points
along a downwave boundary.
In the wavecurrent interaction, the determination of wave angle is very important
because phase speed, group velocity and intrinsic frequency in the wave equation
are determined through the product of wavenumber and current in the dispersion
equation. The wave angles are calculated at the center of each grid location by the
approximation
0 = tan' L e(V /)J (5.5)
[Ke(VYI/,)J
and the wave height is calculated by
H = 2 /e{} + Im{r} (5.6)
The following relationship between the space interval Ax and Ay and time interval
At must be satisfied to perform stable calculations. The stability condition is given
as
T
At < (L Ax)2 (5.7)
[(Lmax/)2 + (Lmx/Ay)2]1/2
where both L/Ax and L/Ay are recommended to be less than 0(10) to obtain the
accurate solution (Watanabe and Maruyama, 1986).
5.1.2 Hyperbolic Model II
On the staggered grid system as used in hyperbolic model I, the wave amplitude
is specified at the center of each grid, while the wave number vector is situated at the
center of the side of the grid.
The upstream difference scheme has been found to be an excellent choice for
the present system of equations (Yoo and O'Connor, 1986a and 1986b). The finite
difference form of the kinematic wave equation (4.21) for x component of wave number
vector is then represented by
Kx + (Cg + u)D KE + (Cg + v)D;Ky + KD*u + KyD*v
At
+ ka Vh a ,[DrD () + ()] = 0 (5.8)
and for the y component of wave number vector,
K"+1 Kn n K
y +1A + (Cg + u)DZK. + (Cg + v)DV Ky + KD u + KD v
At k
km CCg a a
+ Dh D,[D.(a) + D,,()] = 0 (5.9)
sinh 2kh 2a 0 a
where all finite difference operators indicated by VD are given in finite difference form in
the Appendix. The K vector along the side boundaries is approximated by the Snell's
48
law. The finite difference form of the wave amplitude equation (4.22) is represented
by
(a/o')n+1  (a/Co)n o D{[(Cgx +U a2] ) + Dy[(Cg + v)} = 0 (5.10)
At a k a k a 0 (5.10)
The stability condition in the small diffraction zone may be determined by
T
at <; (5.11)
A< [(naL,,/AX)2 + (nL,/Ay)2]/2 (511)
where na = Cga/Ca with subscript a indicating the absolute.
5.1.3 Elliptic Model I
Equation (4.24) can be treated as an ordinary differential equation for as given
below, so that the ydirectional difference operator, D, is explicitly approximated by
using either a finite difference method or a fast Fourier transform method. In this
section, only the finite difference method is employed.
(u2 CCg).S. + {(2iwu + 2uu, + uv + uvy (CCg),}x
+2uvVy,(.) + {(2iwv + 2vvx + uv, + uv (CCg)y}Dy,( )
+(v2 CCg)Dy,(y) + {iw(u + v,) + ,2 2 k2CCg}O = 0 (5.12)
where, the x axis is suggested to be selected for the main direction of wave propagation
between x and y axis. We convert the above equation to a pair of firstorder equations
by the simple expedient of defining the derivative as a second function.
tx =
i = C [{2iwu + 2uu. + uyv + uvy (CCg),}i1i
CCg U2
+2uvDy,(1) + +()]
where
F(O) = {(2iwv + 2vv, + uv. + u v (CCg)y}VDy + (v2 CCg)DVyy
+{iw(u, + v,) + 2 W2 k'CCg}
sub grids
Ill
i i
i I I
w
a
v
e
s
4
Figure 5.2: Subgrid system for elliptic model I.
In this study, the ordinary differential equations are numerically solved by Gragg's
method whose main algorithm for a differential equation /I(x) = f(x, O(x)) is given
as
Yl = i + hf(xil, i1)
Yj+i = yj1 + 2hf(xi_ + jh, yj)
i (yn + Yn1 + hf(zi, y.))/2
,j = 1,2,..,n 1
where h is a subgrid space defined as h = Ax/n as shown in Figure 5.2.
The upwave boundary condition is merely the specified complex determined
by the incident wave amplitude and wave angle. The lateral boundary conditions
are either nonreflective or reflective. The nonreflective boundary condition can be
expressed as
,y = ikby,
where ky = k sin 0 = ko sin 0,
(5.13)
according to Snell's law in the absence of diffraction effects. While the reflective
boundary condition is expressed as
0 = 0 i.e., ky = 0 (5.14)
50
Using the finite central differences on row i yields
ik, Ay
2+1  (j+1 + i) (5.15)
The lateral boundaries are assumed to be located between the grid columns (1,2) and
(NY1,NY). Therefore, < values at the lower lateral boundary can be given as
=(1 ikAy/2)
d' = (' (5.16)
S (1 + ikAy/2) (56)
on row i. The same procedure is used for the upper boundary.
~. (1 + ikAy/2) ,
'Y = (1 i y/2) 1 (5.17)
Equation (5.13) can also be approximated in the RungeKutta method (Dalrymple
et al., 1988) by
6 = ikfyAy (kA iy)) + 1(kyAy)4] s
l1kk~y)) +5.18)
A = 1 + ikyy 2(k yy)2 iz(ky)3 + 4(k y)4] (5.18)
If there is any reflective structure posed in the ydirection, the direction of the
reflected waves is the mirror image of that of the incident wave. Since the unknown
in this model is the complex surface potential, the reflected wave field can be easily
specified as the conjugate by tracing the computation backward.
The wave angle is calculated by
0 = tan1 ( ) (5.19)
where
=m K= S (5.20)
where,S = K x = tan'[m()/Ze()].
The wave height is calculated easily by
H = 2t ,e{}2 +Im{f }2 (5.21)
9 9
51
The following condition for stability is roughly obtained by taking the von Neu
mann stability analysis on the Helmholtz equation as
Ay = Lmaz/7r (5.22)
where L is a wave length (Panchang et al., 1988).
5.1.4 Elliptic Model II
By introducing wave angle, 0, Eq. (4.30) can be expressed as
AK2 C = 0
where
A = CCg C = DJ[CCgD.( )] + D,[CCgD,()] + k2CCg
a 0 0 a
The solution of K is simply
K=
Using the Generalized LaxFriedrich (LF) scheme which is identical to that used
by Perlin and Dean (1983), 0 can be found by the following expression of Eq. (4.32),
Dt [K sin 0] yD[K cos0] = 0 (5.23)
where
rt[] (1 r)[],+.j + 0.5r([],+l,j_ + [],i+1j+,) []i,j
x Ax
where r is a value between 0 and 1 as a kind of weighting parameter. This param
eter is used to enhance the stability of the numerical scheme. However, 7 = 0 is
recommended for the most accurate results. The above equation for the solution of
0 is solved row by row using implicit FDM. The 0 value along the side boundaries is
approximated by Snell's law.
52
The wave amplitudes can also be found by applying the generalized LF scheme
to Eq. (4.29),
2 a2
D a [ (CCgK cos 0 + ua) + D, [(CCgK sin 0 + vo) =0 (5.24)
The above equation for the solution of a is also solved row by row using implicit FDM
until a certain accuracy of the wave amplitude is achieved.
5.1.5 Parabolic Model
The parabolictype wave equation (4.42) is written in the finite difference form
using the CrankNicolson scheme.
a(Cg + u)D.A' + i(k, k.)(CgS + u)A' + 2D[oa(Cg + u)]A' =
+Dy(CCgDyA') ,yvA' wv~ ,A' (5.25)
Equation (5.24) is now expressed in the tridiagonal matrix form which can be easily
solved by the double sweep method. The first sweep is only for determining the group
velocity of x component.
Presuming Snell's law to hold along the nonreflective lateral boundary
A' = ikA' where ky = k sin 0 = ko sin 0o (5.26)
For the lateral reflective boundary ky can also be given according to the ratio of
reflection, while for the xdirectional reflective boundary this wave equation is in
escapable. Using the finite central difference on row i, the expressions for A' and
A'Ny are obtained as similarly as done in Section 5.1.3.
The wave angle is calculated by
0 = tan1 K) (5.27)
where
K,, k 1 (K/k)2, K, = VDS (5.28)
53
where S = f Kdx kox + Ky = tanr[Im(A')/Re(A')].
The wave height is calculated easily by
H = 2+Ze{ A}2 +Im{A' (5.29)
g g
5.2 Comparison of Wave Models
All wave models will be evaluated through the comparison of computational diffi
culty and capability of handling wave characteristics such as wave shoaling, refraction,
diffraction, reflection, wavecurrent interaction, etc. These characteristics are mainly
due to the bottom variation and current, and occur simultaneously. We restrict all
comparisons to ideal cases in order to give the mathematical solution or expression if
possible.
5.2.1 Computational Difficulty
The degree of computational difficulty is measured in terms of stability as previ
ously provided, as well as CPU time. The stability criteria given below are obtained
for ideal cases only. Therefore, they are not general as well as not vigorous.
Table 5.1
Model Stability
HM I At < T/[(Lm,,/Ax)2 + (Lma/Ay)2]1/2
HM II At < T/[(n,Lmax/Ax)2 + (nLmax,/Ay)2]/2
EM I Ay > Lma/Ir for central difference method
EM II stable but convergable when Ay > Lma,/r
PM stable
n, = Cga/Ca (with subscript a indicating the absolute)
The comparison for the computational time is also not general. Rather, a specific
configuration as shown in Figure 5.3 is used as the test bench mark. This configuration
is a circular shoal used by Ito and Tanimoto (1972) in their laboratory experiment
to study combined diffraction and refraction. This configuration has been cited by
many authors for verification purposes. Here, the same grids and same accuracy
criteria are used in all models. Wave heights along three crosssections as shown in
54
Figure 5.4 are compared with the laboratory data of Ito and Tanimoto. It is seen that
the wave models of hyperbolic type II and elliptic type II, which took the averaging
procedure, don't reproduce local minima and diffraction lobes differently from the
other models. Kirby (1986) showed that the absence of the minima and lobes may
result from the governing equations of models obtained by taking the time average
on periodic motion. The CPU time on a VAX8350 computer and the values of the
agreement parameter are given below.
Table 5.2
Model CPU time d(Sec.1) d(Sec.2) d(Sec.3)
HM I 15 min 0.98 0.97 0.95
HM II 12 min 0.97 0.97 0.96
EM I 24 sec 0.96 0.97 0.94
EM II 5 min 0.98 0.97 0.96
PM 17 sec 0.97 0.97 0.95
The models of ordinary (EM I) and parabolic (PM) types provide the fast solutions
showing close agreement when compared with the other models. The agreement is
based on an index, d, given here as an agreement parameter (Willmott, 1981):
EN(P, 0.)2
d = 1 E IP OS + 0O 0)2
where Pi is the numerical value, 0, is the theoretical or observed value and 0 is the
mean of the variates Oi. The values for d vary between 0 and 1.0, with 1.0 indicating
perfect agreement.
5.2.2 Wave Shoaling and Refraction
The phenomenon of wave shoaling and refraction occurs over variable topography
or current field similar to that which occurs in optics and acoustics. However, the
comparison of these phenomena is restricted over the variable topography. Numerical
results were compared with the analytical solution based on the energy conservation
equation and Snell's law for waves propagating over a uniform slope. The input data
are uniformly given to each model as follows:
SEC.2 SEC.3
D 1 2 3 4 5 a 7 a
i/Li
Figure 5.3: Shoal configuration for comparison of
tours of hiLi).
CPU time (concentric circular con
2 H"11 I
2 H 11r
2C EM It
P"
0
X/Li SEC. I
2 H 11
2 EM11 If
EM 11
T/Ij SEC.2
0 1 2 3 q 5
T/Ll SEC.3
Figure 5.4: Comparison with the laboratory data of Ito and Tanimoto (1972).
4
2
SEC.1
NX NY Ax(m) Ay(m) T(sec)
101 21 0.02 0.14 0.8
The time step in the hyperbolic models is fixed at 0.01 sec.
As shown in figures 5.55.7, all models except hyperbolic model I produce results
of good agreement. Hyperbolic model I, on the other hand, induces periodical fluc
tuations. The numerical error appears to be related to the ratio of grid size to wave
length. As the wave length shortens towards shoreline the error becomes larger and
also propagates upwave as time progresses. The numerical results were taken along a
center grid line in x axis.
5.2.3 Wave Diffraction
Wave diffraction is significant in zones causing the height variation, such as behind
a semiinfinite breakwater or a groin facing obliquely to the wave approaching. In this
section, the model capability of handling wave diffraction was evaluated by comparing
wave height with the analytical solution given by Wiegel (1962) for a semiinfinite
breakwater. The input data are uniformly given as NX=91, NY=75, Az=0.04 m
(=0.1 L), Ay=0.08 m for T=0.511 sec except elliptic model II where NY=38 and
Ay=0.16 m were used to avoid numerical instability. The time step is 0.01 sec in the
hyperbolic models.
Figure 5.8 shows the comparisons for waves approaching normal to the breakwater
axis. All models appear to agree well with the analytical result. For the 300 angle to
the normal, however, only hyperbolic model I and elliptic model II perform adequately
(Fig.5.9). The performance in general can be improved by reducing Ay, with the
exception of elliptic model II which is almost stable regardless of the size of Ay.
5.2.4 Wave Reflection
Wave reflection was tested for the case of waves approaching a seawall at 00 and
30* in constant deep water depth of 3 m, using the listed input conditions,
57
SHOALING
1.6
1.4 Theory
o 000 00 00 HM I
001.2
1.0 0
0.8
1.4 Theory
.1M0 HM II & EM II
1.2
1.0
0.8
1.4 L ^Theory
1.4 Theory
o EM I
. 1.2
1.0
! O0PMI
1.2
1.0 
0.8
0.1 0.2 0.3 0.4 0.5
ko h
Figure 5.5: Comparison of wave shoaling.
58
REFRACTION (height)
1.6
1.4 Theory
o OHMI
1.2
O 0
1.0 .
0.8
1.4 Theory
o C HM II & EM 11II
S1.2
1.0
0.8
1.4 Theory
o0 EM I
1.2 
1.0
0.8
1.4 Theory
o 10 PM I
1.2
1.0 
0.8
0.1 0.2 0.3 0.4 0.5
ko h
Figure 5.6: Comparison of wave refraction (wave height).
59
REFRACTION (angle)
40.0
30.0
S20.0
J g Theory
10.0 O HM I
0.0
30.0
S20.0
S Theory
Z
10.0 0 HM II & EM II
0.0
30.0
20.0
3 Theory
Z
10.0 0 OEM I
0.0
S10.0 ( PM I
0.0 1 1 i
0.1 0.2 0.3 0.4 0.5
ko h
Figure 5.7: Comparison of wave refraction (wave angle).
Hyperbolic Model I
8
a
0 T
4 3 2 1 0 1 2 3 4
8 1 c Y'
32 1 0 1 2 3 4
y/L
Elliptic Model II
I /
\ l /
I/
4 3 2 1 0 1 2 3
y/L
Elliptic Model I
4 3 2 1 0 1 2 3 4
y/L
Parabolic Model
a ,
4 \ *mI I.
' 4 I I
2 I
0
4 3 2 1 0 1 2 3 4
y/L
Figure 5.8: Comparison of wave diffraction for semiinfinite breakwater (00) between
analytic solutions (dotted line) and numerical solutions (solid contour line of 0.8, 0.6,
0.4 and 0.2 diffraction coeff. from left).
a
7
a
5
4
1
3
2
0
Hyperbolic Model I
4
3
2
0
i
o
4 3 2 1 0 1 2 3 4
y/L
Elliptic Model II
4 3 2 1 0 1 2 3 4
y/L
4 3 2 1 0 1 2 3 4
y/L
Parabolic Model
4 3 2 1 0 1 2 3 4
y/L
Figure 5.9: Comparison of wave diffraction for semiinfinite breakwater (300).
Elliptic Model I
input data NX NY Ax(m) Ay(m) At(sec) T(sec)
HM 61 21 0.08 0.08 0.02 1.0
EM I 61 21 0.08 0.50 1.0
Owing to the finite grid size and time step a numerical error is also expected.
Figure 5.10 shows that the numerical results, on the whole, agree well with theory
for both 00 and 300 wave angles. The hyperbolic model tends to yield slightly larger
error in wave height, whereas the elliptic model I produces slightly larger phase error.
The wave reflection against the bottom slope was also compared with the 3
dimensional numerical solution represented by Booij (1983). Both models run in this
study reflect reasonable agreement as shown in Figure 5.11.
5.2.5 WaveCurrent Interaction
Wavecurrent interaction is compared for cases of collinear current (Figure 5.12)
and wave refraction due to the shearing current (Figure 5.13), both in constant deep
water depth of 3 m. The analytic solution for the shearing current is given by Longuet
Higgins and Stewart (1961). The given wave conditions are Hi=0.1 m at the upwave
boundary and T=1 sec. Waves are allowed to freely pass through the downwave
boundaries. The input data are uniform with NX=101, NY=21, Ax=0.1 m Ay=0.6
m for the elliptic models and Ay=0.1 m for the rest. At, whenever applicable, is
taken as 0.01 sec. The comparisons with analytical solutions are given in Figures
5.1416. For the collinear case, all except hyperbolic model I performed adequately.
For noncollinear cases, hyperbolic model II and elliptic model II yield good results;
the rest all produce varying degrees of inconsistency.
5.2.6 Summary
Each model was evaluated or run on a number of bench mark cases. The final
evaluations with assigned rankings are given in the following matrix:
63
REFLECTION AGAINST SEAWALL
4.0
3.0
2.0
1.0
0.0
3.0
2.0
1.0
0.0
3.0
2.0
1.0
0.0
3.0
0.5 1.0
1.5 2.0
2.5 3.0
x/L
Figure 5.10: Wave reflection tests against wall.
0.01
Slope
Figure 5.11: Wave reflection tests against bottom slope.
z
77 v
Figure 5.12: Condition of collinear wavecurrent interaction.
x
x I 
yo
waves i
Figure 5.13: Condition of shearing wavecurrent interaction.
66
COLLINEAR
3.0
2.5
2.0
1.5
1.0
0.5
0.0
2.5
2.0
* 1.5
1.0
0.5
0.0
2.5
2.0
1.5
1.0
0.5
0.0
2.5
2.0
1.5
1.0
0.5
0.0
0.3 0.2
0.1 0.0 0.1 0.2 0.3 0.4 0.5
u/Co
Figure 5.14: Comparison of collinear wavecurrent interaction.
67
SHEARING (height)
3.0
2.5
2.5  Theory
2.0
O 0 HM I
1.5
1.0
0.5
0.0
2.5 Theory
2.0
o OD HM II & EM II
1.5
1.0
0.5
0.0
2.5 Theory
2.0
o 0 EM I
1.5
1.0
0.5
0.0
2.5 Theory
2.0
o 0 PM I
1.5
1.0
0.5
0.0
0.3 0.2 0.1 0.0 0.1 0.2 0.3
v/Co
Figure 5.15: Comparison of shearing wavecurrent interaction (height).
68
SHEARING (angle)
50.0
40.0
30.0
20.0
40.0
30.0
20.0
40.0
30.0
20.0
40.0
30.0
20.0
0.1 0.0
v/Co
0.1 0.2
Figure 5.16: Comparison of shearing wavecurrent interaction (angle).
0.3
0.2
69
Table 5.3
Case HMI HM II EM I EM II PM
Governing equation M O O O M
Programming ease O M O M O
Numerical stability M X X X O
Computational time X X O M O
Shoaling O O O O O
Refraction M O M O M
Diffraction (normal) O O O O O
Diffraction (oblique) O O X O M
Reflection (vertical) O O 
Reflection (slope) O O 
Current collinearr) M 0 0 0 0
Current (refraction) X O M O X
0: good M: marginal X: bad : not applicable
CHAPTER 6
SURF ZONE MODEL
The flow properties in surf zone are utmost complex owing to the strong in
teractions among motions induced by waves, currents, and turbulence. The present
knowledge on surf zone dynmaics is still inadequate and most of the models are rather
rudimentary. Most numerous developments in the study of wave breaking have been
made by approximation of the wave energy dissipation. These models can be classi
fied into two groups: one is based on the similarity of the wave breaking process with
other hydraulic phenomena such as a hydraulic jump (Dally et al., 1984), a tidal bore
(Battjes and Janssen, 1978), etc., and the other is based on estimation of the eddy
viscosity (Mizuguchi, 1980) or turbulence (Izumiya and Horikawa, 1984).
In this chapter a simple surf zone model is presented. This model is based on the
consideration of wave energy balance and wave action conservation so that the wave
current interaction is fully taken into account. The model is capable of predicting
wave decay and yields quasithree dimensional mean current profiles both in the
acrossshore and longshore directions. The model is presented in analytical form for
the case of two dimensional graduallysloped bottoms. Numerical scheme for general
topographic applications including the wave diffraction effects is presented in latter
chapters in conjunction with the nearshore circulation model.
Section 6.1 modifies the wave energy equation given in Chapter 3 for surf zone
application with the addition of an energy dissipation term. In Section 6.2, the
validity of wave action equation in surf zone is established. Subsequent sections are
devoted to solve decays of wave heights, mean surface currents, and mean setups in
surf zone. Whenever possible, the results are compared with available experimental
data.
6.1 TimeAveraged Wave Energy Equation in Surf Zone
It is assumed here that the surf zone is coherent in that the essential wavelike
periodic motion is retained and is quasistationary when timeaveraged over wave
period. Furthermore, the turbulent motions are of much smaller time scale that it
effect on the flow of interest can be simply treated as a dissipative force menifested
by an eddy viscosity term. Thus, the familiar Navier Stokes' equation of the following
form can be applied:
U qv ( + g+ = (Vx U) x U+ V2U (6.1)
at 2 p
where / is the total viscosity coefficient including the eddy viscosity due to turbulence.
Let the surface displacement and the velocity vector, U(u, v, w), be decomposed into
mean value, wave and turbulent fluctuations, which are distinguished by subscripts c
and w, and prime, respectively; thus,
7 = j+ ''= + + 7' (6.2)
U = U + U'=Uc + U,+U' (6.3)
where the superscript'is used to denote turbulent averaging. After turbulentaveraging
Eq. (6.1) becomes
+ V ( + + g = (V x U) x U+ VV2 (6.4)
it 2 p
The superscript is omitted hereafter.
Taking the scalar product of U(u, v, w) with the respective terms in the Navier
Stokes equation and summing the products give the energy equation with dissipation
ah t 2 2 p 
Jh{[2]+V. [U( 2+ + gz)] dz = vtU V2Udz (6.5)
by applying the kinematic boundary conditions at the free surface and the bottom,
wl7 U.Vh7 = 0
wlh+U.Vhh = 0
72
and letting p equal to zero at the surface, we obtain,
[ dz + g + v [U( + + gz)]dz + tU V2Udz = 0 (6.6)
Jh 2 at + J+h p fh
again, utilizing the Leibnitz' rule of integration. This equation is similar to the
wave energy equation given by Eq. (3.26) with the additional disspative term. This
disspative term is in the form of a energy flux and can be treated as a head loss term
in the context of Bernoulli equation, i.e.,
) = h VU V2Udz = Va J" glUdz (6.7)
h h
here I is the head loss due to turbulence. Equation (6.6) can then be expressed as
a [ 7 ]dz + g + V 2 [U( + +gz + gl)dz= (6.8)
Following the same procedures in Chapter 3 by taking time average over wave
period, an energy equation similar to that of Eq. (3.38) can be obtained,
aE
S+ Vh (UhE) = 0 (6.9)
at
Here the transport velocity Uh is the counter part of (Cg + i dclk/a) in Eq.
(3.38). Here, the last term can be dropped because of the quasisteady assumption
made here. Clearly, this transport velocity is different from the nondissipative case
and can be represented by
Uh = Cg + + CgD (6.10)
where the first two terms constitute the transport velocity due to nondissipative
forces much the same as in Eq. (3.38) whereas the last term manifests the effect due
to dissipative force. This term, in general, should be negative indicating a reduced
energy flux due to dissipation. In theory, it can be estimated from the timeaveraged
energy dissipation term can be estimated following the definition,
Vh (CgDE) = D
(6.11)
73
where D is the timeaveraged dissipation given by
D UtU V2Udz = Vh glUdz
h h
6.2 Wave Action Equation in Surf Zone
In this section, the wave action equation given in Eq. (3.41) is shown to be also
valid even in the presence of strong turbulence. It is assumed here that the dynamic
free surface boundary condition given in Eq. (3.10) is still valid with the inclusion
of a head loss term. Following the approach by Kirby (1983), a virtual work term
is proportional to WD is introduced to represent the head loss, where W is an
positive undefined coefficient indicating the strength of the dissipation. We examine
this approach. The dynamic free surface boundary condition then becomes
(1 + W)
1 W) [(P)t + Vh" . VhW] (6.12)
9
which is divided into two secular components that should be separately satisfied,
a
A = g a(6.13)
OA
+ 8,. VA = 0 (6.14)
where OD = (1 + W)(w U, k). The subscript, s, denotes the mean water surface
level. The kinematic free surface boundary condition also yields
0= = (1 + W)gk tanh k(h + ic) (6.15)
aa
a + V (O,a) = 0 (6.16)
Combining Eq. (6.14) and Eq. (6.16), we obtain the following wave action equation
which is valid in the surf zone:
a pgH2 Pg H2
(pg + V g ( H~~D ) = 0 (6.17)
( )+o ( UW )= O (6.17)
9t 8 OD 8 ao
This equation enables us to estimate the surface velocity in the surf zone once the
wave height decay rate can be established.
74
6.3 Wave Height Transformation in Surf Zone
Waves break when their height reaches a certain limiting value relative to their
length or water depth as a result of wave shoaling on a slope. The broken waves
normally keep breaking as the water depth decreases, finally reaching the shoreline.
Svensen et al. (1978) divided the breaking zone into inner and outer regions: From
the breaking point and for some distance shoreward, it is the outer region where a
violent transition of the wave slope takes place large scale vortices are formed in this
region. After outer region, the inner region begins as the wave becomes very similar
to a tidal bore or a hydraulic jump. However, this wave breaking process not yet
been fully clarified since the strong currents and turbulence are generated by the
broken waves, and interacted with wave breaking. In this section, we suggest the
new approach which takes account of the wavecurrent interaction in the surf zone.
Differently from most of the existing wave breaking models, this new model provides
the analytical expression of wave height over the wave breaking zone.
The wave energy equation given in Eq. (6.9), when expressed in terms of wave
height, can be written as
O w H2 w H2
(pg) + Vh [(Cg + fI + CgD)pg ] = 0 (6.18)
,t 0D 8 oD 8
Now, again we assume that the surf zone retains a quasisteady state when integrated
over wave period, then the slowly varying flow properties become time independent,
and the absolute frequency becomes a constant. Accordingly, Eq. (6.18) becomes
Vh [(Cg + U + CgD) ] = 0 (6.19)
8 oD
and applying the wave action equation in the steady state, the wave energy equation
in the surf zone is reduced to
Pg Hz2
Vh [(Cg + CgD) ] = 0 (6.20)
8 UD
75
The quantity corresponding to the dissipative force is assumed here to be proportional
to group velocity at the breaking point, Cgb, since wave height is known to decrease
steadily within the surf zone; therefore,
CgD = /3Cg (6.21)
where / is a positive coefficient. Equation (6.20) now becomes
Vh (Cg*DH2rD) = 0 (6.22)
8
where the real relative group velocity in the turbulent surf zone is estimated by
Cg* = Cg Cgb (6.23)
Equation (6.22) is the final form of the proposed energy transformation model. This
model has only unknown coefficient, namely, the dissipation coefficient, P and is
applicable to the general threedimensional topography and any arbitraly incident
wave angles.
An analytical expression can be obtained for twodimensional beaches of uniform
slopes. Equation (6.22) becomes
S(Cg Cgb) ]= 0 (6.24)
ax UD
where x axis is directed onshore. The cross shore component of the above equation
gives,
H2 C aD (6.25)
cos 0(/CgS Cg)
Applying the dynamic free surface boundary condition given in Eq. (6.14) in 2D
steady state condition, we have 'D proportional to the wave height in surf zone.
Equation (6.25) can be written as
H = )H (6.26)
cos 0(PCgb C9)
76
with P/3 determined later. 0 can be determined by Snell's law as
0 = sinl(C, sin 00/C0) (6.27)
where C, = w/k is an absolute phase speed. Equation (6.26) is nondimesionized as
H 1 (6.28)
Hb HbCgb cos 0(3 Cg')
where Cg' = Cg/Cgb. Now we apply two boundary conditions; H/Hb = 1 at a
breaking point d/db = 1, and H/Hb = H' where Cg becomes zero. We then obtain
H ____(6.29)
Hb cos 0[1 Cg'(1 H'/ cos Ob)]
Applying shallow approximation to Cg' gives
H 1 H'
H I H (6.30)
Hb cos 0 1 (1 H'cos Ob)
where d' = d/db and d is defined as a total water depth, (h + 0c). Of the three
parameters H', 1fH and / only one is independent If the value of H' is determined by
experiment the other two are solved by,
cos ObH'
S= bH, HbCgb (6.31)
cos 0b H,
cos Ob
P = (6.32)
cos Ob HI
Figure 6.1 plots the dimensionless wave height in the surf zone for different 3 val
ues. Figure 6.2 shows the comparison between the present theory and the laboratory
data by Horikawa and Kuo (1966). The dimensionless values of the wave height at
the shoreline H', are determined from the data; they are 0.22, 0.18, 0.14 and 0.14,
respectively, for slopes of 1/20, 1/30, 1/65 and 1/80. It was found that the experi
mental values of H' can be closely approximated by Vt/ana with a being the slope
of the beach. Extensive laboratory experiments have shown that the pattern of wave
height decay across the surf zone is strongly a function of the beach slope.
77
Wave Breaking
0.6 0.8
h/hb
Figure 6.1: Effect of a parameter 3 on the dimensionless wave height in the surf zone.
1.0
0.8
0.6
0.4
Slope 1/20
0.0 0.2 0.4 0.6 0.8 1.0
d/db
Slope 1/65
0.0 0.2 0.4 0.6 0.8 1.0
d/db
0.0 0.2 0.4 0.6 0.8 1.0
d/db
Slope 1/80
0.0 0.2 0.4 0.6 0.8 1.0
d/db
Figure 6.2: Comparison with laboratory experiments presented by Horikawa and Kuo
(1966).
I
Slope 1/30
Hyperbolic Model
Elliptic Model I
1.0
0.8
0.6
S 0.4
0.2
0.0
0.0 0.2 0.4 0.8 0.8 1.0
d/db
Elliptic Model II
0.0 0.2 0.4 0.8 0.8 1.0
d/db
Parabolic Model
1.0
0.8
0.6
a 0.4
0.2
0.0
0.0 0.2 0.4 0.6 0.8 1.0
d/db
0.0 0.2 0.4 0.6 0.8 1.0
d/db
Figure 6.3: Comparison between theoretical and numerical results.
1.0
0.8
0.8
:3 0.4
0.2
0.0
Pz
I
80
In the full 3D model, Eq. (6.22) has to be employed. It is solved numerically
either implicitly or explicitly depending upon whether one treats CgfCgb as a group
velocity (implicit) or solves D independently by the explicit expression of Eq. (6.10).
It is suggested that the implicit scheme be used for uniform slope cases and the explicit
scheme for more complicated topographies.
Figure 6.3 compares the numerical results with the theoretical curves for the
case with P=1.28 (or H' = 0.22) corresponding to a uniform slope of 1/20. Since the
dimensionless wave weights from the numerical compuatation are expressed as
H 2 1 H'2
( = 1 (6.33)
Hb cos 0 1 v/(1 H'2/ cos Ob)
and from the theoretical model based on Eq. (6.30) the following relation has to be
used to compare them,
(,2 vX Hb
Hb \Hb x
where the subsripts m and t denote numerical and theoretical values, respectively.
6.4 Surface Currents in Surf Zone
The main mechanism responsible for current generation inside the surf zone is
suggested to be due to the excess waveinduced momentum also known as the radia
tion stresses. Vertical circulations are known to exist as a consequence mass balance
to maintain a quasisteady state. Solution conerning longshore current and its dis
tribution was originally obtained by LongguetHiggins (1970) and later modified and
refined by numerous other investigators. The model was based on the balance be
tween the friction forces and gradients in the radiation stress. Crossshore current
modeling is more recent effort based on such ideas as wave setup, undertow, etc. In
this section, the surface current vectors in the surf zone containing both cross shore
and longshore components are solved by the applications of wave action equation and
the steady state wave energy equation:
_H2
V (U,) = 0 (6.34)
aD
I
81
K H2
V [(Cgb Cg) ] =0 (6.35)
Elimilating from the above equations we obtain a pair of simple equations for
surface currents,
u, = Po cos 0(3Cgb Cg) (6.36)
v. = fL sin 0(3Cgb Cg) (6.37)
where u, and v, denote, respectively, the crossshore and longshore components of
surface current vector at the mean water level, U,; f/o and 3P are constants of pro
portionarity. When expressed in terms of wave height the above pair of equations
become,
u, = floH (6.38)
sin 0 aD
v, = AH 0 (6.39)
cos 0 H2
which shows that while both onshore and longshore current components are inversely
proportional to the wave height sqaured, only longshore component is a function of
wave angle. The surface current equations can be more conveniently expressed in
nondimensional forms,
u/ = o cos 0(3 Cg') (6.40)
and
v, = fLsin0o( Cg') (6.41)
where u' = u,/Cgb and Cg' = . Applying the shallow water condition, we have,
3 Cgb
u', = I3ocos0(p Vh)
I I
(6.42)
Surface
Onshore
Current
0.4 0.6
h/hb
Figure 6.4: Effect of a parameter 3 on the dimensionless surface onshore current in
the surf zone.
\
BETA= 1.3
BETA=1.25 ..
\\
m\.
*.. \%
T\\ 13 .
BETA= 1.2
 BETA= 1.15 ........
... BETA= 1.
 BETA=1.1I
83
Surface Longshore Current
0.0 0.2
0.4 0.6
h/hb
Figure 6.5: Effect of a parameter P on the dimensionless surface longshore current in
the surf zone.
Slope 1/10
0.2 0.4 0.6 0.8
Slope 1/20
0.2 0.4 0.6 0.8
d/db
Figure 6.6: Comparison of longshore current with laboratory experiments presented
by Visser (1991).
85
v L sin O0 Vh( j) (6.43)
Figure 6.4 illustrates the effect of / on the dimensionless surface onshore current
given by Eq. (6.42). Unfortunately, no experimental data are available at present.
Figure 6.5 illustrates the effect of 3 on the dimensionless surface longshore current
given by Eq. (6.43). Figure 6.6 compares the theory with the laboratory longshore
current data measured by Visser (1991). It should be pointed out here that the
Visser's data are depthaveraged and theory is for surface current. In the longshore
direction, however, one expects the vertical distribution to be rather uniform as will
be shown in Chpater 7. The theory appears to fit the data remarkably well. One
of the major advantages of the present model is that it requires only one empirical
coefficient to control the magnitude and eliminates the troublesome mixing coefficient
appeared in most of the existing theories. In order to fit the data one often has to
assume large mixing strength without justification.
6.5 SetUp and SetDown
The wave setup in the surf zone is solved here by the equation of 'the kinematic
conservation of intrinsic frequency' and the steady state continuity,
V. (C,H) = 0
V [C(h + r77)] = 0
It is assumed here that the magnitude of the depthaveraged return current beneath
the mean water level, U is proportional to the onshore surface current, U,, then
H = c(x)(h + +c) (6.44)
where Ki(x) is, in general, a spatial dependent coefficient. If r.(x) is a constant across
the surf zone, the above equation is simply the frequently adopted extension of Miche's
criterion. From Eq. (6.44) the set up is solved as
H
77 () h (6.45)
tc(x)
or in nondimensional form,
S K(x=Xb)H' h' (6.46)
where r' = rc/hb and za indicates the breaking point. Here the definition of the mean
water level is limited to the where the bed is at all times covered by water. The setup
or setdown can be computed at once when the wave height within the surf zone is
determined. Figure 6.7 illustrates the effect of # on the dimensionless setup given
by Eq. (6.46) for a constant I over the surf zone.
There is one notable feature of the model that is generally lacking in most of
the existing models. The model predicts rather mild, sometimes near constant, set
down in the transition zone immediately after breaking point where the wave height
drops sharply. This phenomenon referred to as transition region 'paradox' has been
noted as a significant feature in the transition region (Basco and Yamashita, 1986;
Thieke,1988). The conventional theory of balancing the momentum due to radiation
stress should produce a jump of setup in the transition region where wave height
reduces sharply. Laboratory data, on the other hand, showed nearly constant setdown
across the transition region where the wave height is reduced nearly proportional to
the reduction of water depth. This phenomenon will be further clarified in Section
7.2.3 through the momentum balance.
87
SetUp
0.0
0.1
0.0 0.2
0.4 0.6
h/hb
Figure 6.7: Effect of a parameter 3 on the dimensionless setup in the surf zone.
CHAPTER 7
MATHEMATICAL MODEL FOR WAVEINDUCED CURRENTS
A prominent feature in the nearshore zone is the wave induced current circulation.
It is commonly accepted that the primary driving force is the wave induced radia
tion stress first introduced by LonguetHiggins and Stewart (1961). Modelling this
circulation has since advanced considerably from the earlier development by Noda et
al. (1974) and Ebersole and Dalrymple (1979). Both of these earlier models were
driven by a wave refraction model with no current feedback. In recent years, Yoo and
O'Connor (1986b) developed a coupled waveinduced circulation model based upon
what could be classified as a hyperbolic type wave equation; Yan (1987) and Winer
(1988) developed their interaction models based upon parabolic approximation of the
wave equation. All these models employed the depthaveraged formulations which
have three major deficiencies 1) the surface effects due to wavecurrent interaction,
which generally is very strong, are being neglected; 2) the bottom friction is expressed
in terms of the mean velocity which makes the model unrealistic in areas where the
current profile is strongly three dimensional but the mean current could be small
such as in the surf zone and 3) the convective acceleration terms are also depth
averaged which has the same problem as (2). Recently, de Vriend and Stive (1987)
improved the nearshore circulation model by employing a quasithree dimensional
technique. This technique is very attractive to accommodate the surf zone in which
the depthaveraged model is no longer valid but the full threedimensional modelling
is currently not attainable. In this chapter, this approach of quasithree dimensional
modelling is adopted by developing a circulation model combining depthintegrated
properties with vertical profiles.
89
In Section 7.1, the fundamental conservation equations of mass and momentum
timeaveraged over turbulent scale are presented. In Section 7.2, the depthintegrated
formulations serving as the basic equations for a quasithree dimensional circulation
model are derived. An amended form of radiation stress is presented and the exis
tence of the surface advection terms is verified through comparisons with wave energy
equation. Section 7.3 develops a new model prescribing turbulenceinduced vertical
current distributions in surf zone. This model employing the surface current bound
ary conditions given in Chapter 6 yields circulation patterns in both crossshore and
longshore directions. Finally, in Section 7.4, a quasithree dimensional model suitable
for the entire nearshore zone is developed by linking the depthintegrated model to
the surf zone model.
7.1 TurbulenceAveraged Governing Equations
The strong presence of turbulence is a prominent feature in surf zone. Con
sequently, the fundamental equations governing the fluid motion should also in
clude the turbulent effects. This is usually accomplished with the introduction of
Reynolds stresses by time averaging over the turbulent fluctuations. Accordingly,
the turbulenceaveraged governing equations are presented here. Again, the basic
equations are: the continuity equation for incompressible fluid,
9u av aw
u+ + =o (7.1)
xz 9y Oz
the horizontal momentum equations,
Ou Ouu Ouv Ouw 1 p 1 (Or, 9ry Or.
+ + + =  + ( + +  ) (7.2)
at O 9 y 9z pax p Ox 9y az
Ov Ovu Ovv Ovw 1 Op 1 aOr Or Or a
S+ + j+ 5z + ~t y ) (7.3)
t xz 9y z p ay p ax Dy Qz
and the vertical momentum equation,
Ow Owu Owv Oww 1 O(p + pgz) 1 Orz Orz Or
T+ + + + ( + Oy +  (7.4)
Ot 9x 9y Pz p Oz p ax ay Oz
90
Let velocity vector, U(u, v, w), be decomposed into mean velocity, wave and turbulent
fluctuations, which will be distinguished by bar, hat and prime, respectively; thus,
U= + U'= U + U" + U' (7.5)
where the superscript is used to denote ensembleaveraging. The characteristic time
scales of the three components are assumed to be quite different, therefore no cor
relation between them is considered. In general, the time scale of the turbulent
fluctuations is considerably shorter than the wave period. Turbulentaveraging Eqs.
(7.17.4) results in
+ + = 0 (7.6)
Ox Oy 9z
+ + + + ( + + (7.7)
at ox ay 8z p dz p az 9y az
9v 9 9 a9w o 1_ l1a9f, a9 a9f.
S+ + +  + ( + + ) (7.8)
9t Ox ay 9z p y p ax ay 9z
O Owii Owi Odwt 1 O(P + pgz) 1 fzz +fyz Of
+ + +  + + + ) (7.9)
ot + x+y z p 9z p Ox ay Oz (9
where f includes the Reynolds stresses due to turbulence flow.
7.2 Horizontal Circulation Model
The governing equations for the horizontal circulation model are obtained after
depth integration and waveaveraging. In order to protect from losing the Eulerian
mean quantities at the mean water level, the depth integration is taken prior to
waveaveraging them.
7.2.1 DepthIntegrated and TimeAveraged Equation of Mass
Integrating Eq. (7.1) over depth, hereafter omitting the tildes, we get
aO 9u, 7 O v O a f_ 9a
S+ dz + dz= udz + vdz
t OJh y 0 t O h O y
S 9 (y 7, h h, 0)
+[ + u a+ v w], + [tU + v + w]h = 0 (7.10)
at 8x ay ax ay
