Numerical modeling of waves in the nearshore zone with permeable structures

UFL/COEL-95/013

NUMERICAL MODELING OF WAVES IN THE
NEARSHORE ZONE WITH PERMEABLE STRUCTURES

by

Santiago Alfageme

Thesis

1995

NUMERICAL MODELING OF WAVES IN THE NEARSHORE ZONE
WITH PERMEABLE STRUCTURES

By

SANTIAGO ALFAGEME

A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF MASTER OF ENGINEERING

UNIVERSITY OF FLORIDA

1995

ACKNOWLEDGMENTS

I would like to express my sincere appreciation and gratitude to Dr. Hsiang Wang,

my advisor, who gave me the opportunity to come to the University of Florida and enroll in

the master's program. He, along with the rest of the professors in this department, specially

the other two members of my committee Dr. Robert G. Dean and Dr. Robert J. Thieke, has

guided me through numerous research and class work problems. Appreciation is also

extended to Dr. Masahiko Isobe for his comments and help during the past few months.

The author would also like to thank the other members of the research group he has

worked for during the past two years, Dr. Li-Hwa Lin, Dr. Taerim Kim, and "future Dr." Xu

Wang, their friendship, help and suggestions made those two years infinitely easier.

Gratitude is extended to the staff of the Coastal and Oceanographic Engineering

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

Oceanographic Engineering who helped the author in many ways. Special acknowledgment

is given to Becky Hudson, Helen T. Twedell and John M. Davis for their patience with me

and the rest of the students.

Thanks also go to my fellow students Mike Krecic, Paul Devine, Rob Sloop, Jie

Zheng, Yingong Li, and all the rest; their company and friendship have made my stay in

Gainesville an memorable experience.

Finally I want to thank the most important people in my life, my parents and two

siblings, who have supported my goal of getting this degree even though they did not like the

idea of me being so far away from home. And to the person that has made that distance

bearable, my "novia", Leah Gambal, my deepest love and gratitude.

TABLE OF CONTENTS

ACKNOWLEDGMENTS ............................................. ii

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

ABSTRACT ........................................................ viii

CHAPTERS

1 INTRODUCTION ................................................ 1

1.1 Literature Review ............................................. 3
1.2 Sum m ary of contents .......................................... 6

2 WAVE ENERGY PROPAGATION THROUGH VARYING CURRENT FIELD 8

2.1 Governing Equations and Boundary Conditions
2.2 Wave Energy Equation ...................
2.2.1 Time-Averaged Wave Energy ........
2.2.2 Time-Averaged Wave Energy Flux ...
2.2.3 Time-Averaged Wave Energy Equation
2.3 W ave Action Equation ...................

3 MATHEMATICAL WAVE MODEL ................................ 20

3.1 M ild Slope Equation ........................................ 21
3.2 Derivation of Governing Equation ................................ 25

4 NUMERICAL SCHEME AND TESTING OF WAVE MODEL ............. 28

4.1 Numerical Scheme of Wave Model ............................... 28
4.2 Stability Analysis .................................. ......... 32
4.3 Boundary Conditions ......................................... 34
4.3.1 Lateral Upwave and Downwave Boundary Conditions .......... 34
4.3.2 Boundary Conditions for Groins ........................ 37
4.4 Description of the Wave Model .................................. 42
4.5 Testing of the W ave M odel ...................................... 45
4.5.1 Wave Shoaling and Refraction ............................. 45

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

4.5.2 W ave Diffraction and Reflection ........................... 49
4.5.3 W ave-Current Interaction ................................. 59

5 SURF ZONE MODEL .......................................... 62

5.1 Time-Averaged Wave Energy Equation in the Surf Zone .............. 62
5.2 Wave Action Equation in the Surf Zone ............................ 65
5.3 Wave Height Transformation in the Surf Zone ...................... 66

6 MODEL APPLICATIONS' ......................................... 72

6.1 W ave Field over a Paraboloidal Shoal ............................. 72
6.2 Wave Field around a Detached Breakwater ......................... 76
6.3 W ave Field around Groins and Jetties ............................. 81
6.4 Wave Field around Permeable Groins ............................. 83
6.4.1 Single Groin in Constant Water Depth ....................... 83
6.4.2 Single Groin on a Plane Beach ............................. 88
6.4.3 W aves around a Field of Groins ............................ 90
6.5 Preliminary Results from a Circulation Model ....................... 91
6.5.1 Circulation Results for a Plane Beach ......................... 95
6.5.2 Circulation Results for a Groin Field ......................... 98

7 CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER STUDY .. 101

7.1 Summary and Conclusions ..................................... 101
7.2 Recommendations for Further Study ............................. 103
7.2.1 W ave Propagation M odel ................................ 103
7.2.2 Boundary Conditions ................................... 104
7.2.3 Breaking Criteria ...................................... 105
7.2.4 Radiation Stresses ...................................... 105

APPENDIX

A REFLECTION AND TRANSMISSION FROM POROUS STRUCTURES
UNDER OBLIQUE WAVE ATTACK .............................. 107

BIBLIOGRAPHY ................................................. 113

BIOGRAPHICAL SKETCH .......................................... 117

LIST OF FIGURES

4.1 Sub-Grid system for wave model. ................................ 30

4.2 Location of grid points where the wave potential
will be solved using a boundary condition. ......................... 35

4.3 Description of wave field around a groin. ........................... 37

4.4 Schematic diagram .......................................... 41

4.5 Computational domain of the model. ............................... 43

4.6 Comparison of wave shoaling ................................... 47

4.7 Comparison of wave refraction, 60= 20. .............................. 47

4.8 Comparison of wave shoaling/refraction, 6 = 20. ...................... 48

4.9 Comparison of wave refraction, 0;=30. ................................ 48

4.10 Coordinate system for analytical diffraction/reflection solution. ........... 50

4.11 Grid and structure layout used for first diffraction test.
The values of x and y refer to model coordinates. ...................... 51

4.12 Comparison of wave diffraction/reflection coefficient.
Top 0=5, bottom 0=15 ........................................ 52

4.13 Comparison of wave diffraction/reflection coefficient
along five different sections, 0=50. .................................. 53

4.14 Comparison of wave diffraction/reflection coefficient
along five different sections, 0= 15. ............................... 54

4.15 Comparison of wave diffraction/reflection coefficient
along five different sections, 0=300. ................................. 56

4.16 Comparison of wave diffraction coefficient, 0=90. .................... 57

4.17 Comparison of wave diffraction coefficient, 0=70 ..................... 58

4.18 Conditions of collinear and shearing current. ......................... 59

4.19 Comparison of collinear wave-current interaction. ..................... 60

4.20 Comparison of shearing wave-current interaction,
amplification of wave height. ...................................... 60

4.21 Comparison of shearing wave-current interaction,
change in wave angle. .......................................... 61

5.1 Comparison between breaking analytical and numerical results.
Test 1: T=3 s, Test 2: T=5 s, Test 3: T=7 s .......................... 71

6.1 Bathymetry for paraboloidal shoal configuration,
after Ito and Tanimoto (1972). ................................... 73

6.2 Wave height comparisons between computed and
experimental results by Ito and Tanimoto (1972). ...................... 74

6.3 Numerical model results for shoal configuration
used by Ito and Tanimoto (1972) .................................. 75

6.4 Numerical model results for detached breakwater test. ................... 78

6.5 Comparison of location of breaker line behind a detached breakwater,
experimental data after Watanabe and Maruyama (1986). ............... 79

6.6 Comparison of distributions of wave height around a detached
breakwater, experimental data after Watanabe and Maruyama (1986) ...... 79

6.7 Numerical model results for jetty test. ............................... 81

6.8 Comparison of location of breaker line around a jetty,
experimental data after Watanabe and Maruyama (1986). ............... 82

6.9 Comparison of cross-shore distributions of wave height around a jetty,
experimental data after Watanabe and Maruyama (1986). .............. 82

6.10 Free surface and amplification factor around one groin, Case 1. ........... 85

6.11 Free surface and amplification factor around one groin, Case 2. ........... 85

6.12 Free surface and amplification factor around one groin, Case 3. ........... 86

6.13 Free surface and amplification factor around one groin, Case 4. ........... 86

6.14 Amplification factor along four sections parallel to groin. ............... 87

6.15 Comparison of model results for the cases of fixed and variable
reflection/transmission coefficients. ................................ 89

6.16 Offshore incident wave angle influence on
reflection and transmission coefficients. ............................. 89

6.17 Model results for groin field #1. .................................. 92

6.18 Model results for groin field #2. .................................. 92

6.19 Wave height along the upwave side of the groins. ...................... 93

6.20 Wave height along the downwave side the groins. ..................... 93

6.21 Circulation in a plane beach: Case. ................................. 96

6.22 Circulation in a plane beach: Case 2. ............................... 97

6.23 Baltic coast bathymetry and groin location. ........................... 98

6.24 Wave-Current field. ............................................. 100

6.25 Mean free surface and Currents .................................. 100

A.1 Schematic diagram of porous structure. .......................... . 107

Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Engineering

NUMERICAL MODELING OF WAVES IN THE NEARSHORE ZONE
WITH PERMEABLE STRUCTURES

By

Santiago Alfageme

August 1995

Chairman: Dr. Hsiang Wang
Major Department: Coastal and Oceanographic Engineering

This study presents a numerical model that predicts the transformation of

monochromatic waves over complex bathymetry and includes refractive, diffractive and

reflective effects induced by structures placed along the wave main direction of propagation.

Particular emphasis is placed upon permeable structures which are common in the nearshore

zone. Since the model has been developed with the idea of describing wave propagation in

the surf zone, and in anticipation of the possibility of combining it with a circulation model,

wave-current interaction effects are included in the derivation of the governing equations,

namely the 'wave energy conservation equation', 'wave action equation' and 'kinematic

conservation of intrinsic angular frequency'. Finite difference approximations are used to

solve the governing equations, and the solution is obtained for a finite number of rectilinear

grid cells that comprise the domain of interest. The solution technique employed, based on

a marching Gragg's method, is direct, and eliminates the computer storage problems

associated with large matrices inherent in other methods. This efficiency is obtained at the

cost of resolution of the wave components normal to the incident wave direction.

Boundary conditions for all sides of the numerical domain are derived, as well as the

necessary conditions to reach a solution on grid points next to permeable structures. The

effects induced by different permeability values of the structure on the wave field are well

simulated through the use of the appropriate boundary conditions that include calculated

reflection and transmission coefficients, which in turn are based on the incident wave field

and the physical characteristics of the structure.

A surf zone model based on the work of Lee and Wang (1993) is derived. The model

takes a different approach to surf zone dynamics by showing that wave action is conserved

in the surf zone. This, together with the wave energy equation, allows for the derivation of

a surf zone model that requires fewer empirical coefficients and is better suited for a general

three dimensional bathymetry.

Model results are compared first with data from analytical solutions. Then, the

capability and utility of the model for real coastal applications are illustrated by application

to three different experimental test layouts, where measured data were available. Finally, to

study how the model performs when solving problems that include permeable structures,

several cases are presented, from a single groin in constant water depth, to a shoreline

protected by a groin field. The preliminary results from a circulation model coupled with the

wave model are shown as well.

CHAPTER 1
INTRODUCTION

Coastal engineers and scientists face an increasing demand for more accurate

solutions to a wide range of coastal problems, many of them involving the transformation of

wind generated waves as they propagate from deep water, over complex bathymetry, into

shallow water where they may be further affected by man made structures. The nearshore

wave field will generally need quantitative evaluation in order to design those coastal and

harbor structures, to study and comprehend coastal processes, and to predict nearshore

currents and the resulting sediment motion. Thus, there is a need for an efficient and

accurate method to estimate nearshore wave conditions with arbitrary shoreline shape and

nearshore bathymetry, as well as the presence of coastal structures.

Several mathematical models have been proposed over the years, as well as their

corresponding numerical solutions. Most of them are capable of producing reliable results

under certain conditions, but with a series of drawbacks as well. That is, most of the models

are restricted to a certain range of imposed conditions.

The main objective of this study was to develop a model that would effectively

simulate the nearshore wave field in the presence of structures aligned approximately in the

main wave direction such as a field of groins with different permeability values. Over the

years, different kinds of groins have been used for coastal protection. There are a great

variety of shapes and materials used for groin construction; traditionally these structures

consist of rubble mound. The permeability of this kind of groin could be very different

2

depending of the size of the stone used, the dimensions of the groin, the use of different

layers with different stone sizes across the groin, etc. Completely impermeable groins have

several disadvantages, including the following:

* They are very expensive.

* They tend to induce deep erosion pits seaward of the groin head because of locally

increased currents.

* They may stimulate rip-currents and seaward loss of sand.

* Extensive lee-side erosion may occur on the adjacent unprotected shoreline.

Other types of groins include the ones composed of metal sheets that are driven into

the beach, or even made out of precast concrete units. An interesting solution often used in

places like the Netherlands and Germany consists of single or double permeable rows of

wooden piles perpendicular to the beach. This system costs only 10 to 25% of the cost of

rubble stone groins used in those countries. In contrast to the impermeable groins, which

seek to form a complete obstruction to the longshore current and sediment transport, highly

permeable groins like pile screens are meant as an artificial hydraulic resistance in order to

gradually reduce the longshore rate of sediment transport. The advantages of this kind of

groin can be summarized as follows:

* Low cost structures.

* Less pronounced current concentration seaward of the head.

* Reduced tendency to form rip-currents.

* The aforementioned ability to reduce longshore transport in a gradual way, thus

reducing the lee-side erosion behind massive structures. An adequate screen field

3

layout will distribute the total sediment deficit over a larger area, thus decreasing the

rate of recession of the shoreline.

Because of the different effects induced by completely impermeable groins, and

groins with different degrees of permeability, it is important to develop a model that includes

the correct boundary conditions at the grid points next to these structures, so that their effects

can be properly accounted for. Wave propagation models reviewed by the author usually do

not allow for reflective boundary conditions, while a few do by imposing a complete

reflection condition at the upwave side and zero transmission on the other side of the

structure. In this thesis, adequate boundary conditions that allow for different values of

reflection and transmission coefficients are derived and implemented in the model.

To solve the wave propagation problem with the inclusion of wave-current interaction

effects, such that the model can eventually couple with a nearshore circulation model, two

equations governing the dynamics are needed, the wave action equation and the wave energy

equation. Another basic equation to be satisfied is the conservation of intrinsic wave

frequency. All of them are formally derived first in this study.

Since the model is to be applied in the nearshore zone, a surf zone model is

necessary. In this work, the model first developed by Lee and Wang (1993) is used. The

model is based on wave action and wave energy conservation and fully includes wave-

current interaction.

1.1 Literature Review

Much of the early work done in the area of wave propagation was based on wave

refraction theory and entailed the construction of wave rays using linear potential wave

theory. This approach is an extension of Snell's law from optical wave refraction to the

4

analogous problem of water waves propagating over straight and parallel contours. Thus,

an important assumption made was that the contours of the given bathymetry were locally

straight. These ray methods do not yield solutions at prescribed points, but rather along the

wave rays only. During the 1960s and early 1970s the linear wave refraction problem was

solved in a more efficient way through the use of computers. Noda et al. (1974) was the first

to devise a numerical scheme that solved the wave energy equation and used the

irrotationality of the wave number to obtain wave height and direction at points on a

rectangular grid. Background currents were also included in the formulation.

A serious limitation of the refraction models is the omission of diffraction effects,

which creates unrealistically large wave heights when the wave rays converge to a small area.

The problem of water waves diffraction was first studied analytically by Penney and Price

(1944), who showed that Sommerfeld's solution of the optical diffraction problem is also a

solution of the constant depth water wave diffraction problem. Again, Penney and Price

(1952) solved the constant depth diffraction problem for breakwaters assuming a semi-

infinitely long, infinitesimally thin barrier. In order to include diffraction by coastal

structures in existing wave propagation models, first a solution to the diffraction problem

was given, using tabulated analytical solutions for constant water depth, and the results from

the ray method were superimposed after a few wave lengths.

Berkhoff (1972) derived an elliptic equation describing the complete wave

transformation process for linear waves, including both refraction and diffraction, from deep

to shallow water in terms of velocity potential function with the assumption of mild slope.

The elliptic equation was first solved using finite element methods with good success but at

high computational effort. Also, the treatment of boundary conditions is generally difficult.

5

Since then many others have solved the wave transformation problem for complex but

idealized bottom configurations. Ito and Tanimoto (1972) proposed a numerical method for

harbor wave analysis that allowed easy specification of boundary conditions, but requires

modification to describe wave shoaling. Both Radder (1979) and Booij (1981) derived

parabolic approximations to the original mild slope equation, the latter extended the scope

of the equation to include the effects of currents. Kirby (1984) made corrections to the

equation derived by Booij. Ebersole (1985) suggested a finite-difference solution to the

elliptic equation in a manner suitable for large coastal areas. Watanabe and Maruyama

(1986) proposed a time dependent version of the mild slope equation consisting of two first

order equations, which are obtained by separating the original mild slope equation in terms

of the water surface elevation and the flow rate. Madsen and Larsen (1987) also developed

a model based on a system of hyperbolic equations similar to those governing horizontal

flow. The model used an ADI algorithm to iteratively converge to a stationary solution. The

difference with the set of equations used by Watanabe and Maruyama, as well as Copeland

(1985), was that the harmonic time variation was extracted from the equations, so that the

method converges to a final solution more rapidly.

Panchang et al. (1988) presented a method to solve the boundary value problem of

the elliptic mild slope equation in large domains based on a marching 'Error Vector

Propagation' method that solved backscattering and propagation in the -x direction, thus

overcoming the major limitation of the parabolic approximation.

Lee and Wang (1993) evaluated the performance of several numerical wave models,

including two models based on a set of hyperbolic equations, two based on an elliptic

equation and finally a fifth using a parabolic approximation. They concluded that no single

6

model outperformed the other and thus the selection of a model for application depends on

the intended purpose.

More recently, in an attempt to model irregular nonlinear waves in very shallow

water, models based on Boussinesq equations have been used (Abbot et al. (1978)). Those

equations include non-linearity as well as frequency dispersion. Furthermore, they operate

in the time domain, so that irregular waves can be simulated.

Theoretical solutions for reflection and transmission coefficients and for the wave

field have been derived for porous structures using different methods of approach (e.g. Sollit

and Cross, 1972; Madsen 1974). Recently Darlymple, Losada and Martin (1991) have

obtained results for transmission and reflection coefficients for obliquely incident waves.

In their studies the equivalent work condition by Lorentz is established yielding a potential

flow problem. Matching pressure and horizontal mass flow at the region interfaces solutions

of reflection and transmission coefficients are then obtained. They also showed how

neglecting the non-propagating evanescent wave modes the solution technique can be

simplified for all practical applications.

1.2 Summary of Contents

It has been mentioned already that the objective of this study is to develop a

numerical model that simulates the nearshore wave field in the presence of porous structures

such as offshore breakwaters, jetties and groins.

In Chapter 2 the equations governing the wave-current interaction problem are

derived. Basic assumptions and fundamental equations are first presented, from them

dynamic equations describing wave energy transport and wave action conservation as well

as the kinematic intrinsic frequency conservation equation are obtained. Chapter 3 derives

7

both the hyperbolic and elliptic mild slope equations that govern the refraction-diffraction

of waves. Based on the equation of hyperbolic type, and eliminating the harmonic time

dependence, an elliptic equation more suitable for a numerical solution is presented. Chapter

4 describes in detail the finite difference form of the governing equation and the numerical

solution by a forward marching Gragg's method. Boundary conditions are established on all

sides of the numerical domain as well as on grid points immediately next to possible

structures. A section in this chapter gives a general description of the wave model, detailing

its components and how they interact until a final solution is reached. The resulting model

is tested against several analytical solutions to basic wave propagation and wave-current

interaction problems. Chapter 5 introduces the theory for wave transformation in the surf

zone based on the work of Lee and Wang (1993). The surf zone model is implemented in

the numerical model and computed values are compared to analytical solutions for a beach

with straight and parallel contours. Chapter 6 demonstrates the applicability of the model

based on the results obtained for several situations in which measured data are available from

the literature. The model is also used to obtain the wave field around permeable groins in

constant water depth and on a planar beach. The case of a groin field is also simulated with

the model. The preliminary work and results leading to the complete incorporation of a

nearshore circulation model are summarized too in this chapter. Chapter 7 discusses the

results, the advantages and drawbacks of the model, and makes several recommendations for

further study. One Appendix, where the theory for a linear wave impinging obliquely on a

vertical sided porous structure is presented, is included.

CHAPTER 2
WAVE ENERGY PROPAGATION THROUGH VARYING CURRENT FIELD

When waves propagate through a region with currents, wave energy propagation

speed is no longer equal to the wave group celerity, but it is modified by the current effect.

Other characteristics of the wave train will also be altered, including the wave height, length

and period. This situation is commonly observed in regions with strong currents, like tidal

inlets. In this chapter, two fundamental dynamic equations governing the behavior of wave-

current interactions are derived. They are the wave energy equation and the wave action

equation.

In Section 2.1, the basic problem formulation and boundary conditions are given. In

Section 2.2, the depth-integrated wave energy equation is derived and the exact forms of

wave energy and wave energy flux are presented in wave-averaged quantities. Section 2.3

introduces the wave action equation and the conservation of intrinsic frequency equation.

It will be shown how the wave energy equation originates from the conservation of energy,

while the wave action equation is derived from the free surface boundary conditions.

2.1. Governing Equations and Boundary Conditions

The first assumption made is that the fluid motion is irrotational and thus a velocity

potential q exists and water particle velocities are given by Vb. The kinematic and dynamic

boundary conditions to be satisfied at the free surface, z= 77, are, respectively,

77t + Vh'Vhrl 0z =0 (2.1)
(2.1)

+ ()2 + gz = C(t) (2.2)
2

where C(t) may depend on t, but not on the space variables. We may take C(t)=O without any

essential loss of generality. The subscripts t and z indicate the differentiations with respect

to time and z-axis, respectively. The symbol Vh representents the horizontal gradient

operator.

The cartesian coordinate system is used with origin at the still water level, x(x,y) in

the horizontal plane and z directed vertically upwards. The velocity vector, U(u,v,w), is

related to q5by

u-a, v- and w
ax ay az

The velocity potential and free surface displacement are assumed to be composed of

current and wave components.

O(x, t) = 0c(x;z, t) + efw(X, z,t) (2.3)

r(x, t) = rc(x; t) + erw(X, t) (2.4)

where e is an undefined factor used to separate the current from the wave part of the velocity

potential. The symbol ';' is used to express that the current component of the velocity

potential as well as the free surface vary slowly over time when the time scale is much larger

that the wave period. Also, q5c can account for small vertical variations of the current.

Equations (2.1) and (2.1) are expanded in a Taylor series around the mean water level z= q,

[it +Vho'Vhi z z=c 0,'+e [a(it + 2 +-..=0

[t + ( )2 + = + [ t +(V 2 + gZ] +... = 0
2[t + 2 ]7,

Substituting equations (2.3) and (2.4) into the above equations gives

[(1,)+Vhc.Vhec-()z]z^=1 wt h c. hw hw wz+ hw c]z= (2.5)

[(c)t + (q) + g ]Z + e (w) +V Vc Bw + wz. =0 (2.6)

The above equations are now separated for wave and current parts and truncated to retain

only the first order terms, 0(e):

(c) +Vhc Vhc (2.7)

1 0A h )2 (2.8)
g 2
2 (2.8)
) + +h h
W -Dt h wc) w c (2.9)

1 D b (2.10)
g Dt

where D/Dt- /lt+Vh qC-h.

Note that the last two terms in Eq. (2.9) are, in general, of higher order than the first

one, only when diffraction occurs they become important. For slowly varying water depth,

the wave part of the velocity potential may be written as

w(x,z,t) =f(z) w(x,t) (2.11)

where f(z) = cosh k(h+z)/cosh k(h+ qt) is a slowly varying function of x, k is a real wave

number. For progressive waves, the velocity potential and free surface displacement can be

written in terms of wave averaged, slowly varying quantities as

w(x,z,t) =f(z)A(x;t)ie i' (2.12)

rw(x,t) =a(x;t) e i (2.13)

where a is commonly defined as the wave amplitude. The phase function is defined as

fr= (K.x-ot), where K is the wave number vector including the diffraction effects, and w is

the absolute frequency. The relation between a and A can be established by the dynamic free

surface boundary condition specified in Eq. (2.10), which, after substituting Eqs. (2.12) and

(2.13) into it, yields

D6
Dt

-gaei=" { +U V} {Aie'i'}
at

= odAeiA+ {A +~ -AA}iei (2.14)

where od is the intrinsic wave frequency including the diffraction effects. Its value is defined

as ad = K, and U as Vh .

Equation (2.14) states that a and A should have a phase difference unless we impose

the condition

A +UVA =0 (2.15)
at

Then, the relation between A and a can be given by the following familiar expression

A = -g- (2.16)

Similarly, substituting Eqs. (2.12) and (2.13) into the kinematic free surface boundary

condition given by Eq. (2.9) yields

d = gktanhk(h + c) gVA Vic (2.17)
A
a +V-(Ua)+AK-V O=0 (2.18)
at

Again, the last term in both the above equations reflects the wave diffraction effect and,

under normal circumstances, is of a higher order.

2.2. Wave Energy Equation

In this section, the Eulerian expression of energy equation, which governs the local

balance between the rate of change of energy and the divergence of energy flux at a point,

will be presented. The Euler equation for incompressible and inviscid flow is

DU
P -V(p + pgz) (2.19)
Dt

taking the scalar product of U(u,v,w) with the respective terms in Eq. (2.19) and then

summing the three components yields

P [I-] =-U V(p+pgz) (2.20)
Dt 2

where q2=(u22+v2+2). With the use of the continuity equation, the mechanical energy

conservation equation becomes

-[ +V [U( +p +pgz)] =0 (2.21)
at 2 2

Integrating Eq. (2.21) over the water depth,

qf/ aW2 + W2 (2.22)
f 2dt 2] 2 z

The first term inside the integral represents the volumetric rate of energy change and the

second term gives the energy flux through the enclosing surface. The depth-integrated wave

energy equation has been derived from equation Eq. (2.22) by many including Longuet-

Higgins and Stewart (1961) and Witham (1962). A brief account is given here.

Using Leibnitz's rule, Eq. (2.22) can be written as

[a 2 ]dz f [U( +p+pgz)dz
-h -h
at2 22
-[U( +p+pgz)]*Vh-L[U( +p+ gz.Vhh (2.23)
2^ pg

2
Substituting the cinematic boundary conditions at the free surface anpgz)d at =-h,

Substituting the kinematic boundary conditions at the free surface and at z=-h,

w -J -U-Vh=77 (2.24)

W _h+UVhh =O (2.25)

and letting p be zero at the water surface, Eq. (2.23) becomes

S ]dz + pg 7 + U( 2 +p+pgz)dz=O (2.26)
at (2 .t 2
-h -h

Defining the total energy, E, and the total energy flux, F, as

E- P [ -dz + pg (2.27)
S2 2
-h

F f U( +p +pgz)dz (2.28)
-h

Equation (2.26) yields the well known energy equation given by Longuet-Higgins and

Stewart (1961) and Witham (1962),

FE
-t +Vh-F=O (2.29)
at

Here, E is the energy density in the water column per unit surface area and F is the energy

flux through the vertical surface enclosing the water column.

2.2.1 Time-Averaged Wave Energy

The energy per unit surface area given in Eq. (2.27) is expanded in a Taylor series

with respect to the mean water level z= lc,

S[(V,0)2 +()2 +-pg[(7 + e7)2 -h2] +e7[(V h )2+()
E = pf (V 2 2 h 2 w h2 2 1c
-h

-p f [(vh c+ +Vh 0)2+ (q5cz+ C wZ)2 + Ipg[(y7 +e,)2 h 2]
-h

+ p-en (V,(O +f eVOc 2+ ( +e wz (2.30)
2 "W['hC h2wI "C17 cwzIJ

Taking time average over the wave period and collecting terms associated with current (0(1))

and wave (O(e2)) separately, we obtain the following pair of equations,

EcP f 1[(V 2 + ( pg( h)2 (2.31)
-h

E p +V (w2 )2 P7 + P pg P w[V c w (2.32)
-h

where Ec and E are, respectively, the mean values of current and wave energy. The mean

wave energy density can be expressed in terms of the slowly varying quantities by

substituting Eqs. (2.12), (2.13) and (2.17) in Eq. (2.31):

E -g H2 (2.33)
8

where H is the wave height defined as twice the wave amplitude a.

A few general remarks regarding the definition of wave energy density are made here:

i. E is the energy density directly associated with the O(e) fluctuation motions only.

It does not take into account contributions associated with mean water level change.

ii. The last term in Eq. (2.32) represents the contribution due to wave-current

16

interaction. Longuet-Higgins and Stewart (1961) and Phillips (1977) all included this

term in the mean flow energy rather than in the mean wave energy. However, this

is truly a O(e2) term from the fluctuation motion.

iii. In the absence of current, E reduces to the conventional definition of energy density

in a wave field.

2.2.2 Time-Averaged Wave Energy Flux

The energy flux expression given in Eq. (2.28) can be given, using the Bernoulli

equation, as

F = -pfU-dz (2.34)
-h

Introducing the current and wave components defined in Eqs. (2.3) and (2.4) and expanding

in a Taylor series with respect to the mean water level z= qJ, we have

a a
F = -p f V c(c + w)-(c )dz -per,w[Vh(c w+w) c w+^ (2.35)
-h

Taking time averaging over the wave period and collecting terms of O(L), we have the mean

wave energy flux,

1 ap aqo a(2
F= -p Vh w dz [Vhq- +Vhw ] (2.36)
t -ht
-h

Substituting Eqs. (2.12) and (2.13) into Eq. (2.36) it can be obtained that,

= Cg+U k E (2.37)

The above equation contains an additional term, q,tk/wd, that does not appear in the

conventional form of this equation. This term will be of importance only in the case of

unsteady currents, and will not be included in the derivation of the governing equations for

this model. As explained above, the definition of E is different than the conventional term

too.

2.2.3 Time-Averaged Wave Energy Equation

Now Eqs. (2.33) and (2.37) can be substituted in Eq. (2.29), obtaining,

+Vh (Cg +U--k)E =0 (2.38)
at 6)

2.3 Wave Action Equation

Subtracting Eq. (2.10)xpgr7, from Eq. (2.9)xpgd4, and ignoring the higher order

diffraction effect, we obtain

aa 2
a +*(Us1(prlw7.~Pw ) -, =0 (2.39)

where U, is the current velocity of the mean flow at the water surface level. Substituting

Eqs. (2.12) and (2.13) into Eq. (2.39), the following equation is obtained:

a (Bie2i0) +V'(iBie + Bei2rgktanhk(h + 77C) + 1 0
(Bie2) v.( sBie ) aBei2I + 1=0-
at

where B is defined as

B g H2 (2.40)
8 a

Expanding and separating the harmonic motions,

ie 1 a+V B) + aBei2 2+ +1gknhk(h =01
Bt 1 2

which yields the dispersion relation, o = gk tanh k(h + ir), and the following wave action

equation:

aB
S+Vh [UsB] =0 (2.41)

If we substitute Eq. (2.16) into Eq. (2.18), the following equation is obtained,

aoA
dA +Vh [UsA]=0 (2.42)

Eliminating A from Eqs. (2.42) and (2.15), we arrive at the final equation that governs the

change of the intrinsic wave frequency in a current field:

S+ Vh [s a] = 0 (2.43)
at

which is termed as the 'the kinematic conservation equation of intrinsic frequency', or simply

the 'conservation equation of intrinsic frequency'. It is possible to derive two other

19

kinematic conservation equations. The phase function was defined as fr= (k.x-wt), where

k is the wave number after dropping the diffraction effects. Thus we can obtain k and 0 as

k=V o=^
at

Then, if we take the curl on k, we find that

Vxk=O

That is, the irrotationality condition on k. It also follows that

known as the 'kinematic conservation equation of the wave number'.

CHAPTER 3
MATHEMATICAL WAVE MODEL

A significant advance in the field of wave modeling was made by Berkhoff (1972,

1976), who derived the combined refraction-diffraction mild slope equation. This equation

eliminated the usual problem encountered in refraction studies caused by caustics. It can be

used for a wide range of ocean wave frequencies, since it converts, in the limit, to the deep

and shallow water equations. The usefulness of Berkhoff's equation in yielding good

simulations of wave behavior in a wide variety of situations has been proven by many

investigators. Berkhoff, Booy and Radder (1982) have used it to study short wave

propagation over arbitrary variations in bottom topography. Berkhoff (1976) and Rottmann-

Sode and Zielke (1984) employed it to study wave propagation into harbors. However, when

attempting to simulate wave conditions in large coastal areas, the numerical solution

becomes a formidable task in terms of computer time and memory.

In view of this problem, it became necessary to use the parabolic approximation

(Radder, 1979) of the mild-slope equation. With this method a solution is rendered feasible

in fairly large open coastal areas, but there are two limitations to it, the waves must have a

dominant direction of propagation, and reflections are neglected.

The extension of the mild slope-equation to cover wave-current interaction has been

developed by Kirby (1984), and the corresponding equation eventually yielded the energy

equation suggested by Bretherton and Garrett (1969):

D2 + (VU) -V'(CCgV) +(o2-k2CCg)k=O
Dt2 Dt

where is the complex potential at the water surface.

In Section 3.1, the mild slope equation revised by Kirby (1984) will be obtained from

the energy equation, assuming that the steady state condition of wave and current field could

be retained for the wave-induced currents which are usually accompanied by the waves. In

Section 3.2, the governing wave equation of the wave-current interaction model will be

derived from the mild slope equation. The solution obtained does not rigorously satisfy the

boundary condition for a sloping bottom, but even the first-order effect of the slope is

negligibly small. Booij (1983) and Smith and Sprinks (1975) numerically examined the

validity of the mild-slope equation and concluded that it is applicable to waves propagating

over slopes as steep as 1/3, and even to waves propagating across as step.

3.1 Mild Slope Equation

The mild slope equation has been derived by the variational principle and Green's

second identity. In this section, however, the mild slope equation will be derived directly

from the energy equation as described below. Introducing the velocity potentials and free

surface displacements composed of current and wave components as given in Eqs. (3.3-3.4),

Eq. (3.26) becomes

a 7 ( V +c2 + ( +
[ e)2 ]dz +g(77 + ew)-( 7+ erw)
2 2t
-h (3.1)

-Vh- f [(V + w c)( + eP dz =
-h

Taking Taylor expansion about z= qc,

a T (V + eq )2 a (v_ + _)_ 2 a
]dz +C [- rw + +g(77 +erw)-(+ ( l+ew)
2 at 2 at
-h

-Vh" f [(Voc + 'eV )()c + qw)t]dz + ew[(Vc + V ,)(c + ] =0
-hI

Collecting the terms of O(i) and ignoring the long-term fluctuation of current field:

S t E2 af (
[titi 2 ]dz +-( [ V' h w z +V = i07 (3.2)
-h -h

Substituting Eq. (2.11) into Eq. (3.2),

a Re2 (Vh ) 2 2 a (A5)a2
dz[ + fdz [ ] + [q V;VP]
-at 2 at 2 at
-h -h

+e Vh, f2dzVit +Ves f 0
-h

where, substituting the expression forf

f2dz -
-h

faz-

(3.3)

(3.4)

(3.5)

CCg

d -k2CCg
g

Therefore, Eq. (3.3) becomes

CCg a (Vh )2
g at 2

at a

+ 72-k2 CCg a( ) A
+ [ ]- [ V c w
g at 2 at

-V[h C w t Vc w+V w 0
IT

23
Hereafter, we omit the subscript denoting the wave mode, w.

CCgQVh-Vht + t( o2 k2CCg) + g(7V'V)t +g2177t
tVht(CCgVO) CCgV- Vh -gVh (Vc 17 ) =

The first and sixth terms offset each other and expanding the third and last terms,

,(o2 -k2CCg)+ g c7,Vc.V + g17Vc.V + g 217
tVh(CCgV ) gh(V 17) -gV h V Ot 0

and also offsetting the third and seventh terms,

t(o2 k 2Cg) +g gj7tV + g 2
tV-(CCgV) -gtVh-(V cq) =0

Substituting the dynamic free surface boundary condition of O(e) given in Eq. (2.10),

t(2 k 2CCg) t + g7tV V g( q + V'c.V) 17t
tVh(CCgV ) g h(V ) =0

t(o2 k 2CCg) g t tVh.(CCgV gtVh(V' q7) =0

Combining the second and fourth terms by use of the kinematic free surface boundary

condition of O(e) given in Eq. (2.9),

(o2 -k2CCg) -Vh(CCgV) g' =0 (3.6)

Now we consider two kinds of expressions of z. The first is obtained by substituting the

dynamic free surface boundary condition, Eq. (2.10), into the kinematic boundary condition,

(3.7)

g Dt2'+(V Dt)
g Dt2 Dt I ;

and the second expression is simply obtained from Eq. (2.11) as

(3.8)

which eliminates the slowly varying motion through the mean surface. With the first

relation, we finally get the mild slope equation as derived by Kirby (1984);

D2 + (VhU) -Vh,(CCgV) +( -k2CCg)=0O
Dt2 Dt

and with the second expression, we get the mild slope equation of elliptic type:

Vh(CCgV) + k2CCg =0O

which has the same form as the case of no current.

Equation (3.9) can be fully expanded, becoming,

+ 2UJ.(Vh)t+ 2uv + (u2 CCg) +(v2 -CCg)- 7
8t2 OxOy Ox2 y 2

where

J-= -(Vh.U)D- +V hV (CCg) -( k2CCg)q
Dt
u=.i+vj

(3.9)

(3.10)

(3.11)

~="i

25

It can be easily seen that Eq. (3.11) is a three dimensional second order partial

differential equation of hyperbolic type. This type of equation is also known as the Klein-

Gordon equation or telegraph equation.

It can be shown that Eq. (3.10) reduces to the well-known Helmholtz equation in the

case of deepwater,

V + k2k =0 (3.12)

and in shallow water Eq. (3.10) becomes the long wave equation,

k2CCg() + V()Vh(CCg) =0 (3.13)

3.2 Derivation of Governing Equation

In this section, the governing wave equation of the numerical model selected is

derived from the linearized mild slope equation Eq. (3.9) for waves interacting with currents.

It is rewritten below with the symbol descriptions

D2q+ (V.U)D -V-(CCgV) +(o2 k 2CCg) =0
Dt2 Dt

where 0 is the two-dimensional complex velocity potential,

C is the relative phase velocity (a/k),

Cg is the relative group velocity (Oaldk),

ais the intrinsic frequency (o2=gktanhkh),

w) is the absolute frequency,

k is the wave number,

h is the water depth,

Uis the steady current velocity vector (u,v),

and V = /alx + al/y, omitting the subscript h.

The values of aand k are determined by the Doppler relation, w = a+ U k.

If we express the two-dimensional velocity potential in the linear stationary field as:

(x,t) = (x)e -it (3.14)

where ( is the surface potential in steady state. Substituting Eq. (3.14) into Eq. (3.9) gives,

W 2iUT-V.V + (U-V)(U-V) + (U-V)(-iw, + U-V)
-V-(CCgV) + (o k 2CCg) =0

Rearranging and factoring out e -wt,

-iw[2U-V0 + (U-V) ] + (U-V)(I-Vq) + (V)(U-V) (3.15)
-V-(CCgVo) +(o -k2CCg) 0 =

The above equation is valid as long as the wave motion is time harmonic. The three-

dimensional (x, y and t) problem has been reduced to a two-dimensional (x and y) one with

an equation of elliptic type.

If we expand Eq. (3.15), the first term will result in,

-iw(25xu +2qyv +ux+ y = -i)(u + Vy)O-i 2uix -i)2vqpy

The second and third terms will give,

(U-V)(u y + Vy) + (u +Vy)(u + Vxy)=
(2uxu+vuy + uv x+v u2 xx+ 2uvxy + (uvx y+xy + V2yy

And the fourth and fifth terms,

-(CCg)xx CCg (CCg)yy (CCg)q + (02 2 -k2CCg)q

Combining all terms we get,

(u 2 CCg) + { -2iwu + 2uux + u y + uvy-(CCg)x}

+ 2uv + +-2iwv + 2Vy + uvx + uv-(CCg)y} ~,

+(v2 CCg)y + { -i(ux + vy) + 2 2-k2CCg} =0

To study the nature of the above equation we can rewrite it as,

aq5x +2b +cqy =F( x y, ,x,y)

where,

a = u2-CCg
b = uv
c = v2-CCg

in general, CCg is >> U2=2+v2, and,

ac b 2 = (CCg)2 CCgU2>0

Thus, Eq. (3.16) is an elliptic partial differential equation. In Chapter 4, a numerical method

to solve this governing wave-current interaction equation is presented.

(3.16)

CHAPTER 4
NUMERICAL SCHEME AND TESTING OF WAVE MODEL

In the last chapter a governing equation of elliptic type was derived from the original

mild slope equation of hyperbolic type. A general solution to this kind of governing equation

is not available for a wave field in the vicinity of structures in the nearshore zone; therefore,

numerical computations are needed. In this chapter, a numerical method is developed based

on the work of Lee (1993). The method falls under the general category of finite difference

schemes. Also, the boundary conditions used for the edges of the computational domain and

the grid points next to a groin are obtained. In the next section a general description of the

wave model is given. The last part of this chapter includes comparisons between the results

obtained from the wave model and a series of analytical solutions to basic wave propagation

and wave-current interaction problems.

4.1 Numerical Scheme of Wave Model

In considering the different possible numerical approaches to Equation (3.16), the

author tried to use one that would combine a minimum computational effort and with results

that would include refraction, diffraction as well as reflection due to structures. Several

mathematical models have been proposed for combined refraction and diffraction of simple

harmonic linear waves. All of them have advantages and drawbacks depending on the kind

of problem that one is trying to solve. The solution to the mild slope equation derived by

Berkhoff (1972) is very accurate and complete, but the treatment of boundary conditions is

generally difficult. The numerical method proposed by Ito and Tanimoto (1972) allows easy

29

specification of boundary conditions but requires modification to describe wave shoaling.

A parabolic approximation to the mild slope equation is a popular approach to alleviate some

of the problems addressed above, but it does not include reflections from structures such as

jetties and groins. The numerical solution adopted by the author will include such reflections

without making the computational effort excessive.

Equation (3.16) can be treated as an ordinary differential equation for 0, as given

below, so that the y-directional difference operator, D, is explicitly approximated by using

finite differences,

(u 2 CCg) + { -2iwu + 2uux + UyV + y-(CCg)x}
+2uvDy(x) + { -2iYv +2vVy +uvx + Uxv-(CCg)y}Dy(P) (4.1)
+(v2- CCg)Dyy() + { -iw(ux + Vy) + 02 -w2 -k2CCg} =0

where, the x axis is selected for the main direction of wave propagation. We convert the

above equation to a pair of first-order equations by simply defining the x derivative as a

second function.

OX = 01
x -- -1 2[{ -2iwu +2uux + uv + u-(Cg)} 1 +2uvDy) + -(]
CCg -u

where

7 = {-2iwv+2vVy + uvx + Ux-(CCg)y}D y( ) +( v2-CCg)Dy( )
+ {-iO(u y) + o2 62 -k2CCg}I

In this study, the ordinary differential equations are numerically solved by Gragg's method,

whose main algorithm for a differential equation Q0(x) =f(x, 5x)) is given as

30

yi1 = _i- +hf(x-1, i-1)

Yj, = yi- +2hf(xi-,+jh,yj) j= 1,2,..,n-1

i (Yn + n- + 2hf(xiy,))/2

where h is a subgrid space defined as h = Ax/n as shown in Figure 4.1.

J WAVES
W V y axis
i- -------------

i+1

I Main Grid
S ub-Grid
Figure 4.1: Sub-Grid system for wave model.

If we obtain the truncation error of the above equations as a power series in h:

Error = ath2i
i=1

That is, it only contains even powers of h. If n is even, the result for n steps would be q5n, and

for n/2 steps n/2. Combining both solutions:

4 0
3

31

where q5 is now fourth order accurate, but requiring about one and a half derivative

evaluations per step h instead of four.

The governing wave equation used, as well as the numerical solution used are capable

of solving the reflected wave field generated by structures posed along the x direction, that

is, the reflected potential propagates in the direction of positive x like the incident does.

However, the numerical solution is limited in the sense that it does not include reflected

waves propagating in the offshore direction, like the ones created by a detached breakwater,

which would be posed along the y direction. The model is capable though, of computing the

reflected potential backwards. This reflected potential would be obtained as the conjugate

of the incident potential at the grid points immediately next to the structure. A solution of

this kind would work well in the case of seawalls or 'infinitely' long detached breakwaters,

that is, as long as the only potential propagating offshore is due to reflections and thus easy

to obtain from the incident potential. In the case of a finite detached breakwater though, the

diffraction effects produced by the structure and propagating offshore from its position are

impossible to simulate with the method used. Nevertheless, since the objective of this model

is to simulate the wave field in the nearshore area around a field of groins, the solution

adopted fits the goal. However, this does not mean that detached breakwaters and T-groins

cannot be included in the model, since the diffraction effects induced by such structures

downwave of their location will be accurately reproduced.

Once the wave potential is obtained the remaining characteristics of the wave field

are obtained as follows. The wave angle is calculated by

0=tan-l(Ky/Kx)

(4.2)

where

K =Im{ / K=DyS (4.3)

where, S = Kx = tan-'[Im()/Re(O)].

The wave height is calculated easily by

H =2 Re { (2}2+Im{- 2 (4.4)
g g

4.2 Stability Analysis

It is possible to roughly obtain a stability condition for the method used by

performing a von Neumann stability analysis on the Helmholtz equation (Panchang et al.,

1988),

V2q+k20 =0

which has a solution for constant k

~- e imxe iny

such that m2+n2=k2. The solution for b can be written as:

o(p,q) = e im(p-1)Ae in(q-1Ay = e in(q-l)Ay

Expressing the Helmholtz equation in finite difference form and using 0(p,q),

P-20 + p- 2cosnAy-2 +k20=
(Ax)2 P (Ay)2

33

Defining the amplification factor as G= p+i/p, the above equation gives

G2+BG+1 =0

where B=k2 Ax (2(Ax/Ay) sin a)2-2, and a=nAy/2. When the amplification factor satisfies

a quadratic equation such as the one above, it is convenient to use the formulas derived in

Cushman-Roisin (1984) for investigating stability. In this case those formulas yield the

following condition for stability (Panchang et al., 1988):

<1 or, -2 2. (Ay)2

The right inequality is always satisfied if kAx < 2, or

Ax L/tx (where L=2ix/k)

The left inequality is always satisfied if

Ay L/7

If the marching Gragg's method described in the preceding section is used with these

stability criteria, a stable scheme results. While the condition obtained for Ax is similar to

commonly obtained stability criteria for many partial differential equations, the condition for

Ay is unusual, and it has both advantages and disadvantages. Most important it explains

why a marching method can be used for a wave equation, but is inherently unstable for most

other commonly occurring elliptic problems. For instance, for the Poisson or Laplace

equations, where k2=0, or L-o, it is obvious that any finite grid size Ay will result in an

34

unstable scheme. In addition to making the scheme stable, this type of condition reduces the

number of computational points required, which is important when trying to model the wave

field in the nearshore area, which the area of interest will be usually larger in the y direction

than in the x direction.

There are two disadvantages as well: the grid points being spaced L/h length units

apart, information about depth variations between grid columns is not well represented in the

model, and the wave components in the y direction are not well resolved. The first

disadvantage is not a major drawback, since in coastal areas we do not expect rapid

variations of bathymetry along the y axis. The second one makes the numerical scheme

approximate, in that, as in the parabolic approximation, the dominant wave direction is still

x. This problem will be described latter when testing the model against analytical solutions.

4.3 Boundary Conditions

Radiation boundary conditions are defined at both upwave and downwave sides of

the computational domain. In addition, each grid row is divided in ng+1 subdomains, where

ng is the number of groins crossing that row (see Figure 4.2). Each one of these subdomains

will need lateral boundary conditions. The upwave side corresponds to either the lee side of

a groin or simply a virtual boundary, and the downwave side to either the front side of a groin

or again a virtual boundary.

4.3.1 Lateral Upwave and Downwave Boundary Conditions

Side boundaries of the computation area are generally virtual open boundaries

established for performing the numerical analysis. The presence of such virtual boundaries

should not affect the solution in the study area. In other numerical models of this kind it is

assumed that there is only one wave train crossing the side of the domain. A simple radiation

35

S Wave Direction o o o o olnternal boundary grid points
a M M a MVirtual boundary grid:points
0 a . ... .. . . . .... .

a ,a : .
a i
10 a

i
0 0 0
0 0 STRUCTURE
S0 300

M. 0 0 :.0 0 0 0 M
M. 0 0 o 0

S0 0 30 0 0 0 M0
i 0 0 0 a 0
S0 : 0 0 : 0 M

50nd: ...c........ .... -0 in t s .. 0 .
0 .0 0 0 0 0
o8 : o o I

W 0 0 0 0 0
M 0 0 :0 0 0 0 X
0 0 0 0 0 0

Column Number (j)
Figure 4.2: Location of grid points where the wave potential
will be solved using a boundary condition.

boundary condition is used in this case
m~- oo- .- o ..

S= ik, where k -ksin ksin (4.5)
ay o

where 0 represents the potential associated with the only wave field considered. However,

for the numerical solution to let reflected waves escape the computational domain we need

to modify this expression. We simply consider q5 the total potential, that is, the sum of

incident and reflected potentials, q5 = gri+ 5pR.

'9 + R = ik q5-ik'Y.R (4.6)
ay ay ay

Assuming k' =k we can rewrite the boundary condition as:

a = ik k, w(h- k,) = ik(2i ki- ) (4.7)
dy

36

This expression is only valid when the absolute value of the wave number in the y

direction is the same for both incident and reflected wave fields, which is the case for

problems involving reflecting structures posed in the x direction and with straight parallel

contours along the y direction. The incident potential will be obtained at the upwave lateral

boundary using Snell's law. The boundary condition can then be expressed in finite

difference form to relate the wave potential at the first and second columns of the grid. The

value of qb on j=1 will be known and thus we can obtain qb for j=0.

Using finite central differences on row i yields

= iky(25 1- )
Ay 2

The lateral boundaries are assumed to be located between the grid columns 0,1 and NY,

NY+1. Therefore, 0 values at the upper lateral boundary can be given as

S-iky2Ayq 1 +ikyAy/2
S+ t (4.8)
1-ikAy/2 1 -iky/2 (4.8)

on row i. The same procedure is used for lower boundary, the only difference is that there

is no reflected wave field escaping the domain through this boundary, thus we can substitute

q5, for q5 in Eq. 4.7 and using finite differences yields

l+ik Ay/2
1NY 1 i--k V (4.9)
1 -ik Ay/2

37

The two boundary conditions can also be further approximated using higher order terms from

the Taylor expansion of the y derivatives.

4.3.2 Boundary Conditions for Groins

We will derive the boundary conditions for the general case when waves are incident

on both sides of the groin and from different directions. The groin itself will separate two

regions, (I) and (II), see Figure 4.3.

B

Groin /
\ /

\ //d y axis

Region (I) \ Region (II)
/ N\
/ /

Figure 4.3: Description of wave field around a groin.

We will consider Region (I) the upwave region, with incident and reflected wave field, and

Region (II) the downwave region, with transmitted and diffracted wave field.

The total wave potentials for those regions are, respectively,

Region (): = i + r +
Region (II): OI1) = + 0r+ t

where the velocity potential components can be expressed as:

i= ae ik(ycosa+xsina)
b = R a e i[k(-ycosa+xsina)+6r ra ik(-ycosa+xsina)
t Ta e i[k(ycosa+xsina)+ 6] ta e ik(ycosa+xsinq)

Assuming the depth in regions (I) and (I) is the same, thus, the wave numbers are the same

too, then:

-i = Q e ik(-ycosa'+xsina')
- r =R'a' ei[k(ycosa'+xsina')+6'r] P 'ra' ik(ycosa'+xsina')
It = T'a'e i[k(-ycosa'+xsina')+65't] '' ik(-ycosa'+xsina')
(/0 t T'1a'e' A0ta e i(yoo'xit'

where:

fr = Red, t = Te ', fr = R'ei,

l' = T'ei ',

Now we can obtain the wave potential and its derivative in the y direction at y=O and y=B.

) |y= = (1 + lr)aeikxsina + fta 'eikxsina'

) I y- = ikcosa(1 -,)ae ikxsin- i kcos a''ea eikxsin
ay ,r a

fI) Iy=B = (1/ y'+ 'r y')a 'e ikxsina'+ rae ikxsina

n0) y = ikcosa'(-l/y'+ 'ry')a'eikxin +ikcosa rtaeikxsina
9y
Where: y= e ikBcosa = e ikBcos'

We can rewrite the above expressions as:

= (1 +P)Y+ tY'

S= ikcosa(1 -Pr)Y-ikcosa'P',Y'

) = (1/y' + f'r')Y' + tY

S= ikcosa'(-1/y'+,',y')Y' +ikcosaypty

Where: Y = aek ixina, Y' =a'e ikxsina'

Now we can combine equations (4.10) and (4.11):

(5.10)xik :
(5.11):

(5.10)xik'y +(5.11):

From equations (4.12) and (4.13):

(5.12)xiky :

(5.13):

(5.12)xiky -(5.13):

ik' ") = ikfryY'+ik (1/y' +P'rY)'
(= ik yY+ik (-l/y' +P'y')Y'

ikyq)l y = [iky(l/r' +P'ry') +ik',(l/y' -fl',r)Y'

Combining (4.14) and (4.15) to obtain:

S- (1 +prk) [i ) + t 0[ik 1(-
( =ik'y(l +Pr)+ikO(l -r) iky(1/y'+p,'ry')+ik' y(1/y'-f'lry') yB

Rearranging terms in the equation above:

(5.10)

(5.11)

(5.12)

(5.13)

ik'I) = ik' (l+r)Y+ik''tY'
= ik(1 -fr)Y-ik' B'tY'

ik' + '= [i(k'(1+r)+iky(l--fr)]Y

(5.14)

(5.15)

__(1+/
iky(1-pr) lp ,) +
ik'Y(l+fr) +iky(l1-r) ik' (1+fr) +iky(1+r)
i k 't .(I) P

(5.16)

AII)

Combining (4.14) and (4.15) to obtain:

Sik'y(1 +fr)+iky(1-fr) [ik (ll/'+P'ry')+ik' '(1y'-f'r') [kB -yB

Rearranging terms:

ik'y tY ,(I)+ _
ik'(l+fr)+iky(l- r) 0 ik' (l +r)+ik (1+Pr) (5
i k(1/ -fl'r') (1/y 'r') =O(5.17)
iky(1/y'+P,'ry') +ik' y(l/y'-',ry') i k(1/y'+fplry) +ik' (1/y'-fl,',') ryB

If we assume a' = 900 for pure diffraction waves in downwave side, then,

k' =0 and f'r ='t=0

and substituting in Eqs. (4.16) and (4.17):

(1-Prl)
) ik 0)
qyO Y, (1 +fr)

() tY (If)
S(1 -r)yO

(5.18)

(5.19)

However, equations (4.16) and (4.17) are to be used in the general case when

a' 900. For a single groin, we can assume that the incident potential in Region (II) is only

ik (1/ +P'rY') + i k(1/ry'-P',)

- ik(1/y'+p', ry) +ik'y(1/'- ', yB =

41

due to the diffracted wave field. The wave fronts corresponding to it are perpendicular to the

groin, which means that we can in fact assume a'=900. In the case of a groin field, this

assumption may create problems, since the reflected waves from a groin 'downwave' may

be part of the incident potential in Region (II) for the previous groin in the upwave direction.

In that case a' will be different than 900.

Equations (4.18) and (4.19) will now be expressed in finite difference form. The

upwave groin boundary condition will be located between grid columns m and m+1, while

Eq. (4.19) will be located between columns n and n+l, the groin width, B, is then equal to

Ayx(n-m), (see Figure 4.4).

y axis

Upwave Downwave
S I yaxis

I I
] j, Groin

Upwave Downwave

I I B I "I
I I I I

m m+1 n n+1
mn m+ n +1
Column Index
Figure 4.4: Schematic diagram.

(1 -fr) Ay
1 +ik
(1+B) 2
n+1 = y(1 +) (5.20)
1 -ik(1 ------) Ay
Y(1 +P) 2

To obtain bf we will use (0,, and 0m, which are already known since the upwave boundary

of the groin has already been solved.

+1n 11 m
Ay (1-fP) Ay

Pr t ml- ^m
6t1 y -0 y- (5.21)
(1-6r) Ay

If the complete boundary equations are to be used, it will be necessary to solve a system of

equations involving unknowns ,m,,l and 0,, as well as ky which has to be assumed to be

equal to ky. This assumption is necessary when solving the wave field numerically. This is

because the numerical solution gives the total wave potential, which in the lee of a groin

could be the combination of transmitted, diffracted, incident and reflected potentials,

consequently, in general we do not know in what direction the wave field represented by qi'

is propagating,

Appendix A includes the theoretical derivation of transmission and reflection

coefficients for a permeable structure. The case solved represents the simple situation in

which there is no other potential in region (II) than the transmitted one. But the coefficients

obtained in that simple situation are also valid when solving the general case that includes

incident potential on both regions. This is true when assuming linear wave theory, the

general case is just the superposition of two simple cases like the one solved on Appendix

A. The coefficients can be calculated independently for each incident wave field.

4.4 Description of Wave Model

This section describes the numerical model that has been constructed to calculate the

wave field within a given domain. The domain will be usually bordered by a shoreline along

one boundary and an offshore region along the opposite boundary, although a completely

43

'open' region is also possible. The set of governing equations is solved using a finite

difference method which requires that a computational grid system be constructed.

Definition of the grid scheme is shown in Figure 4.5. A grid mesh comprised of

(NX+l)x(NY+1) rectangular grid cells of constant size (Ax by Ay) will be considered.

i= 0 i I i i

i=2
wave Direction

i=NX
i= 3 . .. .. ... . . .. . .. . .. .. .. .. . ... . . . ... .

i ...... .. . .. ...... .... ... ....... . . ....

i x. .. . .. ... . '. : .. .. . .. .. .. .. .

i = N X i i i .. . . . . .. . .. . .. .
i=NX+1
j=- j=l j=2 j=3 j=NY-1 j=NY j=NY+1
y axis
Figure 4.5: Computational domain of the model.

Specification of the wave parameters at the offshore grid row, height (H), angle (0)

and period (T) is the only wave input required for the model when no structures are present.

In the case when one or more structures are situated in the area of interest, their dimensions,

location as well as physical characteristics, porosity (e) inertial coefficient (S) and friction

or damping coefficient (f) will have to be given. The model will calculate the corresponding

reflection (R) and transmission (T) coefficients depending on the characteristics of the

impinging wave train. It is also possible to directly assign the values of the reflection and

transmission coefficients if known in advance.

Based on the given wave period and the specified depth matrix describing the

bathymetry in the area of interest, the wave number, k, for each grid point is computed using

the dispersion relationship. The dispersion relationship already includes the effects of a

possible current, which can be also given in a matrix form for every grid point. With the

44

computed wave number and period we can calculate the phase speed (C) and wave group

speed (Cg). Once all these values are known, the subroutine containing the numerical

solution to the finite difference equation governing wave propagation can be called. The

subroutine contains an algorithm that ultimately calculates the complex steady free surface

potential (qb) at each grid point and row by row, starting on row one and 'marching' towards

row NX+1. Each row will have two boundary grid points at j=0 and j=NX+1, the

appropriate boundary conditions are then solved relating the wave potential at those grid

points to the already solved at j= 1 andj=NX.

When the domain to be solved includes structures the method of solution varies

slightly, new boundary conditions in the numerical 'subdomains' are necessary, the

expressions governing those boundary conditions as well as the external ones were presented

earlier in this Chapter. The so called subdomains are simply the two divisions that a structure

crossing a row will create in that row. The problem then becomes a matter of bookkeeping

and adequate indexing. When the row being solved is limited by a structure on one or both

of its sides, a call to a subroutine that computes the reflection and transmission coefficients

used in the boundary conditions is made. The subroutine not only calculates R and Tbut also

gives the wave number inside the porous structure, the real part of that wave number

corresponds to the propagating wave mode inside the structure, from it we can obtain a phase

speed (Cpor) and wave group celerity (Cgp,) inside the structure which will be also necessary

when solving the governing wave propagation equation on the boundary grid points.

During the first call to the wave propagation subroutine the model is run over every

groin in the domain but without considering breaking, once the values of the wave height (H)

and angle (0) are obtained at every grid point from the complex free surface potential, the

45

breaking and surf zone subroutine is called. This subroutine numerically applies the surfzone

model derived in Chapter 5. The subroutine first checks every grid point for a simple

breaking condition relating water depth and wave height, when the condition is fulfilled the

model keeps track of the grid point were it occurred, thus obtaining a 'breaking line' that

connects all those grid points. After the model has located the grid points were the waves

break, for every column the model obtains a modified group celerity at every grid point with

an index greater than the index where the waves break for that column. The modified group

celerity (Cg*) includes the dissipative effects present in the surf zone.

Once the breaker line and the modified wave group celerity are determined the wave

propagation subroutine is called once again, the waves are propagated in exactly the same

fashion as in the first call, with the only difference that when the grid point being solved is

beyond the breaker line, the modified group celerity is used instead of the one obtained from

the dispersion relationship.

After the second run of the wave propagation algorithm is over, wave angles and

wave heights are obtained at every grid point.

4.5 Testing of Wave Model

The model will be tested against theoretical and experimental results for problems

that include wave shoaling, refraction, diffraction, reflection, wave-current-interaction, etc.

4.5.1 Wave shoaling and refraction

In shallow water waves transform over a sloping bottom, and if incident waves are

normal to a beach with straight and parallel contours, the change in wave height is solely due

to the change in water depth. This transformation is called wave shoaling. The wave celerity

also changes due to the change in water depth or due to a current field. The gradient of wave

46

celerity along the wave crest results in a modification in wave direction. This process is

known as wave refraction. To solve the wave shoaling process we make use of the energy

conservation equation, while to calculate the refraction we assume a gradually sloping

bottom, constant period and thus irrotationality of the wave number. In the particular case

of a simple bottom topography in which the contours are parallel to the y axis the

irrotationality of the wave number yields Snell's Law.

Numerical results from the wave model were compared to the analytical solution

based on energy conservation and Snell's law. The input data are given as follows:

NX NY Ax (m) Ay (m) Hi (m) T (sec) Slope

60 20 .1524 .3048 .04 0.8 1:20

In order to test how the model reproduces shoaling without refraction, the wave

direction is given as perpendicular to the beach. Next, the waves are run at an angle with

respect to the depth contours, and from the results we can obtain, dividing by the previously

calculated shoaling coefficient, a refraction coefficient.

As shown on Figures 4.6 to 4.8, the model produces results that agree well with the

theory for both the refraction and shoaling coefficients. However, if the incident angle is

increased too much, due to the nature of the numerical method applied, the refraction

coefficient differs considerably from the exact solution: this problem is shown on Figure 4.8.

This might induce errors at the boundaries as well. In the numerical scheme, the boundary

condition is based on Snell's law, which will yield the correct incident wave angle along the

boundary. However, due to numerical error, the wave potential escaping the domain might

propagates in a slightly different direction when the angle becomes large; this will induce

unrealistic values of the wave height and direction at the boundary.

47

SHOALING TEST

Theory

00000 Numerical Model
0

- ,
-^ o^^ ^ ^

0U.8
0.6 0.8 1 1.2 1.4
kh
Figure 4.6: Comparison of wave shoaling.

1.6 1.8 2

REFRACTION TEST

1.15

1.1

1.05

1

0.95

0.9

0.85

0.6 0.8 1 1.2 1.4 1.6
kh
Figure 4.7: Comparison of wave refraction, 60= 20.

1.15

1.1

1.05

0.95

0.9

0.85

SHOALING plus REFRACTION

1.15 Theory

1.1 00000 Numerical Model

1.05

0.95

0.9

0.85

0.8 '
0.6 0.8 1 1.2 1.4 1.6
kh
Figure 4.8: Comparison of wave shoaling/refraction, 6= 20.

1

0.95

0.9

0.85

0.8

REFRACTION TEST

Theory

S00 00 Numerical Model

000000000

0.6 0.8 1 1.2 1.4
kh
Figure 4.9: Comparison of wave refraction, 0 = 300.

4.5.2 Wave Diffraction and Reflection

Diffraction is a phenomenon by which wave energy diffuses or flows transversely to

the main direction of wave propagation. This allows the waves to propagate into sheltered

regions like the lee of a breakwater, jetty or groin.

Penney and Price (1952) derived a solution for constant water depth based on the

Helmholtz equation and similar to the Sommerfeld theory of the diffraction of light. The

expression they obtained governs diffraction of small amplitude waves around a semi-infinite

thin breakwater and propagating through water of uniform depth. It has been already

confirmed that this analytical solution agrees well with hydraulic model experimental

results. The diffraction coefficient, H/Hi, that is, the ratio of diffracted to incident wave

amplitudes, is given as follows in terms of the polar coordinates r and a (Figure 4.10):

Kd=H/H= -I r sin-a e ikrcos(-)I r sin -ikrcos(8+a) (5.22)
=H/ I sin e sin
d 7t 2 7U 2

where

I() i e 2 dA (4.23)

Equation (4.22) includes both the effects of diffraction and reflection from the structure

posed along the x direction.

Equation. (4.23) can also be expressed as

I 1 +C(A) +S(A) C() -S() (4.24)
2 2

50

in which the Fresnel integrals, C(/) ans S(A) are defined as

2 22
C() fcos 2 dA, S(A) =sin2 dA (4.25)
0 0

/

/ sStructure
r

(x,y) /

/Waves Y

Figure 4.10: Coordinate system for
analytical solution.

The model capability of handling wave diffraction was evaluated by comparing the

numerical results with the above analytical solution. Since the model is limited in terms of

the wave incident angle, two different structure layouts where used. In the first one (see

Figure 4.11) the given conditions were H,=l m, depth=5 m andoL =30.0 m, which

corresponds to a wave period, T, of 4.96 seconds. The grid used had dimensions NX= 60,

NY=60, Ax= 10.0 m and Ay= 15.0 m. The structure was arranged along the x axis, that is

the main axis of wave propagation. Therefore, the numerical solution included not only the

diffracted wave field behind the structure, but the combination of reflected, diffracted and

incident wave fields in front of it.

6-
Sec. 3
S9-

Sec. 4
12 ............ .
-STRUCTURE
15 -

Sec. 5

0 3 6 9 12 15 18 21 24 27
y/Lo
Figure 4.11: Grid and structure layout used for first diffraction
test. The values of x and y refer to model coordinates.

Figure 4.12 shows the results for two different incident angles, 5' and 150. Figures

4.13 and 4.14 include comparisons between theory and the numerical solution along five

different cross sections. In all the aforementioned figures x and y as well as angle 0 refer to

the coordinate system for the analytical solution. The results seem to compare well for these

small angles of incidence, of special significance is the agreement obtained in the region

upwave of the structure were reflected and diffracted wave fields interact, since the final

objective of this model is to simulate the wave field around a field of groins. Such groins

would be posed along the x axis and the reflected wave field that they generate would greatly

affect the resulting nearshore circulation.

In the upwave region (Figure 4.14, Sec. 1) the amplification factor (H/Ho) actually

becomes larger than two, this is due to the diffraction effects that are superimposed over the

reflected and incident waves. Such effects will tend to disappear away from the tip of the

breakwater.

8 7 6 5 4 3 2 1 0 -1 -2 -3 -4 -5 -6 -7 -8
y/Lo

-18
8 7 6 5 4 3 2 1 0 -1 -2 -3 -4 -5 -6 -7 -8
y/Lo
Figure 4.12: Comparison of wave diffraction/reflection coefficient. Top 0=
5, bottom 0= 15.

2.5
2 L

-Thery Sec.
0.5 C. OO0000 Numerical Model

-18 -16 -14 -12 -10 -8 -6 -4 -2 0 1
x/Lo

2.5
2
S 1.5

0.5

-18 -16 -14 -12 -10 -8 -6 -4 -2
x/Lo

8 7 6 5 4 3 2 1 0-1-2-3-4-5-6-7-8
y/Lo

2.5
2

0.5
n

2.5
2
1.5

0.5

8 7 6 5 4 3 2 1 0-1-2-3-4-5-6-7-8
y/Lo

8 7 6 5 4 3 2 1 0-1-2-3-4-5-6-7-8
y/Lo

Figure 4.13: Comparison of wave diffraction/reflection coefficient along five different
sections, O=5.

Sec. 2

_______ ^^^^aaWRffffiSBBEe^8
30oarcx~oniociX~oaaijjjLu u^

Sec. 4

I 1

I

eeeBaeeBBBea~PE

S1.5-
I 1 Theory
0.5 0 Numerical Model ec 1
0.
-18 -16 -14 -12 -10 -8 -6 -4 -2 0 1
x/Lo

2.5
2-
1.5 Sec. 2

0.5
0
-18 -16 -14 -12 -10 -8 -6 -4 -2 0 1
x/Lo

2.5
2 Sec. 3
1.5-

0.5
0
8 7 6 5 4 3 2 1 0-1-2-3-4-5-6-7-8
y/Lo

2.5
2 Sec. 4

1.5

0.
8 7 6 5 4 3 2 1 0-1-2-3-4-5-6-7-8
y/Lo

2.5 ... . .
2 0 Sec. 5

0.5
0.5- 0 Aoaoo- 5^

8 7 6 5 4 3 2 1 0-1-2-3-4-5-6-7-8
y/Lo

Figure 4.14: Comparison of wave diffraction/reflection coefficient along five different
sections, 0= 15.

55

For larger angles of incidence the model starts to give inaccurate results, this problem

is illustrated in Figure 4.15. The extension of the diffraction zone seems to fall short of the

theoretical solution, the same can be said about the reflection region on the other side of the

groin. The variations of wave height in the latter region along a section perpendicular to the

groin increase with increasing angles of incidence. The model is not capable of simulating

these variations since there is a limit to the size of the grid the y direction, Ay.

In order to simulate cases with large angles a different layout is used in which the

breakwater is posed along the y axis. The input conditions are the same as for the first layout

with the exception of Ax which is equal to 5 meters now. Figure 4.16 shows the comparison

for waves approaching a breakwater at 900, the model seems to agree reasonably well with

the theory for this case. Figure 4.17 shows the results when the waves approach the

breakwater at a 70 degree angle and directed away from the breakwater. The results for this

case are not as good, most likely due to the limitations of the numerical model for higher

angles of incidence. The results can be improved by reducing Ay, but as explained earlier

in this chapter, there is a numerical limit to that reduction. It is noticeable how the

diffraction bar that starts at the tip of the breakwater is higher than what the theory predicts.

This effect is not observed when the structure is posed in the y direction as in the first layout

used. As explained in the first part of this chapter, the numerical solution is based on a

marching method along the x axis, this means that diffraction and reflection effects that

propagate in the -x direction from the body and tip of the breakwater are not obtained.

2.5
2

1.5

0.5
I

2.

0.

-18 -16 -14 -12 -10 -8 -6 -4 -2 0 1
x/Lo

5
2
5 Sec. 2

1
5
A ... .,

-18 -16 -14 -12 -10 -8 -6 -4 -2 0 1
x/Lo

1.5 -
i.5 ,

0.5-

8 7 6 5 4 3 2 1 0-1-2-3-4-5-6-7-8
y/Lo

8 7 6 5 4 3 2 1 0-1-2-3-4-5-6-7-8
y/Lo

2.5
2
1.5

0.5

8 7 6 5 4 3 2 1 0-1-2-3-4-5-6-7-8
y/Lo

Figure 4.15: Comparison of wave diffraction/reflection coefficient along five different
sections, 0=30.

STheory Md
0 Numerical Model ec. 1
-- .

Sec. 3

0 000 0

10

57

DIFFRACTION COEFFICIENT

5 4 3 2 1 0 -1 -2 -3 -4 -5
y/L

WAVE FRONTS

5 4 3 2 1 0 -1 -2 -3 -4 -5
y/L

Figure 4.16: Comparison of wave diffraction coefficient, 0= 90.

DIFFRACTION COEFFICIENT

5 4 3 2 1 0 -1 -2 -3 -4 -5
y/L

WAVE FRONTS

5 4 3 2 1 0 -1 -2 -3 -4 -5
y/L

Figure 4.17: Comparison of wave diffraction coefficient, 0=70.

4.5.3 Wave-Current interaction

When surface waves of any kind propagate over the surface of a medium in steady

but non-uniform motion, they tend to undergo changes in length, direction and amplitude.

In the case of water waves propagating over a current field the common assumption was, at

first, that the wave energy is simply propagated with a velocity equal to (U + Cg), where Cg

is the vector group velocity and U is the total stream velocity, and that no coupling between

the waves and the current took place. However, as described in Chapter 2, waves are

modified to a much greater extent than would be predicted using that assumption when riding

through a non-uniform wave field.

Collinear Current Shearing Current
y axis y axis
Waves Waves

I

I .

Current Current

Figure 4.18: Conditions of collinear and shearing current.

The performance of the model in terms of wave-current interaction is tested through

two different cases, one with collinear wave and current directions, and the other with a

shearing current that produces refraction of the wave field. Figure 4.18 describes the two

60

situations tested. The analytic solution to both cases is given by Longuet-Higgins and

Stewart (1961). The comparisons with their solutions are shown on Figures 4.19-21. The

given conditions are Ho =0.1 m, depth=3 m and T= 1 sec. In the case of the shearing current

the initial angle between the wave fronts and the current is 30 degrees.

The grid used had dimensions NX= 101, NY=21, Ax=0.1 m and Ay=0.6 m. The

model solution for the problem of collinear wave-current compares very well with the theory.

COLLINEAR CURRENT(Amplification)
3

2.5 T=1 sec
depth=3 m
2
S- Theory
1.5 oooo Numerical Model

0.5

0
0 -----------------------
-0.2 0 0.2 0.4 0.6 0.8
U/Co
Figure 4.19: Comparison of collinear wave-current

interaction.

SHEARING CURRENT (Amplification)

1.

0

Theory
0.5 o o o o Numerical Model -

-0.2 -0.1 0 0.1 0.2 0.3
U/Co
Figure 4.20: Comparison of shearing wave-current
interaction, amplification of wave height.

T=1 sec
5 depth=3 m
Initial Angle=30 degrees

I _______________ENO_

61

SHEARING CURRENT (Angle change)
50

45- T=1 sec
depth=3 m
40 Initial Angle=30 degrees

CD 35 - Theory
oo 0 0 Num. Model

25

20
2 0 .....-' ------ --- --- -
-0.2 -0.1 0 0.1 0.2 0.3
U/Co
Figure 21: Comparison of shearing wave-current interaction,
change in wave angle.

In the case of the shearing current, the model predicts adequately the amplification

factor, but the change in wave angle deviates considerably from the theory, especially for

high positive currents.

CHAPTER 5
SURF ZONE MODEL

The flow properties in the surf zone are extremely complex due to the strong

interactions among motions induced by waves, currents, and turbulence. The present

knowledge on surf zone dynamics is still inadequate and most of the models are rather

rudimentary. Most of the available models that describe wave breaking are based on an

approximation to wave energy conservation. These models can be classified into two groups:

one based on the similarity of the wave breaking process with other hydraulic phenomena

such as a hydraulic jump (Dally et al., 1984), a tidal bore (Battjes and Janssen, 1978), etc.,

and the other is based on estimation of the eddy viscosity (Mizuguchi, 1980) or turbulence

(Izumiya and Horikawa, 1984).

In this chapter a simple surf zone model proposed by Lee and Wang (1993) is

presented. The model is based on the consideration of wave energy balance and wave action

conservation so that the wave-current interaction is fully taken into account. The model is

presented in analytical form for the case of two dimensional gradually sloped bottoms.

5.1 Time Averaged Wave Energy Equation in the Surf Zone

It is assumed here that the surf zone maintains a wave-like periodic motion that is

quasi-stationary when time-averaged over a wave period. The turbulent motions occur on a

much smaller time scale and thus can be treated as dissipative forces represented by turbulent

friction terms which include the eddy viscosity. Therefore, we can use the Navier Stokes'

equation to describe the flow in the surf zone:

au+V -q2+ +gz =(VxU)xU+pV2xU (5.1)
at 2 p

where p is the kinematic viscosity coefficient. Let the surface displacement and the velocity

vector, U(u,v,w), be decomposed into mean value, wave and turbulent fluctuations, which

are distinguished by subscripts c and w, and prime respectively, thus,

7 = + 77' = 77, + +77' (5.2)

U= U+U'= Uc +U +U' (5.3)

where the superscript ~ is used to denote turbulent averaging. After turbulent averaging, Eq.

(5.1) becomes

f+ V i-+ + g 4 =(VxU)xU + V2xU (5.4)
at 2 p

where v is the total viscosity including the eddy viscosity due to turbulence. The subscript

~is omitted thereafter.

Taking the scalar product of U(u,v,w) and the respective terms in Eq. (5.4), and

summing the products gives the energy equation with a dissipation term:

Sf 9 +V'[U(q- + +gz)] dz= -U dz (5.5)
-h t -h

by applying the kinematic boundary conditions at the free surface and the bottom,

at
wr U'Vr =0O

w h + UVhh = 0

and setting p equal to zero at the surface, we obtain,

dr 2 17 2
S )dz + gh + [U( + +gz)]dz- fvUV2Udz (5.6)
tJ 2 at 2 p
-h -h h

using the Leibnitz' rule of integration. This equation is similar to the wave energy equation

given by Eq. (2.26) with the additional dissipative term. This volumetric dissipative term is

replaced here by equivalent energy flux term introducing a head loss term in the context of

Bernoulli sum, i.e.,

D= vUV2Udz = Vhf glUdz (5.7)
-h -h

where I is defined as the head loss due to turbulence. Eq. (5.6) can then be expressed as

S(q )dz + gh + Vhf [U( + +gz +gl)]dz =O (5.8)
t 2 t h 2 p
-h -h

Following the same procedure used in Chapter 3 by taking time average over the

wave period, an energy equation similar to that of Eq. (3.38) can be obtained,

6E+ Vh-(UhE) = 0 (5.9)
(t (5.9)
dt

65

in which the transport velocity, Uh, is the counterpart of (Cg+U) in the non dissipative case

and can be represented by

Uh=Cg+U+CgD (5.10)

The first two terms in (5.10) constitute the transport velocity due to non-dissipative forces,

whereas the last term manifests the effect due to the dissipative forces. The latter, in general,

should be negative indicating a reduced energy flux due to dissipation. In theory, it can be

estimated from the time-averaged energy dissipation term as follows:

Vh (CgDE)=D (5.11)

where D is equal to the time average of D.

5.2 Wave Action Equation in the Surf Zone

In this section, the wave action equation given in Eq. (2.41) is shown to be also valid

even in the presence of strong turbulence. It is assumed here that the dynamic free surface

boundary condition given in Eq. (2.10) is still valid with the inclusion of a head loss term.

Following the approach by Kirby (1983), a virtual work term proportional to W(DJ/Dt) is

introduced to represent the head loss, where Wis a positive undefined coefficient indicating

the strength of the dissipation. The dynamic free surface boundary condition then becomes

= (1 + W) [(w)t +Vh/.Vh (5.12)
g

which, after substituting Eqs. (2.12) and (2.13) into it yields

aei (1+ W) [AeiP+iei' {A +U sVA=0}]
g at

66

where ad is again the intrinsic frequency including the diffraction effects. Following the

same reasoning as in Chapter 2 we obtain

A=-g-- where Oa,=(1+W)(w<-U,-k). (5.13)
D
aA
A+ U+sVA =0 (5.14)
9t s

The subscript, s, denotes the mean water surface level. The kinematic free surface boundary

condition also yields

o =(1 + W)gk tanhk(h + 7) (5.15)

Ba
+V-(Usa)=0 (5.16)
dt

Combining Eq.(5.14) and Eq. (5.16), we obtain the following wave action equation which

is valid in the surf zone:

(pg H pg H2 (5.17)
at 8 OD 8 D

This equation enables us to estimate the surface velocity in the surf zone once the wave

height decay rate can be established.

5.3 Wave Height Transformation in the surf zone

In this section an approach developed by Lee and Wang (1993) to wave breaking and

wave-current interaction in the surf zone is presented.

The wave energy equation given in Eq. (5.9), when expressed in terms of wave

height, can be written as,

d H) 2 -(g--- ) ++ h(CgCo)pgg+--] =0 (5.18)
at OD 8 OD 8

Assuming that the surf zone retains a quasi-steady state when integrated over a wave period,

then the slowly varying flow properties become time independent, and the absolute frequency

becomes constant. Therefore, Eq. (5.18) becomes

Vh [(Cg+J +Cgo)- 2] (5.19)
8 oD

Now we apply the wave action equation in the steady state, and the wave energy equation in

the surf zone is reduced to

,og H2
Vh [(Cg +CgD)- f ] =0 (5.20)
8 oD

We will assume that the dissipative term in equation (5.21) is proportional to the group

velocity at the breaking point, Cgb,

CgD = -Cgb (5.21)

where /f is a positive coefficient. Substituting in Eq. (5.20)

h [Cg* H]=O (5.22)
8 D

where the Cg* is the group velocity in the surf zone and is calculated by

Cg*=Cg Cgb (5.23)

This energy transformation model given in its final form by Eq. (5.22) has only one unknown

coefficient, /p. The model is applicable to three dimensional bathymetry and any incident

wave angles.

In the case of two-dimensional bathymetry of uniform slope, it is possible to obtain

an analytical expression that describes the wave height transformation in the surf zone.

Considering that x is normal to the beach and directed onshore, Eq. (5.22) becomes

(C Cg- Cgb) =0 (5.24)
ax OD

The x component of the above equation results in,

H 2 (5.25)
cos (f/Cg, Cg)

Using the dynamic free surface boundary condition given in Eq. (5.13),which states that oD

is proportional to the wave height in the surf zone, equation (5.25) can be written as

gH
H H (5.26)
cosO(fpCg, Cg)

with pf determined later. The wave angle 0 can be obtained using Snell's law as

0= sin-'(Ca sin 0o/Co)

(5.27)

69

where Ca,= /k is an absolute phase speed. Equation (5.26) is adimensionalized as

H 1 PH
Hb HbCgb cosO( -Cg)

where Cg'=Cg/Cgb. After applying boundary conditions at breaking point, H/fl =1 and

d/db=l, and H/Hb=H', where Cg becomes zero, we obtain

H H's
Ss (5.29)
Hb cos[1 -Cg'(1 -H's/cosOb)]

Now, using a shallow water approximation to Cg',Eq. (5.29) yields

H H'S
=- (5.30)
Hb cosO[1 -P(l -H,/cosb)] (5.30)

Where d'=d/db. Only one of the three parameters H's, PH and P is independent. If the value

of H', is obtained from an experiment, the other two are solved by,

cos 0
cosOb -H', (5.31)

H COOb H's HbC (5.32)
ScosOb -H'

Lee and Wang (1993) showed good agreement between the present theory and the laboratory

experiments by Horikawa and Kuo (1966). The values of H' are evaluated from the data,

they depend on the slope of the beach face and can be closely approximated by (tan a)2

70

where a is the slope of the beach. This idea coincides with the fact that wave height decay

across the surf zone is strongly influenced by the beach slope.

In the numerical wave propagation model in which this theory is to be included Eq.

(5.22) will be used since full three dimensional cases will be solved. The model treats the

value of Cg-Cgb as the new group velocity inside the surf zone and solves the wave

propagation equation using it. Figures 5.2 to 5.3 compare the numerical results with the

theoretical curves obtained with Eq. (5.30) for a 1 to 20 slope with H,'=0.22, that is, /f= 1.28.

The figures show how the numerical results oscillate around the analytical solution,

and the amplitude of the oscillations increases with increasing wave period. The author has

compared the solutions from other numerical wave propagation models that include this surf

zone theory and the deviations between the analytical and numerical solutions were smaller.

Therefore it is reasonable to conclude that the fluctuation of the wave height after breaking

is entirely due to the numerical scheme used.

BREAKING TEST 1

d/db

BREAKING TEST 2

d/db

BREAKING TEST 3

0 0.2 0.4 0.6 0.8 1
d/db
Figure 5.1: Comparison between breaking analytical and
numerical results. Test 1: T=3 s, Test 2: T=5 s, Test 3: T=7 s.

CHAPTER 6
MODEL APPLICATIONS

In the previous chapters a numerical model describing wave propagation in the

nearshore zone has been developed and tested against several analytical solutions. The

model will be now used to simulate the wave field in various cases where experimental data

are available in the literature.

The first section of this chapter evaluates the performance of the model for the case

when waves propagate over a paraboloidal shoal surrounded by a region of constant water

depth. In the second section, the model is compared with measurements of wave height

obtained in a laboratory experiment simulating the nearshore wave field around a detached

breakwater on a sloping beach. The third section concerns the wave field in the vicinity of

ajetty placed perpendicular to the shoreline, again against data from a laboratory experiment.

In the last section, the model application is illustrated by a field of groins in constant

and sloping bottoms. It will be shown how the model can simulate the wave field around

such structures allowing for different dimensions and permeabilities.

6.1 Wave Field over a Paraboloidal Shoal

Since the conventional wave refraction theory is based on the geometrical optic

approximation, it fails to predict wave height at and near caustics where ray intersections

occur. The effect of diffraction should be included in the analysis of waves in the vicinity

of ray convergence. The basic equation for the present model includes this effect. As a

demonstration, the model is applied to wave propagation over a submerged shoal with

73

concentric circular contours where the conventional refraction theory indicates the formation

of caustics. Hydraulic model experiments for this situation were conducted by Ito and

Tanimoto (1972). All tests were conducted for non breaking waves. The arrangement of the

shoal in the numerical computation is shown in Figure 6.1.

Cross Section

y/Lo
Plan View

8
0 1 2 3 4 5 6
y/Lo
Figure 6.1: Bathymetry for paraboloidal shoal
configuration, after Ito and Tanimoto (1972).

The water depth around the shoal is constant with a value of 15 cm, the top of the shoal is

at a depth of 5 cm. The initial wave length on the upwave row of the numerical is 40 cm, the

rest of the input conditions are given as follows

NX NY Ax (m) Ay (m) Hi (m) T (sec) 0 ()

63 17 .05 .15 .02 .5107 0

The results of the numerical calculation and the hydraulic model tests are presented

in Figure 6.2. Figure 6.3 shows three dimensional plots of both free surface and wave height.

It clearly shows how the wave fronts bend and wrap around the shoal. The maximum wave

height is generated right behind the shoal, where the wave rays would be crossing according

to conventional wave refraction theory.

3
2.5
2
1.5
S1"5

0.5
0

2 3 4 5 6
x/Lo

7 8

3
y/Lo

0 1 2 3 4 5 6
y/Lo
Figure 6.2: Wave height comparisons between computed and
experimental results by Ito and Tanimoto (1972).

SNumerical Model
000OMeasured 6 0 O 0-~ -~

Sec. 1

3
.5-
2 Sec. 30

0 0O
n C
.5- 00
0. 1

-
-

2.

0.

75

FREE SURFACE (ETA/ETAi)

x/Lo

WAVE HEIGHT (H/Hi)

6 0

x/Lo

Figure 6.3: Numerical
(1972).

model results for shoal configuration used by Ito and Tanimoto

y/Lo

y/Lo

76

6.2 Wave Field around a Detached Breakwater

A detached breakwater is also a common coastal structure. The numerical model is

tested here to predict the wave field, breaking and surf-zone decay around one of these

structures.

Watanabe and Maruyama (1986) performed a hydraulic model test that included a

detached breakwater made of a steel plate 2.67 m long and placed at a water depth of 6.0 cm,

and parallel to the shoreline on a plane beach made of mortar with a slope of 1/50. The

breakwater was designed to almost perfectly reflect the incident waves. The waves were

incident normal to the breakwater and the contour lines, and the deepwater wave height and

period were 2.0 cm and 1.2 sec., respectively. The depth at the offshore grid row of the

numerical domain was 12 cm, the rest of the input data were given as

NX NY Ax (m) Ay (m) Hi (m) T (sec) ()

60 60 0.10 .3814 .0202 1.2 0

The results of the numerical calculation for both free surface and wave height are

presented in Figure 6.4. The location of the computed breaker line is compared with that

obtained in the experiment in Figure 6.5. The overall agreement seems to be fairly good,

with the exception of the area located in the middle, behind the breakwater, where the

predicted breaker line is closer to shore than the measurement. Two explanations can be

given for this difference. One related to the fact that the computed results do not include

nearshore circulation, which in this area would be directed offshore against the incoming

waves, thus, the waves break earlier. A second reason could be related to the model's partial

inability to fully simulate waves that do not propagate mainly in the x direction.

Comparisons of cross-shore distributions of wave height are shown in Figure 6.6. The

77

results compare reasonably well, except the area affected by the reflected waves from the

breakwater. As explained in the last chapter, the model is not capable of simulating

diffracted waves propagating in the negative x direction. Although it would be possible to

include the reflected wave field from the breakwater, in the present model is not included

since usually our main interest for a detached breakwater is the wave condition on the down-

wave side.

6.3 Wave Field around a Groins and Jetties

Groins and jetties, in general, are shore perpendicular structures used on the coast for

multiple purposes. Sometimes jetties are constructed to keep a tidal inlet in place so as to

stabilize the inlet entrance. An important use of groins is as terminal structures for beach

nourishment projects, with the idea of maintaining a wide beach for a longer time by holding

the sand in place. They might also be used to simply slow down erosion along a natural

shoreline. In any of these cases, to accurately predict the effects of these structures on the

shoreline, the first step is to predict the wave field around them.

Watanabe and Maruyama (1986) constructed a hydraulic model to test this situation,

the basic layout of the test was the same as the one described in section 6.2, but instead of

a detached breakwater, the experiment set up included a jetty made of steel plate extended

to an offshore distance of 4 m, where the water depth was 8 cm. The deep water incident

wave angle is also different, 60 for this case. After applying Snell's law and calculating the

refraction and shoaling coefficients between deep water and the offshore grid row at a water

depth of 12 cm, the rest of the input conditions were given as

NX NY Ax (m) Ay (m) Hi (m) T (sec) (o)

60 60 0.1 0.4 .0143 1.2 28.3

78

FREE SURFACE ETA

-3 -2 -1 0 1 2 3 4
Distance from center of breakwater y (m)

WAVE HEIGHT H (cm)

-4 -3 -2 -1 0 1 2 3 4
Distance from center of breakwater y (m)

Figure 6.4: Numerical model results for detached breakwater test.

-4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4
Distance from center of breakwater y (m)
Figure 6.5: Comparison of location of breaker line behind a
detached breakwater, experimental data after Watanabe and
Maruyama (1986).

4- x=4 m Numerical Model
3 0 0 00 Measured

,I Q ^ , 0

1____

-4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2
Distance from center of breakwater y (m)

2.5 3 3.5 4

3 3.5
Distance onshore x (m)

5
4 y=2 m
g3
o3 0000 0

0
1 1.5 2 2.5 3 3.5 4 4.5 5 5.5
Distance onshore x (m)
Figure 6.6: Comparison of distributions of wave height around a detached
breakwater, experimental data after Watanabe and Maruyama (1986).

80

The computed results for both free surface and wave height are presented in Figure

6.7. The location of the computed breaker line is compared with that obtained in the

experiment in Figure 6.8. The overall agreement is acceptable, although there is some

discrepancy, in particular, near the jetty on both sides. Part of this deviation is attributable

to the effect of the nearshore current which has not been incorporated here. Figure 6.9 shows

comparisons of cross-shore distributions of wave height. They reveal that the model predicts

with sufficient accuracy the wave field around the jetty, especially taking into account that

the model is based on linear theory.

One of the limitations of the model can be identified studying Figure 6.9 in more

detail; the wave height along y=4.8 m does not seem to decrease as much as the data when

approaching the shoreline at x=6 m. This is thought to be due to two different reasons: (1)

the breaking criteria is simple and probably inaccurate when more than one wave field is

superimposed. Watanabe, Hara, and Horikawa (1984) presented a new breaking condition

and diagrams in which the main governing parameter is the ratio of orbital velocity at the

crest to the phase velocity, in place of the conventional depth to height ratio. Under this kind

of breaking criteria, waves would break in cases where the conventional depth limiting

breaking criteria has not been reached yet; (2) the surf zone model used, which affects the

wave height reduction rate after breaking might not be adequate for the case when incident

and reflected waves are combined.

81

FREE SURFACE (ETA)

0-

1

2

3

4
or

5

6
0

8 9 10

WAVE HEIGHT H (cm)

6 1 I I I I I I I
0 1 2 3 4 5 6 7 8 9 10
Distance alongshore y (m)

Figure 6.7: Numerical model results for jetty test.

1 2 3 4 5 6 7
Distance alongshore y (m)

82

Wave Direction
1.5
.- Numerical Model
2- 0 0 0 0 Measured

2.5
a
3-

S3.5 c

4- 0
0 0 O0 0 0O
4.50o0o o0 0 0 000 ooo0 000
5-0 0

5.5

1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9
Distance alongshore y (m)
Figure 6.8: Comparison of location of breaker line around a jetty,
experimental data after Watanabe and Maruyama (1986).

5 i 1 1 1 1 1 1 1

4- y=5.2 m Numerical Model
0 0 0 0 Measured
3-

S2 000000 0on
1 oon000000000000no o

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6
Distance onshore x (m)

4- y=4.8 m

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6
Distance onshore x (m)

4 y=7.0 m

E 3 00ooO0
2 nOnoo00000000

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6
Distance onshore x (m)

Figure 6.9: Comparison of cross-shore distributions of wave height
around a jetty, experimental data after Watanabe and Maruyama (1986).

83

6.4 Wave Field around Permeable Groins

The main objective of this thesis was to obtain a wave propagation model that could

simulate how a field of permeable groins alters the wave propagation in the nearshore region.

The basic numerical solution and appropriate boundary conditions were described in Chapter

5. In this section it will be shown how it is possible to assign different physical

characteristics to a groin and obtain the corresponding wave field

6.4.1 Single Groin in Constant Water Depth

To first study how the characteristics of the groin affect the amount of energy that is

reflected and transmitted, and thus the wave field around the structure, a simple case is

presented. The model was run over a domain with dimensions 600x900 meters and a 400

m long groin located at y=450 m. The rest of the necessary input conditions are given as,

NX NY Ax (m) Ay (m) Depth (m) Hi (m) T (sec) 0 (o)

60 60 10.0 15.0 5.0 1.0 7.0 20.0

Several cases were run for this configuration with different values assigned to the

physical properties of the groin. The following table includes the values used as well as the

reflection and transmission coefficients calculated by the model and employed in the

boundary conditions around the structure,

Case 1 Case 2 Case 3 Case 4

Inertial Coefficient, S 1.0 1.0 1.0 1.0

Friction Coefficient,f 0.0 1.0 0.5 0.01

Permeability, e 0.0 0.1 0.5 0.95

Reflection Coefficient, R 1.00 0.66 0.14 0.05

Transmission Coefficient, T 0.00 0.12 0.36 0.96

84

The model results are presented in Figures 6.10 to 6.13, the contour lines represent the

instant water free surface and show the bending of the wave crests behind the groin as well

as the formation of short crested waves perpendicular to the groin in the upwave area. The

relationship between the calculated and the incident wave height (H/Ho) is shown in shades

of gray with a scale on the right hand side of each figure. Figure 6.14 presents the wave

height change along four sections parallel to the groin, the first two at 6.5 and -6.5 meters and

the others at 67.5 and -67.5 meters. It can be seen that H/Ho always increases for increasing

R values in the upwave side of the groin (negative y). However in the lee of the groin H/Ho

does not seem to always decrease with decreasing transmission coefficient. For Case 3,

which has a larger transmission coefficient than Case 1, HIHo values at y=6.5 are greater

initially than the values from Case 1 near the head of the groin, but become smaller towards

the root of the groin. This effect is thought to be due to the superposition of transmitted and

diffracted wave fields on the lee side of the groin. At y=67.5 m, further downwave where

the diffraction effect becomes less important, H/Ho for Case 3 is clearly greater than that for

Case 1. Case 4 shows the results when the groin is close to being almost 'transparent' to

wave propagation.

The special case of groins composed of pile screens can also be simulated

approximately by choosing the right physical characteristics for the groin. Here, an inertial

coefficient, S, of around 2 would seem to be more suitable. The permeability will also be

considerably higher than that of a rubble mound groin, with values usually around 0.8 to 0.9.

The friction factor on the other hand will be reduced. Ideallyf should be a value dependent

upon the wave field also, thus, successive iterations would lead to the rightfvalue based on

the Lorentz's equivalent work concept..' In this model, for the sake of computational time,

such dependency is not included, andfis a given input.

H/Ho
100

-200 -150 -100 -50 0 50 100 150 202.2

1 -20 0 -150 -100 -50 0 50 100 150 200

1.8

1.6

0.8

S 0.4

S0.2

-1 -, ---------... ..... L 0- o
-200 -150 -100 -50 0 50 100 150 200
Distance from center of jetty y (m)

Figure 6.10: Free surface and amplification factor around one groin, Case 1.

H/Ho

2.2

0 1.8

1.6

3001.4

0.8

0.6

0.4

-200 -150 -100 -50 0 50 100 150 200
Distance from center of jetty y (m)

Figure 6.11: Free surface and amplification factor around one groin, Case 2.

H/Ho
In1

2.2

-102

1.8
30
1.6

0.8

-200 -150 -100 -50 0 50 100 150 200
0.6

0.4

0.2

fi...ll ----_- 0

-200 -150 -100 -50 0 50 100 150 200
Distance from center of jetty y (m)

Figure 6.13: Free surface and amplification factor around one groin, Case 4.
H/Ho

2.2

1.8

91.4

0.6

0.4

-200 -150 -100 -50 0 50 100 150 200
Distance from center of jetty y (m)

Figure 6.13: Free surface and amplification factor around one groin, Case 4.

2.5
2.25
2
1.75
1.5
1.25
1
0.75
0.5
0.25

100 150 200 250 300 350 400 450 500 550 6(
x (m)

2.25
2
1.75
S1.5
1.25
1
0.75
0.5
0.25
0
1(

2.5
2.25
2
1.75
o 1.5
1.25
1
0.75
0.5

y=-7.5 m

)0 150 200 250 300 350 400 450 500 550 600
x (m)

100 150 200 250 300 350 400 450 500 550 600
x (m)

1.75

1.25

0.75
0.5 y=67.5 m
0.25

100 150 200 250 300 350 400 450 500 550 600
x (m)

Figure 6.14: Amplification factor along four sections
parallel to groin.

CaseI
S Case 2

- --- Case. 3

y=-67.5 m

- Case 1
- Case 2
- - Case 3
-- - Case 4

y=7.5 m --
, . . .. ." . . .

6.4.2 Single Groin on a Plane Beach

For the case of a varying bottom topography the computation of reflection and

transmission coefficients for a porous groin becomes more complicated, as transmission and

reflection coefficients are functions of wave properties which, in this case, change along the

length of the structure. At first it was decided not to include these changes in the model.

However, after further consideration it was decided to include the complete solution to

compute reflection and transmission coefficients to examine the effect. A case was run with

constant reflection and transmission coefficients based on the groin characteristics and angle

of incidence at the offshore tip of the groin. The groin was located in plane beach with a 1

to 50 slope. Then, the full wave-dependent solutions of transmission and reflection

coefficients were used. The results are compared in Figure 6.15. The table below gives the

computed TIR values at the tip and the base of the groin for the case of = 20.

S f E 0 tip ~ base tip tip Rbase Tbase

2.0 1.0 0.8 20 16.55 2.68 0.24 0.07 0.80 0.00

It can be seen that the reflection coefficient increases from 0.26 at the offshore end of the

groin to 0.80 at the its base; this effect will always occur when groins are placed across a

surfzone, thus, it seems important to compute R and T coefficients based on a porous flow

model rather than just using constant values.

Using the full solution to compute R and T coefficients not only accounts for changes

in their value along the groin, but also simplifies the input given to the model. Figure 6.16

compares the cases when same R and T coefficients are used for two different offshore wave

angles and when the model computes those coefficients depending on the value of 0,. The

89

fixed coefficients correspond to an offshore angle, 06, of 20 degrees, and the variable values

to 0/=5.

1.75 -- Fixed R,T
1.5 Variabe R,T
1.25

0.75
0.5 Front side of groin
0.25

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x/Lo

2
1.75 Fixed R,T
1.5 Variabe R,T
1.25

0.75 ... ... .....................

075
0.5 Lee side of groin
0.25

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x/Lo

Figure 6.15: Comparison of model results for the cases of fixed and
variable reflection/transmission coefficients.

2
1.75 Fixed R,T
1.5 ...... Variabe R,T
1.25 -
1-
0.75
0.5 Front side of groin
0.25 -
0 -1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x/Lo

2
1.75 Fixed R,T
1.5 .... VariabeR,T
1.25
1
0.75
0.5 Lee side of groin
0.25
0-
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x/Lo

Figure 6.16: Offshore incident wave angle influence on reflection and
transmission coefficients.

6.4.3 Waves around a Groin Field

Groin fields are used frequently as coastal defense systems. Their performance will

depend very much on their ability to disrupt both wave energy propagation and longshore

currents. The applicability of the model to this kind of situation will be shown in the

following examples. Here, three groins with arbitrary dimensions are introduced. An

unlimited number of groins can be introduced in the computational domain as long as

physically permissible. The input wave conditions, as well as the computational domain

dimensions are given as follows,

NX NY Ax (m) Ay (m) Slope Hi (m) T (sec) 0 ()

60 70 5.0 15.0 1/60 1.0 7.0 20.0

The three groins are located at y=225, 525 and 825 meters, with lengths 100, 150 and 125

meters respectively. The dimensions are definitely larger than what is usually constructed

on real coastlines, such that the effects of groins are exaggerated to aid in visual appreciation

in the figures. Anyway, structures up to 300 meters long have been used in places like the

Mediterranean Spanish coast. Figure 6.17 shows the results when the groins are completely

impermeable, referred as Field 1. On Figure 6.18 the case when the three groins have

different physical characteristics, Field 2, is presented. The characteristics of the groins for

this case are given as:

S f e

Groin at y=225 m 1.0 1.0 0.4

Groin at y=525 m 1.0 1.0 0.2

Groin at y=825 m 1.0 2.0 0.1

MILLISECOND	CLASS.METHOD	MESSAGE
0	sobekcm_page_globals.constructor
0	sobekcm_page_globals.constructor	Application State validated or built
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.constructor	Navigation Object created from URI query string
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.display_item	Retrieving item or group information
0	sobekcm_page_globals.get_entire_collection_hierarchy	Retrieving hierarchy information
0	sobekcm_assistant.get_entire_collection_hierarchy
0	cached_data_manager.retrieve_item_aggregation
0	cached_data_manager.retrieve_item_aggregation	Found item aggregation on local cache
0	item_aggregation_builder.get_item_aggregation	Found 'all' item aggregation in cache
0	system.web.ui.page.page_load (ufdc.page_load)
0	sobekcm_page_globals.constructor.on_page_load
0	html_echo_mainwriter.add_style_references	Adding style references to HTML
0	html_echo_mainwriter.add_text_to_page	Reading the text from the file and echoing back to the output stream
27	html_echo_mainwriter.add_text_to_page	Finished reading and writing the file