Implementation of Wave Effects in The Unstructured Delft3D Suite Final
Implementation of Wave Effects in The Unstructured Delft3D Suite Final
Implementation of Wave Effects in The Unstructured Delft3D Suite Final
Author: B. Stengs
Email: basstengs@gmail.com
Committee:
Prof.dr.ir. M.J.F Stive Delft University of Technology, HE-CE
ir. A.P. Luijendijk Delft University of Technology /Deltares, HCO
dr.ir. M. Zijlema Delft University of Technology, HE-EFM
Prof.dr.ir. J.A. Roelvink Deltares/UNESCO-IHE
Prof.dr.ir. G.S. Stelling Delft University of Technology, HE-EFM
MSc Thesis: May 28, 2013
ii
MSc Thesis: May 28, 2013
Abstract
Recently an unstructured version of Delft3D, named D-Flow FM, has been developed at Deltares. This
tool shows high potentials for coastal modelling due to the flexibility in domain resolution and
unstructured grids. This MSc thesis focussed on implementing wave effects in D-Flow FM with
satisfying results. In addition, a coupling system is built, enabling a flexible coupling between D-Flow
FM and the wave module SWAN, for modelling the wave-current interaction. An innovative program-
structure is presented which provides a dynamic coupling tool enabling interactive modelling.
This coupled tool shows large opportunities for coastal models where flexible resolution and
geometry fitting are desired. Applications like the Sand Engine or the Wadden Sea modelling tidal
velocity patterns at offshore regions and meandering gullies and rip channels at nearshore regions,
all in one domain. These kinds of models are usually very computationally intensive, especially when
morphodynamics are also included. With the introduction of unstructured modelling, coastal
applications like the Sand Engine or the Wadden Sea, are modelled more efficiently.
The wave formulations implemented in D-Flow FM are state-of-the-art and correspond to the
implementation as in Delft3D, which have been tested and validated with success. In the validation
of wave modelling by D-Flow FM, four different test cases are considered. First, a simplified open
channel flow is analysed, where the general performance of D-Flow FM is verified to the analytical
solution. Second, an M2-tide for an alongshore uniform beach is modelled. Next, oblique incident
waves are modelled on a planar beach. And finally, a non-uniform barred beach is modelled (Egmond
case) with rip channels induced by waves. By means of these test cases the wave formulations in D-
Flow FM are validated.
From the test cases, it is concluded that modelling waves in D-Flow FM generally gives good
agreement with validated Delft3D models. The velocity patterns, both in magnitude and direction,
show a good match, as confirmed by the Egmond case. The water levels modelled by the planar
beach case match the Delft3D simulation as well. However, despite a good performance inside the
surf zone, small model deviations are found outside the surf zone, at intermediate water depths.
These model deviations should be reduced by correcting the bed shear stress for Stokes drift.
Past studies have proven the benefits of unstructured modelling over structured techniques, which
shows high potential for D-Flow FM. With the inclusion of wave effects in D-Flow FM, improvement
to nearshore modelling is made, but an extensive validation needs to be done. The dynamic software
coupler provides an interactive modelling experience, and serves as a comprehensive tool for general
model understanding. These developments have contributed significantly to improve nearshore
modelling by D-Flow FM.
iii
MSc Thesis: May 28, 2013
iv
MSc Thesis: May 28, 2013
Acknowledgements
This MSc Thesis is made possible by the support and consult of many experts from Deltares. Thanks
to the internship of eight months at Deltares prior to the graduation work I was able to push the
limits of my MSc work. I would like to give special acknowledge to those who made this possible for
me, especially Arjen Luijendijk.
The Delft3D software has been an inspiration for the development of the wave formulation in D-Flow
FM. I would like to give special thanks to Jan van Kester for his assist and knowledge contribution of
the Delft3D applications. I want to thank Adri Mourtis, Jan Mooiman, Sander van der Pijl and Dano
Roelvink for their support and consult during glitches and fatal errors.
Most of all I want to thank Fedor Baart for his assistance in programming and his large contribution
to coupling D-Flow FM to SWAN, with innovative technics. Despite a doctoral thesis and becoming
dad, he always managed to make time for (Fortran-, Python- or MATLAB-) programming both ends
together. I have always enjoyed the endless sessions with Fedor Baart which were very pragmatic
and enlightening. Many thanks.
Furthermore, I would like to thank my committee for their support and guidance throughout my MSc
thesis. Arjen Luijendijk, my daily supervisor, who established the assignment and always managed to
find time for bridging obstacles. Prof. Dano Roelvink, my expert at Deltares, for showing interest and
giving consult in numerical modelling. Prof. Marcel Stive and Prof. Guus Stelling for supporting me on
behalf of the TU Delft. Marcel Zijlema, expert at the TU Delft, for his review and assistance on behalf
of the TU Delft.
Finally I thank my dear fellow student, Pieter Visser, for his review and his mental support
throughout the curriculum. Last but not least, I thank my fellow students of the ‘Feestruimte’ for
enjoying me with fascinating topics over a cup of coffee.
v
MSc Thesis: May 28, 2013
vi
MSc Thesis: May 28, 2013
Contents
1. Introduction.......................................................................................................................... 1
1.1 Background.............................................................................................................................. 1
1.2 Research question and objectives ........................................................................................... 1
1.3 Approach ................................................................................................................................. 2
1.4 Outline of thesis ...................................................................................................................... 2
vii
MSc Thesis: May 28, 2013
viii
MSc Thesis: May 28, 2013
References ................................................................................................................................. 59
Appendix A. Applied numerics of nearshore modelling with Delft3D............................................ 65
A.1 Introduction ........................................................................................................................... 65
A.2 Integration of FLOW – WAVE module ................................................................................... 65
A.3 Equations FLOW .................................................................................................................... 65
A.4 Time integration FLOW ......................................................................................................... 67
A.5 Spatial integration FLOW....................................................................................................... 67
A.6 Mesh (structured).................................................................................................................. 69
A.7 Delft3D WAVE........................................................................................................................ 70
A.8 Interpolations Delft3D + SWAN ............................................................................................. 73
Appendix B. Conceptual wave coupling ....................................................................................... 75
B.1 Introduction ........................................................................................................................... 75
B.2 Surfbeat model ...................................................................................................................... 75
B.3 Boussinesq and non-hydrostatic models .............................................................................. 76
B.4 Integrated wave modelling ................................................................................................... 76
B.5 Discussion .............................................................................................................................. 77
Appendix C. Sensitivity analysis planar beach case ...................................................................... 79
C.1 Modelling approach .............................................................................................................. 79
C.2 Model results ......................................................................................................................... 83
ix
MSc Thesis: May 28, 2013
x
MSc Thesis: May 28, 2013
Chapter 1
1. Introduction
1.1 Background
In the study of hydrodynamic systems, the analytical solution of fluid motion is often unknown. This
is where numerical models come into play by means of CFD modelling (Computational Fluid
Dynamic). CFD models became invasive with the rise of computational power (e.g. Evans and Harlow,
1957). With modern technology the demand for high resolution model keeps rising. Structured grids
are often used, based on quadrangles only. High resolution modelling is possible, but is less relevant
at offshore regions. A rather new development is CFD modelling by unstructured grids, made from all
kinds of polygons. This allows high flexibility for domain resolution which is much more efficient for
CFD modelling (MacWilliams et al., 2007).
Deltares has recently developed the unstructured version of Delft3D, named D-Flow FM (D-Flow
Flexible Mesh), for integrated modelling of rivers, coastal areas and estuaries. This program covers
the modelling of flow patterns, morphological developments and ecology in an (unstructured)
discretized way. Complex bottom patterns, such as shallow sandbanks, gullies and rip channels, can
form difficulties for modelling with curvilinear grids due to the orthogonality requirements of the
model. Whereas unstructured grids can be applied with much more flexibility and overcome these
adverse effects relative easily (Kernkamp et al., 2011).
There is a strong demand in evaluating this expertise to a higher level of hydraulic engineering, as
stated by Casulli and Walters (2000), Kernkamp et al. (2011) and Verwey et al. (2011). Coastal
modelling belongs to one of these categories. For improving coastal models, an integration of wave
modelling is required. Realizing this piece of engineering requires a validated coupling of D-Flow FM
with SWAN (Simulating WAves Nearshore).
“What developments could improve modelling nearshore hydrodynamics of complex coastal systems
by D-Flow Flexible Mesh?”
1
MSc Thesis: May 28, 2013
1.3 Approach
Essential elements for comprehensive numerical modelling are the physical interpretations of
systems and the model response. A systematic approach is followed in this research which
distinguishes four phases. At first, relevant nearshore processes are investigated. Next the general
operation of Delft3D is analyzed and the implemented physical processes are studied. The same
analysis is carried out for D-Flow FM and extended for numerical shortcomings. Finally, the extended
formulations are tested relative to Delft3D and the analytical solution when possible.
The coherence of numerical modelling and physical modelling are a key element throughout the MSc
work. These two lines of approach connect and collaborate via the model development and model
validation. A fine collaboration of these two approaches is fundamental throughout the process.
After this research, D-Flow FM is extended with nearshore processes that were not included in the
program. The robustness and performance of this piece of coding is tested by means of different test
cases. Finally, conclusions and recommendations are drawn based on the performance of the test
cases.
In the appendix of this report an introduction is given to the operation of Delft3D for nearshore
modelling. Next a conceptual wave-coupling technique is discussed for full spectral modelling. Finally,
a sensitivity analysis is carried out for modelling wave propagation in the nearshore.
2
MSc Thesis: May 28, 2013
Chapter 2
Here for, comprehensive knowledge is required of physical processes in the system before a model is
set up. The relevant contribution of certain system processes determines whether or not such
processes should be accounted for in the model. This schematization is an important procedure to
achieve a sufficient model.
The response of physical systems can properly be described by mathematical expressions, especially
in system dynamics. These systems are described by a rate of change through space or time. In
mathematics this is described as a small unit change per small step through the considered
dimension. By taking the limit to zero for this small step, gives a mathematical expression which is
called a differential equation.
A relative simple physical system is described by a relative simple mathematical expression. For such
systems it is sometimes possible to solve the expression in an analytical way. These analytical models
are divided into elemental parts or basic principles. This allows a direct understanding of the systems
behaviour in terms of dependency of parameters and processes. If a system is influenced by multiple
(complex) relevant processes, it gets rather hard to find an analytical expression for the system. That
is where the numerical models come into play.
A numerical model is a very handy tool for solving complex mathematical processes, but it does not
give a comprehensive knowledge of the system at first sight. There is no one-to-one mapping
possible of the system dependency. Another drawback is the numerical inaccuracies that arise due to
truncation errors and convergence errors to the analytical solution (Zijlema, 2011). These
inaccuracies require a certain tolerance of acceptance before making the numerical model
operational. These drawbacks show the weakness and reliability of numerical models.
3
MSc Thesis: May 28, 2013
The relation to the frequency of the moon-earth system, the frequency of the earth-sun system and
the frequency of the earth’s rotation, determines the astronomical tide. When these frequencies
align spring-tide will occur, and opposite when these frequencies are out of phase neap-tide occurs.
There is a short phase shift of these spring-neap frequencies in relation to the astronomical
frequencies due to the traveling time of the tidal wave from the free ‘forcing tidal wave’ near the
Antarctic.
The rise and fall of the astronomical tide is called the vertical tide. The tide-induced hydrodynamics
are called the horizontal tide. Considering a small tidal amplitude compared to the water depth and a
more or less horizontal bottom slope leads to the following reduced shallow water equations:
Eq. 2-1
Eq. 2-2
This expression describes a balance between inertia, Coriolis, the pressure gradient and bottom
friction. The tidal wave deflects, in relative shallow basins, due to the Coriolis effect. In the Northern
hemisphere this deflection is to the right due to the difference in the earth’s orbital velocity as a
function of latitude. Considering the North Sea basin, this results in a tidal wave originating from the
Atlantic and traveling along the coastlines of the basin. This kind of wave is called a Kelvin wave.
Along the Dutch coast the water depths are considerably shallow, i.e. small depths compared to the
length of the tidal wave. The phase velocity of the tidal wave becomes a direct function of the water
4
MSc Thesis: May 28, 2013
depth and bottom friction, according Eq. 2-1 and Eq. 2-2. Neglecting bottom friction for simplicity,
the expression for wave propagation speed is as follows:
√ Eq. 2-3
The wind exerts a shear stress on the water and induces a setup of the water level. This is expressed
in the equilibrium situation as:
̅
Eq. 2-4
with:
Eq. 2-5
In this equation represents the wind velocity and the wind-stress coefficient which is
determined by Wu (1982). Eq. 2-4 shows that the water level setup due to wind shear stress is
inversely proportional to the water depth, which gives larger water level setups in shallower regions
(the nearshore). This water level difference is compensated by a return flow in the deeper parts of
the water column.
Ocean waves formed by offshore storms and approaching the coast are influenced by many different
processes. In this thesis the process of shoaling, refraction and wave breaking is highlighted in short
for an understanding of the mean mechanism behind wave-induced hydrodynamics in the nearshore.
The nearshore is considered to be the region from the beach to the offshore depth where the orbital
wave motion of the water is not influenced by the bottom anymore. After a theoretical introduction
of ocean waves, the nearshore hydrodynamics induced by wind-waves are divided into cross-shore
and alongshore components.
5
MSc Thesis: May 28, 2013
Physics of waves
Ocean waves generated by a storm can be described by a statistical representation of the sea state
by means of spectral analysis. Assuming an offshore stationary wave record (typically of 15-30min) of
waves which are not too steep, allows a stochastic description of the surface elevation ( ). This is
done by means of Fourier series which takes the sum of a large number of harmonic wave
components each propagating with different amplitudes, frequencies and phases. From this random-
phase/amplitude model the variance density spectrum for frequency and direction can be found,
representing a statistical analysis of the sea state:
⁄ ⁄ Eq. 2-6
{ } [ ]
Integrating the two-dimensional spectrum over all frequencies and over the full
directional domain, gives a the time-averaged surface elevation squared which is also known as the
total variance:
Eq. 2-7
̅̅̅ ∫ ∫
Figure 2-1: Measured wave spectrum near IJmuiden expressed in and time. (Source:
http://www.meetadviesdienst.nl)
The orbital motion of the water particles under a propagating deep water (linear) wave can be
considered as closed circles from a Eulerian perspective. This circular motion decays exponentially
downward over the water column. When waves approach the shoreface, the orbital motion of water
particles under waves start to touch the seafloor and become closed ellipses instead of closed circles.
Time averaging these orbital motions gives a net current velocity of zero.
6
MSc Thesis: May 28, 2013
The propagation speed of a surface wave slows down as it approaches the nearshore. This
propagation speed is called the phase speed of a wave and is best described by the dispersion
relationship for arbitrary depths:
Eq. 2-8
with,
Eq. 2-9
where represents the radian frequency and k the wave number. From Eq. 2-8 one can concluded
that the phase speed of a surface wave depends on the wave length and the frequency. For deep
water Eq. 2-8 becomes:
Eq. 2-10
√ Eq. 2-11
In according to Eq. 2-10 the phase speed of a deep water wave is a function of the wave period only.
By means of spectral analysis, a sea state can be considered as a large sum of harmonic waves, each
with a different wave amplitude, wave period and phase. By Eq. 2-8, waves with a longer wave
period travel faster than waves with shorter wave periods. In space, the longer waves (swell waves)
filter from the short waves and are much more regular and long crested.
Thus, in according to the dispersion relationship, the wave frequencies sort and group together over
space. Summating two of these waves with slightly different frequencies will reinforce at a certain
moment and somewhat later dampen out (e.g. two people clapping non-synchronic). This forms a
longer wave with its crest at the reinforced moments, and a trough when they cancel out. These long
waves are called infra-gravity waves and are bound by the wave groups (phase locked), see Figure
2-2. Bound long waves travel with the speed of the group velocity, which is defined as:
Eq. 2-12
where is the phase speed of the individual waves defined by Eq. 2-8 and,
( ) Eq. 2-13
which is considered to be a relative deepness parameter of the incoming waves with wave number .
7
MSc Thesis: May 28, 2013
Inside the surf zone, where the short waves break, the bound long wave gets released from the
group and is no longer phase locked. The bound long wave becomes a free long wave (Van Dongeren
et al., 2003).
The wave period of infra-gravity waves vary from 30sec – 5min. Under its crest the wave energy at
the coast is much larger than under the trough. This varying wave forcing over time in the nearshore
induces a difference in water level setup, called surfbeat. The wave-induced forces are by definition
related to the radiation stresses under the waves by:
Eq. 2-14a
Eq. 2-14b
Radiation stress is the wave-induced momentum flux, by the particle bodily motions and wave-
induced pressures. Wave forces arise thus from a change in the transfer of momentum due to waves
from one location to another (Longuet-Higgins and Stewart, 1964). The radiation stresses in Eq. 2-14)
are defined as:
( )
( ) Eq. 2-15
Eq. 2-16
8
MSc Thesis: May 28, 2013
Cross-shore hydrodynamics
When waves travel from the deep waters to the nearshore they start to deform due to the influence
of the bottom topography. The waves start to shoal in decreasing water depths, and non-linear
effects force a skewness of the wave. Approaching even shallower waters, the crest of the wave
starts to take over the trough of the wave, in according to the dispersion relationship of Eq. 2-8. This
results in a pitched forward shape of the wave until the wave breaks. The kinematic and potential
energy of the wave is dissipated by the roller energy and turbulence.
The wave-induced energy is directly related to the wave forces. One can expect that spatial varying
wave energy in the cross-shore induces a difference in wave forcing. In the cross-shore, this spatial
varying wave force finds equilibrium with the pressure gradient. This results in a water level set down
just outside the surf zone, and a water level setup inside the surf zone.
To comply with mass conservation during the different stages of wave deformation, hydrodynamics
in the nearshore are triggered in the water column (2DV effects). Outside the breaker zone wave
skewness drives an onshore current near the bed. After wave breaking the mass of the rolling face is
compensated by the undertow, which is an offshore directed current. Inside the breaker zone the
infra-gravity waves drive an onshore current near the bed. Outside the breaker zone this current is
offshore directed (Roelvink and Stive, 1989).
Alongshore hydrodynamics
Ocean waves that approach the coast under an angle (oblique incident waves), are an important
driving force of alongshore hydrodynamics. Tides and wind-induced hydrodynamics also have their
contribution, especially in the outlet of lagoons and shallow waters respectively. For the sake of
simplicity an alongshore uniform coast is considered with oblique incident waves – the terms in
Eq. 2-14) become zero.
Oblique incident waves refract when entering shallow waters, in according to Snell’s law. The wave
rays bend towards shallow bottom contours. The angle of incidence thus converges to zero when
waves enter shallow waters.
From Eq. 2-15 is found that the radiation stresses due to waves are a function of the angle of
incidence, the relative depth and the wave height squared. These parameters become an important
measure for alongshore wave forcing. This alongshore force finds equilibrium with the bed shear
stress outside the surf zone for an infinitely long coastline. Inside the surf zone can be found, by
means of Snel’s Law and conservation of energy principle, that the alongshore driving force becomes
a function of the wave energy dissipation, or wave breaking. By i.a. (Reniers, 1999) was found that
the maxima alongshore current velocities for oblique incident waves are found at the breaker line(s),
i.e. near the bars crest and another maxima close to the shore. These results are shown in Figure 2-3.
9
MSc Thesis: May 28, 2013
Figure 2-3: Alongshore forcing over a uniform barred beach, with local pressure gradient (solid line) and without (dashed
line). Measured data is denoted by ‘x’ (Reniers, 1999)
Modelling wave-driven alongshore currents on a uniform sloping beach, generally show good
agreements with laboratory data and field data. For a barred beach this is not generally the case, i.e.
shear instabilities of the alongshore currents show difficulties in modelling (Reniers, 1999). This is
especially the case with alongshore inhomogeneities like coastal structures, bars, rip channels, and
gullies. Due to these 3D features, complex current patterns may arise which are challenging to
translate into the laws of physics.
Eq. 2-17
This expression is proportional to the velocity squared, which makes the differential equations Eq.
2-1 and Eq. 2-2 non-linear. This makes solving Eq. 2-1 and Eq. 2-2 analytically nearly impossible
without simplifications. Solving non-linear equations like Eq. 2-1 and Eq. 2-2 efficiently requires
numerical analysis.
10
MSc Thesis: May 28, 2013
Figure 2-4: Example of a flexible mesh from the Continental Shelf Model
A rather new development is CFD modelling of unstructured grid, also called a flexible mesh (FM),
made from polygons instead of only quadrangles. This enables much more flexibility for refining the
resolution within one domain (Casulli and Walters, 2000). If desired, locally high resolution can be
applied in modelling vast areas, as shown in Figure 2-4.
These advantages of unstructured models where recognized by many researchers whom contributed
to the development of unstructured CFD modelling. Examples of unstructured models are ADCIRC
(Dietrich et al., 2011), TELEMAC (Jones and Davies, 2006), UnTRIM (Casulli and Zanolli, 2002;
MacWilliams et al., 2007) and MIKE Flexible Mesh (Sørensen, 2004). Deltares has recently been
working on an unstructured model of Delft3D, named D-Flow FM. Primary test cases found positive
results with a need to further development (Kernkamp et al., 2011).
The objective of this thesis is to extend flexible modelling, by means of D-Flow FM, to coastal regions
with wave-dominated systems.
11
MSc Thesis: May 28, 2013
12
MSc Thesis: May 28, 2013
Chapter 3
3.1 Introduction
Since 2010, Deltares has been working on a new computational solver for unstructured grids. This
numerical engine is called D-Flow FM and is mainly based on the original Delft3D module. Primary
tests show promising results for the development of computational modelling. In according to Casulli
and Walters (2000), Kernkamp et al. (2011) and Verwey et al. (2011), there is a strong demand in
evaluating the expertise of unstructured modelling to higher levels of hydraulic engineering.
⃗ Eq. 3-1
⃗
⃗ 〈⃗⃗⃗ 〉 ⃗ Eq. 3-2
In which [ ] is the horizontal gradient operator, q and S representing the source term and
external forcing respectively, ⃗ is the advection term, Ω is the angular velocity of rotation of
the earth and H is defined as,
Eq. 3-3
with ζ the water level relative to a reference plane and d the bottom level positively upward, in
contrast to Delft3D (see Eq. A-14). The SWEs of Eq. 3-1 and Eq. 3-2 express conservation of mass and
momentum (Kernkamp et al., 2011).
13
MSc Thesis: May 28, 2013
The same staggered grid principle as Delft3D is used, but now in terms of flow nodes and net nodes
(Casulli and Zanolli, 2002). The nodes are connected by links. The net nodes are defined on the
corners of the grid cell. The flow nodes are defined at the cell circumcenter; the intersection of all
perpendicular bisectors of the grid cell. A perpendicular bisector is found by drawing two circles with
the net nodes as center and connecting the two intersections of the circles (see Figure 3-1).
{( )( ) ( )( ) ( )( )}
Eq. 3-4a
{( ) ( ) ( ) }
Eq. 3-4b
with;
{ ( ) ( ) ( )}
Eq. 3-4c
At the flow nodes the water levels are defined and at the net nodes the depth points are imposed.
The velocity points are defined at the cell faces, along the flow links. This gives a similar staggered
grid principle as Delft3D but now the spatial integration is in one-dimension, namely along the flow
links.
In Delft3D, the water level in a single cell and depth values at the cell’s interface are determined by
the DPSOPT and DPUOPT respectively. In D-Flow FM this option is determined by choosing a
⃗
BotLevType. The bottom level at the velocity points ( ) are found by the minimum-, mean- or
maximum depth value of the net nodes ( ). And then (by default) the bottom level at the water
⃗
level points ( ) are determined by the lowest connected link of .
14
MSc Thesis: May 28, 2013
Figure 3-2: Example of 2DH cells defined on an unstructured grid with a sloping bathymetry
To assure numerical stability, the Courant number (CFL condition) may not be chosen too large due
to the explicit part of the numerical solver. Defining the numerical time step by means of the CFL-
criteria, shows to be efficient for unstructured grids due to the varying flow links (Kernkamp et al.,
2011). With λ as a dimension of space, the Courant number is defined as:
√ Eq. 3-5
After discretizing Eq. 3-2, can be expressed in terms of the water level and substituted back
into Eq. 3-1, leaving a system of equations with only the water levels as unknown. Eq. 3-1 transforms
into a matrix vector system of the type , with representing the unknown water level vector
at the new time level and matrix a sparse symmetric diagonally dominant matrix (Casulli and
Walters, 2000). The full discretization of Eq. 3-1 and Eq. 3-2 to a system of equations in matrix form,
can be found in Kernkamp et al. (2011).
Per time step the computation is divided into sections and solved parallel, through the use of
OpenMP. In this sense the computational workload is divided over multiple cores of the processor,
providing high efficiency to the available computational capacity. Conventional models like Delft3D
are not supported with OpenMP and are thus based on single core computations.
15
MSc Thesis: May 28, 2013
A flexible mesh is made of polygons like triangles, pentagons, etc. Rather than conventional models
whom are limited by rectangles. Due to the flexibility in resolution, nested models are no longer
required. Here for, points of interest with different length scales are simulated with one single mesh.
This is economically interesting due to a reduction in computation time. Another benefit of FM
modelling is the ability to merge several grids (3D, 2D and 1D) to one mesh. For example a delta
network with many branches and junctions as shown in Figure 3-3.
Figure 3-3: Example of an unstructured grid for the Mahakam Delta (Borneo)
There are two design requirements for creating an unstructured grid, similar as Delft3D. It needs to
be smooth and orthogonal. The smoothness of a cell is defined by the relative cell surface area to
adjacent cells. This ratio ideally equals unity from which an example is shown in Figure 3-1. The
orthogonality is the degree of perpendicularity the flow link intersects the net link (cell face). This
measure is defined by the cosine of the angle ϕ between both links, ideally being perpendicular
(Casulli and Zanolli, 2002). To comply with these design criteria, D-Flow FM has an
‘Orthogonalise/Smooth’ feature which replaces the net nodes and flow nodes for finding a better
network criteria.
16
MSc Thesis: May 28, 2013
To guaranty numerical stability the length of the flow links needs to be sufficient to comply with the
CFL criteria of Eq. 3-5 and to avoid small denominators. This means that the cell circumcenter needs
to lie within the grid cell. Here for grid cells need to have acute angles. Figure 3-4 gives an example of
a wrong computational cell in red, but with still an acceptable flow link due to the geometry of the
adjacent cell. It would be better to correct this cell by avoiding the obtuse angle and matching the
surface areas to the surrounding cells (see green cell).
Figure 3-4: Example of a bad computational cell in red due to the obtuse angle. The green cell corrected the obtuse angle
which positions the cell circumcenter inside the cell.
Although a flexible mesh allows all kinds of polygons in the domain, mostly curvilinear grid cells are
used. Unstructured cells are mainly used in transition areas with different grid resolution, grid
connections and meandering channels. As demonstrated by Verwey et al. (2011), a curvilinear grid
more or less in line with the main flow pattern is preferred over a fully unstructured grid, due to the
reduced velocity fluxes in space. Relative large velocity gradients make the non-linear advection
terms of Eq. A-6 until Eq. A-8 relative dominant, introducing relative large numerical errors. This
theory suggests an application of grid cells that are elongated and aligned in the flow direction, as
illustrated in Figure 3-3. Models like ADCIRC (Dietrich et al., 2011), TELEMAC (Jones and Davies, 2006)
and MIKE Flexible Mesh (Sørensen, 2004), are less flexible in that perspective due to their fully
unstructured grid restrictions.
Water level:
Discharge:
17
MSc Thesis: May 28, 2013
Neumann:
Riemann invariant: √
The implementation of the harmonic water level boundaries in D-Flow FM are verified to the
WAQUA model by Kernkamp et al. (2011). An application of the Continental Shelf Model (CSM) is
used for a tidal simulation with astronomical boundary components. A uniform grid with spherical
coordinates is applied in the WAQUA model. D-Flow FM used the unstructured grid illustrated in
Figure 2-4. In according to Kernkamp et al. (2011), the computed water levels of both models match
well, which verifies the implementation of the harmonic boundary condition in D-Flow FM.
∑̂ ( ) ∑̂ ( ) Eq. 3-6
where ̂ represents the amplitude of the harmonic signal in [m], and the radian frequency and
wave number respectively determined by Eq. 2-9 and the phase difference over a distance of
is determined by:
Eq. 3-7
The numerical stability is enhanced by adding some numerical damping. Stelling (1983) added the
time-derivative of the Riemann invariant multiplied with a reflection coefficient , to Eq. 3-6. This
gives the following water level boundary:
( √ ) Eq. 3-8
The implementation of this formulation in Delft3D is validated by Gerritsen et al. (2007). In D-Flow
FM the formulation of Eq. 3-6 has been verified by Kernkamp et al. (2011) in terms of astronomical
components. However, the formulation of Eq. 3-8 has yet (2013) to be implemented in D-Flow FM.
18
MSc Thesis: May 28, 2013
Neumann boundary:
A Neumann boundary type is often used for relative small nearshore models but should always be
applied with (at least) one water level boundary to make the system numerically well-posed (Zijlema,
2011). For example, a tidal forcing is simulated by one harmonic offshore water level boundary and
two harmonic lateral Neumann boundaries. In this sense, the Neumann boundaries represent an
alongshore water level gradient, enabling a discharge and a change in water level at the boundary.
Roelvink and Walstra (2004) showed the relevance of implementing the Neumann formulation in
Delft3D.
The Neumann boundary of the type is found by taking the spatial derivative of Eq. 3-6,
which gives the following expression:
∑ ̂ ( ) ∑ ̂ ( ( )) Eq. 3-9
The amplitude of the Neumann boundary is thus a function of the wave number times the wave
amplitude. Due to the cosine formulation an extra phase shift of ⁄ is required. The water level
gradient is implicitly implemented in de numerical solver. Expression Eq. 3-9 is first multiplied with
and then placed on the right hand side of the matrix solver.
For both Delft3D as D-Flow FM the post-processing of the water level gradient have great similarities.
Only D-Flow FM requires a positive slope for model inflow. An application of the Neumann boundary
is given in the test case: Application to a tidal model (Duck schematisation).
19
MSc Thesis: May 28, 2013
3.6.1 Introduction
Relative simple but well understood systems help us to understand the behaviour of numerical
models. In this case a river flow debouching in a large open lake is considered with a constant river
discharge, based on the validation case of Gerritsen et al. (2007). The backwater curve finds
equilibrium between the hydrostatic pressure and the bed shear stress. This river case is well
understood in terms of analytical modelling. This allows a direct understanding of system behaviour
in terms of dependency of parameters and processes. This kind of comprehensive modelling is
preferred over many other kinds of modelling techniques.
( ) Eq. 3-10
√
Solving Eq. 3-10 for the parameters of Table 3-1 gives an equilibrium depth of:
20
MSc Thesis: May 28, 2013
Eq. 3-11
Figure 3-5: Open channel sketch including the water level definitions
Consider a domain of 20 rectangular grid cells. This returns a staggered grid of mxn equals 21x2 grid
points. The grid enclosure, in terms of the grid point (array) number, equals 22x3 as clarified in Figure
3-6. Applying no boundary condition coincides with a ‘closed boundary’ where the normal velocities
are zero. For the river case this coincides with the river banks. Both upstream and downstream open
boundaries are applied with respectively a constant discharge and a fixed water level. The discharge
boundary is positioned on the cell face. The water level boundary is located outside the
computational control volume; in Figure 3-6 this is positioned on array m=22. The observation point
10.000m upstream is at the array m=2.
Figure 3-6: Open channel definition sketch of staggered grid with grid enclosure in grey
With a subtraction of and again the 0.025m due to the staggered grid solver,
one finds:
21
MSc Thesis: May 28, 2013
Eq. 3-12
Eq. 3-13
Figure 3-7: Open channel definition with flow nodes in white and bottom level dotted on network nodes
The condition of the upstream discharge boundary at the cells face, determines the upstream
information of the ghost cell. This is where numerical errors arise. For calibration matters, the D-Flow
FM model is provided with a new feature which does one forward step in space, only at the
discharge boundary. The performance of D-Flow FM with this new feature seems to be stable for this
case study, although the robustness of this solution remains questionable.
Another D-Flow FM model is built, by taking this matter into account. All water levels are now
directly determined on the water level points (the cell circumcenter) instead of the cell face, so no
corrections need to be applied for the so-called ‘staggered grid solver’. The is determined
by the condition of the second last flow node of Figure 3-7. One assumes that the condition in this
22
MSc Thesis: May 28, 2013
particular node equalizes the equilibrium depth. The depth information in the flow node comes from
the right network node:
Eq. 3-14
In the second last flow node to be considered, the depth of Eq. 3-14 equals zero. Due to this matter
the downstream boundary conditions equals the equilibrium depth minus the vertical difference
over the length of one grid cell:
Eq. 3-15
The observation point is found in the second first flow node of Figure 3-7, that is 20 flow links
upstream of . Applying the correction of the bed slope and adding the vertical difference
over one grid cell to obtain the computed equilibrium depth over the full 20*500=10.000m, gives the
normative computational output:
Eq. 3-16
As one may notice, the accuracy compared to the analytical solution is much higher compared to the
outcome of Eq. 3-13. Solving Eq. 3-10 with an accuracy of 9 decimals for the parameters of Table 3-1
gives an equilibrium depth of:
which makes the numerical outcome of Eq. 3-16 accurate by O(e-09) compared to the analytical
solution.
23
MSc Thesis: May 28, 2013
The steady state solution of both models show acceptable accurate results. Although both model
settings are kept the same, the Delft3D model gives a smaller relative difference to the analytical
solution than D-Flow FM. The ratio of the order of magnitude of this relative difference is O(e-08) to
O(e-06). In D-Flow FM the largest errors were found upstream of the flow model, near the discharge
boundary, due to the upwind solver.
A feature is built for D-Flow FM which does one forward step in space only at the boundary. Although
the robustness of this solution still needs to be investigated, the numerical outcome is much better.
The relative error for this specific case is reduced to O(e-13).
These findings are of interest for further model comparisons to more complex case studies.
24
MSc Thesis: May 28, 2013
3.7.1 Introduction
A model simulation is carried out in D-Flow FM to verify the nearshore hydrodynamics induced by a
tidal forcing for the Duck model. Special attention is given to the flooding/drying principle near the
waterline. The emphasis of this test case is on the implementation of the boundary conditions for
nearshore models. For the sake of simplicity, a stationary moment in time is modelled. These
performances are tested against the Delft3D model. Both model settings are kept similar.
In the simulation, wave effects and morphodynamics are excluded to achieve the main objective of
the test case. The obstruction of the pier is removed from the model in order to make the model
quasi alongshore uniform. The model bathymetry is based on the field data from SandyDuck ’97,
which shows more or less an alongshore uniform bathymetry.
25
MSc Thesis: May 28, 2013
The threshold depth, above which a grid cell is considered to be wet, is set from e-01m to e-04m.
These settings are of main interest near the waterline, where the flooding/drying principle plays a
dominant role. For a rising tide (flooding of the cells), locally high velocities at the shoreline could
appear due to sudden flooding of the cells. By lowering the threshold depth, this instant flooding is
prevented. Flooding of the cells happen much smoother in this sense and prevents numerical
instabilities near the waterline.
√
̂ 6.23e-006 [rad]
Phase difference:
√
3°
26
MSc Thesis: May 28, 2013
∑ ̂ ( ( ))
̂ ( ( ))
6.23e-006 ( ( ))
6.22146e-006 [-]
Similar the water level boundary and the Northern Neumann boundary are deduced. These
amplitudes are used as model input for both models. In D-Flow FM, the sign of the water level
gradient determines model inflow or outflow. A positive specified Neumann amplitude leads to
inflow. The ebb velocities of the Duck model are towards the north, which holds for a model outflow
at the Northern boundary and thus a negative sign. The applied boundary settings for the D-Flow FM
model are presented below.
Neumann north Water level north Water level south Neumann south
Period (min) 0.0 0.0 0.0 0.0
Amplitude (ISO) -6.23000e-006 0.0 0.02617 6.22146e -006
Phase (deg) 0.0 0.0 0.0 0.0
The model settings of both Delft3D and D-Flow FM are kept similar for a one-to-one model
comparison. The computational time step is fixated and kept similar to Duck Delft3D ( ). As
well as the initial conditions, friction definition ( ), threshold depth ( ), grid
resolution and depth file.
The relative performance of the model forcing is tested by means of the Root Mean Square Error
(RMSE). This parameter shows the overall model performance relative to the Delft3D model. The
model parameters to be considered are the water level and current velocity:
27
MSc Thesis: May 28, 2013
Both the water levels and the current velocities show a relative good match to the Delft3D model.
For a further analysis, the current velocities of both models are investigated. The results are plotted
in Figure 3-9.
Figure 3-9: Depth-averaged velocities of the Duck model with 2x Neumann and one offshore water level boundary.
These results show qualitative similarities, which verify stationary Neumann boundaries in D-Flow
FM. The model differences are analysed in more detail by subtracting both models from one another.
These results are shown in Figure 3-10, by means of a plan view and a cross section at y = 2000m.
Two locations of interest are found from the difference plot; near the waterline and at the offshore
boundary. Despite tuning the flooding/drying settings, numerical differences are found at the
waterline. Offshore, the differences between both models are rather significant. Near the Southern
boundary (at y < -1000m) this difference is hardly noticeable but reinforces moving north. Both these
model deviations should be taken into account for further developments.
28
MSc Thesis: May 28, 2013
Figure 3-10: Velocity differences of the Duck model with Delft3D minus D-Flow FM
Delft3D D-Flow FM
Duration simulation [hours] 24h 24h
Averaged time step [seconds] 6s 5.999s
Spin-up time [hours] 4h 4h
Computation time [minutes] 15min 25min
The spin-up time for both models shows to be more or less the same which based on the velocity
fields. The significant difference in computation time could mainly be induced by the visualizer of the
graphical user interface of D-Flow FM. In according to Kernkamp et al. (2011), computation times of
D-Flow FM are in the same order of magnitude as WAQUA, TRIWAQ and Delft3D-FLOW.
29
MSc Thesis: May 28, 2013
The model forcing is simplified to an offshore sloping water level boundary and two lateral Neumann
boundaries. The settings are derived from the validated Delft3D model and tuned for a stationary
case. Attention is given to the numerical implementation of the boundary settings for both models.
Other model settings are kept similar as well, with special attention to the flooding/drying settings.
Based on the model comparison, the results of the water levels and velocities show realistic and
accurate results. This concludes a qualitative correct formulation of stationary Neumann boundaries
in D-Flow FM. Though, predominated model differences are found at the waterline and at the
offshore region. These differences should be taken into account for using Neumann boundaries for
further developments.
In addition, further research is required for validating harmonic Neumann boundaries in D-Flow FM.
30
MSc Thesis: May 28, 2013
Chapter 4
4.1 Introduction
Modelling the evolution of wind-generated waves in D-Flow FM requires a coupling to a wave
module. In Delft3D the standardized wave module is the third-generation SWAN model (Ris et al.,
1999; Booij et al., 1999). For the online coupling of Delft3D + SWAN, the communication file is used
for data transfer between both modules. By means of this report a new couple technique is
presented for D-Flow FM + SWAN. With the data transfer through the internal memory, a much
faster coupling tool is established which allows interactive modelling.
The layout of this chapter starts with an overview of the applied software methodology. The dynamic
coupling tool of D-Flow FM + SWAN is highlighted with the use of innovative platforms. To establish
this coupling system, the encoding of D-Flow FM is adapted by the inclusion of wave effects. This is
done by means of this MSc work and mainly based on the Delft3D encoding. These adaptions are
tested and validated by means of two nearshore schematisations.
An overview of the software methodology used for the D-Flow FM + SWAN coupling is given in
Diagram 1. The layout of this structure will be discussed in the following paragraphs. On the top, the
main program (framework) is positioned as an executable. Next the ESMF (Earth System Modelling
Framework) platform is built, where the main model loop through time is tracked and data
31
MSc Thesis: May 28, 2013
regridding is performed. Next, the Basic Model Interface (or BMI) is used to enable model
communication between D-Flow FM and SWAN, but also with DeltaShell for example. And finally, by
the use of an API (Application Programming Interface), the dynamic-link libraries are run. Note in the
diagram the colored arrows, which represent the allocated variables through the online
communication.
DFlowFM_SWAN.exe
DFLOWFM
DFLOWFM BMI SWAN BMI
DeltaShell
DFLOWFM SWAN
engine engine
By means of this architecture multiple programs (or modules) can ‘communicate’ with one another,
which could be interesting for modelling other physical processes like ecology for example. In
addition, this architecture enables dynamic modelling. The engines do not run the model from start
to end, but can be run per time step with intermediate model adjustments to the bathymetry or
wave conditions for example. Here for the model is only initialization once, namely at the start of the
time loop. Within the time loop the results are directly visualized to the user.
Another advantage of this dynamic layout is the data transfer via the internal memory. In contrast to
data transfer via disk (or file), data transfer by the internal memory is much faster, but somewhat
more complicated to realize. Encoded variable-positions are allocated after each computational
32
MSc Thesis: May 28, 2013
performance of the module. With modern technology (2013) this speeds up the process with a factor
of 250 to more or less 50 GB/s. With relative short coupling intervals, this difference becomes
significant.
This makes the D-Flow FM + SWAN coupling efficient from a computational perspective, but also
from the users’ perspective. By the work of Donchyts et al. (2013) the user can run the model,
visualize the model and edit the model all at the same time. This interactive way of modelling gives a
whole new nature to coastal modelling which promotes general model understanding.
These improvements of the dynamic coupler show high potentials for the field of coastal engineering
and oceanography. The interactive coupling tool provides a much shorter feedback loop of the model
response, promoting the model understanding. This is illustrated by Figure 4-1 which shows the
interactive improvements.
Figure 4-1: Left the conventional run by an executable. Right an ‘interactive’ run by means of a dynamic-link library.
The choice for using the ESMF framework is based on the following aspects:
Good support for coupling structured and unstructured models
Written in a clear Fortran interface; (Delft3D and D-Flow FM are both written in this same
programming language).
Supports MPI’s (Message Passing Interfaces), e.g. the Linux cluster of Deltares
In- and output files of ESMF regridding application are also in netCDF format.
33
MSc Thesis: May 28, 2013
Main functionalities of the ESMF framework are tracking the model time-loop, managing when and
how individual models are run and data interpolation. Both D-Flow FM and SWAN are managed at
once on this ESMF platform. Both modules can run parallel with different time steps. When both
flow- and wave computations are finished, the relevant data is interpolated by the ESMF regridding
and set as model input for the new computations. In this sense the online coupling is managed by the
ESMF framework. Other applications of the ESMF framework are discussed by the work of Collins et
al. (2005) and Killeen et al. (2006), and are used for climate models like NOAA GFDL model and NASA
GEO-5 model.
The development of the BMI application is mainly the work of Peckham et al. (2012), which
distinguishes three main parts in a model run. At the start of the time loop the model is initialized,
next a model update is performed within the time loop and finishes with a finalize step. An example
of a BMI application processed with MATLAB is given below:
This example describes the essence of the BMI functionality. The dynamic-link library of D-Flow FM
(unstruc.dll) is first initialized by using the model input file (input.mdu). The model is updated every
time step (dt) and retrieves (in this example) the allocated water levels (s1) for that time step. In the
end of the time loop the model is finalized. This is a classic example of the Hollywood principle. The
dynamic-link library of SWAN is called in a similar way.
Note that the model initialization is only performed once; namely in the beginning of the
computation. Now that the model initialization of both D-Flow FM and SWAN is only performed
34
MSc Thesis: May 28, 2013
once, the computational efficiency is significantly increased. Rather than the conventional Delft3D
routine, where an initialization is performed every time a module is started.
In the BMI-example the water level parameter is retrieved (bmi_var_get) after each model update.
Similar, parameters are set in the model (bmi_var_set). This setting- and getting of model
parameters is shown by the colored arrows in Diagram 1, and is the fundamental essence in the
online coupling technique.
By means of this BMI platform a dynamic communication with SWAN is realized. With this interface,
coupling to other (open source) software packages like DELWAQ or XBeach is made possible. The
modelling network can be managed with a graphical user interface (GUI) like DeltaShell (Donchyts
and Jagers, 2010).
4.3.1 Introduction
The wave parameters computed by SWAN give rise to a change in the nearshore hydrodynamics.
Under the 2DH assumption, the bed shear stresses are enhanced and the wave forces trigger
alongshore currents and a water level change in the nearshore. These effects can be significant,
depending on the amount of wave energy. For a planar beach, current velocities of 0.3m/s arise with
waves of (only) 1.0m high. For a beach with rip channels these velocities can be much higher.
The wave formulations implemented in D-Flow FM are based on the Delft3D formulations. Main
differences are the method of integration and the data interpolation. This section discusses the
numerical implementation of the wave formulation in D-Flow FM.
A similar procedure is used for the depth-averaged model of D-Flow FM. Considering the wave-
averaged momentum equation of D-Flow FM and substituting the wave-induced forces as a source
term, gives the following expression:
〈⃗ 〉
〈⃗ 〉 〈 〉 〈⃗⃗⃗ 〉 〈⃗ 〉 Eq. 4-1
where 〈 〉 denotes a time averaging over the wave motion, ⃗⃗⃗ the vectorial bottom shear stress and
[ ] representing the substituted wave forces. The bed shear stresses (⃗⃗⃗ ) from Eq. 4-1
35
MSc Thesis: May 28, 2013
consists of a wave related part and a current related part. Various wave-current interaction models
are built, to account for this phenomenon. Soulsby et al. (1993) developed a parameterization of
these models with one standard formula, each model having its own fitting coefficients. This
formulation is adopted in the D-Flow FM code, based on the Delft3D formulations.
The bed shear stress due to the current related part (| |) is a main function of the current velocity
squared and some kind of friction coefficient, depending on the kind of formulation. The expression
for the wave related part (| |) is proportional to the orbital velocity squared and to a wave friction
factor given by Swart (1974). The orbital velocity is computed from the linear wave theory and is
given by [Deltares, 2011, Delft3D-FLOW]:
√ Eq. 4-2
where represents the Root-Mean-Square wave height, and the radian frequency and wave
number respectively, and the water depth. The radian frequency is defined according to Eq. 2-9, by
over the wave period. In this definition the smoothed peak wave period ( ) is used from SWAN
instead of the standard peak wave period ( ). The is continuous in the frequency domain, rather
than the discrete . There is also an option to retrieve orbital velocities by SWAN, which is not in the
scope of this research. Summarized, the following parameters from SWAN are used to determine the
enhanced bed shear stresses due to waves (| |):
The bed shear stresses are computed in a separate OpenMP section to distribute the overall
computational workload over multiple cores of the processor. First the orbital velocities are
computed by Eq. 4-2. Next the enhanced bed shear stresses are computed based on the Soulsby
formulation. And finally, the friction term of Eq. 4-1 is updated. The wave forces are interpolated to
the velocity points and substituted as explicit source terms in Eq. 4-1. Next Eq. 4-1 is solved for the
new computational time step.
Delft3D takes account for these phenomena by the online coupling method to SWAN. Updated flow
conditions are sent to SWAN on a certain interval. With these updated flow conditions a new SWAN
computation is performed. The following parameters are updated in SWAN:
36
MSc Thesis: May 28, 2013
Water level
Bottom level
GLM velocity in x-direction
GLM velocity in y-direction
Wind velocity in x-direction
Wind velocity in y-direction
Validation and verification of SWAN is extensively carried out by Delft University of Technology.
Different test cases are validated by the work of Ris et al. (1999) and Booij et al. (1999), which are in
agreement with Gautier (2010). The following aspects are considered:
Verification of the proper implementation of formulations for physical processes that are
implemented in SWAN.
Determine the numerical accuracy and robustness of the SWAN model.
Verification of the input and output commands of SWAN.
Validation and verification of the performance of the SWAN model in idealized-, laboratory-
and field cases (with added winds and currents).
Eq. 4-3
Through the wave action density spectrum, wave action is conserved. In SWAN the action balance
equation is solved which for Cartesian coordinates is (Hasselmann et al., 1973):
Eq. 4-4
37
MSc Thesis: May 28, 2013
The first term on the left-hand side represents the inertia or ‘local rate of change’. The second and
third terms represent the rate of change by transportation or advection. The fourth term represents
the frequency shifting of the waves due to ambient currents. Here σ is defined as the radian
frequency in a system moving with the current (relative radian frequency) and is a function of depth
and ambient currents (Holthuijsen, 2007). The fifth term represents depth-induced and current-
induced refraction. The source term on the right of Eq. 4-4 consists of the contribution of winds, non-
linear wave-wave interaction and dissipation.
SWAN can be applied in a stationary mode so time has been omitted from Eq. 4-4. A wave record is
assumed to be stationary (Gaussian) over a period of 15–30minutes (Holthuijsen, 2007).The interval
of FLOW + WAVE coupling by the communication file is, for stationary assumptions, consequently set
to 30minutes. Non-stationary situations, i.e. coupling intervals smaller than 30minutes, are simulated
as quasi-stationary with repeated model runs. This requires a relative large computational domain,
i.e. large compared to the phase velocity of the waves through the computational domain (Deltares,
2011, Delft3D-WAVE).
The SWAN simulation is run on the computational grid for the x-, y-space, and a separate grid for the
θ-, σ-space. In geographic space an implicit upwind scheme is used as numerical solver with an
iterative four-sweep technique. A spatial discretization is used with the information on the vertices.
In spectral space the numerical schemes are implicit mixed upwind/central order schemes. These
implicit schemes in the four-dimensional propagation space ensure unconditionally stable wave
propagation (Booij et al., 1999).
The model initialization and model computation of SWAN are separated by means of SWAN BMI. The
SWAN model is only initialization once, namely at the start of the time loop. This is where the input
files of SWAN are defined, which are used throughout the main time loop. By initializing the model
only once, the computational efficiency is increased, compared to the conventional Delft3D-WAVE
simulation.
The wave parameters in SWAN are stored in one large output array (SWAN team, 2013). This variable
is adopted in the internal memory of the whole framework. This data is internally transferred via the
ESMF regridding to the internal memory of D-Flow FM. This online communication is schematized in
Diagram 1 by the colored arrows. An interesting concept with this configuration is the routine to
which D-Flow FM receives wave information. A SWAN computation is not called (as with the
executable configuration of Delft3D) but received when the main program says so. This is a classic
application of the Hollywood principle. These commands are managed on the ESMF platform.
38
MSc Thesis: May 28, 2013
The relevant hydrodynamic conditions from D-Flow FM are all known on the water level points of the
grid, apart from the wind velocities. In addition, the grid cell administration of D-Flow FM is kept at
the water level points (flow nodes). This is an essential difference compared to Delft3D, due to the
different polygonal cell shapes of the unstructured grid. Communicating the model parameters at the
water level points to the ESMF regridding, is in that sense obvious. By the ESMF regridding the data is
mapped to the vertices of the computational grid used by SWAN.
In general the SWAN domain is much larger than the flow domain due to implications from the wave
boundaries in SWAN. For a successful SWAN computation the flow conditions in the non-overlapping
areas must be known by data extrapolation. A robust and straightforward method is to copy and
extent the condition at the boundaries of the inner domain to the non-overlapping area.
Eq. 4-5
39
MSc Thesis: May 28, 2013
Bilinear interpolation:
The default interpolation method of the ESMF regridding application is bilinear interpolation. The
weighted value of the destination grid point is determined by a linear interpolation in one direction,
and then again in the other direction. Rewriting these terms for a structured grid gives the following
algorithm to generate the bilinear weight matrix:
Eq. 4-6
40
MSc Thesis: May 28, 2013
4.6.1 Introduction
The formulation of wave-induced bed shear stresses in according to Soulsby et al. (1993) is built in
the D-Flow FM code. The wave forces are included in the code by implementing an extra source term
in the momentum equation Eq. 4-1. The preliminary tests of this Fortran code is executed by the use
of spatial varying wave parameters from a Delft3D run with oblique incident waves. The wave
parameters are read from the communication file, converted to the correct grid administration of D-
Flow FM and placed as arrays in the code. In this sense the implemented wave formulation in D-Flow
FM is tested for primary purposes. The results are compared to Delft3D and a 1D wave energy
balance model.
The 1D wave energy balance model solves the equation in the form of Eq. B-2. This model is
MATLAB-based and is established by Deltares for primary purposes. For the sake of accuracy the
model is extended with the formulation of Soulsby et al. (1993) and filled by the Fredsøe 1984 fitting
coefficients.
The Delft3D model built for the specific case contains a rectangular grid with 50x50 cells on a
resolution of 50x50m. A linear sloping bathymetry is used similar to the 1D model, as schematised in
Figure 4-4. In order to make the flow stationary, the SWAN domain is 4 times larger than the FLOW
domain but has a similar resolution. In this sense the numerical disturbances from the lateral wave
boundaries do not influence the FLOW model. The directional spreading in SWAN is put to and
the bottom friction definition of JONSWAP is used with a default coefficient of 0.067m2/s3. An
offshore stationary water level condition is imposed with a reflection parameter alpha set to zero. In
the cross-shore two Neumann boundaries are imposed, with a water level gradient set to zero. In this
sense a stationary steady state solution is found for wave-induced hydrodynamics in the nearshore.
Figure 4-4: Planar beach schematised by 50x50 grid cells, 2 lateral Neumann conditions, one stationary offshore water
level condition and oblique incident waves.
41
MSc Thesis: May 28, 2013
The wave parameters are read from the communication file and pre-processed in MATLAB in
according to the grid administration of D-Flow FM. For this routine the feature ‘renumber (internal)
flow nodes’ is switched off, which normally serves for a faster matrix solver. The wave parameters on
the boundaries are defined by a copy of the nearest neighbour, similar as in Delft3D.
The same grid, bathymetry and boundary conditions are used from the Delf3D simulation. The
alongshore water level gradients on the lateral boundaries are put to zero due to the alongshore
uniform flow of the steady state solution. Stationary Neumann boundaries in D-Flow FM show
realistic results as was shown by the test case: Application to a tidal model (Duck schematisation).
The post-processing of the model results are in GLM-velocities (General Lagrangian Mean), which are
Eulerian velocities compensated by Stokes drift. The GLM-velocity is a measure for the velocities
observed at the surface, also called the total velocity, and is related to the Eulerian velocity by:
⃗ ⃗ ⃗ Eq. 4-7
where ⃗ is the GLM-velocity vector, ⃗ is the Eulerian velocity vector and ⃗ is the Stokes drift
vector (Bosboom and Stive, 2011).
Eq. 4-8
Eq. 4-9
The bed shear stress (⃗⃗⃗ ) is in general an important model parameter. The enhanced bed shear stress
due to waves modelled by D-Flow FM is validated for the stationary alongshore uniform current.
42
MSc Thesis: May 28, 2013
Figure 4-5: Alongshore velocities due to oblique incident waves on a planar beach. Upper panel are the Delft3D results
and lower panel the D-Flow FM results.
43
MSc Thesis: May 28, 2013
Figure 4-6: Cross-section of oblique incident waves approaching a planar beach. Delft3D model in blue and D-Flow FM
model in red.
The velocity patterns of Figure 4-5 show a general good agreement, which gives a primary confirm of
a correct implementation of the 2DH wave formulation in the D-Flow FM code. The cross-sectional
plot of Figure 4-6 shows a good match of the alongshore velocity profile, considering the primary
purposes of this test case. The alongshore velocities of Delft3D and the 1D model are of the same
order of magnitude.
The wave propagation modelled by the 1D model and Delft3D show striking results. The difference in
wave shoaling induces an offshore shift of wave energy dissipation. An extensive analysis is carried
out in Appendix C, but with no better results. To verify the Delft3D model, a stand-alone SWAN
computation is performed with the exact same input variables. The wave height and wave direction
of this SWAN model is also plotted in Figure 4-6, and coincides with the Delft3D model.
The wave forcing in cross-shore direction is verified by the water level plot of Figure 4-6, in according
to Eq. 4-8. The wave-induced water level setup in the nearshore of D-Flow FM show qualitative
similarities to the Delft3D simulation. The alongshore velocity profile of D-Flow FM shows a general
good match to the Delft3D velocity profile. These results are most accurate in shallow waters. At the
offshore regions differences between both models arise. At this region the effect of wave forces
decrease and bed shear stresses become dominant. From [Deltares, 2011, Delft3D-FLOW] is given
that in Delft3D the bed shear stresses are computed with the Eulerian velocities. In D-Flow FM the
GLM velocities are used, which is a plausible cause for the model differences. This Stokes drift
correction by means of Eq. 4-7 is yet (2013) to be implemented in D-Flow FM.
44
MSc Thesis: May 28, 2013
Table 4-1: Velocity values between Delft3D and D-Flow FM with different formulations of Soulsby et al. (1993):
In the post-processing of the Delft3D communication-files, the wave parameters are identical due to
the non-extension of the hydrodynamic FLOW results to WAVE. The wave parameters used as model
input for the D-Flow FM simulations, have thus not been changed. The results of Table 4-1 show a
correlation between the RMS values and the model deviation. A linear up scaling of the model
deviation can be recognized with respect to the velocity magnitude.
The water level setup and velocity profile of D-Flow FM shows a qualitative good match to the
Delft3D simulation and the 1D model. This suggests a correct implementation of the wave
formulations in the D-Flow FM.
From a closer perspective is found that the effect of Stokes drift is not neglectable at offshore
regions. The enhanced bed shear stress computed by D-Flow FM is yet (2013) to be extended with
the Stokes drift correction. Despite this effect, the overall model performance relative to Delft3D is
considered good.
45
MSc Thesis: May 28, 2013
4.7.1 Introduction
The wave formulations in D-Flow FM are tested for the Egmond case. This barred beach has been a
study-ground for years now (Elias, 1999; Luijendijk et al., 2010; Winter, 2012). The coastline of
Egmond attracts many bathing people in summer time but is also prone to strong rip currents. This
combination makes this coastline potentially dangerous for swimming, which is a strong motivation
for coastal research.
Models are built to give insight in the hydrodynamic patterns in the nearshore and to predict the
possible occurrence of rip currents. Different model packages are used, but a widely applied tool is
Delft3D. These settings are calibrated and validated with different field data applications (i.e. wave
buoys, drifter data, Argus images) by i.a. Elias (1999) and Winter (2012).
A replica of the hindcast model from Winter (2012) is built in Delft3D, for the field survey period of
22-8-2011 until 26-8-2011. The model includes nearshore hydrodynamics induced by wind, waves
and tides for a non-uniform barred coastline. The model results are compared to field measurements
by means of drifter data. This validated Delft3D model forms the basis of this test case. Next this
Delft3D model is simplified by considering wave-induced hydrodynamics only. A similar D-Flow FM
model is built and is tested to the Delft3D model.
46
MSc Thesis: May 28, 2013
Clearly developed rip channels have been reported and captured on Thursday 25-8-2011 around
18.00GMT. This was during falling tide, 1hour prior low tide. The offshore wave conditions obtained
are a significant wave height of 0.59m from the south-west, a peak wave period of 5.1sec, an angle of
incidence of 48° and neglectable winds.
This stationary moment in time is used to hindcast the field measurements and to validate a D-Flow
FM wave model of the Egmond coastline. A time independent model is used where tidal influences
are neglected. The lateral boundary conditions are two zero-Neumann boundaries and an offshore
water level boundary of -0.50m. The wave condition implemented in Delft3D is a stationary 2D
spectrum file (*.sp2) for the wave conditions of 18.00GMT. The spectrum file is obtained by a nesting
routine from a Dutch coastline model, nested in a continental shelf model, nested in a global model.
These models are forced by astronomical tidal forcing and wind fields, based on metrological files.
The wave forcing in the D-Flow FM model is built with the same routine as the test case: Application
of wave-driven alongshore currents (planar beach). The wave parameters are retrieved from the
communication-file of a previous Delft3D simulation, and put as parameter-arrays in the D-Flow FM
code. In the online coupling of Delft3D-FLOW with SWAN, the hydrodynamics results from Delft3D-
FLOW are not used as model input for SWAN. In this sense it represents the D-Flow FM simulation
with the stationary wave parameters as model input.
The same model settings as the Delft3D simulation are used in D-Flow FM, which gives a one-to-one
model comparison. The simulations are run on the same rectangular grids. Though, the powerful
feature of D-Flow FM is to apply high resolutions in areas of interest. Unstructured grids show high
potentials for models like the Egmond case, especially for capturing rip currents on high resolution.
The rectangular grid (which is used) has a grid resolution of 5x10m in the rip channels and 20x10m
offshore. Here for 131x138 grid cells are used, 17 938 cells in total. This is a computational extensive
job, given that the high resolution at offshore regions is not necessary.
47
MSc Thesis: May 28, 2013
Figure 4-8: Delft3D Egmond model compared to drifter measurements. Inclusion of wave-, tide- and wind effects.
The nearshore velocity patterns of the Delf3D model shows a general good match to the drifter data.
The velocity field of the rip channels match the trajectories of the drifters rather well. On the
sandbanks the increase in wave action is clearly modelled by an increase in the velocity field, feeding
the rip channel. Local velocity differences are likely caused by effects on a much smaller scale than a
5x10m resolution or 3D effects. One can conclude a realistic model representation of the 2DH
Delft3D Egmond model.
48
MSc Thesis: May 28, 2013
for a one-to-one model comparison. The SWAN computation is performed on the same
computational grid as the reference case. Other flow-model settings are kept similar, with special
attention to the following settings:
boundary conditions
initial conditions
friction definition
horizontal eddy viscosity
depth file
depth definition
threshold depth
Soulsby formulation
49
MSc Thesis: May 28, 2013
Figure 4-9: Wave simulation of the Egmond model. Nearshore GLM-velocities of the Delft3D model (left) compared to the
D-Flow FM model (right). Points of interest are the rip channel and the nearshore velocity extreme.
50
MSc Thesis: May 28, 2013
Figure 4-9 highlights the nearshore hydrodynamics of both models, solely induced by a stationary
wave field. A clear rip channel is shown which is driven by a high alongshore velocity in the
nearshore. The location of this rip channel coincides with the Southern rip channel of Figure 4-8. A
time-series plot is made of the absolute velocities at this location (102 818; 514 920) and at the
nearshore extreme (102 850; 514 630).
Figure 4-10: Time-series plot of the depth-averaged velocity; inside the rip channel and at the nearshore extreme.
Despite the fact that this computation is stationary in time, this time-series plot shows some
interesting results. The spin-up between Delft3D and D-Flow FM is more or less the same. After the
spin-up time, D-Flow FM still shows distortions in the velocity signal. The lower panel of Figure 4-10
clearly shows these distortions, while the Delft3D model shows a smooth stationary solution. This
clarifies the constraint for using the CFL-criteria for the D-Flow FM simulation to make the model
numerically stable. These instabilities are most likely induced by the explicit solver of the advection
term.
From the lower panel of Figure 4-10 is shown that the nearshore velocities are slightly overestimated
by D-Flow FM, relative to Delft3D. These nearshore velocities ‘feed’ the rip channel by means of
advected flow. This is clearly visual in Figure 4-9. Despite the overestimate in the nearshore extreme,
the velocities inside the rip channel are underestimated by 10-15%, as shown in the upper panel of
Figure 4-10. These model differences are likely induced by the different advection solver of D-Flow
FM.
51
MSc Thesis: May 28, 2013
Figure 4-11: Velocity difference of the Egmond model; Delft3D minus D-Flow FM.
These results validate the results of the test case: Application of wave-driven alongshore currents
(planar beach). At the offshore region effects of wave-induced forces are small. The enhanced bed
shear stresses are dominant at these regions. These differences are likely caused due to the exclusion
of Stokes drift in D-Flow FM. In Delft3D the enhanced bed shear stresses are computed by means of
the Eulerian velocities, which are corrected for Stokes drift by Eq. 4-7 (Deltares, 2011,Delft3D-FLOW).
52
MSc Thesis: May 28, 2013
Other striking results from Figure 4-11 are the local velocity differences in the breaker zone and at
the shoreline. Model differences at the shoreline were also found in the test case: Application to a
tidal model (Duck schematisation), where this is explained by the flooding/drying principle. These
numerical differences are in the same order of magnitude as Figure 4-11. In addition, the velocity
differences in the breaker zone are likely caused by the advection solver of D-Flow FM which is
different compared to Delft3D. This effect is clearly noticeable inside the rip channel, where
advection is dominant.
The overall model performance of D-Flow FM compared to Delft3D is assessed by means of the Root
Mean Square Error (RMSE). The model parameters to be considered are the water level and the
absolute current velocity. The Root Mean Square (RMS) of the model parameters is also shown which
is aggregated over space:
The water level differences are some order of magnitude smaller than the velocity differences. The
underestimated velocities of D-Flow FM are predominated at intermediate water depths, as shown
in Figure 4-11. The model differences inside the surf zone are much smaller.
From the model comparison is concluded that the velocities of both models (in magnitude and
direction) show good similarities inside the surf zone. The spatial varying velocities in the nearshore
are in agreement and the rip channel modelled by D-Flow FM is captured well with velocity fields
that match the Delft3D simulation accurate by 85-90%.
The results presented are in agreement with the test case: Application of wave-driven alongshore
currents (planar beach). The effect of excluding Stokes drift in D-Flow FM is not neglectable at
offshore regions where predominated model differences are found. In the nearshore, the effect of
the different advection solver is clearly captured. D-Flow FM gives critical velocities at the nearshore
due to the explicit scheme. These nearshore velocities (of the advected flow) are slightly
underestimated by D-Flow FM, which is noticeable inside the rip channels.
From this qualitative assessment is concluded that the implemented wave formulations in D-Flow FM
give accurate results for primary purposes.
53
MSc Thesis: May 28, 2013
54
MSc Thesis: May 28, 2013
Chapter 5
5.1 Introduction
The fundaments of an innovative coupling tool are built to run a D-Flow FM + SWAN computation
efficiently. In the development of this coupling tool, the code of D-Flow FM is adapted by the
inclusion of wave effects. This is done by means of this MSc work and mainly based on the Delft3D
encoding. The adapted D-Flow FM model is tested and validated by means of two nearshore models.
Improvements are made for the D-Flow FM + SWAN coupling, compared to the conventional
coupling of Delft3D + SWAN. Parallel computing, combined with communication via internal memory
and a single initialization of the modules, makes the D-Flow FM + SWAN coupling significantly faster.
In addition, the coupling tool is made interactive for the user with a much shorter feedback loop of
the model response. The user can run the model, visualize the model and edit the model all at the
same time. These improvements allow great opportunities in the fields of coastal engineering and
oceanography.
Four test cases are considered. The first two test cases validate solely current related hydrodynamics
for stationary Chézy flow and tidal forcing, respectively. Next the wave formulations in D-Flow FM
are validated for a planar beach case and the Egmond case, respectively. Oblique incident waves are
considered for an alongshore uniform coast and for a non-uniform barred beach. The performance of
D-Flow FM is analysed by these test cases and based on this analysis, recommendations are
presented for improving nearshore modelling with D-Flow FM.
The second test case simulates the tidal forcing for an alongshore uniform beach. The M2-tide is
modelled by means of a stationary water level boundary and two stationary Neumann boundaries.
Based on the model comparison, the resulting water levels and velocities show realistic and accurate
results, leading to the conclusion that the formulation of stationary Neumann boundaries in D-Flow
55
MSc Thesis: May 28, 2013
FM is qualitatively correct. However, further research is required for validating harmonic Neumann
boundaries in D-Flow FM. Here for, only stationary boundary conditions are used for the wave
simulations.
The wave formulations implemented in D-Flow FM are validated by considering two wave
simulations: an alongshore uniform model (planar beach case) and a non-uniform barred beach
(Egmond case). Both models are forced by a stationary wave field, which triggers a stationary velocity
field and water level difference in the nearshore. The model performances are tested against
validated Delft3D models.
Wave parameters are retrieved from a Delft3D + SWAN model and used as model input for the D-
Flow FM simulation. The wave parameters are read from the communication file and are mapped to
the D-Flow FM model by studying the grid administrations. Next, the parameter arrays are copied
into the D-Flow FM code to mimic the SWAN coupling. The following wave parameters are used as
model input:
The model settings and grid resolution of D-Flow FM and Delft3D are kept similar. In this sense a one-
to-one model comparison is made. The model response is tested by analyzing the following
parameters:
Based on the results of both wave cases, it can be concluded that wave-driven hydrodynamics
modelled by D-Flow FM are accurate within 85-90% with respect to nearshore Delft3D models. The
water levels show a good match, as confirmed by the planar beach case. The wave-induced velocity
patterns, both in magnitude and direction, are in agreement with Delft3D, as demonstrated by the
Egmond case. The wave-driven hydrodynamics in the nearshore are captured well, including the rip
channel. Based on these conclusions is found that the implemented wave formulations in D-Flow FM
are successfully validated.
Based on the results, it is concluded that predominant model differences are found at intermediate
water depths. From the planar beach case it is found that the effect of Stokes drift is not neglectable
and may give a significant contribution to these model deviations. In Delft3D, the enhanced bed
shear stresses are computed with Eulerian velocities, rather than GLM-velocities. This essential
difference is confirmed by the Egmond model. In addition, the difference between both advection
solvers is evident in nearshore regions for the Egmond case, where the absolute velocities inside the
rip channel were slightly underestimated by D-Flow FM.
56
MSc Thesis: May 28, 2013
5.3 Recommendations
Wave modelling by D-Flow FM:
Include depth-averaged Stokes drift. The enhanced bed shear stresses are computed with
GLM-velocities in D-Flow FM, rather than Eulerian velocities as in Delft3D.
Validate the model response for time-varying wave fields. In the test cases, only stationary
wave fields are considered where inertial effects of flow are neglected.
Verify SWAN BMI in the D-Flow FM + SWAN coupling. The model initialization and model
computation of SWAN are separated by means of the BMI application. In addition, SWAN is
called by means of a dynamic library rather than an executable, as in Delft3D. Variables are
set to- and retrieved from the internal memory at each SWAN computation. These
modifications to SWAN need to be verified.
Build a feature to run multiple nested SWAN models. D-Flow FM is able to model large
coastal areas with a single mesh, but with the inclusion of waves, a nesting routine for the
SWAN domains is required.
Couple D-Flow FM to UNSWAN. This is the unstructured module of SWAN (Zijlema, 2010),
which shows high potential for coastal modelling in combination with D-Flow FM.
Verify the ESMF regridding functionality with the different interpolation methods.
3D implementation of the wave formulations:
o Vertical mass flux due to Stokes drift.
o Streaming near the bottom, as a consequence of a wave-induced currents in the
wave boundary layer.
o Turbulent effects.
o Adjustment of bottom shear stress formulation.
57
MSc Thesis: May 28, 2013
Verify implemented viscosity parameter in D-Flow FM. Model runs performed in debug mode
with tuned viscosity parameters caused program errors; thus, the viscosity parameter is set
to zero for the wave simulations.
58
MSc Thesis: May 28, 2013
References
Booij N., Ris, R.C., Holthuijsen, L.H., (1999). A third-generation wave model for coastal regions. Part I:
Model description and validation. Journal of Geophysical Research: Oceans (1978-2012), Volume 104,
No. C4, pp. 7649-7666.
Borsboom, M.J.A., Doorn, N., Groeneweg, J., van Gent, M.R.A., (2000). A Boussinesq-type wave model
that conserves both mass and momentum. Proc. 27th, International Conference on Coastal
Engineering, Sydney, Australia.
Broomans, P., (2003). Numerical accuracy in solutions of the shallow-water equations. MSc Thesis,
Delft University of Technology, EWI, Numerical Analysis.
Casulli, V. and Walters, R.A., (2000). An unstructured grid, three-dimensional model based on the
shallow water equations. Int. J. Numer. Meth. Fluids, 32:331–348.
Casulli, V. and Zanolli, P., (2002). Semi-implicit numerical modeling of nonhydrostatic free-surface
flows for environmental problems. Mathematical and Computer Modelling 36, 1131-1149.
Collins, N., Theurich, G. DeLuca, C., Trayaonv, A. Li, P., Yang, W., Hill, C. (2005). Design and
implementation of Earth System Modeling Framework components. International Journal of High
Performance Computing Applications.
Deltares (2011). Delft3D-QUICKPLOT Visualisation and animation program for analysis of simulation
results. Delft, Deltares.
Deltares (2011). Delft3D-WAVE Simulation of short-crested waves with SWAN. Delft, Deltares.
Dietrich, J.C., Zijlema, M., Westerink, J.J., Holthuijsen, L.H., Dawson, C., Luettich Jr., R.A., Jensen, R.E.,
Smith, J.M., Stelling, G.S., Stone, G.W., (2011). Modeling hurricane waves and storm surge using
integrally-coupled, scalable computations. Coastal Engineering, 58, 45-65.
59
MSc Thesis: May 28, 2013
Dingemans, M.W., Radder, A.C., de Vriend, H.J., (1987). Computation of the driving forces of wave-
induced currents. Coastal Engineering 11, pp. 539-563.
Donchyts, G.V., Baart, F., van Dam, A., Jagers, H.R.A., (2013). The joy of interactive modelling.
European Geoscience Union, Vienna.
Donchyts, G.V. and Jagers, H.R.A., (2010). DeltaShell – an open modelling environment. International
Congress on Environmental Modelling and Software (iEMSs), Ottawa, Intario, Canada.
Dykes, J.D., Hsu, Y.L., Kaihatu, J.M., (2003). Application of Delft3D in the nearshore zone. Proc. 5th
AMS Coastal Conference, Seattle, Washington, 57-61.
Elias, E.P.L., (1999). The Egmond mode – Calibration, validation and evaluation of Delft3D-MOR with
field measurements. Delft University of Technology, Institutional Repository, Z2394.20/Z2396.20.
Evans, M.W., Harlow, F.H., (1957). The Particle-In-Cell method for hydrodynamic calculations. Los
Alamos Scientific Laboratory report, LA-2139, University of California.
Gautier, C., (2010). SWAN calibration and validation for HBC2011. Delft, Deltares.
Gerritsen, H., de Goede, E.D., Platzek, F.W., Genseberger, M., van Kester, J.A. Th.M, Uittenbogaard,
R.E., (2007). Validation document Delft3D-FLOW – a software system for 3D flow simulations. Delft
Hydraulics [s.I.], 266 pp.
Goede, E. de, (2011). Validatie van Villemonte overlaatformulering in WAQUA met praktijkmetingen.
Rijkswaterstaat – Deltares, Institutional Repository, Delft.
Hasselmann, K., Bernett, T.P., Bouws, E., Carlson, H., Cartwright, D.E., Enke, K., Ewing, J., Gienapp, H.,
Hasselmann, D.E., Kruseman, P., Meerburg, A., Müller, P., Olbers, D.J., Richter, K., Sell, W., Walden,
H., (1973). Measurements of wind wave growth and swell decay during the Joint North Sea Wave
Project (JONSWAP). Deutsches Hydrographisches Institut Hamburg, Ergänzungsheft zur Deutschen
Hydrographischen Zeitschrift Reihe A (8), Nr. 12.
Havinga, H., Visser, P.J., de Vriend, H.J., Wang, Z.B., (2006). River Engineering. Reader, Delft
University of Technology, Institutional Repository.
Hill, C., DeLuca, C., Balaji, V., Suarez, M., da Silva, A., (2004). Architecture of the Earth System
Modeling Framework. Computing in Science and Engineering, Vol. 6, Nr. 1, pp. 18-28.
Holthuijsen, L.H. (2007). Waves in Oceanic and Coastal Waters. Cambridge University Press.
Hsu, Y.L., Dykes, J.D., Allard, R.A., Kaihatu, J.M., (2006). Evaluation of Delft3D performance in
nearshore flows. NRL Memorandum Report, Naval Research Laboratory, Stennis Space Center, MS
39529-5004.
60
MSc Thesis: May 28, 2013
Jones, J.E. and Davies, A.M., (2006). Application of a finite element model (TELEMAC) to computing
the wind induced response of the Irish Sea. Continental Shelf Research, 26, 1519-1541.
Kernkamp, H.W.J., van Dam, A., Stelling, G.S., de Goede, E.D., (2011). Efficient scheme for the shallow
water equations on unstructured grids with applications to the Continental Shelf. Ocean Dynamics,
61:1175-1188.
Killeen, T., DeLuca, C., Gombosi, T., Toth, G., Stout, Q., Goodrich, C., Sussman, A., Hesse, M., (2006).
Integrated Frameworks for Earth and Space Weather Simulation. American Meteorological Society
Meeting, Atlanta, GA.
Kirkpatrick, M.P., Armfield, S.W., Kent, J.H., (2003). A representation of curved boundaries for the
solution of the Navier-Stokes equations on a staggered three-dimensional Cartesian grid. Journal of
Computational Physics 184, 1-36.
Kramer, S.C., Stelling, G.S., (2008). A conservative unstructured scheme for rapidly varying flows. J
Numer Meth Fluids 58:183–1212.
Longuet-Higgins, M. S. and R.W. Stewart, (1964). Radiation stresses in water waves; a physical
discussion, with applications. Deep-Sea Research, 11, 4, 529–562.
Luijendijk, A.P., Henrotte, J., Walstra, D.J.R., van Ormondt, M., (2010). Quasi-3D modelling of surf
zone dynamics. Delft University of Technology, Institutional Repository.
MacWilliams, M.L., Gross, E.S., DeGeorge, J.F., Rachiele, R.R., (2007). Three-dimensional
hydrodynamic modeling of the San Francisco Estuary on an unstructured grid. IAHR, 32nd Congress,
Venice, Italy.
Peckham, S.D., Hutton, E.W.H., Norris, B., (2012). A component-based approach to integrated
modeling in the geosciences: The design of CSDMS. Computers & Geosciences: Modeling for
Environmental Change, 53:3 – 12.
Reniers, A.J.H.M., (1999). Longshore current dynamics. Delft University of Technology, Institutional
Repository.
Rijn, L.C. van, (2007a). Unified view of sediment transport by currents and waves. I: initiation of
motion, bed roughness and bed-load transport. Journal of Hydraulic Engineering, vol. 133, no. 6.
Ris, R. C., (1997). Spectral Modelling of Wind Waves in Coastal Areas. Delft University of Technology,
Institutional Repository.
Ris, R.C., Holthuijsen, L.H., Booij, N., (1999). A third generation wave model for coastal regions. Part
II: Verification. Journal of Geophysical Research: Oceans (1978-2012), Volume 104, No. C4, pp. 7667-
7681.
61
MSc Thesis: May 28, 2013
Roelvink, J.A., Petit, H.A.H., Kostense, J.K., (1992). Verification of a one-dimensional surfbeat model
against laboratory data. Coastal Engineering. Part II: Long Period Waves, Storm Surges and Wave
Groups, pp. 960-973.
Roelvink, J.A., Stelling, G.S., Hoonhout, B.M., Risandi, J., Jacobs, W., (2012). Development and field
validation of a 2DH curvilinear storm impact model. ICCE 2012: Proceedings of the 33rd International
Conference on Coastal Engineering, Santander, Spain, 1-6 July 2012.
Roelvink, J.A. and Stive, M.J.F., (1989). Bar-generating cross-shore flow mechanisms on a beach.
Journal Geophysical. Research., Vol. 94, pp. 4785-4800.
Roelvink, J.A. and Walstra, D.J.R., (2004) Keeping it simple by using complex models. In Proceedings of
the 6th International Conference on Hydro-Science and Engineering. Advances in Hydro-Science and
Engineering, vol. VI, page p. 12. Brisbane, Australia. 217-219.
Sørensen, O., (2004). Floodplain modelling using unstructured finite volume technique. DHI Technical
Note, January.
Soulsby, R.L., (1997). Dynamics of marine sands: A manual for practical applications. Thomas Telford
Publications, London.
Soulsby, R.L., Hamm, L., Klopman, G., Myrhaug, D., Simons, R.R., Thomas, G. P., (1993). Wave-current
interaction within and outside the bottom boundary layer. Coastal Engineering 21: 41-69. 243, 245,
246, 247.
Stelling, G.S., (1983). On the construction of computational methods for shallow water flow problems.
Delft University of Technology, Institutional Repository.
SWAN team, (2013). SWAN Cycle III version 40.91AB. User manual, Delft University of Technology,
Faculty of Civil Engineering and Geosciences Environmental Fluid Mechanics Section.
Swart, D.H., (1974). Offshore sediment transport and equilibrium beach profiles. Delft University of
Technology, Institutional Repository. Delft Hydraulics Publ. 131. 248, 366.
Van Dongeren, A., Reniers, A., Battjes, J., and Svendsen, I., (2003). Numerical modeling of infragravity
wave response during DELILAH. J. Geophys. Res., 108(C9), 3288, doi:10.1029/2002JC001332.
Verwey, A., Kernkamp, H.W.J., Stelling, G.S., Tse, M.L., Leung, W.C., (2011). Potential and application
of hydrodynamic modelling on unstructured grids. Asian and Pacific Coasts 2011: pp. 1254-1261.
Winter, G., (2012). Rip current characteristics at the Dutch coast: Egmond aan Zee. MSc Thesis, Delft
University of Technology, Institutional Repository, 1204386-000-HYE-0004.
Wu, J., (1982). Wind-stress coefficients over sea surface from breeze to hurricane, J. Geophys. Res.,
87, C12, 9704–9706.
62
MSc Thesis: May 28, 2013
Zijlema, M., Stelling, G.S., Smit, P., (2011). SWASH: An operational public domain code for simulating
wave fields and rapidly varied flows in coastal waters. Coastal Engineering, 58: 992-1012.
Zijlema, M., (2010). Computation of wind–wave spectra in coastal waters with SWAN on unstructured
grids. Coastal Engineering 57, 267-277.
Zijlema, M., (2011). Computational modelling of flow and transport. Reader, Delft University of
Technology.
63
MSc Thesis: May 28, 2013
64
MSc Thesis: May 28, 2013
Appendix A
A.1 Introduction
Deltares has developed Delft3D which is a 3D computational model for rivers, estuaries and coastal
regions. It can simulate hydrodynamic flows, sediment transport, waves, water quality,
morphological developments and ecology. This program has several modules with amongst others:
open channel flows (FLOW), wave-driven hydrodynamics (WAVE) and water quality (WAQ). These
modules cover an integrated 3D modelling of specific flow patterns, morphological developments
and ecological changes.
Eq. A-1
Eq. A-2
Eq. A-3
65
MSc Thesis: May 28, 2013
where u, v and w denote the velocity components in x-, y- and z-direction respectively, the density,
the initial density, p the pressure, the kinematic viscosity and , and represent the
components of the Coriolis force. The equations can be recognized by the momentum balance
equations on the left hand side and the forcing on the right hand side. These can be derived based on
conservation of momentum under the assumption of an incompressible fluid:
Similarly, the continuity equation can be derived under the assumption of conservation of mass:
Eq. A-4
In the 2DH approach, the assumption is made that the vertical accelerations are small enough to be
neglected. Then Eq. A-3 describing the rate of change in the vertical, reduces to the hydrostatic
pressure distribution over the vertical:
Eq. A-5
Substituting Eq. A-5 into Eq. A-1 and Eq. A-2, and elaborating the viscous stresses yields the 2DH
SWEs:
Eq. A-6
( ) Eq. A-7
( ) Eq. A-8
These three equations contain three unknown parameters, namely the water level (ζ) and the
velocities in two directions (u and v). However, the analytical solution of these equations in a more
simplified form can only be calculated in a very limited number of cases, which is very unlikely to
occur in nature. A numerical approximation, by means of e.g. Delft3D-FLOW, is therefore required.
66
MSc Thesis: May 28, 2013
For example the numerical scheme of the two-dimensional heat equation (diffusion equation):
( ) Eq. A-9
( ) Eq. A-10a
( ) Eq. A-10b
Note that in Eq. A-10a the second order derivative in x-direction is evaluated implicitly and the
second order derivative in y-direction is evaluated explicitly. For the next time stage Eq. A-10b) this is
vica versa. The SWEs in Delft3D-FLOW are solved in a similar way. A full description of the discretized
SWEs is given by Broomans (2003).
To assure numerical stability, the Courant number (CFL condition) may not be too large due to the
ADI-effect. The numerical solver turns out to be robust for a CFL upper bound of √ (Deltares,
2011, Delft3D-FLOW). With λ as a dimension of space, the Courant number is defined as:
√ √ Eq. A-11
67
MSc Thesis: May 28, 2013
Legend:
The continuity equation of Eq. A-6 is solved on the water level points and the momentum equations
Eq. A-7 and Eq. A-8 on the velocity points. The choice of the spatial discretization of the advection
terms had great influence on the accuracy, monotony and efficiency of the computational method,
as shown in the river case. The discretisation of the advection term of Eq. A-6 in x-direction at (m,n,k)
for example is (Broomans, 2003):
( ) Eq. A-12
with:
( ) Eq. A-13
Eq. A-12 is solved for both time stages. Note that in the first stage the expression is implicit and it is
explicit in the second stage. For the explicit expression, the velocity on a cell face is partly
determined by upstream information. In case of modelling a Chezy flow this effect may have large
consequences near the upstream boundary, as found in the river case.
In Eq. A-13 the water levels are defined regarding a certain reference level, here MSL or z=0.00m. In
Delft3D, the coordinate system of the vertical is positively upwards for the water level and for the
bottom level downward (see Figure A-3). Eq. A-13 can also be rewritten in a more general form:
Eq. A-14
68
MSc Thesis: May 28, 2013
In 2DH the total water level in a single cell ( ) is determined by the four surrounding depth points
( ). In Delft3D this is defined with the option DPSOPT; taking the minimum-, mean- or
maximum value of the surrounding depth points. In a similar way depth values need to be
determined at the cell interfaces (e.g. ), as the velocities are determined by the discharge
times the interface area. This is done with the option DPUOPT; taking the minimum-, mean- or
upwind value of the depth defined at the water level points ( and ).
Figure A-4: Example of 2DH cells defined on a staggered grid with a sloping bathymetry
69
MSc Thesis: May 28, 2013
The flexibility of modelling with curvilinear grids in Delft3D is limited by the orthogonality
requirement. To reduce the numerical error to a minimum, the grid cells need to satisfy a certain
order of orthogonality in order to comply with the formulations of Eq. A-6 until Eq. A-8. Building
complex models (e.g. deltas or estuaries with different length scales) is challenging and could take
weeks due to these limiting factors.
Wave forces
( ) com-file FLOW
Shear
stresses
water level
com- currents
WAVE file bathymetry
(wind)
70
MSc Thesis: May 28, 2013
The wave forces are by definition related to the radiation stresses under the waves, in according to
Eq. 2-14). This is an important physical mechanism for wave-driven hydrodynamics in the nearshore.
According to Dingemans et al. (1987), the wave forces formulated in terms of gradients of the
radiation stresses could give rise to spurious currents. A more robust formulation is an expression in
terms of wave energy dissipation. Dingemans et al. (1987) found a dissipation formulation with
qualitatively correct current patterns:
Eq. A-15a
Eq. A-15b
where k and represent the wave number and radian frequency respectively. The rate of wave
energy dissipation (D) is computed in SWAN, from the right hand side of Eq. 4-4, and consists of a
contribution due to wave breaking and bottom friction.
Table A-2 describes which parameters are stored in the com-file, with the specific units and the
location on the grid. The abbreviations used for the locations on the grid are clarified in Table A-1
(Deltares, 2011, Delft3D-QUICKPLOT). For a clarification of the staggered grid one is referred to
Figure A-2.
Abbreviation Meaning
D depth points, corners of grid cells
Z water level points, cell centers of grid
UV velocity points
Z (UV) data averaged to water level points (data on file
defined at velocity points)
c defined at cell centers in vertical direction
i defined at cell interfaces in vertical direction
71
MSc Thesis: May 28, 2013
(1) In the recent version of the program, the hydrodynamic grid connects the 3D co-ordinates of the
cell centers of the computation grid points. This may change in a future release to represent the
bounding boxes of the grid cells.
(2) The names and units of these fields are adjusted according to the roughness formulation used in
the simulation. Units: Manning n in s/m1/3, White-Colebrook/Nikuradse k in m, Chezy in m1/2/s, z0 in
m.
72
MSc Thesis: May 28, 2013
(3) New 2D and 3D fields of the same dimension as the grid are automatically detected. They are
assumed to be located at water level points.
The input grid equals the output grid and is directly used in the FLOW module. The computational
grid is used for the SWAN computation. In according to Deltares (2011, Delft3D-WAVE) and SWAN
team (2013), there are no interpolation errors when Input grid = computational grid = output grid.
73
MSc Thesis: May 28, 2013
74
MSc Thesis: May 28, 2013
Appendix B
B.1 Introduction
From a research of Luijendijk et al. (2010) has been concluded that the conventional wave modelling
of Delft3D shows a good performance and is computational efficient. Also for the vertical velocity
distribution with hydrodynamics induced by wave breaking, dissipation due to bottom friction, wind
shear stresses acting on the surface and the tidal forcing. Modelling these nearshore hydrodynamics
are based on the Delft3D FLOW + WAVE coupling. The coupling between both modules is a function
of:
( ⃗ ) Eq. B-1
One desire is to include the modelling of both short-wave (gravity wave) and long wave (infra-gravity
wave) effects. The stationary JONSWAP spectrum showed to be a realistic representation of a young
sea state but is not applicable to swell (Holthuijsen, 2007). Nevertheless, the low frequency waves
introduce the phenomena of surfbeat, as discussed in chapter 2. This is an important mechanism for
sediment transport in the nearshore, as was found by Roelvink and Stive (1989).
Different models already exist which account for this phenomena. A short introduction for three of
these models is given as a study to find a possible integrated wave model for D-Flow FM.
( ) ( ) Eq. B-2
Here represents the short wave energy defined by Eq. 2-16 and represents the dissipation of
wave energy. The model is phase-averaged and contains a phase-randomizer with a uniform
distribution on [ ] (Deltares, 2011, Delft3D-FLOW).
75
MSc Thesis: May 28, 2013
A limitation of this surfbeat model is the restrict application for narrow-banded wave spectra, both in
frequency and directional space. A separate wave model is needed to compute the wave directions
and wave periods of the short waves. This could be the phase-averaged wave model SWAN for
example.
An interesting configuration is by making the input spectrum for SWAN time-dependent. After each
wave computation the action density spectrum can externally be stored, becoming time-dependent
. This can be transferred into a wave-record which is quasi-phase-resolved by imposing a
similar phase-randomizer as for the surfbeat model (Deltares, 2011, Delft3D-FLOW). Only to take
account for the low frequency waves, the phase-randomizer needs to be correlated to the historical
wave record in order to find a smooth overlap of e.g. the infra-gravity waves. This non-stationary
wave record is then stored and expands over time, with the phases of the bound long waves being
locked. An example of a time expanding wave record is shown in the figure below.
Figure B-1: Example a non-stationary wave record with a statistical analysis over the grey region
For the next wave computation, a spectral analysis is done over the last 30minutes of the wave
record in order to make the SWAN computation stationary. In this spectral analysis the wave phases
76
MSc Thesis: May 28, 2013
are again made uniformly distributed. Now the amplitudes and frequencies of the infra-gravity waves
are included for the new wave run.
The inclusion of the wave-induced hydrodynamics with the ‘flow’ induced hydrodynamics is now a
function of:
( ⃗ ) Eq. B-3
where represents the action density spectrum defined as:
Eq. B-4
The time variable in this expression is related to the wave record of Figure B-1. This wave record
varies over space. Translating this into a model would contain a space varying spectral condition in
the form of Eq. B-4. The input file for SWAN would then be a 2D spectrum file (*.sp2) for example.
B.5 Discussion
Although this conceptual wave model takes account for the full frequency domain, it will most likely
be computational intensive. The relevance of developing this methodology in a most efficient way,
outweighed to the conventional wave modelling packages, depends on the model developers. A
presumable extension to modelling infra-gravity waves will be the inclusion of the surfbeat model in
D-Flow FM. This model has been verified in according to Roelvink et al. (1992), and is (also)
successfully integrated in the Delft3D-FLOW module.
Another development within Deltares, for morphological modelling purposes, is the research of
implementing the XBeach formulation directly in D-Flow FM. This formulation uses the action
balance equation in the form of Eq. 4-4. Recent developments of XBeach are curvilinear grid models
(see Roelvink et al., 2012) which showed good morphological predictions to field data.
A plausible configuration of modelling full spectrum waves in D-Flow FM will be a selection of both of
these options. For hydrodynamic purposes SWAN combined with the surfbeat model, and for
morphological purposes the XBeach model.
77
MSc Thesis: May 28, 2013
78
MSc Thesis: May 28, 2013
Appendix C
The original settings of the planar beach case are presented at first. Followed by a sensitivity analysis
to the following parameters:
Directional spreading
Horizontal eddy viscosity
Horizontal eddy diffusivity
Breaker index
Bottom friction formulation SWAN
Bottom friction coefficient SWAN
For the sack of accuracy the 1D model is extended with the Soulsby fitting coefficients of Fredsøe
1984 (Soulsby et al., 1993). The breaker index parameter of the 1D model is set to 0.73. The wave
parameters of the 1D model are kept constant. The results are shown in Figure C-1 through Figure
C-6.
79
MSc Thesis: May 28, 2013
Case 0:
Directional spreading 4
Horizontal eddy viscosity 0
Horizontal eddy diffusivity 0
Soulsby fitting coefficients Fredsøe 1984
Breaker index 0.73
Bottom friction formulation SWAN JONSWAP
Bottom friction coefficient SWAN 0.067
Case 1:
Directional spreading 4
Horizontal eddy viscosity 1
Horizontal eddy diffusivity 10
Soulsby fitting coefficients Fredsøe 1984
Breaker index 0.73
Bottom friction formulation SWAN JONSWAP
Bottom friction coefficient SWAN 0.067
80
MSc Thesis: May 28, 2013
Case 2:
Directional spreading 60
Horizontal eddy viscosity 0
Horizontal eddy diffusivity 0
Soulsby fitting coefficients Fredsøe 1984
Breaker index 0.73
Bottom friction formulation SWAN JONSWAP
Bottom friction coefficient SWAN 0.067
Case 3:
Directional spreading 4
Horizontal eddy viscosity 0
Horizontal eddy diffusivity 0
Soulsby fitting coefficients Fredsøe 1984
Breaker index 0.73
Bottom friction formulation SWAN None
Bottom friction coefficient SWAN -
81
MSc Thesis: May 28, 2013
Case 4:
Directional spreading 4
Horizontal eddy viscosity 0
Horizontal eddy diffusivity 0
Soulsby fitting coefficients Fredsøe 1984
Breaker index 0.60
Bottom friction formulation SWAN JONSWAP
Bottom friction coefficient SWAN 0.067
Case 5:
Directional spreading 10
Horizontal eddy viscosity 0
Horizontal eddy diffusivity 0
Soulsby fitting coefficients Fredsøe 1984
Breaker index 0.60
Bottom friction formulation SWAN JONSWAP
Bottom friction coefficient SWAN 0.08
Figure C-6: Case 5 with a modified directional wave spreading, breaker index and bottom friction coefficient
82
MSc Thesis: May 28, 2013
Tuning the breaker index parameter, results in an offshore shift of wave breaking. This matches the
1D velocity profile closely, but mismatches the wave propagation. Combining all three wave
parameters (directional spreading, bottom friction, and the breaker index) does not make things
better. The shoaling effect is incorrect and the offshore velocities are too high.
The model settings of case 0 are a compromise between realistic absolute velocities and reasonable
wave propagation in the nearshore. An interesting notice from these test cases is the effect of
horizontal eddy viscosity, between case 0 and case 1. In case 1 this effect smoothes the velocity
profile at the offshore region.
83