• TABLE OF CONTENTS
HIDE
 Title Page
 Acknowledgement
 Table of Contents
 List of Figures
 Abstract
 Introduction
 Description of nearshore hydrodynamic...
 Energy propagation through varying...
 Mathematical wave models
 Numerical scheme and result of...
 Surf zone model
 Mathematical model for wave-induced...
 Numerical scheme and result of...
 Conclusion and further study
 Appendix A: Wave-averaged energy...
 Appendix B: Wave-averaged energy...
 Appendix C: Green's identities...
 Appendix D: Numerical difference...






Group Title: Technical report – University of Florida. Coastal and Oceanographic Engineering Program ; 95
Title: Wave-current interaction and quasi-three-dimensional modeling in nearshore zone
CITATION PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00075322/00001
 Material Information
Title: Wave-current interaction and quasi-three-dimensional modeling in nearshore zone
Physical Description: x, 158 leaves : ill. ; 29 cm.
Language: English
Creator: Lee, Jung Lyul, 1959-
Publication Date: 1993
 Subjects
Subject: Coastal and Oceanographic Engineering thesis Ph. D
Dissertations, Academic -- Coastal and Oceanographic Engineering -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
 Notes
Thesis: Thesis (Ph. D.)--University of Florida, 1993.
Bibliography: Includes bibliographical references (leaves 153-157).
Statement of Responsibility: by Jung Lyul Lee.
General Note: Typescript.
General Note: Vita.
Funding: Technical report (University of Florida. Coastal and Oceanographic Engineering Dept.) ;
 Record Information
Bibliographic ID: UF00075322
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved, Board of Trustees of the University of Florida
Resource Identifier: aleph - 001934048
oclc - 30811712
notis - AKB0130

Table of Contents
    Title Page
        Title Page
    Acknowledgement
        Acknowledgement
    Table of Contents
        Table of Contents 1
        Table of Contents 2
        Table of Contents 3
    List of Figures
        List of Figures 1
        List of Figures 2
        List of Figures 3
    Abstract
        Abstract 1
        Abstract 2
    Introduction
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
    Description of nearshore hydrodynamic model
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
    Energy propagation through varying current field
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
    Mathematical wave models
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
    Numerical scheme and result of circulation models
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
    Surf zone model
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
    Mathematical model for wave-induced currents
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
    Numerical scheme and result of circulation model
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
    Conclusion and further study
        Page 144
        Page 145
        Page 146
    Appendix A: Wave-averaged energy density
        Page 147
        Page 148
    Appendix B: Wave-averaged energy flux
        Page 149
        Page 150
    Appendix C: Green's identities consistent with energy equation
        Page 151
    Appendix D: Numerical difference operators
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
Full Text












WAVE-CURRENT INTERACTION AND
QUASI-THREE-DIMENSIONAL 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. Li-Hwa 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 Wave-Current 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 Quasi-3D 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 Time-Averaged Wave Energy .................. 26
3.2.2 Time-Averaged Wave Energy Flux ............... 27
3.2.3 Time-Averaged Wave Energy Equation ............. 28
3.3 Wave Action Equation .......................... 28
3.4 Excitation Due to Wave-Current 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 Hyperbolic-Type Model I (HM I) ............. ... 36
4.2.2 Hyperbolic-Type Model II (HM II) ............... 38
4.2.3 Elliptic-Type Model I (EM I) .................. 39
4.2.4 Elliptic-Type Model II (EM II) ................. 40
4.2.5 Parabolic-Type 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 Wave-Current Interaction . . . . ... 62
5.2.6 Summary ............................. 62

6 SURF ZONE MODEL ........................... .. 70

6.1 Time-Averaged 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 Set-Up and Set-Down ........................... 85

7 MATHEMATICAL MODEL FOR WAVE-INDUCED CURRENTS .... 88

7.1 Turbulence-Averaged Governing Equations . . . .. 89
7.2 Horizontal Circulation Model .. ............... .. 90
7.2.1 Depth-Integrated and Time-Averaged Equation of Mass .. 90
7.2.2 Depth-Integrated and Time-Averaged 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 Three-Dimensional Topography 112
7.4 Quasi-3D M odel .............................. 116
7.4.1 Modification of 2D Depth-Integrated 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 Quasi-3D Circulation Model . . ... 122
8.1.1 Fractional Step Method . . . . ... 123
8.1.2 Finite-Volume 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 WAVE-AVERAGED ENERGY DENSITY ................. 147

B WAVE-AVERAGED 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 Sub-grid 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 semi-infinite 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 semi-infinite 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 wave-current interaction. . ... 65

5.13 Condition of shearing wave-current interaction. . ... 65

5.14 Comparison of collinear wave-current interaction. . ... 66

5.15 Comparison of shearing wave-current interaction (height). . 67









5.16 Comparison of shearing wave-current 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 set-up in the surf zone. 87

7.1 Vertical profile of cross-shore currents measured by Hwung and
Lin (1990) ............ ....... .. . . 99

7.2 Vertical profiles of cross-shore 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 three-dimensional profiles. . . . ... 115

7.8 Vertical profiles of cross-shore 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 set-up contours represented by Gourlay (1974). 136

8.5 Measured current pattern. . . .... . . 137









8.6 Comparisons of wave height contours. . . . ... 138

8.7 Comparisons of set-up contours. . . . ... 139

8.8 Comparisons of depth-mean current pattern resulted from quasi-
3D m odel ............................ 140

8.9 Comparisons of depth-mean 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

WAVE-CURRENT INTERACTION AND
QUASI-THREE-DIMENSIONAL 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 wave-current interaction, surf-zone hydro-

dynamics, and wave-induced three dimensional current pattern. Based on this study,

a quasi-three 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 wave-current 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 two-dimensional gradually-sloped bottoms. This

surf zone model provides the surface current pattern of the vertical circulation model,

and consequently, significantly contributes to solving the three-dimensional 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

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

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

set-up and the wave-induced current such as the longshore, undertow, etc.

To solve the wave-current 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 three-dimensional 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 three-dimensional 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 two-dimensional profiles. Most of the proposed surf zone models whether the-

oritical or numerical were also restricted to two-dimensional. The model developed

here allows for applications to more general three-dimensional topographies.


I







3

The surf-zone 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 cross-shore component in the surf zone contains a shore-

ward mass transport and an offshore-directed 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 three-dimensional 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 Wave-Current Interaction Model

Wave-current interaction problems have interested a considerable number of math-

ematicians beginning apparently with Longuet-Higgins and Stewart (1960-1961), and

Whitham (1962). Longuet-Higgins 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 wave-current

interactions. The problems of wave-current 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 time-dependent, depth-integrated 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 refraction-diffraction effects of wave-current 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 wave-current 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 wave-current 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 cross-shore 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 cross-shore surf zone modeling, the usual elements include the fixation of

breaking criterion, the energy decay rate across the surf zone, the wave set-up and

finally the cross-shore 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 two-parameter

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 cross-shore current pattern has been studied during the

past ten years for the simple two-dimensional 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.

Wave-induced 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 Longuet-Higgins (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 wave-induced 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 wave-induced 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 vertically-averaged longshore current. Recently, however, de Vriend

and Stive (1987) improved the nearshore circulation model by a quasi-three dimen-

sional approach, which employed a combined depth-integrated 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

cross-shore and longshore direction. Svensen and Lorenz (1989) assumed that the

equations in the cross-shore and longshore directions could be decoupled and solved

the cross-shore 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 wave-current 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 refraction-diffraction. 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 cross-shore and longshore circulation model is established

by solving the time averaged mass and momentum equations. Finally, a quasi-3D

model is constructed by combining the horizontal circulation model and the vertical

circulation model.

In Chapter 8, The numerical scheme solving the quasi-3D 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 waves-induced 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 sub-model 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 surf-zone 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 wave-induced

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

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

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 quasi-three 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+ (Cg-K + 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 method-FDM

EM II a, K and 0 Generalized Lax-Friedrich FDM

PM complex A' Crank-Nicholson FDM

The wave breaking model nested in to each wave model can be applied by replacing
the Cg-PCgb 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 cross-shore 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, turbulent-induced 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)
S-24(1 7)/L (P Cg/Cgb)
P V. (KH/k)
F C
4rPL (P Cg/Cgb)








where I|orbl is defined as gH/2C.
2.2.2 Quasi-3D Model

The quasi-3D depth-varying 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 three-dimensional model.
The continuity equation
c Q gH(2k 8L gH2k
7+ + (Q + )+a -(Q,+ 8 )=
Wt TX 8o Ty 8o
The x-directional 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 y-directional 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 Longuet-Higgins (1970).

E, = NlxI\Vgd


It was suggested that a dimensionless constant N, should be less than 0.016. The
y-directional 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 time-dependent, three-dimensional incompressible Navier-Stokes 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 finite-volume method allows one to handle arbitrary grids

and helps to avoid problems with metric singularities and non-conservative nature

that are usually associated with finite-difference 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

wave-current 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 depth-integrated wave energy equation is derived and the
exact forms of wave energy and wave energy flux are presented in wave-averaged 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
wave-current interaction models both inside and outside surf-zones 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 3-dimensional 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 z-axis, 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,

wave-induced 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) + Enon-propagating 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 wave-averaged, 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 depth-integrated







25
energy equation has been derived from Eq. (3.22) by many including Longuet-Higgins
and Stewart (1961) Witham (1962). A brief account is given here.
Using Leibnitz's rule (f2h D fdz = D f h fdz fl,Dr + fl-hDh), Eq. (3.22) can
be written as

7 +[ 2 j]dz + V. 7 [UQ + p + pgz)]dz
-h 2, J-h[ 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
WI-h + 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 J-h[ ^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 well-known energy equation given by Longuet-Higgins 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 Time-Averaged 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 wave-current
interaction. Longuet-Higgins 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 Time-Averaged 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 pe-77[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 Time-Averaged Wave Energy Equation

By substituting Eq. (3.37) into Eq. (3.29) the time-averaged 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 Wave-Current Interaction

Three basic equations governing wave-current 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+ = c-k) =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 no-current 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 (Longuet-Higgins 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. (Cg-ck)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 x-direction,

(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 long-term 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= +U-k


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 wave-averaged quantities, and the parabolic type is expressed in terms of
wave-averaged quantities for the main propagating axis, and wave mode for the the
other axis.

The extension of the mild slope equation to the wave-current 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 wave-current
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 wave-induced 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.3-3.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
J-H[ ]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 long-term 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 f-h 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 semi-steady 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 Hyperbolic-Type Model I (HM I)

The first governing equation of the hyperbolic-type is given as a pair of first-order

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 time-dependent 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.13-15) 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
0-at aJ c
This inconsistency is caused by Eq. (3.15) which satisfies a slowly varying property
of the dynamic motion.









4.2.2 Hyperbolic-Type 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 wave-period

and wave-length 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 wave-averaged

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 'Elliptic-type 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+U-K= 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
time-dependent dynamic conservation equation is generally given by

Oa o a2
+ 2V [(Cg + ,)a] =0 (4.22)
0t 2a a

4.2.3 Elliptic-Type Model I (EM I)

The governing equation for the first elliptic-type 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 time-harmonic.









4.2.4 Elliptic-Type Model II (EM II)

If the two-dimensional velocity potential is assumed through use of the wave-

period and wave-length 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-) + CCg-V K = 0 (4.25)
U9dd Od

Imaginary part

CCga-K' V (CCgV ) oda + (a2 k2CCg)a- = 0 (4.26)
O-d 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 Parabolic-Type Model (PM)

The parabolic approximation to the elliptic-type 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 x-direction 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(u-y) -
Ox oy

+-[(CCg-v2) + 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 pseudo-operators 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 (CCg-u2)


(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 y-directional part of the phase function, e f kydy, and the x-directional
phase modulation due to the difference between the local and averaged wave num-
bers, ei(k-fk'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.13-4.15) are solved explicitly for

the x-y 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) = t-r(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 x-axis.

The non-reflective 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 wave-current 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 y-directional 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 first-order 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: Sub-grid 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(xi-l, i-1)

Yj+i = yj-1 + 2hf(xi-_ + jh, yj)

i (yn + Yn-1 + 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
(NY-1,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 Runge-Kutta 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 y-direction, 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 = tan-1 ( ) (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 Lax-Friedrich (L-F) 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 L-F 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 parabolic-type wave equation (4.42) is written in the finite difference form

using the Crank-Nicolson 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 tri-diagonal 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 non-reflective 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 x-directional 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 = tan-1 K) (5.27)


where

K,, k 1 (K/k)2, K, = VDS (5.28)







53

where S = f Kdx kox + Ky = tan-r[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, wave-current 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 cross-sections 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 VAX-8350 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.5-5.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 semi-infinite 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 semi-infinite

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'













--3-2 -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 semi-infinite 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 semi-infinite 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 Wave-Current Interaction

Wave-current 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.14-16. For the collinear case, all except hyperbolic model I performed adequately.

For non-collinear 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 wave-current interaction.







x












x I -

yo
waves i


Figure 5.13: Condition of shearing wave-current 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 wave-current 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 wave-current 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 wave-current 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 quasi-three dimensional mean current profiles both in the

across-shore and longshore directions. The model is presented in analytical form for

the case of two dimensional gradually-sloped 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 set-ups in

surf zone. Whenever possible, the results are compared with available experimental








data.
6.1 Time-Averaged Wave Energy Equation in Surf Zone

It is assumed here that the surf zone is coherent in that the essential wave-like
periodic motion is retained and is quasi-stationary when time-averaged 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 turbulent-averaging
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

a-h 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,

wl-7 -U.Vh7 = 0
wl-h+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)
J-h 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 quasi-steady assumption
made here. Clearly, this transport velocity is different from the non-dissipative case
and can be represented by

Uh = Cg + + CgD (6.10)

where the first two terms constitute the transport velocity due to non-dissipative
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 time-averaged
energy dissipation term can be estimated following the definition,


Vh -(CgDE) = D


(6.11)







73
where D is the time-averaged 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 wave-current 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 quasi-steady 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 three-dimensional topography and any arbitraly incident

wave angles.
An analytical expression can be obtained for two-dimensional 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 2-D
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 = sin-l(C, sin 00/C0) (6.27)


where C, = w/k is an absolute phase speed. Equation (6.26) is non-dimesionized 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 3-D model, Eq. (6.22) has to be employed. It is solved numerically

either implicitly or explicitly depending upon whether one treats Cg-fCgb 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 wave-induced momentum also known as the radia-

tion stresses. Vertical circulations are known to exist as a consequence mass balance

to maintain a quasi-steady state. Solution conerning longshore current and its dis-

tribution was originally obtained by Longguet-Higgins (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. Cross-shore current

modeling is more recent effort based on such ideas as wave set-up, 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 cross-shore 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

non-dimensional 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 depth-averaged 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 Set-Up and Set-Down

The wave set-up 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 depth-averaged 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 non-dimensional 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 set-up

or set-down 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 set-up 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 set-up in the transition region where wave height

reduces sharply. Laboratory data, on the other hand, showed nearly constant set-down

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











Set-Up


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 set-up in the surf zone.














CHAPTER 7
MATHEMATICAL MODEL FOR WAVE-INDUCED 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 Longuet-Higgins 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 wave-induced 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 depth-averaged formulations which

have three major deficiencies 1) the surface effects due to wave-current 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 quasi-three dimensional

technique. This technique is very attractive to accommodate the surf zone in which

the depth-averaged model is no longer valid but the full three-dimensional modelling

is currently not attainable. In this chapter, this approach of quasi-three dimensional

modelling is adopted by developing a circulation model combining depth-integrated
properties with vertical profiles.







89
In Section 7.1, the fundamental conservation equations of mass and momentum

time-averaged over turbulent scale are presented. In Section 7.2, the depth-integrated

formulations serving as the basic equations for a quasi-three 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 turbulence-induced vertical

current distributions in surf zone. This model employing the surface current bound-

ary conditions given in Chapter 6 yields circulation patterns in both cross-shore and

longshore directions. Finally, in Section 7.4, a quasi-three dimensional model suitable

for the entire nearshore zone is developed by linking the depth-integrated model to
the surf zone model.

7.1 Turbulence-Averaged 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 turbulence-averaged 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+ + --+ + ( + O-y + -- (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 ensemble-averaging. 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. Turbulent-averaging Eqs.
(7.1-7.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 wave-averaging. In order to protect from losing the Eulerian

mean quantities at the mean water level, the depth integration is taken prior to
wave-averaging them.
7.2.1 Depth-Integrated and Time-Averaged 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 OJ-h 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




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs