# Numerical Modeling of Irregular Water Wave Transformation

Authors: K. M. Fassieh; O. Fahmy; M. A. Zaki
Journal: International Journal of Oceanography
 (2014)

---

## Abstract

Propagation of irregular water wave from deep water to a shoreline has been numerically modeled. Linear and irregular waves are considered. Model equations govern effects of shoaling, refraction, and diffraction over a varying bathymetry. The model requires the input of the incoming directional random sea at the offshore boundary. Statistical energy dissipation model is incorporated to predict realistically energy losses due to wave breaking in surf zone. Unlike most of the previous models, this model can predict wave transformation in surf zone where energy dissipation and bottom friction must be taken into consideration. The model does not have the limitation of parabolic approximation models (PAM) that are valid only in case of weak refraction. Finite difference approximations have been used to solve the governing equation. The model results are compared with experimental data for directional random wave propagation over a submerged shoal. Good agreements between the model results and experimental data are shown. Applicability of the model to real coastal areas is shown by application to coastal areas along the Nile Delta Coast, Egypt.


---

## Body

## 1. Introduction

The propagation of surface gravity waves over varying topography is an important problem in coastal engineering. Hence there was a definite need for an efficient and effective method of the estimation of wave conditions over a large expanse of shoreline with arbitrary shape and bathymetry. Several engineering wave models for calculating wave characteristics for coastal engineering applications have been developed in recent years. In particular, the introduction of the mild-slope equation (MSE) has enabled the development of a number of models which are usable for practical purposes [1]. In 1990 Liu gave a derivation of mild-slope equation which includes the effect of currents [2]. MSE is an elliptic partial differential equation that can model many processes such as refraction, diffraction, shoaling, and wave breaking. The numerical difficulties associated with solving an elliptic partial differential equation are considerable. To overcome these problems the method of the parabolic equation has been introduced to the study of surface gravity waves propagation and it has been an effective model for dealing rapidly and accurately with propagation problems in coastal areas, for example, [3–5]. This allows efficient numerical solution over large coastal regions. The major drawback with such a model is the restriction on the directional variation of the propagating wave field. This method is based on the assumption that waves propagate within a limited range of angles about an assumed propagation direction. Although many numerical models have been developed to study wave transformation over an irregular bathymetry based on MSE or its parabolic approximation, for example, [6–8], most of these models did not include the energy dissipation and hence it is not applicable in surf zone.To model the offshore irregular sea state many researchers used representative parameters of its directional spectrum, such as significant wave height,H
s, peak energy period, T
p, and dominant wave direction. By using such a description a monochromatic propagation model may be used to transform a single representative wave train inshore. From this the nearshore wave climate is then assessed. Wave tank experiments have demonstrated that monochromatic models can significantly overpredict the maximum amplification of irregular waves due to bathymetric variations [9].In this paper a numerical model is developed that uses the mild-slope equation to study random wave transformation from deep water to shoreline. This model is based on discretization of the offshore spectrum into individual monochromatic directional components. Then a transformation of each component is obtained and the statistical characteristic at any point in the study area is obtained using linear superposition. A statistical energy dissipation model is used to predict energy loss due to wave breaking. The model results are compared with experimental data, Vincent and Briggs 1989, for directional random wave propagation over a submerged shoal. Good agreements with experimental data are shown. Applicability of the model to real coastal areas is shown by application to coastal areas along the Nile Delta Coast, Egypt.


## 2. Wave Model for Individual Component

The mild-slope equation is given by(1)
∇
·
(
c
c
g
∇
ϕ
)
+
ω
2
c
g
c
ϕ
=
0
,
where ϕ
=
ϕ
(
x
,
y
) = complex velocity potential.∇
=
∂
/
∂
x
,
∂
/
∂
y is the horizontal gradient operator.c
(
x
,
y
) is the wave celerity, c
g
(
x
,
y
) is the group velocity, ω is the angular frequency, and
(2)
c
=
g
k
tanh
k
h
,
c
g
=
n
c
,
n
=
1
2
+
k
h
sinh
⁡
2
k
h
in which k
=
k
(
x
,
y
)
=
ω
/
c is the wave number, h
=
h
(
x
,
y
) is the local water depth, and g is the gravity acceleration.ϕ can be written as [10]
(3)
ϕ
=
a
e
-
i
S
,
in which a
=
a
(
x
,
y
) = wave amplitude = H
/
2
;
H
(
x
,
y
) = wave height; and S
(
x
,
y
) = phase function.Substituting the above expression of velocity potential in (3) we obtain (4) and (5), respectively, from the real and imaginary parts. Consider
(4)
1
a
c
c
g
∇
·
(
c
c
g
∇
a
)
+
k
2
-
∇
S
·
∇
S
=
0
,
(5)
∇
·
(
c
c
g
a
2
∇
S
)
=
0
.
The first term in (4) represents the diffraction effects. Neglecting this term gives the refraction equations (eikonal and transport equations)
(6)
∇
S
=
k
_
=
k
cos
⁡
θ
i
+
k
sin
θ
j
,
∇
·
(
ω
a
2
c
g
)
=
0
,
where θ
=
θ
(
x
,
y
) = the wave angle and i
,
j = unit vectors in x and y directions, respectively.The gradient of the phase function takes the form(7)
∇
S
=
|
∇
S
|
cos
⁡
(
θ
)
i
+
|
∇
S
|
sin
(
θ
)
j
.
Substituting (7) in (4) and (5) gives
(8)
|
∇
S
|
2
=
k
2
+
1
H
(
∂
2
H
∂
x
2
+
∂
2
H
∂
y
2
+
1
c
c
g
(
∂
H
∂
x
∂
c
c
g
∂
x
+
∂
H
∂
y
∂
c
c
g
∂
y
)
)
,
(9)
∂
∂
x
(
H
2
c
c
g
|
∇
S
|
cos
⁡
(
θ
)
)
+
∂
∂
y
(
H
2
c
c
g
|
∇
S
|
sin
(
θ
)
)
=
0
.
Irrotationality condition gives ∇
×
∇
S
=
0; that is,
(10)
∂
∂
x
(
|
∇
S
|
sin
(
θ
)
)
-
∂
∂
y
(
|
∇
S
|
cos
⁡
(
θ
)
)
=
0
.
The corresponding dispersion relation is as follows:
(11)
ω
2
=
g
k
tanh
(
k
h
)
.
Equations (8)–(11) constitute a complete set of equations that describe the monochromatic water wave transformation due to change in bathymetry without energy dissipation. These equations will be solved to get the transformed wave characteristics θ, H, and |
∇
S
|. Equation (10) describes the kinematical aspect of wave transformation (wave angle θ), while (9) describes the dynamical aspect (wave height H). Equation (8) describes the change in phase function due to diffraction effects.


## 3. Energy Dissipation

In the previous section the model was developed based on the assumption that no energy dissipation occurs during the wave transformation process. However, in most coastal problems the effects of energy dissipation, such as bottom friction and wave breaking, may become significant, Liu 1990 [2]. As a wave proceeds into shallow water, it continues to shoal and increase in wave height. However, at some depth, a wave will become unstable and break. Seafloor and wave characteristics determine how a wave will break. The mild-slope equation may be modified in a simple manner to accommodate these phenomena by including an energy dissipation function describing the rate of change of wave energy.Following Liu 1990 the modified mild-slope equation becomes(12)
∇
·
(
c
c
g
∇
ϕ
)
+
ω
2
c
g
c
ϕ
=
i
ω
F
,
in which F is the energy dissipation function.The energy dissipation functions are usually defined empirically according to different dissipative processes, for example, Dalrymple et al. 1984 [11]. To consider wave breaking in the parabolic approximation of the mild-slope equation for spectral wave conditions, it is possible to introduce a wave breaking term. Here we use the breaking model by Thornton and Guza 1983 [12] where energy dissipation function is taken to be
(13)
F
=
α
ϕ
,
where α is the dissipation coefficient and is defined as
(14)
α
=
3
π
4
ω
B
3
γ
4
h
5
H
5
,
where h = local water depth; the constants B and γ are chosen to be equal to 1 and 0.6, respectively.The dissipation coefficient is very small when breaking does not occur. However, once breaking starts it begins to take on significant values and energy is dissipated from the wave field. Now (9) becomes
(15)
∂
∂
x
(
H
2
c
c
g
|
∇
S
|
cos
⁡
(
θ
)
)
+
∂
∂
y
(
H
2
c
c
g
|
∇
S
|
sin
(
θ
)
)
=
ω
α
H
2
.
Now (8), (15), (10), and (11) constitute the complete set of equations that describe the one component of the irregular monochromatic water wave transformation due to change in bathymetry with energy dissipation.


## 4. Boundary Conditions

The system of equations to be solved is of elliptic type, needing a boundary condition along the whole contour of the solution domain. Possible boundary conditions are as follows.(a)
At the offshore driving boundary, the incident directional spectrum is completely specified and discretized into components.
(b)
At fixed boundaries, partial or full reflection is considered.
(c)
At open boundaries, radiation condition is assumed.


## 5. Numerical Solution

The incident directional spectrum is discretized into components. Each component has its own frequency and direction. A discrete representation of the directional spectrum, for the direction and frequency domains(
θ
min
⁡
,
θ
max
⁡
), (
f
min
⁡
,
f
max
⁡
), is considered as follows:
(16)
f
m
=
f
min
⁡
+
(
m
-
1
)
Δ
f
for
m
=
1
,
…
,
M
,
θ
n
=
θ
min
⁡
+
(
n
-
1
)
Δ
θ
for
n
=
1
,
…
,
N
,
where Δ
f
=
(
f
max
⁡
-
f
min
⁡
)
/
(
M
-
1
) and Δ
θ
=
(
θ
max
⁡
-
θ
min
⁡
)
/
(
N
-
1
).For each component, (8), (10), (11), and (15) are solved using finite difference technique. The x direction is chosen to be the offshore direction coordinate, while the y-coordinate is chosen to be the along-shore direction.The study area is divided into grid mesh ofM
×
N rectangular grid cell, each having lengths of Δ
x and Δ
y in the x and y directions, respectively. The mesh sizes Δ
x and Δ
y can be different from each other. The Crank-Nicolson method is used to rewrite the governing differential equation in a finite difference form. Both the x and y derivatives are written using central differences.A summary of the computation process is described by the order of execution; deep water parameters are given as wave period (T), incident wave height (H
o), and incident wave direction (θ
o). Using the wave period (T) the dispersion relation equation (11) is solved to compute the wave number (k). Wave celerity (c) and group velocity (c
g) over the computational grid points are given as follows:
(17)
c
=
ω
k
,
c
g
=
1
2
(
1
+
2
k
h
sinh
⁡
2
k
h
)
c
.
To set up the model, the diffraction effects are neglected at the first iteration; that is, |
∇
S
|
=
k. Then (10) and (15) are solved to compute the wave direction θ and wave height H over all grid points. In solving (10) and (15), forward differences are taken in x direction and central difference in y direction. Lax technique is used to achieve unconditional stable scheme [13]. For the second iteration wave phase function |
∇
S
| is computed using (8) and then (10) and (15) are solved again using the obtained value of |
∇
S
| to account for diffractive effects. These steps are repeated until the required accuracy is achieved.The transformed spectrum can be evaluated at any grid point by the superposition of the different wave components and the significant wave height can be computed as(18)
H
s
(
x
,
y
)
=
2
∑
m
=
1
M
∑
n
=
1
N
|
H
m
n
(
x
,
y
)
|
2
.


## 6. Model Verification

To test the model for irregular wave propagation an experiment of Vincent and Briggs 1989 is reproduced. The experimental setup is shown in Figure1.Figure 1
The experimental setup for Vincent and Briggs, [9].The directional spectral wave basin of 35 m wide by 29 m long has a constant water depth of 0.457 m. The elliptic shoal has a major axis of 3.96 m and minor axis of 3.05 m and a maximum height of 0.3048 m. The water depth above the shoal is given as(19)
h
(
X
,
Y
)
=
0.914
-
0.762
1
-
(
X
3.81
)
2
-
(
Y
4.95
)
2
,
where X and Y are local coordinates aligned along the minor and major axes, respectively, and the shoal boundary is defined by (
X
/
3.05
)
2
+
(
Y
/
3.96
)
2
=
1.A frequency spectrum of the incident wave is defined by the TMA model where the energy densityS
(
f
) is given as
(20)
S
(
f
)
=
α
g
2
(
2
π
)
4
f
5
exp
⁡
(
-
1.25
(
f
m
f
)
4
mmmmmmmmmm
+
ln
⁡
(
γ
)
exp
⁡
[
-
(
f
-
f
m
)
2
2
σ
2
f
m
2
]
)
φ
(
f
,
h
)
,
where α is a constant, f
m is peak frequency, γ is peak enhancement factor, σ is shape parameter, and φ
(
f
,
h
) is a nondimensional function that take into account the effect of finite water depths.For the experiments,φ
(
f
,
h
) was allowed to vary as
(21)
φ
=
{
0.5
ω
h
2
ω
h
<
1
1
-
0.5
(
1
-
2
ω
h
)
2
1
≤
ω
h
≤
2
1
ω
h
≥
2
.
The directional spreading function is given by Borgman [14]:
(22)
G
(
θ
)
=
1
2
π
+
1
π
∑
l
=
1
L
exp
⁡
[
-
(
l
σ
m
)
2
2
]
cos
⁡
(
l
(
θ
-
θ
m
)
)
,
where σ
m is the spreading parameter, θ
m is the mean wave direction, and L
=
20 is the total number of wave components.The following four test cases are considered to verify the model results:N3: narrow directional spreading and narrow frequency spreading,
α
=
0.00155; γ
=
2; and σ
m
=
10;
N4: narrow directional spreading and broad frequency spreading,
α
=
0.00047; γ
=
20; and σ
m
=
10;
B3: broad directional spreading and narrow frequency spreading,
α
=
0.00155; γ
=
2; and σ
m
=
30;
B4: broad directional spreading and broad frequency spreading,
α
=
0.00047; γ
=
20; and σ
m
=
30.Figures2, 3, 4, and 5 show a comparison of the model results for narrow and broad directional spectra against the experimental data along the transect number 4 which lies behind the shoal. For the four test cases, wave period is 1.3 sec and representative wave height is 2.54 cm. The comparisons show a good agreement between the model results and the experimental data.Figure 2
Normalized wave heights along transect 4 for case N3.Figure 3
Normalized wave heights along transect 4 for case N4.Figure 4
Normalized wave heights along transect 4 for case B3.Figure 5
Normalized wave heights along transect 4 for case B4.


## 7. Wave Transformation along the Nile Delta Coast

To study the transformation of wave characteristics along the Nile Delta Coast, the area was covered by a mesh of grid 2 km × 2 km. The bathymetry of the study area is shown in Figure6.Figure 6
Bathymetry of the study area.The Mediterranean weather is highly seasonal in nature. Three seasons with respect to atmospheric circulation patterns have been distinguished by Hamed, 1983 [15], which are winter, spring, and summer. Spring season extends from April to May, while the winter season extends from November to March and the summer season covers the period from June to September.Two stations for the measurement of wave energy and direction were installed off the Nile Delta at Abu Quir and Ras el Bar and reported by Naffaa 1995 [16]. Measurements at Ras el Bar are used in the current study. The predominant wave direction is from the NW and WNW whereas waves from the NNE and NE sectors are limited in magnitude and occur primarily during the summer months.The incident directional distributions of the wave heights during the three seasons are extracted from the available wave roses as shown in Figure7 [16]. Wave roses of November, May, and August are considered as representatives for winter, spring, and summer seasons, respectively. Wave rose is a graph that shows the directional distribution of wave heights during a certain period of time. The probability of different wave heights in each direction is measured from the wave roses for winter, spring, and summer seasons.Figure 7
Monthly wave roses along the Egyptian northern coast 1985.Histograms of wave period during the different seasons are also available in [16] and are shown in Figure 8. A wave period of about 7-8 sec characterized the winter and summer seasons reducing to 5-6 sec during the spring season.Figure 8
Histograms of wave period during the different seasons.The incident directional spectrum is discretized into components. Each component has its own frequency (or wave period) and direction. The solution of each component is carried on. The transformed spectrum can be evaluated at any grid point by the superposition of the different wave components and the significant wave height can be computed as(23)
H
s
(
x
,
y
)
=
2
∑
m
=
1
N
f
∑
n
=
1
N
θ
|
H
m
n
(
x
,
y
)
|
2
.
Figures 9, 10, and 11 show the significant wave height distribution during the different seasons. Comparing Figures 9 and 10 it is shown that the wave height during winter season is much more severe than the summer season. The maximum significant wave heights were found to be 3.2 m, 2.9 m, and 2.1 m for winter, spring, and summer, respectively. Figure 10 shows stronger refraction effects in summer. A large spatial variation in the significant wave height is observed around the head of Ras El Bar and this can be explained by the large variation in depths in this area.Figure 9
Wave height distribution during winter along the Nile Delta Coast of Egypt.Figure 10
Wave height distribution during summer along the Nile Delta Coast of Egypt.Figure 11
Wave height distribution during spring along the Nile Delta Coast of Egypt.


## 8. Discussion and Conclusions

A numerical model to study random wave propagation has been presented. The model is based on the superposition of results from a number of spectral components. An iterative procedure to solve the elliptic type MSE is introduced. Unlike most of previous work, the energy dissipation due to wave breaking is taken into consideration and the applicability of the model extended to the shoreline. The presented results show good agreement with experimental data. The model is applied to study the random wave propagation along the Nile Della Coast of Egypt.



---

*Source: 563467-2014-08-12.xml*