UFL/COEL95/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. LiHwa 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 TimeAveraged Wave Energy ........
2.2.2 TimeAveraged Wave Energy Flux ...
2.2.3 TimeAveraged 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 aveCurrent Interaction ................................. 59
5 SURF ZONE MODEL .......................................... 62
5.1 TimeAveraged 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 SubGrid 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 wavecurrent interaction. ..................... 60
4.20 Comparison of shearing wavecurrent interaction,
amplification of wave height. ...................................... 60
4.21 Comparison of shearing wavecurrent 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 crossshore 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 WaveCurrent 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,
wavecurrent 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 ripcurrents and seaward loss of sand.
* Extensive leeside 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 ripcurrents.
* The aforementioned ability to reduce longshore transport in a gradual way, thus
reducing the leeside 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 wavecurrent 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 finitedifference 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 nonlinearity 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 nonpropagating 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 wavecurrent 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 refractiondiffraction
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 wavecurrent
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 depthintegrated wave energy equation is derived and the exact forms of
wave energy and wave energy flux are presented in waveaveraged 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 zaxis, 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
ua, 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 qCh.
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.xot), 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)+AKV 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 depthintegrated 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)]*VhL[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 UVh=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 LonguetHiggins and
Stewart (1961) and Witham (1962),
FE
t +VhF=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 TimeAveraged 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
+ pen (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 wavecurrent
16
interaction. LonguetHiggins 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 TimeAveraged Wave Energy Flux
The energy flux expression given in Eq. (2.28) can be given, using the Bernoulli
equation, as
F = pfUdz (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 TimeAveraged Wave Energy Equation
Now Eqs. (2.33) and (2.37) can be substituted in Eq. (2.29), obtaining,
+Vh (Cg +Uk)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.xwt), 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 refractiondiffraction 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 mildslope 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 slopeequation to cover wavecurrent 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) +(o2k2CCg)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 waveinduced currents which are usually accompanied by the waves. In
Section 3.2, the governing wave equation of the wavecurrent 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 firstorder effect of the slope is
negligibly small. Booij (1983) and Smith and Sprinks (1975) numerically examined the
validity of the mildslope 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.33.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 longterm 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
+ 72k2 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.
CCgQVhVht + 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 wellknown 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 twodimensional 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 twodimensional 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 2iUTV.V + (UV)(UV) + (UV)(iw, + UV)
V(CCgV) + (o k 2CCg) =0
Rearranging and factoring out e wt,
iw[2UV0 + (UV) ] + (UV)(IVq) + (V)(UV) (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 twodimensional (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)Oi 2uix i)2vqpy
The second and third terms will give,
(UV)(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 2k2CCg} =0
To study the nature of the above equation we can rewrite it as,
aq5x +2b +cqy =F( x y, ,x,y)
where,
a = u2CCg
b = uv
c = v2CCg
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 wavecurrent 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 wavecurrent 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 ydirectional 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 firstorder 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( ) +( v2CCg)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(x1, i1)
Yj, = yi +2hf(xi,+jh,yj) j= 1,2,..,n1
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 ubGrid
Figure 4.1: SubGrid 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 Tgroins
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=tanl(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(p1)Ae in(q1Ay = e in(ql)Ay
Expressing the Helmholtz equation in finite difference form and using 0(p,q),
P20 + p 2cosnAy2 +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)22, 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
CushmanRoisin (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 Lo, 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 q5ik'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
Siky2Ayq 1 +ikyAy/2
S+ t (4.8)
1ikAy/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 ik 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)Yikcosa'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)Yik' B'tY'
ik' + '= [i(k'(1+r)+iky(lfr)]Y
(5.14)
(5.15)
__(1+/
iky(1pr) lp ,) +
ik'Y(l+fr) +iky(l1r) 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(1fr) [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):
(1Prl)
) 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(nm), (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 (1fP) Ay
Pr t ml ^m
6t1 y 0 y (5.21)
(16r) 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=NY1 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, wavecurrentinteraction, 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 semiinfinite
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 sina 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 012345678
y/Lo
2.5
2
0.5
n
2.5
2
1.5
0.5
8 7 6 5 4 3 2 1 012345678
y/Lo
8 7 6 5 4 3 2 1 012345678
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 012345678
y/Lo
2.5
2 Sec. 4
1.5
0.
8 7 6 5 4 3 2 1 012345678
y/Lo
2.5 ... . .
2 0 Sec. 5
0.5
0.5 0 Aoaoo 5^
8 7 6 5 4 3 2 1 012345678
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 012345678
y/Lo
8 7 6 5 4 3 2 1 012345678
y/Lo
2.5
2
1.5
0.5
8 7 6 5 4 3 2 1 012345678
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 WaveCurrent interaction
When surface waves of any kind propagate over the surface of a medium in steady
but nonuniform 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 nonuniform 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 wavecurrent 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 LonguetHiggins and
Stewart (1961). The comparisons with their solutions are shown on Figures 4.1921. 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 wavecurrent 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 wavecurrent
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 wavecurrent
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 wavecurrent 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 wavecurrent 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 wavelike periodic motion that is
quasistationary when timeaveraged 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 nondissipative 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 timeaveraged 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
wavecurrent 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 quasisteady 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 twodimensional 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 CgCgb 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 surfzone 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 crossshore 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 crossshore 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
50 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 crossshore 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 wavedependent 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
