Numerical modeling of current and wave interactions of an inlet-beach system

Material Information

Numerical modeling of current and wave interactions of an inlet-beach system
Series Title:
Yan, Yixin, 1949-
University of Florida -- Coastal and Oceanographic Engineering Dept
Place of Publication:
Gainesville Fla
Coastal and Oceanographic Engineering Dept.
Publication Date:
Physical Description:
xi, 163 p. : ill. ; 28 cm.


Subjects / Keywords:
Water waves -- Mathematical models ( lcsh )
Hydrodynamics -- Mathematical models ( lcsh )
Inlets ( lcsh )
Sediment transport ( lcsh )
Engineering Sciences thesis Ph. D
Dissertations, Academic -- Engineering Sciences -- UF
bibliography ( marcgt )
non-fiction ( marcgt )


Bibliography: leaves 159-163.
General Note:
Originally presented as the author's thesis (Ph. D.)--University of Florida.
This publication is being made available as part of the report series written by the faculty, staff, and students of the Coastal and Oceanographic Program of the Department of Civil and Coastal Engineering.
Statement of Responsibility:
by Yixin Yan.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright Yan Yixin. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
19278601 ( OCLC )

Full Text







The author would like to express his deepest appreciation to his adviser, Dr. Hsiang Wang, Professor of Coastal and Oceanographic Engineering, for his continuing advice and encouragement throughout this study, also to co-adviser Dr. James T.Kirby, Assistant Professor of Coastal and Oceanographic Engineering, for his guidance and support. Appreciation is extended to Dr. Robert G. Dean, Graduate Research Professor of Coastal and Oceanographic Engineering, Dr. D. Max Sheppard, Professor of Coastal and Oceanographic Engineering, Dr. Ulrich H. Kurzweg, Professor of Engineering Science and Dr. James E. Keesling, Professor of Mathematics, for their valuable suggestions in this research.
The author also wants to thank Dr. Y. Peter Sheng and Dr. Joseph L. Hammack for their Suggestions and help.
The author is indebted to all staff of the Coastal Engineering Laboratory at University of Florida, specially to Marc Perlin, Sidney Schofield, Vernon Sparkman, Charles Broward and Jimmy E. Joiner, for their cooperation and assistance with the experiments.
Thanks also go to Subarna B. Malakar for the helping with programming and Ms. Lillean A. Pieter for the drafting of figures.
Finally, the author would like to thank his wife Ran Liu for her support, encouragement and patience.

ACKNOWLEDGEMENTS ................................ ii
LIST OF FIGURES .................................... vi
LIST OF TABLES .................................... x
ABSTRACT ........................................ xi
1 INTRODUCTION ................................... 1
1.1 Statement of Problem .............................. 1
1.2 Past Study .. .. .. .... .. .... .. .. ... .. ... ... ... 2
1.3 Scope of Study .................................. 3
2.1 Introduction .. .. ... .. ... .... .. ... ... .. ... ... 5
2.2 Governing Equations ............................... 5
2.2.1 Depth-Averaged Equation of Continuity ................ 7
2.2.2 Depth-Averaged Equations of Momentum ............... 10
2.2.3 Wave Energy Conservation Equation .................. 14
2.2.4 Dispersion Relation ........................... 30
2.2.5 Radiation Stresses, Lateral Mixing Stresses and Bottom Stresses 35 2.2.6 Wave Breaking Criterion ........................ 42
2.3 Sum m ary ..................................... 44
3 NUMERICAL SCHEME . . . . . . . . . . . . . . . 46
3.1 Grid System and Definition of Variables . . . . . . . . . . 46
3.2 Finite Difference Scheme for Combined Continuity- Momentum Equation 47

. . 137

3.2.1 Linearized Implicit Model . . . . . . . . . . . .
3.2.2 Numerical Scheme for the Full Equations . . . . ... . .
3.2.3 Method of Calculation and Boundary Conditions
3.2.4 Stability Criteria .
3.3 Finite Difference Scheme for the Wave Equation . . . . . . . .
3.4 Method of Solution . . . . . . . . . . . . . . . .
4 PERFORMANCE OF THE MODEL . . . . . . . . . . . .
4.1 Wave Setup on Plane Beach . . . . . . . . . . . . . .
4.2 Open Lateral Boundary Conditions and Longshore Current Distributions on
a Plane Beach .
4.3 Nearshore Circulation . . . . . . . . . . . . . . .
4.4 Combined Refraction, Diffraction and Current Interaction Case . . . .

4.4.1 The Circular Shoal on a Flat Bottom . . . . . . . . .
4.4.2 Circular Symmetric Shoal on Plane Beach . . . . . . . .
5 EXPERIMENTAL COMPARISONS . . . . . . . . . . . .
5.1 Purpose of Experiments . . . . . . . . . . . . . . .
5.2 Facilities and Instruments . . . . . . . . . . . . . .
5.2.1 The Laboratory Facilities . . . . . . . . . . . .
5.2.2 Instrumentation . . .
5.3 Experimental Arrangement and Measurements . . . . . . . .
5.4 Experimental Results and Numerical Comparisons . . . . . . .
5.4.1 Wave Shoaling on Plane Beach . . . . . . . . . . .
5.4.2 Flow Field due to Inlet on Sloping Beach with no Waves . . . .
5.4.3 Inlet on Sloping Beach with Waves . . . . . . . . . .

66 70 81 81 82 82 86 89 90 91
94 95

64 66


B NON-SLIP BOTTOM BOUNDARY CONDITION APPLIED IN WAVE EQUATION .... .......................................... 149
BIOGRAPHICAL SKETCH ............................... 164


2.1 Outline of Studied Area ........................... 6
2.2 Variations of Dispersion Relation with U/co = -0.1 and c ...... Linear;
-.- Hedges; --- Whitham; Composite ..................... 33
2.3 Variations of Dispersion Relation with U/co = -0.1 and c, ...... Linear;
-.- Hedges; --- Whitham; Composite ..................... 34
2.4 Velocity Profiles for different friction models (2.107); - (2.117) 39
2.5 Longuet-Higgen's Longshore Current Profiles ................... 41
2.6 Variations of Eddy Vicosity Coefficients. = 0. ; ... =0.001;
=0.01; - Von Khrmhn's Model with upper-limit =0.001. unit= m2/s 43 3.1 Grid Scheme ........ ................................. 48
3.2 Definitions of Xnand (Ya)! ........................ 49
4.1 Set-up on a Plane Beach. Eq.(2.138); - Linear Wave; L Data 62
4.2 Longshore Current Generated by Obliquely Incident Waves. --- no mixing; breaking mixing coeff.=0.01m2/8; - breaking mixing coeff.
=0.01m2/8 and eddy vis. coeff.=0.005m /s ..................... 65
4.3 Velocity Field of Periodic Bottom Topography .................. 67
4.4 Wave Amplitude Variations along Centerline A0k0 = 0.16; -
A0k0 = 0.32 ........ .................................. 69
4.5 Wave Amplitude Variations across Centerline A0ko = 0.16; -
Aoko = 0.32 ........ .................................. 71
4.6 Wave Amplitude Variations across Centerline A0ko = 0.16; -
A0k0 = 0.32 ........ .................................. 72
4.7 Bottom Topography : a Circular Shoal with Parabolic Configuration on
the Plane Beach ........ ............................... 73
4.8 Description of Water Depth ..... ........................ ..74

4.9 Velocity Vector Field for Circular Shoal Located on the Plane Beach (no
inlet) ........ .....................................
4.10 Longshore Current Distributions at Lowest Boundary. No Inlet;
- with Inlet ........ ............................. ....
4.11 Velocity Vector field for Circular Shoal Located on the Plane Beach with'
Current from Inlet...............................

4.12 Contour Plot of Wave Amplitudes on the Circular Shoal Located on the
Plane Beach (no inlet) ...... ............................
4.13 Contour Plot of Wave Amplitudes on the Circular Shoal Located on the
Plane Beach with Current from Inlet ........................
5.1 Plan View of Wave Tank ...............................
5.2 Cross-Section Views of Wave Tank .... ....................
5.3 Cross-Sections of Measurement ...........................
5.4 Velocity Distribution at I=Idry .... ......................
5.5 Instrument Calibrations ................................
5.6 Instrument Connections ................................
5.7 Velocity Variations along the Beach, Y7(1,j) = 0; 0 17(1,j) =
a2/(2 sinh kh) ....... ................................
5.8 Wave Height Variations along the Beach, i7(1,j) = 0; 0 v7(1,j) =
a2/(2 sinh kh) ....... ................................
5.9 Comparisons between m = 43 with Ax = 0.12 m, -; and m = 83 with
Ax= 0.06 m, - at Il = 15,30 for m = 43 and 12 = 29,59 for m = 83

5.10 Wave Height Variations along the Beach (no current), - Energy Eq.
(2.138); - Linear Shoaling; x Measured Data ..............
5.11 Velocity profiles for Qo = 0.0056m3/sec. (no waves) .............
5.12 Velocity profiles for Qo = 0.0076m3/sec. (no waves) .............
5.13 Variation of x-Component Velocity. Case a: Q0 = 0.0056m3/s; T = 0.87
sec.; H0 = 0.97 cm. Numerical Model; o Measured Data ........
5.14 Comparisons of Calculating Velocities between with and without Waves.
Case a: Q0 = 0.0056m3/s, T = 0.87 sec.; H0 = 0.97 cm. with Waves; -- no Waves .......................................
5.15 Variation of x-Component Velocity. Case b: Qo = 0.0056m3/s; T = 1.16
sec.; H0 = 0.88 cm. Numerical Model; Measured Data ........

80 83
84 85 87 88 89
92 92 93

5.16 Comparisons of Calculating Velocities between with and without Waves.
Case b: Q0 = 0.0056m/s, T = 1.16 sec.; H0 = 0.88 cm. with Waves;
--- no Waves ........ ................................ 107
5.17 Variation of x-Component Velocity. Case c: Qo = 0.0056ms/8; T = 1.36
sec.; Ho = 0.97 cm. Numerical Model; o Measured Data.. ..... ..108
5.18 Comparisons of Calculating Velocities between with and without Waves.
Case c: Qo = 0.0056mS/s, T = 1.36 sec.; Ho = 0.97 cm. with Waves;
--- no Waves ....... ................................. 109
5.19 Variation of x-Component Velocity. Case d: Qo = 0.0076m3/s; T = 0.87
sec.; Ho = 1.23 cm. Numerical Model; o Measured Data ....... .110
5.20 Comparisons of Calculating Velocities between with and without Waves.
Case d: Q0 = 0.0076m3/s, T = 0.87 sec.; HO = 1.23 cm. with Waves;
--- no Waves ....... ................................. 111
5.21 Variation of x-Component Velocity. Case e: Qo = 0.0076M3/8; T = 1.16
sec.; H0 = 0.95 cm. Numerical Model; Measured Data ....... .112
5.22 Comparisons of Calculating Velocities between with and without waves.
Case e: Q0 = 0.0076m3/s, T = 1.16 sec.; HO = 0.95 cm. with Waves;
- - no Waves ....... ................................. 113
5.23 Variation of x-Component Velocity. Case f: Q0 = 0.0076m3/8; T = 1,36
sec.; Ho = 0.97 cm. Numerical Model; Measured Data ....... .114
5.24 Comparisons of Calculating Velocities between with and without Waves.
Case f: Qo = 0.0076m3/s, T = 1.36 sec.; H0 = 0.97 cm. with Waves;
--- no Waves ....... ................................. 115
5.25 Velocity Field and Contour Lines in Flat Bottom. f = 0 ......... ..116
5.26 Velocity Field and Contour Lines in Flat Bottom. f = 0.01 ....... ..117
5.27 Contour Plot of Velocity. Case i: Qo = 0.0056m3/s (no waves) ..... .118
5.28 Contour Plot of Velocity. Case a: Q0 = 0.0056m3/s; T = 0.87; Ho =
0.97 cm ........ .................................... 119
5.29 Contour Plot of Velocity. Case b: Qo = 0.0056m3/s; T = 1.16; Ho =
0.88 cm ........ .................................... 120
5.30 Contour Plot of Velocity. Case c: Qo = 0.0056m3/s; T = 1.36; Ho =
0.97 cm ........ .................................... 121
5.31 Contour Plot of Velocity. Case ii: Qo = 0.0076m3/s (no waves) . . 122
5.32 Contour Plot of Velocity. Case d: Qo = 0.0076m3/s; T = 0.87; Ho =
1.23 cm. ........ .................................... 123

5.33 Contour Plot of Velocity. Case e: Q0 = 0.0076m$/s; T = 1.16; H0 =
0.95 cm ........ .................................... 124
5.34 Contour Plot of Velocity. Case f: Q0 = 0.0076m3/s; T = 1'36; Ho =
0.97 cm ...... .......................... ........... 125
5.35 Bottom Friction Coefficients f as Function of u../U ... .......... : 126
5.36 Wave Amplitude Variations. Case a: Qo = 0.0056m3/s; T = 0.87 sec.;
Ho = 0.97 cm. Numerical Model; A Data ................... 127
5.37 Wave Amplitude Variations. Case b: Qo = 0.0056m3/s; T = 1.16 sec.;
Ho = 0.88 cm. Numerical Model; A Data ................... 128
5.38 Wave Amplitude Variations. Case c: Qo = 0.0056m3/s; T = 1.36 sec.;
Ho = 0.97 cm. Numerical Model; A Data ................... 129
5.39 Wave Amplitude Variations. Case d: Qo = 0.0076m3/s; T = 0.87 sec.;
Ho = 1.23 cm. Numerical Model; A Data ................... 130
5.40 Wave Amplitude Variations. Case e: Qo = 0.0076mS/s; T = 1.16 sec.;
Ho = 0.95 cm. Numerical Model; A Data ................... 131
5.41 Wave Amplitude Variations. Case f: Q0 = 0.0076m3/s; T = 1.36 sec.;
Ho = 0.97 cm. Numerical Model; A Data ................... 132
5.42 Wave Amplitude Variations along Centerline .................. 133


Test Condition...........................
The values of c. and f in the Tests.. .. .. .. .. ..
Bottom Friction Coefficients Used in the Model . .. Incident Wave Height.... .. .. .. .. .. .. .. .. .

. . .. 90
. .. .. 96
. . .. 98

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
August 1987
Chairman: Dr. Hsiang Wang
Major Department: Engineering Sciences
A numerical model for the prediction of nearshore hydrodynamics of an inlet-beach system is developed in this study. The theoretical formulation is based on the time- and depth-averaged equations of continuity and momentum, and a 3rd-order wave energy equation that retains the important nearshore effects of refraction, diffraction, current-wave interaction and energy dissipation. Variational principle as well as a parallel method by using perturbation technique is applied to arrive at the same wave equation. The application of 3rd-order wave equation is restricted to small incident wave angle and slowly varying topography.
A numerical scheme is developed to solve the above governing equations. The doublesweep Crank-Nicolson finite difference scheme is adopted for the model. In each spatial sweeping, a recently developed two time-level implicit method by Sheng is applied. The numerical model is tested against existing numerical models of simpler cases and is found to yield reasonable comparisons for all cases.
Laboratory experiments are conducted to verify the numerical model and it is found that the model performs well in predicting current field but grossly underpredicts the wave amplitude in shallow water, particular, in the center region of the inlet jet.

1.1 Statement of Problem
Tidal and river inlets are prominent coastal features. They play a vital role in the evolution process of coastlines. Along the east coast and gulf coast of the United States, inlets are numerous because of the extensive barrier island system. In these regions, the rules of inlets are particularly significant owing to the heavy development pressure and the competing interests of various human activities. To maintain navigation and water commerce, for instance, often requires inlets to be stabilized with man-made structures. The disruption of littoral processes owing to the presence of Inlet structures, on the other hand, poses a major threat to down drift beaches, thus diminishing the value of recreational beaches and endangering waterfront structures. Many coastal communities spent millions of dollars installing coastal protective structures or simply transporting sand from one place to another. However, these measures do not always render satisfactory results owing partly to our poor understanding on complicated interactions of an inlet-beach system.
The hydrographic and shoreline changes of an inlet-beach system are the consequence of sediment movement which, in turn, is caused by the hydrodynamic forces exerted on the system. Therefore, to better predict the behavior of inlet-beach interaction one must first be able to predict the hydrodynamics.
Physically, the problem in hand is a complicated one. Both wave and current field are modified by the nearshore topography and the inlet geometry. In addition, the spatially non-uniform current will interact with the unsteady flow field induced by waves.
Mathematically, there is no unified wave theory in the first place capable of describing waves traveling through large areas of varied depth. Furthermore, the current field that we

are dealing with has a different scale of variation from that of the local wave field. The nonlinear interaction between the two systems is expected to be important as well-as the non-linear behavior of the wave field.
1.2 Past Study
In the absence of an inlet, the water motion in the nearshore zone is primarily induced by waves from offshore as modified by the bottom topography and/or shore structure. The influence of bottom topography on wave direction, height and length has long been recognized and as early as in the 1940s analytical methods were already being developed to solved the so called wave refraction problem (Johnson et al. 1948; Arther et al. 1952) by tracing wave rays based on the principle of geometrical optics. Since then, significant progress has been made in treating the shallow water wave transformation phenomenon. Arthur (1950) included the effect of current into the ray meth, ,, Dobson (1967) developed one of the first computer codes to trace wave ray over an irregular bottom. Skovgaard et al. (1975) extended the ray method to include the effects of wave energy dissipation due to bottom friction. Noda et al. (1974) developed a finite difference scheme based on the conservations of wave number and wave energy, that compute spatial distributions of wave number and wave energy in shallow water. His method was extended by Wang and Yang (1981) and Chen and Wang (1983) to wave spectral transformation and the inclusion of non-stationarity, thus enabling the computation of spatial and temporal distributions of random waves. Birkemier and Dalrymple (1976) and Ebersole and Dalrymple (1979), on the other hand, expanded Noda's method by including the depth averaged continuity and momentum equations to computer wave induced nearshore circulations. In the presence of an inlet, the nearshore wave field is further modified due to the interactions with a spatially non-uniform current and the diffraction effect of inlet jetty structures.
Traditionally, the problems of wave refraction and diffraction were treated separately. Penny and Price (1944,1952) were among the first to apply Sommerfeld's (1896) solution to wave diffraction around a sharp obstacle such as a semi-infinite thin breakwater. Other

solutions of similar nature were obtained for different geometries (Wiegel,1962; Memos,1976; Monfetusco,1968; Goda et al.,1971).
Only recently, Berkhoff (1972) developed what is now referred to as the mild slope equation that governs the combined refraction- diffraction phenomenon. Various solutions became available on wave field modified by different beach-structure configurations such as waves in the vicinity of a thin groin (Liu and Lozano,1979; Liu,1982), edge waves by barredbeach topography (Kirby et al.,1981), and waves over circular shoal (Radder, 1979), etc. Booij (1981) further introduced the ambient current element into the basic wave transformation equation and derived the governing equation on variational principle. The resulting equation is a hyperbolic type and is, in general, difficult to solve. Booij's equation reduces to a parabolic type similar to Berkhoff's if the mild slope assumption is invoked. A number of simple examples were given by him. Kirby (1983) systematically studied the current-wave interaction problem and proposed a weekly non-linear model appropriate to thr Stokes wave expansion.
The present study extends the Booij's and the Kirby's model to the problem of interaction between water waves and inlet flows.
1.3 Scope of Study
The primary purpose of this study is to develop a numerical model for the prediction of wave characteristics and current field of an inlet beach system. Specifically, the mathematical formulation of the problem was performed first. This task entailed the establishment of the governing equations sufficient to solve the variables of interest. Out of necessity, simplified assumption are made such that the problem is amenable to the present computational capability.
First of all, the fluid is treated as incompressible. Second, the surface shear stress as well as part of the turbulent shear stress (normal component) is also neglected. Third, also one of the most major ones is that the irrotationality is assumed in the wave equation allowing the existence of a wave velocity potential and a current velocity potential, although

such an assumption is in contradiction to both physical reality and mathematical results of the combined velocity field. The equations governing the mean flow -including both continuity and momentum equations are depth averaged. The equation governing the wave motion with due considerations of refraction, diffraction and current interaction isquasithree dimensional owing to the velocity potential assumption. Two parallel methods were used to derived the wave energy equation: one based on variational principle and the other on perturbation method. By introducing proper scaling factors, the wave energy equation is further simplified by containing terms of only the same order of magnitude. Other relationships including dispersion relation, radiation stresses, bottom shears and lateral mixing compatible to the system under consideration are also formulated.
A numerical scheme is developed to solve the governing equations by double sweep marching method. The scheme is implicit which allows for a large time step, especially the two time level fully alternating direction implicit scheme utilized by Sheng (1981).
Attempts to verify the numerical model are made by comparing the present model with known analytical and numerical solutions of simple geometries and with laboratory experiments carried out in the Coastal and Oceanographic Engineering Laboratory. Due to wave basin limitations, the experiments were carried out only for plane beach with an inlet normal to the contours.

2.1 Introduction
In order to solve the flow field in a specified region, the governing equations with proper initial and boundary conditions must be provided. Furthermore whether a unique solution exists in the domain of interest or whether mathematical singularities will appear that will cause the equations failing to give sensible solutions should also be examined.
In this chapter, the theoretical development of the governing equations for both the mean flow field and the wave field is given in detail.
There are five basic equations employed in the computation including four flow equations: one continuity equation, two horizontal momentum equations and one wave equation, and one equation relating wave number and wave period.
2.2 Governing Equations
The basic equations, as mentioned before, are based upon the principle of conservation of mass, conservation of momentum, and conservation of wave energy and upon the gravitational wave dispersion characteristics. The wave equation under consideration includes the effects of wave refraction, diffraction and shoaling during transformation from deep water to shallow water regions. Thes e equations, together with initial and boundary conditions, give a mathematical representation of the physical problem of interest.
Since the problem of interest entails waves propagating towards shore over a complex bathymetry, wave refraction, diffraction and breaking are expected to occur as well as current-wave interaction due to inlet discharge. Figure 2.1 depicts the physical system being studied. The problem has two apparent time scales: a short time scale associated with the

Inlet 11-:// .9/
- -Shoal

Figure 2.1: Outline of Studied Area

oscillation of the incident waves, and a larger time scale over which the characteristics of the incident wave and mean flow, such as wave height, incident wave angle and current are perceptibly modified. The larger time scale may also include the effect of changes in wind, tides and topography. Since our attention here is towards mean quantities,.-nany of.them may be averaged over a wave period to remove the direct effect of oscillation. In addition, since we are mainly interested in the net transport quantities the knowledge in the detailed structure of the velocity distribution over depth is not needed and the equations may be averaged over depth, reducing the entire problem to two dimensional in the horizontal plane.
In the derivation of wave energy conservation equation, due to the extreme complexity in mathematics, some restrictions have to be imposed on the orders of magnitude of wave height, bottom slope and current strength.
In the model, the quantities to be determined are as follows: the horizontal components of the mean current including wave-induced current and divergent current from the inlet; the local wave amplitude under the influence of refraction, diffraction and breaking ; the local wave angle and wave number and the mean water surface level.
2.2.1 Depth-Averaged Equation of Continuity
A Cartesian coordinate system is defined with x in the onshore direction, normal to the coastal line, y in the longshore direction and z vertically upward, as shown in Fig.2.1.
The summation convention adopted here states that in Cartesian coordinates whenever the same letter subscript occurs twice in a term,that subscript is to be given all possible values and the results in horizontal space add together.
The continuity equation reads
ap + apu +apW0 (2.1)
where i = 1, 2 represents xy direction, respectively. p is the water density and ui = {u, v} represents the {x, y} components of velocity,respectively, w is the velocity component in the z direction.

In application, water density which only varies slightly is considered as constant. The continuity equation reduces to att 8w
i+ = 0 : (2.2)
ax 8z
Integrating Eq.(2.2) over depth, we obtain
h ( X + -)dz = 0 (2.3)
where r is free surface elevation, and h is water depth. Using Leibniz's rule which states that
a a(s) a(s) 8f Ba
f (x, z)dz = -dz + f (x, c) f (x, P)
8T pc) ff ') p(z) 8z z x
Eq.(2.3) can be written as
a fu ar ah
a xifudz u |,7 ai- uj I-, + W O -W |-h= 0 (2.4)
in which the subscripts after vertical stroke are at which the value being calculated. The bottorn boundary condition (BBC) and the kinematic free surface boundary condition (KFSBC) are given as
u8 + w = 0 (BBC) at z = -h(x,y) (2.5)
- w uy- (KFSBC) at z= t(x,y,t) (2.6)
Substituting (2.5) and (2.6) into (2.4) results in a?) a rI
- + a udz = 0 (2.7)
t hz
We now introduce
ui= Ut + u:


where UiJ are time and depth independent uniform current velocities, j7 is the mean surface elevation, uti and 17' are the fluctuating components. u'i can be separated further into two parts:
u's=u +S _. (2.9)
where uF and u! are the slowly varying wave motion and the fast varying turbulent fluctuations, respectively. The properties of these quantities, after taking average over wave period, are expressed as
U' = 0= o
where the upper bar represents the time average over one wave period f= f dt
Therefore, (2.7) becomes, after taking time average ai7 a (2.10)
a...ffh u'dz Ph uvdz
(h + !) (h +I)
the final form of the mean continuity equation is arrived at ay + a (UD) = 0 (2.11)
where D = h+1 is total depth and Ui = 7+fii is total mean current velocity averaged over one wave period with U = {U,V}. Equation (2.11) represents conservation of mass per unit surface area under the assumption of water density being constant in space and time, and the bottom unchanged with time.

2.2.2 Depth-Averaged Equations of Momentum
The horizontal depth-averaged momentum equations used here were originally derived by Ebersole and Dalrymple(1979). A brief summary was given here. The Navier-Stokes equation in the horizontal (x,y) directions takes the following form: Bu a Bug Bui 1 Op 18rij 18rzi
-g; + Uj-y + woT I p,+ + (2.12)
Bt y 8 p zi p Byx p 8z
and in vertical (z) direction it is
Ow aw aw 1 ap 1 8rj, 1Orzz
51+ u- +w z + "+ (2.13)
i = 1,2, j= 1,2
where p is the total pressure including dynamic and static pressures and riq is the shear stress, where the first subscript is the surface that the stress acts on and the second subscript is the direction of shear stress. Introducing the continuity equation (2.2) into (2.12) and (2.13) we have
Bu, Ouiu Ouiw 1 Op 1 Ori 1 8r,
- + + + p 1 + -+ (2.14)
7t 8xi 8z px p Oxi p 8z
8w Owua aww 1ap 1 rzi 1 rzz
Ft + =+ + + - (2.15)
Bx i z P p xj x p( z
Integrating over depth and utilizing the Leibeniz's rule, the following term by term results are obtained:
17 au Uj a 7 9
Ou,u, 8 8sB Oh
i) a-dz = -] uydz- uu1 49 Uis -h
h 3_] $Zj -h 82w
m -) W dz = uiw j1, -u W I-h
Thus the left hand side of Eq.(2.14) becomes O O a17
LHS = T uidz +u iujdz "u- (F+uj 117 1w17)
f a~jfh a3

-i 1-h (Uj 1-h + W I-h)
Recognizing the last two brackets of LHS are the KFSBC and BBC, respectively, we have LHS = uidz + uiudz
Again introducing the definition of m t
Ui = Ui + Ue + u!
into the above equation and time averaging over a wave period yields term by term i) s_ d = (7, + up + e )dz = -(UD)
where Ui and D are total current velocity and total water depth respectively as defined earlier, and
i)indz = (Fe + up + U!)(Fi + "7 + uj)dt
ii) "f-h+
Noting that there is no interaction between wave orbital velocity and turbulent velocity, Phillips (1977), the above term can be expanded as follows
a ,a -a dz+
7 UUidz + (Ugu + Uifu)dz + u udz +
8zi f-h 8z xy f-h
= -UUi)+ a(Ui'LiD)+ (UiD)+--[ ufu.'dz +
axi -h=a( +iaD) + (fij D) + a(- fij D) + a 7uF udz +
a_ o f dz
Introducing Ui = Ui + fi into the above expression, the LHS of Eq.(2.14) becomes
a (UD) + az (UgUD) 8ai (2 D) + 8a. J7 u-u-dz +

- utudz (2.16)
The right hand side of (2.14) is, after integrating over depth and averaging over time,
RHS = f '' dz+-J --idz+ 1faLsdz
P _h _h-sO
p a 8i P -4 8i P -a z
here water density is considered as constant. After neglecting the second term which is the lateral shear stress due to viscosity, we obtain
1 8 1 8 r; 1 8h 1 1
RHf pdz + -p In +-P -, +-G, I -Gi -,
RHS= a 181- all -+ 1PI-h ~+ I~z 117 -z
P i -h p 8 p 8 p p
The datum of pressure is chosen as the atmosphere pressure, so that at free surface, p 1, is zero. The mean pressure at the bottom p I-h consists of two parts: hydrostatic pressure and wave induced dynamic pressure, viz.
P I-h = -pg(h+4) + Pd |= -pgD +Pd l-h
where Pd is the dynamic pressure. Multiplying the above equation by I h we have
1 ah 1 a B 1 Oh
-p I-h --(gD ) -gD + P I-h
p 8:4 2 8i 82; p 8zi
The RHS can now be written as
1 a f 1 a 2) y 1 ah 1
-- O- pdz + (gD) gD'- + -Pd I- + -I-h
p 82i -h 2 8zy 8i p By p
--Yb (2.17)
Yei = rTi is time averaged surface shear stress and
U = ri -h is time averaged bottom shear stress.

Finally, combining (2.16) with (2.17), the time averaged and depth integrated horizontal momentum equation is given by
a(UD) +ai(UgUiD) i aud aud
X-(1D + ij(UiU)D+ -- uvuvs~dz+ +--j h tdz
18 18 1 ah 1 1
Jpdz + -(gD -gD + -p + -ri -rbi
P 8zi -h 2 i 84 p 84 p P
or rearranged as
a(UgD) + azi(UgUiD) = -9D~z p-7_ [z i d4z+ p ?'uIpupdz
a, D aU D 1 8
1 2 1 Ah 1 1 a8 s t t
--P9D6, pf4 s, D] + -h + -f + uujdz (2.18)
2 p 8z P p p az J-h
Collectively, the second to the fifth terms on the RHS are known as the radiation stress or the excess momentum flux due to the presence of waves, i.e.,
Sij= f pidz + p upupdz pgDSi, pii jD (2.19)
-h -h2
1 i=j
Considering that the product of bottom slope and the mean dynamic pressure at the bottom is negligible, we obtain
8 8 BjT 18S-- D 8
-(UiD) + (UiUjD) = -gD- (a -(puVu) +
at B8z px 8 PX
1 1
-ri -rbi (2.20)
after assuming uit is independent of depth. The term -puu is Reynolds stresses and their normal component (when i = j) can usually be neglected.

Writing in component form yields X-Direction
a a U2 a I as." I_ __(UD) + (UD) + -(UVD) = _ y +
D BY; 1 1
S0-+ -res -res (2.21)
p y pz p b
a a a (V2D g~ I_1 as.'Y 1 asyy
(VD)+ -(UVD)+ a = D -Da "
at Tz y ay p 8x p ay
D aT, 1 1
Sz + -reS pre (2.22)
where re = -puu, is the lateral friction due to Reynolds stress. In the model, ri will not be considered. The possible formations of Tv and ft will discuss later in this chapter. They are flow related variables.
The average continuity and momentum equations can be linearized. Noda et al.(1974) gave a simple form of those equations after neglecting the nonlinear and time dependent terms in (2.11) and (2.20). They are a8 (UgD) = 0 (2.23)
0 = -gD p -a, (2.24)
8zi p ax p
The surface shear stresses and lateral shear stresses are also neglected.
2.2.3 Wave Energy Conservation Equation
The equation considering both refraction and diffraction with current has been developed by Boiij (1981) and corrections were made later by Kirby (1984). Both used variational principle. In this study, a method that modified Kirby's derivation is developed. In

parallel, the perturbation technique is also applied. Same results are obtained and are given in Appendix 1. Here only the derivation from variational principle is presented.
The variational principle, first applied in the irrotational flow motion by Luke (1967), states that the variation of a certain conservative quantity, often an integral, vanishes if the function that describes the evolution of the physical process is subjected to small variations. It is represented as
6ff Ldxdt] = 0 (2.25)
where x = xf + y is in horizontal space. For irrotational flow, the energy is conserved and L can be chosen as the Bernoulli sum: L J t + (V4)2 + gz]dz (2.26)
V is the three dimensional gradient operator,
a" a__a,
+ay3 + zk
," and k are unit vectors in x, y and z directions, respectively.
The potential function consists of two parts
, = + 0' (2.27)
where 0.. and 0, are wave potential and current potential functions,respectively. A weak point about the expression of 0, is that the mean current field may be rotational, since we do not know a priori whether the mean flow field derived from momentum equations satisfies the condition of irrotationality everywhere, i.e., curl U = 0

However, the irrotational assumption will be used here.

For a slowly varying bottom topography 4 can be expressed as (Kirby, 1983):
O(x,y,z,t) = 60'('UszUY)+ E 1(z,1s/2Y"UYApZzt)+
C24 w2 (X, A1/2y, py, pX, z, t) (2.28)
where p is a scale parameter with an order of magnitude O(Vhh/kh) representing the spatial change over the space of wave length in which Vh is horizontal spatial gradient and k = 27r/L is wave number where L is wave length. 6 is the strength parameter of current defined as Vh~c/V'g.Tc is the Stokes wave steepness parameter c = 0(ak), a is the wave amplitude.
The spatial scale parameter p is introduced to differentiate two different spatial scales in the model: one is related to the fast wave-like functions, and the other one is connected to the slow spatial modulations, such as bottom topography, the changes of current speed, wave length etc. Therefore, we have
X = x+IIx-X+X2 (2.29)
Y = p11/2y +p Y + Y2 (2.30)
and their derivatives are replaced by
a a a
T-X "- gX + PaX2 (.1
a 1 1/2 a a
P +' (2.32)
The scales as selected imply that spatial derivative with respect to phase function occurs only in the x-direction which is true when the angle of wave propagating is small with respect to the x-axis, i.e., it is assumed that no fast wave-like variations occur in the ydirection. It is consistent with the forward scattering approach that the complex amplitude A is allowed to vary as 0(M1/2) in the y-direction. The problem is further considered to be quasi-steady; time derivatives only appear in wave term.
There are three magnitude parameters 6, p and c in the formulation. Different ratios among those parameters will yield different forms of the energy equation. In the present

study, we adopt Kirby's (1983) model given the relationship as follows
6 =1 (2.33)
which represent slow varying topography, strong mean current, and small incident wave angles.
In equation (2.28) 4, and 4w2 are first and second order potential function, respectively. 4.l and 0w2 are expressed as follows: .1 = Am(x, p'l/2y, px, py, t) fm(px, py, z) (2.34)
0w2 = A'(x, p/2y, px, py, t)fm(p py, z) (2.35)
The Aj can be split into two parts: a modulated amplitude function and a phase function:
A(, p, py, t) =A (p1/2y, pX, py)ei" + c.c (2.36)
where AL is the complex amplitude, 4 is the phase function and c.c is the complex conjugate. The phase function is defined as S= kldx wt (2.37)
where k, is the wave number in main (or x) direction and w = 2r/T is the angular wave frequency in which T is the wave period. The real wave direction which is different from the main direction is absorbed in the amplitude function. The subscripts and superscripts in the equations represent the order and harmonics, respectively. However the absolute superscript values should be no bigger than the subscript values. Substituting (2.34) and (2.35) into (2.28), yields
4 = C-2,Oc(,C2z,C2y) + cA(z, cy,e,c2y,t)fm(e, y, ) +
e2Am(x, Cy, c2x, e2y, t)fm(e2, ey,z) (2.38)
where 6 and p have been replaced by C-2 and c2, respectively, by virtue of Eq.(2.33).
Substituting (2.38) into (2.26), we obtain with (2.37), term by term

Oe = cA'fm + C2Ag f (


VO = Vh4O + effA' I + E2f-A- i + C f VhA- + c3AVhf
+e'C 2 +e ffm It + esffVAAF +e'AC4 f
+ f?7Ac2X' + YfIA',j + C f2VhA' + eA'Vhf
+(cAmfJ + c2A'ff)k (2.40)
Vh is the horizontal gradient operator, Vh= TS +
Therefore, by definition we have
VO = U + Vj (2.41)
where velocities include wave induced mean current.
The first two terms in Eq.(2.40) represents linear wave superimposed upon uniform current over flat bottom. In Eq.(2.26), the (V)2 term when computed to fourth order becomes
(V = (Vh c)2 + eC2ff'n ARA +c 4f fn A An +4 ffA n A n
+e'fgfz-IA~An + nf2nfAFA" + e2Uf-Am +C 2VfA-Yx
+E 32fmVhc, VhA + eC32A-VhA Vhfm + C22UffAx + c32VffAmy
+C42ffVh, VhAm + C42AmVhOc Vhff + C 2fffnA A4 m fn n $r n +a C3 n f mafn
+e 2fffX2AA + c32fmffnAm A + 2AmAzf~f (2.42)
The subscripts X, Y and z are partial derivatives. The superscripts m and n represent interaction between two harmonics.
Free surface envelope can be considered as two parts: one is related to mean current as t7o and another one is driven by waves as r/'. Later we will see to has the same order as IU I.

q(X, V, = 7o(x, y) + 1'(x, y, t)


or in the stretched coordinates
1(x E, cy, 2, 2Y,t) = (6 2x, 2Y) + (z, Y, c2X, C yt) (2.44) where 7' has the order 0(c) at most and it can be expressed as 17' = 07 + + 0(c ) (2.45)
where t7i and 72 are first and second order surface elevations, respectively. By virtue of (2.45), (2.43) can be rewritten as
7 = vo + C'i + C 27 (2.46)
Expanding (2.26) in Taylor series about z = iGo, results in
L = [,t + (V)2]dz + g7 ~gh
-h2 2
= lot + 1(VO)']dz + '[ lot + I(VO)'].=no + 1 17 +, (+ 1 ]=
1 37 1 1 1
6 8-z ot + 2(V) ]'=no + 2 gh+*** (2.47)
Introducing a shifting coordinate defined as z = z o, h' = h + ro
and substituting into (2.47), we get, after dropping the primes on h' and z' for convenience,
L = [#t + 1(V )']dz + 1'[, + (VO)2]=o+ 1 7 2[t + I (VO)2]z=o
2-2 Z 2
1 1(V )2. + 17 2 gho (2.48)
where we denote h0 as still water depth. The expansion is carried out to fourth term only as the higher order terms will not affect the third order wave equation.

Certain integrations in (2.48) that involve fie of Eqs.(2.34) and (2.35) are defined as follows for future use:

ff"dz = Q'j


f"fjdz = Py"

f = Ri,



where i, j, m and n are free indexes and subscript z represents partial derivatives. From Stokes wave theory, the well-known expressions for the fm and f," are

f -1 cosh k(h + z) =1 cosh kh
f2 g_2= cosh 2k(h + z)
cosh 2kh

The typical results in (2.49)-(2.51) that will be used later are given as follows:

cosh kh
11 k sinh3 h
1+2cosh kh
S3k cosh kh sinh kh
p22 1 h sinh 4kh
sinh8 kh (-2 + -8k:
Rif -(e' k'cc,)
I 2 4k
R2 = coshkhsinhkh
123 cosh kh sinh kch

(2.52) (2.53)

= k2 sinh 4kh
2sinh kh k -h)
where c is the wave celerity, c = LT c is the group velocity, cg = c n with n = (1 2kh
2(1 +,-fh$)
a = /gk tanh kh is the local angular wave frequency, it will be discussed in section (2.2.4).
Substituting (2.39) and (2.42) into (2.48),with the definition (2.49)-(2.51), we have
I1pmnI4 mn n 1 2 mn n
L = AcQA' + C'P A An + C ePITAIY AnY +e R AiAn
I~~ ~ it2R i y11~AA
+eQ UAx + c2 QmVAmy, + E31Vh" e V* A7 + c'Pz"AmAn
+ P2 A hxAx + R A ARAA + jg(c217o0q + cE + c32,7mn) + it
+Esl f2A n+ C'm nA n + n(Vh)2 3 clfmfAm A~nx
+ 2 z 2 1 1 nx A + Vfl A-" i
1 k f m f 27n Uf A n am c3fVfA nn +2Vt A n
+ c 2Um1nAL + eg 2mnAn + n + e fk 111 A AN
14,ra ,nn ~snUf Ax+l4r n + C47f.fnU nrA~
+c j~m + A + c3imf~A +Uf?'jf~~A x + S7VAA 1
+ 4 I fizzA A + fA ~ + mV f i, A + nA
24 ~ ~ ~ n n 12 4,nyg An 2L 1 2 17
14 n
it~lzl ~2tJz1 2F it 2 iJz'Y i 2
+ek 14 Uf"A + mA InzzA" + 3 IzA (2.54)
+c lI UfjnAx + CE z + C3vJz
where the condition of Vhfn = 0 has been involved. All terms are evaluated at z = 0. For simplicity, we only keep terms containing ,l and A,, since as will be shown L varies with respect to i and Al only. In (2.54), the superscripts are free indexes (nonrepeating) which may not assume the same value for different terms, but they give same value in one term. The terms of (" are defined as
Sk = mng7, m+ n = k (2.55)

S= S mn7i'1t l71 k+m+n =j


when m = n or k = m = n, 6mn and Skmn equal to unit, otherwise they assume the values of 2 and 3, respectively. Those give
L= L o cL+e L e 4+--(257
17 22%7102 (e)
21 1
-- 71171
1 1 1
-1 =31f~l1171 + 3~0f?1
After introducing the following
L = Lo + cL1 e2L2 + 3 m+ E4L4 +2 (2.57)
Equation (2.54) can be separated by grouping the terms of the same order:
L, =QmAl + Qm'UAx + gto7 + 171m(Vh'O)2 (2.58)
1mn Am An mn m n + Q1VAm m + n
+q, mU nAN X(2.59)
L3 = QmVhc VhA 1 + m 1A A + RmnAmA1n + mn + 2nA
n 1 k nn fmn n .n
+t72mf,"Al + ffmfAmAx + f IfizAf, Aj + iqfVfAnAY nU+1 -n + +qfU fnA + tmU fnA~ + 192m 11Alt 1sUfl.A (2.60)
+1imf~AX 2i 2 i

L 1 = 2 P y' A r~" lar z n a " A n p iDn r Anx +n f f* V hAn + V f 2 nA 2
+ L4"A" + Afm A A + 7mVf "A
+ of n AmA.x n 1m 7 n.fgi + 2n f AnA", (fn 1 n n. 7m xnA
1 I1 1 2 1A 1 2 V I
+2k1ms 1AAn 2 + m n + AYl 2 m 1 1xA U
+ j1 2 ,n n m A + !j n 2 +t 2 +t 1
1 nA A + Vfnn 1 +n + 17n
1 1
+-'m s"f,, + 6 mU fzAx (2.61)
There is no L0 term, as it contains neither 17l nor A, term. For each order, varying with respect to 71 and A,, we will obtain a series of equations. The variational principle gives
8L, a0Li a-L
B- a(a-) Vh[(-)] = De (2.62)
B B 7t a B, a(VAB) where B represents either 171 or A1, D, is the dissipation term that is caused by bottom friction and turbulence that eventually are converted into internal energy, which is irreversible. (One possible formation of De is given in Appendix A by using perturbation method.) Equation (2.62) can be written in detail as
[a- ( a ) (a) a ) ( )]Li = D, (2.63)
B aX2 (aX,2 1Yi 29Y2 8 I
The following discussions are based on the orders of L. 1) 0(C)
Varying (2.58) with respect to 17r yields
170 = -g (Vh ) (2.64)
which states that the surface elevation by mean current is simply balanced by kinetic pressure. The same result is obtained if one uses Bernoulli equation at order 0(1). 2) 0(e2)
Equation (2.59) varied with respect to j7, results in
n" = --fn(Ai' + UA nx)

_ I, DA (2.65)
g Dt
is the total derivative related to the x- direction D a a
5-t t a
For n = 0 or the zeroth harmonic, A' is independent of X and t, which results in n} =O for n=0 (2.66)
From Stokes linear waves, we know that Al = -,g a & (2.67)
AlI= a----(2.68)
where a is a complex amplitude and a, is defined as a1 = w Ukl (2.69)
which is the apparent local angular frequency that represents the waves riding on unidirectional flow. According to linear wave theory, the dispersion relationship is a = w Uk, Vk2 (2.70)
where kj = k cos 0 and k2 = k sin 0 are the wave numbers in x and y direction, respectively and 0 is the incident wave angle as shown in Fig.2.1. From (2.69) and (2.70), we find 01 = a + Vk2 = a + ON
since for small incident wave angle, sin 0 0 = 0(c). Hence for n = 1 and n = 1,f = 4I = 1 at z = 0, we have
1 1DA l a e (2.71)
71 = g Dt 2

1 1DA-1 = ae-i (2.72)
g Dt 2
both are the first harmonic solutions and the same as given by linear wave theory. ,7 and 1j' have the same amplitudes, but traveling in opposite directions. As mentioned before, the complex amplitude "a" absorbs the difference between real wave direction and its main direction, since phase function 0 is only defined in x the direction. Next, varying L2 with respect to Af, (j = m or n), we have Rn"A f"m- PA A~xx Ufl"mx = 0 (2.73)
a) let m = 0 (the zeroth harmonic), it gives ROA = 0
11 1=l--0
A= 0 or fo_= 0 from (2.51) (2.74)
Later we can prove that if fiz = 0, then fo = 0. b) Considering the first harmonic, we obtain [cc9(k2 k2)] A = 0 (2.75)
It is important to note that (2.75) is of the order of 0(E3), since k2 k = k = k sin2 0 = 0(e0). This term later will be combined with other third order terms to form the wave equation which is of the order of 0(cs3) 3) 0o(c)
Varying L3 with respect to ,i and Ai (j = n, m or k) respectively, we obtain
1 1
g'y" + f,"Afe + 2 + 2 A +Vf"A"Y1 + A
+UfenAf + i71UfnzAlx = 0 (2.76)
R + 1f fA1 Vh (QTVhc) 2 PAxx ft2 ff, (A1x)x
-V f -t. UfWrfx fn, k 0 (2.77)-- /lX Z r/1 r1 t -- ] l 1 17 1X - 0 ( 2 .7 7 )

we emphasize again that the superscripts only take the same values in one term; there is no relationship among different terms. There are only two unknowns, 72 and A', in the above two equations and so that they can be solved. a) Collecting all the terms with zeroth harmonics in (2.76), we get
0 .. {f12A +- 12 1A 1 + 1 --1 + + -[1
-= 1 zA + fAA f+1 [At]+fyIAx
Here the condition of f' = fi-1 has been applied. The values on the right hand side of above equation are known, so that gk 9k tanh kh k tanh kh
2 1[gkl gk tanh2kh ktahkh] (2.78)
4a 4o-22
Since we only consider up to 0(0), kxand a1 can be replaced by k and o, and vice versa. Equation (2.78) becomes 0o = a k ) (2.79)
2sinh 2kh + () (2.79)
This is the set-down of mean free surface with reference to the shifted coordinates. From (2.77), we obtain
A0 = Vh- (QVhoc) (2.80)
where m can assume any arbitary value. Hence, A can not be determined as f remains unknown. However we will show later that A0f2 is non-zero. From (2.80), the previous statement in case 2 is confirmed, that is, let m = 0, we must have flo = 0 for fr\ = 0. Otherwise, A will approach infinity. b) The first order harmonics in (2.76) yield the following solution I = 1[VAy if2aiA'] (2.81)
From (2.77), together with (2.81), we have 1~ ~ pmt 1_ 22m lfifm)
A2122 1
A)(Rjms + kfPm1 1gf2 ff) = 0 (2.82)

For m = 1, if we assume that f2 has the same form as fl, the terms in the bracket are identically equal zero. But, if m is any value less than or equal to 1, i.e., m < 1 |, we may select A' = 0 such that
S= -VA}y = iVaye (2.83)
which states that the surface elevation of the first harmonic of the second order of the wave is balanced by the variation of wave amplitude modulation in longshore direction. c) For the first negative harmonic, (2.76) and (2.77) give
A2' = 0 (2.84)
171 = -i-V ay, e-i (2.85)
d) Considering the second harmonic, after lower harmonic terms have been dropped, we obtain from (2.76) and (2.77), respectively
2rl 2, +A 1 .12 1 1a2 12 2 2 11- 1 1 1
g9 + 2A +A + fl1.A1 + f2UA2x + fI't79At + flzUtA'x = 0
2 2 2
R 12 2 -12 141 12 A2 2 1 X_ 12_2t.117
R A2 + fli71A P2Axx 12 (Aix1)x -' Ux 1 i
-f i1 1 = 0
Solving two equations for two unknowns A2 and 9 yield
A2 = -i a Iei2 +0(E) (2.86)
j Ial'kcosh kh
2 = a cos h (2 + cosh 2kh)e"2 + 0(cs) (2.87)
Sr 8 sinh kh
They are the same as the second order Stokes waves.

4) 0 (c4)
Varying L4 with respect to ,i, gives
fVh e VAn + f Am A" + fnfA AnAy + nf nAn + t f An An 17 2 1X 2X2I2t 172 ]lz lt
n +1 n
+Tfx1fzAfxA"x + f1finzzA"1A" +'nU fzAnx + 7Uf"Anx +2 2izA
+ ff~-nAx +VfAYn + VfAny +VfAY = 0
For the first harmonic, we have D A'
Vhqc-VhA1+ f1 1D 2 a
VA~cYA~iflzlizA1 A2 +ffA AlJztZX + fl'f 2-Ax22l+f-, + /z' + f1f/+(2A xAx fl + AGA71) + flfIf
I -All 1 1 oDAI +1 2DAr1
+f~ ~ i 1+fz~D 1 DtT 120-)+ ~Z
(2A(Az '7 + A) 1 ) 2 + f1.zg" = 0
DA' g
Dt 2
DA, g jo
= ae ;
Dt 2
D A 2 3 a e i z v 2 3 O2j a 2, 2'P Dt 8
The only unknown in the equation is f2zA so that
f0A = -[fVhc Vh() + g VhcVha + -j a 12] (2.88)
1 01 oF1
in which
F = gk' (3 2 tanh2 kh + 4 sinh2 kh sinh2 kh cosh2 kh)
16sinh' kh
Next, varying L4 with'respect to Ai, also considering first harmonic, we get
fzf22z12 + 1 jflzft71A2 + f(iA1 + 2i71) + fzf 1(A

+2A-I) + -~c~ir1 1 0) .- t i +0 1 ll,)1
-(P'zA x)x2 P11vAxx -Vh (9Vh1 ) f -0 1'Ax)x (A x
--I IlX2X22A
-12 1X 2x -1 2
-q )x f1(A22A + #Aix) Vj~y = 0 where again f'1 = f-1 = 1 is applied. Each term in the above equation has been calculated before. After some tedious algebraic manipulation, we obtain
_ i(Cg V + 02,. C~kl + UOI)X
2(CCk1 + oU)ax2 + 2oalVay2 icc9 V2)ayy +1 [(cck + U )x2
+(-)y,]a + ikcc,k'l a 12a = 0 (2.89)
where k' = k3 c cosh 4kh + 8 2 tanh2 kh
C9 8sinh4kh
kj and al in (2.89) can be replaced by k and a without introducing errors larger than 0(E4). Equation (2.89) contains information of wave refraction and diffraction, as well as wave amplitude dispersion and wave-current interaction. Finally, as we have noted before the first harmonic for L, varying with respect to A1 as expressed by (2.75) is third order and, therefore, should be combined with (2.89) to reach the wave euqation. Hence, we obtain in the (z, y) coordinates
,,ccgkl + U ( V
2(cc,ki + oIU)as + 2oaVay i(cc, V2)avv +of [( + ). + (-)v
1 1
+iccg (k2 k2)]a + ikcc,k'l a 12a = 0 (2.90)
This is the final form of the wave equation without dissipation. If dissipation is to be considered, then a dissipation term D, could be included in the left hand side. One possible formation of De is given in Appendix 1. Following Kirby and Dalrymple (1983), a phaseshifted wave amplitude is introduced to simplify the wave number computation a = Aei(ko-fz kidz) (2.91)
where ko is the reference wave number. Equation (2.91) implies that

A" = a" ei(kx-koz)


where a" is the real wave amplitude. The final wave equation, after dropping the double primes, takes the following form: 2ccokl-Uoa)
2(cckl + oU)A. + 2alVAy i(cc, V')Ave + 2[(ck-U )
+(-), + icc,(k2 k2) + 2i(ccki + ouU)(ko ki)]A + i(4o 02)A
(i 1) gkI
_ kh )A + 2(ccgki + oU)WA = 0 (2.93)
cosh2 kh 23.
where the cubic amplitude term in Eq.(2.90) has been replaced by i(aO0 02)A (Kirby and Dalrymple,(1986)), in which
a = gk(1 + fle2Dk)tanh(kh + fze) it will be discussed in next section. The last term on the left hand side is the decay of energy due to wave breaking and the details of which will be discussed in section 2.2.6, and W is defined as
W = 0 before wave breaking
W = 2(1 )24A after wave breaking, -y = 0.39 and q = 0.4
The dissipation term De is given in appendix A. Equation (2.93) is the final form of the wave equation. It is similar to the wave equation given by Kirby (1984) except the additional imaginary term in the second last term. (see Liu 1986). For plane beach and in the absence of current, linearized (2.90) reduces to the conventional one-dimensional wave energy conservation equation. It reads a = ao(cgoo/co)(bo/b)
where quantities with subscripts 0 refer to known condition and b is the width between adjacent wave orthogonals.
2.2.4 Dispersion Relation
When waves travel over a varying topography or, are embraced by current, the wave number vector will change its value. The four equations established in the previous section

enable us to solve the four unknowns: ij, U, V and H, provided the wave number vector is known. Since wave number is a horizontal vector, in general two equations are required to solve its magnitude and direction. The wave number irrotationality and the wave dispersion relationship derived from dynamic free surface condition are commonli used. In the presence of a current field the mathematics could become complicated if the nonlinear behavior is fully accounted for. In the present study, a semi-nonlinear dispersion relationship is considered and is incorporated directly in the wave equation; thus the linear dispersion relationship is used in the model. Also, by assuming small incident wave angle, the directional information on wave number is incorporated into the complex wave amplitude as presented in the previous section, thus, avoiding the inclusion of the wave number irrotational equation explicitly. The linear dispersion relation with current is given as (Dean and Dalrymple, 1984)
w=a Uk + Vk2 (2.94)
a= Vgk tanh kh
k, = kcosO
k2 = k sin 0
In shallow water, tanh kh --+ kh approximately, Eq.(2.94) becomes w = kx/ + Uki + Vk2 (2.95)
Several modifications have been proposed. Whitham(1967) gave a Stokes nonlinear dispersion relation for intermediate depth region in the absence of current: 0o = oV11 + C2Dk (2.96)
c = ak
cosh 4kh + 8 2 tanh2 kh
8 sinh4 kh

As waves enter shallow water, (2.96) becomes invalid as Dk approaches (kh)-4
In the shallow water region. Hedges (1976) proposed a different formula
O'o = Nfgk tanh k(h + a) (2.97)
for the case of no current in which a is wave amplitude. For the same wave frequency, it yields longer wave length than that given by the linear dispersion relation.
More recently, by combining both cases discussed above, Kirby and Dalrymple (1986) constructed a composite model which satisfies both intermediate water depth, Eq.(2.96), and shallow water depth, Eq.(2.97),
an0 = [gk(1 + fic2Dk) tanh(kh + f2f)]1/2 (2.98)
f tanh5 kh
f2 = (kh/ sinh kh)'
In deep water region, (2.98) approaches (2.96) as f2 -- 0; conversely, in shallow water region, it becomes (2.97) as f, -- 0. Now assuming there is no nonlinear current effect on wave dispersion, the presence of current causes a Doppler shift and Eq.(2.98) becomes w = ato + Uk1 + Vk2 (2.99)
Figures 2.2 and 2.3 shows the results of w for a constant non-dimensional velocity U/co =
-0.1 and varying of E from 0.01 to 0.4, where co = \/k is the deep water wave celerity. The dispersion equations given by (2.94), (2.96) and (2.97) are also plotted in Figs.2.2 and
2.3 for comparison.
The emprical nonlinear dispersion (2.98) can be used to modify the wave equation as we did in Eq.(2.93), in which the Eq. (2.98) is incorporated directly in the governing equation in the form of a term which directly produces the required amplitude dispersion (Kirby and Dalrymple, 1986). Consequently, the linear dispersion relation (2.94) will be used in this model.

E =0.01






1 .00

Figure 2.2: Variations of Dispersion Relation Hedges; - Witham; Composite







with U/co = -0.1 and c, ..... Linear; -.-



o |

CD 0:3

.00 0.50 1.00 1.50 2.00 2.50
o (1"
o.. ............
.00oo 0.50 1.00 1.50 2.00 2.50
o .
.00 0.50 1.00 1.50 2.00 2.50
Figure 2.3: Variations of Dispersion Relation with U/co = -0.1 and <, ..... Linear; Hedges; - Whitham; Composite

Once the dispersion equation (2.99) is brought into the model, the cubic amplitude term in Eq.(2.90) should be dropped, since the amplitude dispersion has already been includeed in the nonlinear dispersion equation (2.98).
A wave number vector which is irrotational can be defined as:
k = Vw (2.100)
where k is a wave number vector, k = {k1, k2} and w corresponds to the scaler phase function. Using the mathematical property that the curl of a gradient is identically zero, it is shown that
V xVW = 0 (2.101)
V x k = 0 (2.102)
In general, Eq.(2.102) together with Eq.(2.99) is necessary to solve the wave number field of a domain.
However, in the present study the wave angle information has already been included in the complex wave amplitude given by Eq.(2.92). Therefore equation (2.102) does not need to be included explicitly. However, the assumption that wave approaching angle is small must be observed.
2.2.5 Radiation Stresses, Lateral Mixing Stresses and Bottom Stresses
In the wave momentum equations (2.21) and (2.22), there are some undetermined terms: radiation stresses, lateral mixing stresses and bottom shear stresses. They are not independent variables and are calculated by the following equations. Radiation Stresses
The radiation stresses are given in (2.19):
_' IpgD 26j piiiaiD
f-h i -h S u3d 2

The form was simplified by Longuet-Higgins and Stewart (1962) to second order for a linear progressive wave system. For normal incident waves, it can be approximately by S S, E(2n 0.5) 0 (2.13)
S-o YoY 0 E(n 0.5)
where E is local wave energy E = pgj a 12, a is complex wave amplitude and n is the ratio of wave group velocity to the wave celerity.
In general, when waves enter beach with an angle, Sqj can be obtained by
[ Si ] = [A ][ Si ] [ A IT (2.104)
where [ A ] is directional matrix, [A] [ cos0 -sin0 (2.105)
sin 0 cos 0
S. SaY E[n(cos2 0 + 1)- ] (Ensin 20
[sj] 1 (2.106)
S., Su, En sin 20 E[n(sin2 0 + 1) ]
where 0 is local wave angle.
Bottom Shear Stresses
The bottom shear stresses are customarily expressed as
FM = pfut Tut I (2.107)
where ut is total velocity composing of both wave orbital and mean current velocities, and uti is its component form. The quantity f is the friction factor. The magnitude of total velocity I ut is equal to Vfu +v2. The u and v velocity components can be expressed as S= U + u = U + uw cos 0 (2.108)

v = V + ug = U + u' sin 0


where U and V are mean current components as defined before. u' and u' are wave orbital velocities in the x and y direction, respectively and uw = [(u')2 + (uw)2]1/2. The total velocity I ut I can be expressed as ut 1= [U2 + V2 + (u,)2 + 2Uu cos 0 + 2Vu- sin 0]1/2 (2.110)
The local wave orbital velocity u' is purely oscillating in time uW = umcosot (2.111)
where umn is the maximum wave orbital velocity at the bottom um = ia (2.112)
sinh kh
The time averaged bottom shear stresses in the z and y directions then become b- = P (U+umcosO0cosot) j ut j dt (2.113)
T 0T
b = p (V + um sin 0 cos rt) I ut dt (2.114)
Alternative, for computational convenience, it can be rewritten as = fp 2(U + um cos 0 cos t) I ut I d(at) (2.115)
- = p-L f2(U + um sin 0 cos t)I ut d(at) (2.116)
These integrals can be calculated by applying Simpson's rule. The serious problem by applied the above bottom friction model is the computational time. A simplified model is introduced with
7= pfI u' I Us (2.117)

I u't 1= V/U2 + V + urn

since there is no interaction term between wave orbital velocity and mean current velocity. In the component form the above equation gives
Tb- = P1 Iu't U (2.118)
= pf U'tIV (2.119)
The difference between the two models (2.107) and (2.117) is minor for cases U >> un as shown by the example in Fig.2.4, in which the same bottom friction coefficient is used (the testing conditions refer to Fig.5.3 in chapter 5), but the saving in computational time is more than 60%. In fact the small difference can be offset by increasing the bottom friction coefficient in (2.117) to account for the absence of urn term. Equation (2.117) will be used in our model.
Another very important factor is selection of the friction coefficient f. Based on field observations and laboratory experiments Longuet-Higgins (1970) suggested a value about 0.01 for long-shore current calculations. The bottom friction coefficient for waves alone presents a more difficult situation. For smooth turbulent flow over fixed bed, Kamphuis(1975) gave f = 1 (2.120)
where Re is Reynolds number, R, = la w and v is the kinematic viscosity.
For rough turbulent flow, friction coefficient no longer depends upon Reynolds number. Investigations were conducted by Jonsson(1966)and Grant & Madsen(1979). The most commonly used relation is
f lexp[5.213(')0'4 5.977] < 0.63
f = (2.121)
0.15 r > 0.63
where r is the roughness of the bed and a' =1 a 1. Equation (2.121) is an approximate explicit solution given by Swart (1974) to an implicit, semi-empirical formula by Jonsson and Carlsen (1976).

o 0 (o-

0 DZo

(' -iI I I
0.00 '4.00 8.00 1
T (M)

Y (M)




0.00 4.00 8.00 12
Y (M)

a Lo
an D (10



q. 00
Y (M)

Figure 2.4: Velocity Profiles for different friction models (2.107); - (2.117)

I=1 5



. 00

C r%







The situation becomes further complicated once combined wave with current, the state of the art of modeling such flows is rather rudimentary. The total velocity in (2.107) becomes much larger than the wave orbital velocity in the present case of very strong current. The actual values of f used in the model will be adjusted as a function of (a; T) based on experimental results to be discussed in chapter 4.
In application, the bottom friction coefficient may be treated differently in x and y directions, in the present study, the bottom friction coefficient in y direction is taken as one half of the value in x (or main) direction.
Lateral Mixing
Lateral mixing due to turbulent diffusion is a complicated flow phenomenon. Rigorous mathematical treament is beyond the state of art. Figure 2.5 shows the effects of lateral mixing on longshore current distribution (Longuet-Higgins, 1970 a.b). The common assumption is that turbulent mixing is proportional to the local mean velocity gradient. Following Longuet-Higgins' formulation on longshore current mixing, the lateral shear stress is written as
= au av (2.122)
where c. and cy are the horizontal eddy viscosities in the z and y direction respectively. We define
i = i + c i = Xy (2.123)
where ei and eci are the eddy viscosities of wave and current, respectively which are assumed to be independent. The inclusion of Eci is based on the fact that even though there is no wave, the flow from the inlet still generates turbulence. As the simplest possible approach, Longuet-Higgins proposed a model for cw, that is proportional to the distance offshore, X, multiplied by velocity scale

=WZ = N(v22


~Without Mixing LU ~With Mxn
Breaker Line
Figure 2.5: Longuet-Higgen's Longshore Current Profiles were N is a dimensionless constant whose limits Longuet-Higgins gave as
0 0 __ 0.016 (2.125)
Following Ebersole and Dalrymple (1979), N is chosen as 0.01 and c.. is allowed to vary linearly with z to the breaking line. From this point seaward the coeffcient remains at this constant value. There is no elaborate consideration of c... Ebersole and Dalrymple assumed a constant value u it is chosen as half the value of throughout the field.
The eddy viscosity theoretically, will change from point to point. One of the formulas was due to Von K~rma'n's similarity hypothesis: C0 "- MN 01 (2.126)
CC 2 -Y 21 dU (2.127)
FolwigEbrol ndDlrml (99) i hoe a .1 n e s lowdtovr

dV d2V dU d2 U
ml = k-x- /-c--2 and 2 = k-/and k = 0.4 is Von Khrmhn's constant. In the numerical calculation, these coefficients go to infinit near a point of inflection. An upper-limit value should be imposed. Alternatively, both eddy viscosities can be treated as constants. Figure 2.6 give the effects of eddy viscosities on the velocity profiles (the testing conditions refer to Fig.5.3 in chapter 5). It shows that they have only minor difference between constant 0.001 m2/sec and Von Krmhn's mixing model with upper-limit equal to 0.001 m2/8ec. Therefore in this investigation, instead of Eqs.(2.126) and (2.127) the constant eddy viscosity is selected.
2.2.6 Wave Breaking Criterion
When waves enter shallow water, they will eventually become unstable and break, dissipating the energy in the form of turbulence and work against bottom friction. For a very mildly sloping beach, spilling breaker occurs and the simplest criterion predicted by solitary wave theory is (McCowan (1894))
y= (-)b = 0.78 (2.128)
where D is total water depth and subscript b denotes the value at breaking.
Dally et al. (1984) proposed that the decay of energy flux with distance in the surf zone is proportional to the excess of energy flux above a stable value, given the relation ,,(Ec.) q -D [Ec (Ec,).] (2.129)
where q is a constant related to the rate of energy decay. The quantity (Ecg), is the stable energy flux for a broken wave in water of total depth D. The last term in the wave equation (2.93) which appears only after wave breaking

2(cckl + oriU)WA

C- 1=31
'000 0.80 1.60 2.40 3.20 4.00
Y (M)
0 1=1 6
'0.00 0.80 11.60 2 .40 3.20 4.00
Y (M)
Figui'e 2.6: Variations of Eddy Vicosity Coefficients. = 0. ; ... =0.001; - =0.01; -
Von Khrm&n's Model with upper-limit =0.001. unit= m2/e

can now be completely defined as W =q 1 2D2
W ( 41A 12 (2.130)
by virtue of Eqs. (2.128) and (2.129). The value of q is chosen as 0.4 following Kirby and Dalrymple (1986).
2.3 Summary
In this chapter, the analytical model of a current-wave-topographic interaction system is formulated under the assumption of strong current and slowly varying bottom. Basic equations together with essential terms are developed. They are summarized as follows: Depth-Integrated Continuity Equation 8Bt BUD 8VD
"- + D + ---- 0 (2.131)
at ax ay
Depth-Integrated Momentum Equation in X-Direction
aU aU vaU aBy 1 asz 1 asz 18aT 1
-T + U + V -g DDa + Y (2.132)
t ax ay az pD az pD ay p -y pD
Depth-Integrated Momentum Equation in Y-Direction
8V 8 V 8 V = Byf 1 BS. 1 BSv 1 8, 1
av + U + V a= -g a. I + 1 (2.133)
-t+ x -z y gy pD ax pD ay p 8 p- (2.133)
in which we have substituted (2.131) into (2.21) and (2.22), where
S, Sf, E[n(cOs 0 + 1)- } (Ensin 20
Si] (2.134)
Sov Su Ensin20 E[n(sin2 0 + 1)
Fe = pE au- + PEav (2.135)
ay az
ei = ei + Cci; ci is a constant and e, is given by e = Nzv-D with N = 0.01;
Cty = 2

Y = pf 'I Uj (2.136)
,, = /U2 + V2 + U. (2.137)
um sinh kh
Wave Equation
i(CCO V2 + 02,- C'gkl +U') 2(ccgki + olU)As + 2a'VA, i(cc V2)Ay, + o(ck + U) a?
+(V), + icq(kl k') + 2i(cckl + alU)(ko ks)]A + i(eo a')A
+( V)y + iccg(k 2 k 2) + 2i(cc~kl + oi U) (ko ki)] A+ i(o2 a2 u)A
(i 1)gk2
- ) -k V-)A + 2(ccgk + o'IU)WA= 0 (2.138)
cosh2 kh 2w
o = [gk(1 + fle2Dk) tanh(kh + fze)11/2
fA = tanhs kh
f2 = (kh/sinhkh)4
Dk = cosh 4kh+ 8 2 tanh kh
8 sinh' kh
=W = 0 before wave breaking (2.139)
2(1 ) after wave breaking, y = 0.39 and q = 0.4 (2.139)
When waves enter the surf zone, the last term in (2.138) will be brought into calculation till waves reach stable condition.
The dispersion relation used in this model is
W = 0+ Uk i+Vk2 (2.140)

and the wave angle is included in the complex wave amplitude A.

In this chapter, the numerical schemes for solving the continuity, momentum and wave equations are presented, along with boundary and initial conditions.
The numerical scheme requires (1) the definition of a grid system with a system which provides a systematic method for identifying variables of interest and (2) the conversion of the governing equations into finite difference forms.
3.1 Grid System and Definition of Variables A rectangular grid mesh was established over the area of interest as shown in Fig.3.1, where x and y denote the onshore and longshore direction, respectively. Two kinds of variables are employed in the computation: the grid-center variable and the grid-side variable. Practically all the quantities of interest are represented by the grid-center form which expresses the value at the center of the grid and has the following notation: X1.
where X is the grid-center variable, the subscripts i, j represent the ith grid side in the x-direction and the jth grid side in the y-direction, respectively, and the superscript n represents the nth time step. The grid-side variable expresses the value of the quantity along any grid side. The only variable that appeared in this form in the present study is the velocity vector; it is denoted by a subscript s to differentiate it from the grid-center value and has the following form:

Y is the vector quantity along a grid side: the subscripts are such that i denotes the ith grid point in the x-direction on the left of the vector and j denotes the jth grid point in the y-direction on the right of the vector.
In the present study, the grid-side velocity vector was applied in the continuity and momentum equations and a grid-center velocity vector was applied in the wave equation. The definitions of grid-side and grid-center variables are shown in Fig.3.2
3.2 Finite Difference Scheme for Combined Continuity- Momentum Equation
The Crank-Nicolson finite difference scheme in double sweep form is applied here to solve the continuity and momentum equations and later also to solve the wave equation.
3.2.1 Linearized Implicit Model
A variety of numerical schemes are presently available to solve the equations for nearshore hydrodynamics similar to those developed in Chapter 2. Ebersole and Dalrymple (1979) presented an explicit model for nearshore circulation which requires small time steps due to the stability criterion, thus relatively long computational time. Vemulakonda (1982) presented an implicit model allowing for much larger time steps. Both models adopt a three time level schemes. Butler (1981) and Sheng(1981)proposed a two time level implicit model. This method not only saves computational time but also shortens the development of the difference scheme for the governing equations. To demonstrate the structure of their model, the simplified linearized equations are discussed here. The depth averaged linearized continuity and momentum equations are
at + Q + o-g = 0 (3.1)
at ax ay(31
axL (3.2)

I Y!
- -4- -- 4!

n+1 n

2 1

Figure 3.1: Grid Scheme

M M-1





Figure 3.2: Definition of X and (Y)



a = -=gDa (3.3)
at ay
where Q, = UD and Qy = VD. U and V in this section denote the grid-side values. In matrix form, linearized equations can be rewritten as Wt + AW. + BW, = 0 (3.4)
W= Q. A= gD 0 0 B= 0 0 0
Q 0 0 0 gD 0 0
Sheng (1981) implemented a two time level fully implicit scheme in computing the external model in his three dimensional mode.
1 +1 -q"
,,j + I r(Q-) n+l ti 1 nl ~
Z,-, ( )(Q) ] + Y [(Q*)i+ (Q)j 0 (3.5)
" =i~i (Q =i + gD r-~ --~
- Qt + As -ni -1i] = 0 (3.6)
(,) j (Q), gD n+ --n+
At + y i 0 (3.7)
or in the matrix form
t (Wn+ Wn) + ( z 6 + 6Y)Wn+l 0 (3.8)
where S and 8u are central spatial difference operators which provide the differentiated quantities at the same points as the time difference quantities. Equation(3.8) involves only the (n+1) time level in the last two terms, i.e., the whole equation is staggered in time.
Rearranging (3.8) by time levels, we obtain
(1 + 2f8 + 2)Wn+1 = Wn (3.9)
At At
fz- 2AzA.; fl = -2yBS
26: 2Ay

By adding the quantity of 4zflzy(Wn+l Wn), which has a second order truncation error as
(At)2 AB x
we have
(1 + 2P,)(1 + 2py)W""+ = (1 + 4pzp,)Wn (3.10)
This factorized finite difference equation is still a second order approximation to the differential equation (3.4). The advantage of using (3.10) lies in the fact that the solution of the factorized form can be split into two separate one dimensional operations. To achieve this, we introduce an intermediate value, W*, such that (1 + 2).,)W* = (1 2py)W (3.11)
(1 + 2pu)W"n+ = W* + 2pyW (3.12)
From now on, quantities without the time level denoted are the values at time level n. Equation (3.10) is recoverable from (3.11) and (3.12) by eliminating W*.
A double sweep method is used to solve (3.11) and (3.12). At the first sweep, the values at time level n are known we solve for W*. The full solution W"+1 then is obtained based upon known W* values. After some manipulations, the first sweep in x direction gives S- + -z Q + -6vQY = 0 (3.13)
Q+ Q. + =At
{ Qu+ Q + gD6.V = 0 and the second sweep in the y direction yields ry+1 T* + b~-,Qn AtQ =0 (3.14)
Qn+1 Qu + A ggD6ur+ = 0 On the first sweep, Un+1 and (* are found, and on next sweep the values of Vn+1 and r+1 are obtained. Therefore, after one loop, the quantities at time level n + 1 are known. The

advantages of this scheme are (a) the finite difference equations are brief, at second sweep (3.14) there is no Q,, involved; (b) a significant savings in computation time is realized with large time steps and less computations.
3.2.2 Numerical Scheme for the Full Equations
Following the above model, the full equations (2.131), (2.132) and (2.133) in chapter 2 can be written in numerical form as (see Appendix C) x-sweep
t z(Un+1D) + St y (VD) = 0
(1+ Atc)Un+l U + g6j* + Atm = 0(3.15)
- + (V+1D) (VD) = 0
(1+ AtcY)V"+1 V + A-gY86,"+ + Atm2 = 0 where c.U and cuV are bottom shear stresses in x and y direction, respectively. By (2.136) in chapter 2, we have
f f
cZ = Ut I: and c.= = t| u|, where I u I is the total velocity defined in (2.137) and here is specified in each direction. ml and m2 are the remaining terms in Eqs. (2.132) and (2.133).
Again, on the x sweep the solutions for i* and U"+1 are obtained; these up-dated values are brought into the y sweep to obtain i-+1 and Vn+'. Writing the above two set of equations in detail, we have a) X-Sweep
Continuity equation
1 1
-At( ,i ii) 2Azx U,+ -(Di+1,y + D1,) Ui '(Di,y + Dj- ,)]
+ 1[V,+(D,+ + V(D + = 0 (3.17)
2Ay- [TVi,i+l(Di,i+1 +r Di,i) V4,y( Did +[ Didj-1)] = 0 (3.17)

X-Direction Momentum Equation
I (u"f.- U,.) + (U,+l, U ,,) + (V i+1 + V + V-I,+I V-1.)
At Ax 8Ay
g _2 (Uj+1 Ui,-1) = (iv, -*,) (Di + Di-) [(s ,)ii (S)Ax p(Di,j + Dj-1,1) Ax
1 1 S
p(Did Di-1,j) 2Ay [(Sy)ij+l (S2y)i,-1 + (Sy)i-lj+l (Sy)i-lj-1
(Di, + D+i-1,) (u1)ii U + (FLy)id + (FLys),i (3.18)

(ut)ii j = [(Um)i,i + (um)i-1,j + 2(ut)i,j um is wave orbital velocity as defined in (1.112) in chapter 2.
(,,),, = {Ulj + [0.25( .i + J+ + Vi, +-1,+1 2l

(FLyy)i,i = (FLuz)ij =

4(Ay)2 {(U4+ Uij[((eY)i,, + (eW,)d+1 + () + (ey)-1+) +2((cv)Uj + (if),+c)] (U,, U E,w-)[((WY),d + ( iJ-1 +(Cmy)i-,i + (esy)i-1.i-1) + 2((cey)ij + ( gMy),.i-)]}
+ (fz~~ + 2((ex~)- +
4Azay {(Vi~i+1 V-1.i+1)[((ews).i+ 1+ ( ),i + (es)i-1d +(cw.)i.-1,+,) + 2((cz)i,i+1 + (cz)i-1.+1)] (Vi,i Vi-1.i)[((ews)i, +(Cw)i-,i + (ews);,-1 + (ews)--11) + 2((en=),i + (cs)i-l)]

b) Y-Sweep
Continuity Equation
(" saly ,i) + Vni (id+ 1Di+ + Did)* Vni ",,( + Dii-1)*
2Ay[Vj+(Dii+ + Dij). Vi,i(Di,i + Did-1).] = 0 (3.19)
Y-Direction Momentum Equation (U(+ l + U U 8 + 1 i + ) (- ,.+ U+ + ,1 J+,j(t -V,,) 8A (V+id V-,d) +

g _-n+1 -n+1 p(~ 1 ~ 1S (V,i+1 Vij-,1) = t(,f +- p(Did Did-1)* 2 (SmY)i+1,y
-AY p(DSj Dli,,1
-(Say),-xy + (S.Y),+1,y-1 (S.A)-,-11 ]- C ((Sy)i,
p(Diji + Did-1)* y
f y+
-( )i-1)- 2(Di+ Di-).I (ut)i,i jVi' + (FL.y)i,i + (FLz.),j '(3.20)
2(D,,3 + Di .
where the bottom shear stress coffficient has taken as half the value of that in x the direction as stated in (2.2.5), and
I (ut)ij 1, = [(um),,i + (um),i-,]] + 2(u,)i,j
(tz),> = v{ + [0.25(U, + + U1" + ,,I1 + U, z.j]2}1/2
(FLzz)id = 4(Az)2 {(V, Vi ew)+,+ ( ),i+ (Ws);l,-1 + (e)+)i,(-1)
+2((e),+1,i + (e))]- (V,1 Vi-l4 ,)[((eCw);, + (Cw,),-1,j
+(ew.)i,- + ( )~-xi-1) + 2((e~.),, + (s)-)
( =4A1A+
+(E+ 2((y)+, + (c)i+,-1 -+ ( )
+(c )i,, + (cwy)i-j + (cwy)i-,1-.) + 2((c.E)i, + (C'Y)i, _1)]}
In (3.17) to (3.20), the quantities with no time level are at time level n; and the subscript
(.) outside the bracket in(3.19) and (3.20) denotes up-dated values, i.e., D. = h + V*. Also we have used the latest x direction velocity, Un+l, in the y sweep to accelerate the convergence. In application, the (en) and (eqy) are considered as constants as dicussed in chapter 2.
3.2.3 Method of Calculation and Boundary Conditions
A double sweep method has been developed in each direction,i.e., in both x-sweep and y-sweep. The advantage of this sweep method lies in the fact that it gives an explicit formula at each sweep to calculate the results.

Collecting the terms from (3.17) to (3.20) by time level, so that the unknown quantities are to the left and the known qantities are to the right, we have X-Sweep
iT ..r prn+l. p n+l
- -U{- + Ri ,Ri'i -1, ijU; = c -(3.21)
A ,"f + RI V'+Ij R1 (.j = Bij Y-Sweep
--q1 n+l n+l = -1 =
'+ QVAy z, ,(3.22) ,.n+ R2 ,1 R2-- +1 R, g V". + R2 rI- R2( i = Spy where the quatities without superscripts are at time level n. The definitions of undefined quantities are given in Appendix D, which also shows the details of the following derivation.
The x sweep equations, after some manipulations, can be reorganized as ( ,i Fli + E1Ug32
Un+i =- F2i_1 + E2-li-Ij where
=Ci Pi,F2
S 1 + P,y E2i Eli = Pi-1,y
1 + Pi jE2i
F Bi+1j R1Fli+1 F24 =
Ai+ij + R1Eli+1 R1
E2i =
Ai+1j + R1Eli+1
The boundary conditions state that at the shoreline (i=idry), the velocities should be equal to given velocities in the inlet section and they disappear elsewhere, and the mean water displacement V, assumes to be identically equal zero. Therefore, the first sweep will solve coefficients Fl, F2, El and E2 starting from i=idry till offshore line (i=1). Once those coefficients are determined, the second sweep will directly yield results of #* and Un+1 starting from i=1 (offshore boundary) till i=iwet, where iwet = idry-1, is the last grid at

which water depth has non-zero positive value. The offshore boundary condition is given at deep water, i=l, by specifying the mean water displacement to be zero. This is approximate assumption to maintain a steady state mean flow in the region of concern with no increase or decrease of water mass. A more rigorous condition is to allow for water level variation at the offshore boundary such that the change in hydrostatic force is balanced by the net momentum flux and bottom shears in the x-direction.
The lateral boundary conditions basically have three kinds for the concerned region: the symmetric free boundary conditions; the no-flow boundary conditions and the nonsymmetric free boundary conditions. They will be specified for specific problems later.
The y sweep is similar to the x sweep and (3.22) can be rewritten as =F3 +E3V.n+l (3.24)
=n~1 F4I. + E4-ljli,, 1
- Zij QijF4j
1 + Qi,E4j
E3j = Q____1 + Q,E4
-F ,i+ R2F3j+I
F4j= R,j+l + R2E3i+l
E4j = R2
Rj+1 + R2E3j+l
The first sweep will determine F3, E3, F4 and E4, then the following sweep will give the values of r+1 and Vn+1. The shoreline boundary condition for V is set equal to zero at i=idry; and offshore boundary condition at i=1 supposes that it satisfies the mass conservation, i.e., (VD),.. = 0. The lateral boundary conditions will be discussed for different situations later.

3.2.4 Stability Criteria
The stability criteria of these schemes can not be established theoretically. In general, the stability criterion for explicit scheme can be guided by the following inequality: A/t (A )2+(A,)2 (3.25)
1t_< I ,+ lgv-D
where I ut I is the magnitude of total velocity. For application, it becomes AtMin.[ (3.26)
for jut << 'D
Actually, the critical time step, At, in our study by using a two time level implicit model is at least one order larger than that given by (3.26).
3.3 Finite Difference Scheme for the Wave Equation
The final form of wave equation (2.138) in chapter 2 is mixed parabolic-hyperbolic partial differential equation. The advantage of parabolic type partial differential equation is that the marching method can be used in step by step calculation. Booij (1981) proposed a parabolic model following on Radder (1979). Then the wave field can be split into a transmitted and a reflected field. In his parabolic approach, the reflected field is completely neglected. Kirby (1983) gave another approximation, after assuming that the complex amplitude can absorb the difference between real and image wave directions for small angle incident waves, and also assuming that wave amplitudes are changing more rapidly in y direction than in x direction for slowly varying topography. Both of those assumptions have been introduced in our derivation of the wave equation.
For the parabolic partial differential equation, the Crank-Nicolson finite difference scheme is applied. For simplicity, (2.138) in chapter 2 can be rewritten as
A.,+ oA -iCiAYV + I i21aC3
+ --A + A + (C4 + iCs)A = 0 (3.27)
C2 2a- C -y

Co 2V
cc V2
CP =
C2 = _:-2Os =
C3 =
C4 = R+W
Cs = I cc,(k k2)] + (ko k) R+ 1(o ,2)
in which
P = 2(ccgkl +oiU)
R gk2
Pcosh2 kh 2w
As the waves enter the surf zone, the second term in C4 will be added into calculation. The
quantity of j A 12 in a2 is the value evaluated at last iteration, unlike Kirby who specified
this value at the same iteration. The values of velocity components U and V are grid-center
values in this calculation.
The Crank-Nicolson finite difference scheme constructs Eq.(3.27) as
A; A 1 C
A+,x ,i + 4y [(Co)i+j(A+i+j+ Ai+1.i-1) + (Co)i,
(Aii+1 -Aij-1)] [(C )+1,i(A4.+j+1 2A,+1,j + A,+1,j-1) + (Cl)j,j
1 ((C2),+l,, (C21,) (A+, (Ai~+1 2Ai + Ag-I)] + 2A ((C2)(+,j + (C2)) (+I + Aj) 2Az ((C2)i+1; + (C2)idj)
1 ((C3)+1,i+1 (C3)i+1,.i-1)A ((Cs),,i+1 (C3),i-1)
4Ay (C2)i+ ,j (C2)i,i
+1 [((C4),+l, + i(Cs),+1i)A.+I,1 + ((C4),. + i(C5),i)] = 0 (3.28)
A similar double sweep method is used to solve (3.28), resulting in
Aj = R2-.1 + R1, Ah4j-1 (3.29)
The definitions of R1, R2 and other details are also given in Appendix C.

The boundary conditions are stated as follows:
a) at offshore, i=1, the wave amplitudes are specified in complex forms; b) at shoreline, i=idry, the wave amplitudes are supposed to die off completely; c) the lateral boundary conditions can be considered as Ay = iksin OA (3.30)
which assumes that the waves travel mainly in the x direction on a slowly varying topography and Snell's law is valid for first order correction. The finite difference form of (3.30) gives Aij+1 AijI + Q (3.31)
Q = (iksinO A-)i,j+1
The local wave propagation angle 0, as discussed in chapter 2, can be found through the complex wave amplitude. By Eq.(2.92) of chapter 2, we have In()A'a i(ki ko)Ax (3.32)
In( A,+1) -- i2k2AY (3.33)
where the double primes in Eq.(2.92) have been dropped. The local wave angles are included both in k1 and k2, so that either of above equations can be used to obtain 0.
3.4 Method of Solution
With the specified initial and boundary conditions, the finite difference equations (3.21), (3.22) and (3.28) can be solved in a loop of iteration.
The initial condition is assumed to be the state of rest. The velocity field, both U and V components as well as the mean free surface displacement, 'q, are initially set equal to zero. A bottom topography is input into the model. The initial wave amplitude field is calculated

by solving the wave energy equation with velocities equal to zero, which corresponds to Snell's law for plane beach before waves break.
To avoid "shock loading" the model, both wave amplitude and velocities of inlet discharge are built up from zero to their termination value over a number of iterations. We choose 100 iterations to set up those values by using the hyperbolic tangent function
F = Fo tanh (IT/25) (3.34)
where F is either wave amplitude or velocity of inlet discharge; F0 is their termination values; IT is the number of iteration. The approach declares that the termination value will be reached as total iteration is up to 100, while tanh 4 approaches the value of unity.
The program will run till the variables U, V, !7 and A reach their steady state. The convergency is also checked by the balance of total discharge in and out of the controlled boundary.
There are two coefficients to be determined in a running case: one is the bottom friction coefficient and the other is the eddy viscosity coefficient. For fixed inlet discharge, the eddy viscosity is assumed to be constant, but the bottom friction coefficient may vary for different combinations of wave and current. Further explanations on how to select these two coefficients are given in the next two chapters.

To evaluate the performance of the numerical model, a number of tests were run, mainly utilizing existing numerical models of simple cases as bench marks. Little attempt was made to compare directly with existing experimental data for two reasons:
1(1) Experimental data that can be used to test the model to its full term, i.e., currentwave interaction due to inlet on sloping beach, does not exist.
(2) The bench mark numerical models selected here usually already did their comparisons with laboratory data, whenever available. Since the purpose of the present model is not to seek improvement of any specific numerical model but rather to extend the capabilities to new applications, it is deemed sufficient for verification purpose to compare with numerical results only.
4.1 Wave Setup on Plane Beach
The model was run for a case of normal incident wave on a plane smooth beach. The conditions were as follows: T=1.14 sec, deep water wave height Ho = 6.45 cm, and beach slope = 0.082. This was the same case used by Ebserole and Dalrymple (1979) to compare with laboratory data by Bowen, Inman and Simmons (1968). Since it is a two dimensional case, we have selected m = 52 in the x- direction (on-offshore) with Ax = 0.25m and n = 4 in the y-direction (longshore) with Ay = 0.4m. The water depth at the deep water end in the numerical model was about 1 m and the mean water level displacement due to set down can be neglected. The results are shown in Fig.4.1. The wave height was computed on the basis of linear wave shoaling and the wave equation as expressed by (2.138) was used for wave height computation. The data from Bowen, Inman and Simmons (1968) were also plotted for comparisons. Apparently, these two approaches are similar except around 61

0 0
r) e
0 A
* IJ
0 0
0.00 80.00 160.00 240.00 320.00 400.00 480.00
X (CM)
Figure 4.1: Set-up on a Plane Beach. Eq.(2.138); - Linear Wave; A Data

breaking line.
4.2 Open Lateral Boundary Conditions and Longshore Current Distributions on a Plane Beach
Two typical open boundary conditions can be specified for the present numerical model. One of them is the periodic lateral boundary conditions, which requires that Q(i, 1)-Q(i, N+ 1)
Q(i, 2) =Q(,N + 2)
Q(i, 3) =Q(i,N + 3)
and so forth, where Q is any quantity of interest. This boundary condition can obviously be applied to topographies with periodicity such as the example given by Ebsersole and Dalrymple (1979), but can also be relaxed to apply to an arbitrary topography that has no periodicity as long as the boundaries are chosen far enough away from the area of interest and there is no source or sink in the area.
In the case where a source or sink exists, such as the present study with an inlet, the periodic lateral boundary condition can no longer be applied as the out-flow from downstream would not be balanced by the in-flow from the upstream. An approximate free boundary condition in this case is to let
at far enough distance on both upstream and downstream. The above boundary condition implies no longshore variations far away from the area of interest.
A typical example for open-boundary-condition computation is the longshore current generated by obliquely incident waves. The problem was studied analytically by LonguetHiggins (1970). The case illustrated here is for the conditions: beach slope = 0.075; incident wave angle = 50; bottom friction coefficient = 0.01; wave height = 0.02 m and wave period = 1 second. The results are shown in Fig.4.2.


If there is no lateral mixing effects at all, the velocity decreases rapidly from maximum to zero at breakline. The apparent existence of velocities outside the breaker zone is due to the fact that the advective acceleration terms in the differenced y momentum equation cause velocities outside the surf zone to be contaminated by velocities inside the breaking zone. The strong mixing effects are demonstrated in the other two cases: one only has the lateral mixing due to breaking with a mixing coefficient = 0.01 m2/sec and the other includes both mixing from breaking waves with coefficient= 0.01 m2/sec and flow turbulence with an eddy viscosity coefficient = 0.005 m2/sec.
4.3 Nearshore Circulation
In this example, currents induced by waves over irregular bottom topography are computed. The bottom topography is slightly modified from that used by Noda et al. (1974) and Ebersole and Dalrymple (1979) and is generated by the following equation: h= a. x[ 1 + Ae-(z120) sin1 0(y Vx/tanf3)] (4.1)
s: beach slope = 0.025
x: distance from shoreline
A: periodic beach length = 104 m
A: amplitude of bottom variation = 10 m f8: = 30*
in which Vf is used instead of x on the right hand side as in the original equation by Noda. This modification is necessary to reduce the variations on the lateral boundary conditions. We choose m=40 in x-direction with Ax = 5.0 m and n=27 in y-direction with Ay = 4.0 m. The other parameters are given as
T = 4.Osec.; ao = 0.75m; 0 = 50; f = 0.012
breaking mixing coefficient = 0.01m2/sec.

a :
S- -. %
10 ----2. 20 0.00 2.20 4. 40 6.60 8.80
Figure 4.2: Lonshore Current Generated by Obliquely Incident Waves. -- no mixing; breaking mixing coeff.=O.01m2/s; - breaking mixing coeff. =0.01m2/s and eddy vis. coeff.--O.005m2/B

The resulting velocity field is shown in Fig.4.3. Although this case offers no direct comparison with the result by Ebersole and Dalrymple, the flow patterns are similar and the current field is physically plausible.
4.4 Combined Refraction, Diffraction and Current Interaction Case
This example deals with the wave refraction problem as illustrated in the previous three cases, but also the combined effects of refraction and diffraction under the influence of topography and the interaction with current. A case of a submerged circular shoal with a parabolic configuration resting on a beach is studied. Similar topography is commonly seen around inlet where an ebb shoal is located outside an inlet. This idealized topography was first investigated by Radder (1979) and then by Kirby and Dalrymple (1983) on a fiat bottom with normal incident waves. Therefore, the numerical model was first run at the same condition used by Radder. Then, the case of obliquely incident waves was considered with the shoal on a plane beach. Finally, the model was treated with the presence of an inlet to the above configuration.
4.4.1 The Circular Shoal on a Flat Bottom
Radder (1979) suggested two configurations of shoals in his study of parabolic wave equation. Configuration I, of a circular shoal on a fiat bottom is represented by the depth profile
{ hhm+er2 r h=h+ r> h =ho r >R
r 2 = (X X.)I + (y_ V
eo = (ho hm)/R2
xm and ytm are the coordinates of the center of the shoal, and hm is the water depth at that point; h0 is the water depth at the flat bottom area. The portion of r < R prescribes a

.t . . . . 1 0 I t 1 t I i I I I t tIf f I I I
..........., tt\\ t t,,t it illittI I t
~~. . . . . . a t t l it\ \ 1t 1I t ,/ t tl l i f fft tt
~~. . . . . . s 1 1 1 1\ t r\ t t l f i t i t t i f
O. .. . . . .1 1 1 1 ,1 / t /l/tI t
~~. . . . . . . \1 1 1 1 1 1 % 1 1 1 1 // / I I1 11 1 t t i
" . . . . ., . 1 1\ \ \ \\xx \ It ti t t it I I t e t t
.. .... . ., . . \ \ \ \ \ \ i t t t ; ;,
o. ....... ,..., ...... \ \\\\\\\\\ lt lt ire.
. .... ...... ..11
. . . . . .. . \ \ \ \ \ \ \ \ \ \ \ \
Ix .. . . . . . , . 1/ / 1 1 1 1 \ \ \ \ \ 1
o. . . . . . ,t tl / //,, .I / / I f/ / / / t fit
,a- . . . . . . t i t I Iti t l / / / / / / / ti f fl r
S . . . . . . \, ,\ 1 1 t i t f t / / / //I t / / ,t t s , .
. . ..W. 1 1 1 1 t / / / / t


parabolic curve. The numerical values selected are
h,,,/R = 0.0625; ho/R = 0.1875; Lo/R = 0.5
where L0 is the wave length in the flat bottom area.
For such a configuration, the wave rays of a linear refraction approximation are focus into cusped caustics. Peregrine (1983) and Kirby and Dalrymple (1983) discussed this case and have shown that two jumps in wave condition fan out in approximately the same manner as the caustics emanating from the region of the cusp. The region between the jump conditions, directly in the lee of the submerged shoal, would then be dominated by waves of approximately uniform amplitude propagating in the incident wave direction.
The grid spacings used by Radder are described as
Ay/A= 0.5; and Lo/Ax = 1,2,4 and 8. and it is noted here that the larger the value of Lo/Az, the better the resolution. For the shoal studied here, L0 = 0.25m is used, the corresponding Lo/Ay = 2.5. The center of the shoal, with the radius R = 1OAz(5Ay) is located along the line of symmetry with respect to the y-axis.
Without consideration of current, the wave energy equation becomes
2ccgkiA, iccgAyV + [ (kiccg). + iccg(k2 k2) + i2(kiccg)(ko ki) ]A
-c0,h ,2 (i 1)gk )A + (2ccgkl)WA 0 (4.3)
cosh 2 kh V(2w
Following Kirby and Dalrymple (1983), the non-dimensional incident wave amplitudes are chosen as k0Ao = 0.16 and 0.32 in which ko = 21r/Lo. The variations of wave amplitudes along the center line (y/R = 0) in the x direction are shown in Fig.4.4 in which x/R = 0 starts from the fore-edge of the shoal with the center of shoal located at z/R = 1.0. For both cases, in the neighborhood of the center of the shoal the wave amplitudes drop significantly below the incident wave amplitude. The peak value of koA0 = 0.32 is less than half the value of k0A0 = 0.16; it indicates the effects of nonlinearity with the increase of wave steepness.

Figure 4.4: Wave Amplitude Variations along Centerline Aoko = 0.16; --- Aoko = 0.32

In other words, the wave energy is much easier to be trapped and harder to decay for the shorter waves at the line of symmetry. Also, there is a slight phase shift in the position of maximum wave amplitudes. Figures 4.5 and 4.6 gives the wave amplitude variations in the longshore direction at the location of z/R = 3.0, 4.0, 5.0 and 6.0, respectively. y/R 0 is the symmetric center line of the shoal. These plots show that as waves travel downstream from the shoal the influence zone becomes wider. The spreading of the influence zone for both cases are almost the same. Since wave amplitudes for koAo = 0.16 have a large amplification at y/R = 0, which can also be shown in Fig.4.4, they drop rapidly to below the incident wave amplitude to keep the wave energy balance.
These results are similar to those described by Kirby and Dalrymple (1983).
4.4.2 Circular Symmetric Shoal on Plane Beach
In the case of the circular shoal located at the center of a uniform sloping beach as shown in Fig.4.7, the water depth, referring to Fig.4.8, is prescribed by the following formulas: h = x tan o r > Rcos(
h = xtanc h' r < Rcosa(
r I= (X X.l+ (Y _-.2
h' = h"/ cos ce,
h -B + V-B2 -4C
B = 2r and
hn tan2 a sin a
C r 2 R 2
sin 2 a tan2 a
The values used in the examples are: R = 0.7521m; h, 0.085m; Ax A y = 0.15m and the center of shoal is located at grid point i=17 and j=13.
The model is first run for a case with no inlet and then for a case with the presence of an inlet. In both cases, the input waves form an angle of 5 with the normal of the shoreline.

a: CI
Sco 1'.30 2.60 3.90 5.20 6.50
In in
00.oo 1.30 2.60 3.90 5.20 6.50
Figure 4.5: Wave Amplitude Variations across Centerline Aoko = 0.16; - Aoko =0.32

\:A \
'0 oo 1.30 2.60 3.90 5.20 6.50
M C.
T.oo 1.30 2.60 3.90 5.20 6.50
Figure 4.6: Wave Amplitude Variations across Centerline Aoko = 0.16; - Aoko = 0.32

Figure 4.7: Bottom Topography a Circular Shoal with Parabolic Configuration on the Plane Beach

h oh
Figure 4.8: Description of Water Depth
The other input conditions are
T = l.Osec.; ao = 0.02m; a = 0.075; f 0.01
breaking mixing coefficient = 0.01m2/sCc.
For the case of no inlet, the eddy viscosity coefficient is equal to zero and for the case with inlet, a value of D.005m2/sec is used. The inlet discharge is equal to 0.007ms/sec.
The velocity field with no inlet is shown in Fig.4.9. Waves break at the top of the shoal. The current induced by the waves is characterized by a strong onshore flow over the shoal consistent with the wave incident angle; at the lee side of the shoal, the flow splits into two directions, one upstream and one downstream. The upstream flow develops into two vortexes rotating in opposite directions; the downstream flow forms a longshore flow system with two portions. This longshore current pattern at the downstream boundary is sketched in Fig.4.10.
The velocity field with the inlet is shown in Fig.4.11. The effects of inlet flow are seen

to be quite pronounced. First of all, the two upstream vortices are strengthened and pushed farther offshore. Secondly, the downstream longshore current system is also strengthened with a flow reversal close to the shoreline, which is caused by the lower mean water surface elevation at the mid-section of the shoreline in the shadow of the shoal and also cause d by the lower pressure due to the higher velocities in the inlet, and the influence is extended further offshore which is due to the fact that the momentum of the inlet discharge pushes the flow into the lower pressure region. The mean water surface elevation is also lowered in this region. The longshore current distribution at the downstream boundary for the case with inlet is also shown in Fig.4.1O. The location of null flow on the lee side of the shoal is now being pushed seaward onto the shoal due to the inlet flow momentum.
Figures 4.12 and 4.13 show the contour lines of wave amplitudes. The unit in these figures is 10-4 meter. For the case of no inlet, the wave refraction effect causes the waves to focus on the shoal and this leads to wave breaking and a low region just leeward of the shoal crest. The diffraction effect, on the other hand causes waves to focus on the leeward of the shoal and creates a local high region there.
The effects of current are most pronounced in the center region between the inlet and the shoal. The local high region behind the shoal as mentioned above is enhanced considerably. The waves at the location of inlet entrance are also steepened and the heights increased.

* .-..-~- S.. ~
- S \ %
- I t / ~ .- *. I t / ~ -.5- 5. .. S
- \ ** I
- / .- I ~* 5- ~
- '- S
- .'
- / /1 /
- 'I I 5
- '~ / a
- I S
* a S .
0. 222~+00

t O~. ?L_2 .E+0
Figure 4.9: Velocity Vector Field for Circular Shoal Located on the Plane Beach (no inlet)


CD 0.00 0 .0 oo0I.00 20.0 oo
Figure 4.10: Longshore Current Distributions at Lowest Boundary. No Inlet; -- with
0 eA
0 0
0-20.00 -10.00 0.00 10.00 20.00 30.00
Figure 4.10: Longshore Current Distributions at Lowest Boundary. No Inlet; - with Inle,t

-"s4.. 11111/
.. j 1 1/ -----. -1 1
._._ 1 1_. 1. t4 1 4 .... .... ,.
. .*. . .
*- *- *. \ , ,. . a ,
** I I *
- / 111' i
* 0 I t* * I I
** * . d 4.
- S I S .
O. 216E+00
Figure 4.11: Velocity Vector field for Circular Shoal Located on the Plane Beach with Current from Inlet

Figure 4.12: Contour Plot of Wave Amplitudes on the Circular Shoal Located on the Plane Beach (no inlet)

Figure 4.13: Contour Plot of Wave Amplitudes on the Circular Shoal Located on the Plane Beach with Current from Inlet

In order to test the validity of the numerical model, laboratory experiments were conducted in a square basin in the laboratory of the Coastal and Oceanographic Department, University of Florida. The experimental set up, procedures and results together with numerical simulations are reported in this chapter.
5.1 Purpose of Experiments
Although there are many theoretical studied in recent years on the subject of wavecurrent interaction, laboratory experiments and field data are scarce. Ismail-Awadalla (1981) conducted laboratory experiments in wave-current interaction on a flat bottom with a nozzle injecting dye against the normal incident waves. The results demonstrated the physical phenomenon of current-wave interaction but were not sufficiently quantitative to verify numerical models.
In the present experiments, we focus on a beach-inlet system, where the phenomena of wave transformation over a gentle slope, coupled with the effects of a non-uniform current from an inlet can be studied. There are four basic quantities to be investigated: mean velocities: U in x direction and V in y direction; the mean water surface displacements iT and the wave amplitudes A. Since W1 is a small quantity and requires a long time to establish a steady value, the effects of wave reflection in the wave basin preclude meaningful measurements. Only wave amplitude and current field were investigated.

5.2 Facilities and Instruments
5.2.1 The Laboratory Facilities
The experiments were conducted in a square 7.22 m by 7.22 mn wavebasin with a depth of 0.8 m as shown in Fig.5.1 and 5.2. A fixed plane beach was constructed of cement with 1:13.3 slope (or s = 0.075). The toe of the beach was 0.85 m from the wave generator. The inlet was located at the center of the beach end. It was 1.44 m wide and the bottom was horizontal. The total length from the leading edge of the inlet to the diffuser was 1.2 m as shown in Fig.5.2. For all cases tested the water depth was kept 0.369 m measured from the bottom of the basin. Thus, the water depth was 0.045 m in the inlet. A grid system was established as shown in Fig.5.3 to facilitate numerical computation and laboratory measurement. The grid size is 0.12 mn x 0.18 m. There are 43 grids in the x-direction and 40 in y-direction. The inlet was located at the center of basin from j = 17 to j = 24. The leading edge of the inlet was at i = 37. The toe of the sloping beach was located at i = 1.
Waves were generated by a flap type wave generator. The bottom of the flap was elevated 0.2 m above the bottom to permit the excess water from the inlet to return freely under it. The wave generator was driven by a D.C. motor. Wave period and amplitude were controlled by adjusting the frequency of the motor and the eccentricity of the linkage. Uniformly long crested waves can be obtained for small amplitude waves. Cross waves would appear as the wave amplitudes became larger. Standing waves also existed under certain circumstance due to the short distance from the toe of the beach to the wave flap.
To prevent the leakage of wave energy from behind the wavemaker, horse hair mats were placed between the wave paddle and the wall, so that the wave motion behind the paddle was greatly reduced. On the other hand, due to space limitation, wave energy absorber could not be placed at the beach end. Therefore, reflection of waves was serious at times.
The current system consisted of a pipe system, a flow filter, a diffuser and a reservoir. As shown in Fig.5.1, a 2.5 horsepower centrifugal pump located between the flow filter and the diffuser was used to transport the water. The flow filter situated behind the wave flap


4 7.22m


Flap Type Wavemaker

Cross-Section A-A


~- FIIter

Cross-Section B-B

Figure 5.2: Cross-Section Views of Wave Tank

0.8m I1-



Figure 5.3: Cross-Sections of Measurement



- I I I II F l !l l l l F -l iT l l l l l

29 27 24

14 12

has ten control valves, allowing the water to flow naturally to the recirculation pipe without disturbing natural flow condition. A diffuser with a length 1.4 m was installed in a trench at the end of the inlet. It was a section of pipe with a series of different sized holes to distribute the flow uniformly. The gap around the diffuser was filled in with pebbles and a filter screen made by two layer a net filter with horse hair in between was set up in front of the trench to promote uniform flow across the inlet. However, the distance between the diffuser and the intersection of the shoreline and the inlet entrance was only 0.6 mn. It was very difficult to attain a uniform flow condition for such a short distance. After many fine tunings, an acceptable uniform flow condition was finally reached. Fig.5.4 shows the velocity distribution in the inlet at the cross-section that intersects the shoreline. The adopted final velocity distribution was 7/25/86 profile.
5.2.2 Instrumentation
The main properties measured included wave characteristics and current. Four resistance type wave gages were used to measured the spatial distributions of wave heights and periods. Current was measured by a Marsh-Mcburney type electro-magnetic current meter. In addition to wave and current, the discharge from the inlet was monitored by a commercial flow meter manufactured by Sitnet Scientific. The wave gage was calibrated statically by varying the gage submergence stepwise. It was found that the gage could not maintain linearity for the full gage length. Therefore, gage calibration was finally done for each location prior to each test. The current meter was supplied with calibration but was further calibrated, verified and adjusted by towing test in a long flume in the same laboratory. The flow meter was also calibrated by measuring the discharge filling in a reservoir. Fig.5.5 shows the calibration curves of the three instruments; the results were consistent and stable.
A 6-channel strip chart recorder was used to display the data. The wave gage signal was first amplified before feeding into recorder. The current meter signal was time-integrated to eliminate short term fluctuations. The layout of the instrument arrangement is shown in Fig.5.6

x 7/9/86
0 7/21/86 0 7/25/86

0 0
0 0 0 0 x
0 0 x 0 x 9 0 0
0 O 00 0 0 0 0



0.2 -


Figure 5.4: Velocity Distribution at I=Idry

0.8 -

0.6 I-





2 4 6 8 10

>: 10 00
a 0
> 10

cc 0


6 8 10
RECORDING (cm/sec)


Figure 5.5: Instrument Calibrations

x X Direction
O -X Direction
-- O Y Direction
6 -Y Direction