# Benchmarking of numerical models for wave overtopping at dikes with shallow mildly sloping foreshores: Accuracy versus speed

Christopher H Lashley, Barbara Zanuttigh, Jeremy D Bricker, Jentsje van der Meer, Corrado Altomare, Tomohiro Suzuki, Volker Roeber, Patrick Oosterlo

► To cite this version:

Christopher H Lashley, Barbara Zanuttigh, Jeremy D Bricker, Jentsje van der Meer, Corrado Altomare, et al.. Benchmarking of numerical models for wave overtopping at dikes with shallow mildly sloping foreshores: Accuracy versus speed. Environmental Modelling and Software, 2020, 130, pp.104740. 10.1016/j.envsoft.2020.104740. hal-04946540

HAL Id: hal-04946540

https://hal.science/hal-04946540v1

Submitted on 13 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Check for
updates

# Benchmarking of numerical models for wave overtopping at dikes with shallow mildly sloping foreshores: Accuracy versus speed

Christopher H. Lashley<sup>a, *</sup>, Barbara Zanuttigh<sup>b</sup>, Jeremy D. Bricker<sup>a, c</sup>, Jentsje van der Meer<sup>d, e</sup>, Corrado Altomare<sup>f, g</sup>, Tomohiro Suzuki<sup>a, h</sup>, Volker Roeber<sup>i</sup>, Patrick Oosterlo<sup>a</sup>

## ARTICLE INFO

### Keywords:

Infragravity wave
OpenFOAM
BOSZ
XBeach
SWASH
SWAN

## ABSTRACT

Practitioners often employ diverse, though not always thoroughly validated, numerical models to directly or indirectly estimate wave overtopping (q) at sloping structures. These models, broadly classified as either phase-resolving or phase-averaged, each have strengths and limitations owing to the physical schematization of processes within them. Models which resolve the vertical flow structure or the full wave spectrum (i.e. sea-swell (SS) and infragravity (IG) waves) are considered more accurate, but more computationally demanding than those with approximations. Here, we assess the speed-accuracy trade-off of six well-known models for estimating q, under shallow foreshore conditions. The results demonstrate that: i) q is underestimated by an order of magnitude when IG waves are neglected; ii) using more computationally-demanding models does not guarantee improved accuracy; and iii) with empirical corrections to incorporate IG waves, phase-averaged models like SWAN can perform on par, if not better than, phase-resolving models but with far less computational effort.

# 1. Introduction

## 1.1. Background

Coastal engineers often employ numerical modelling in the design, assessment and rehabilitation of coastal structures to accurately forecast nearshore waves and currents, sometimes including the consequences (Akbar and Aliabadi, 2013; Sierra et al., 2010; Smith et al., 2012; Suzuki et al., 2017). Of particular interest is the extent to which waves reach and pass over the crest of a structure, referred to as wave overtopping. Extreme overtopping events are characterized by considerable flow velocities which impose serious hazards to both people and infrastructure; with flooding or coastal inundation as the most critical consequence.

The integration of numerical modelling in estimating wave over-

topping and the design of coastal structures is becoming increasingly more attractive given the progress in available computing power and the ability of numerical models to incorporate both unconventional structure configurations and extreme environmental conditions. Furthermore, as many of the empirical overtopping models (e.g. EurOtop (2018)) require the incident significant wave height ($H_{m0,T, toe,in}$) and spectral wave period ($T_{m-1,0, toe,in}$) at the toe of the structure as input, numerical models are often needed to accurately capture the nonlinear effects associated with the shoaling and breaking of high-frequency sea-swell (SS) waves in shallow water (Altomare et al., 2016; Mase et al., 2013). Such effects include a rise in mean water level—known as wave-induced setup—and the growth of low-frequency infragravity (IG) waves (Fig. 1) which not only contribute to $H_{m0,T, toe,in}$ but also result in higher values of $T_{m-1,0, toe,in}$ (Hofland et al., 2017).

A variety of numerical models, which may be broadly classified as

<footnote><sup>a</sup> Dept. of Hydraulic Engineering, Delft University of Technology, Stevinweg 1, 2628 CN, Delft, the Netherlands</footnote>

<footnote><sup>b</sup> Dept. Civil, Chemical, Environmental and Materials Engineering, University of Bologna, Via Zamboni, 33, 40126, Bologna, BO, Italy</footnote>

<footnote><sup>c</sup> Dept. of Civil and Environmental Engineering, University of Michigan, 2350 Hayward St., Ann Arbor, MI, 48109, USA</footnote>

<footnote><sup>d</sup> Water Science and Engineering Department, IHE Delft, Westvest 7, 2611 AX, Delft, the Netherlands</footnote>

<footnote><sup>e</sup> Van der Meer Consulting, P.O. Box 11, 8490 AA, Akkrum, the Netherlands</footnote>

<footnote><sup>f</sup> Universitat Politècnica de Catalunya – BarcelonaTech, carrer Jordi Girona 1-3, 08034, Barcelona, Spain</footnote>

<footnote><sup>g</sup> Ghent University, Department of Civil Engineering, Technologiepark 60, 9052, Gent, Belgium</footnote>

<footnote><sup>h</sup> Flanders Hydraulics Research, Berchelei 1m15, 2140, Antwerp, Belgium</footnote>

<footnote><sup>i</sup> E2S chair HPC-Waves Laboratoire SIAME, Univ Pau & Pays Adour, Anglet campus, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: c.h.lashley@tudelft.nl (C.H. Lashley).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Schematic representation of wave transformation over a shallow fore-shore (from an XBeach model simulation), showing a) the growth of IG waves; b) the increase MWL at the dike toe; and c) the shift in the peak in energy density to lower frequencies from offshore (P1) to the dike toe (P2). Vertical line in panel 'c' indicates the separation between SS and IG frequencies.</figcaption></figure>

phase-resolving or phase-averaged, have been developed for such applications; each with strengths and limitations owing to the physical parameterization of processes and the numerical schemes incorporated within them (Cavaleri et al., 2007; Vyzikas and Greaves, 2018). Models which attempt to resolve the vertical flow structure and those that consider the full frequency range of nearshore waves (i.e. both SS and IG waves) are considered not only more accurate, but also more computationally demanding than those which make use of approximations.

Within the phase-resolving class of wave models, those that resolve the vertical flow structure and solve the fully nonlinear, time-averaged Navier-Stokes (NS) equations—often referred to as Computational Fluid Dynamics (CFD) or depth-resolving models—have the least theoretical limitations and are generally considered the most accurate. CFD models, such as the mesh-based Eulerian approach OpenFOAM (Jasak et al., 2007) or mesh-less Lagrangian approach DualSPHysics (Crespo et al., 2015), are able to simulate complex wave problems, such as: nonbreaking and breaking waves, wave-current interaction and wave-structure interaction from deep to shallow water conditions, including the overturning (Lowe et al., 2019) and roller formation of breaking waves (Higuera et al., 2013). However, these models require a significant amount of computational effort (unless a coupling method is applied (Altomare et al., 2015; Altomare et al., 2018; Verbrugghe et al., 2018)); thus, limiting their application so far to very local phenomena—for example, wave overtopping.

As depth-resolved (fully 3D or 2DV) models are generally considered too computationally expensive for operational use, the problem may be further simplified by depth-averaging. These models, in which the vertical structure is not directly resolved but only modelled parametrically, are referred to as two-dimensional in the horizontal (2DH), or 1DH where only a cross-shore transect is simulated (Brocchini and Dodd, 2008). As a result of depth-averaging, processes such as wave overturning, air-entrainment and wave generated turbulence are not directly solved. Those that simulate the amplitude and phase variation of SS waves are often referred to as phase-resolving. Within this type of model, there are generally two main sets of governing equations: i) the Non-linear shallow water (NLSW) equations; and ii) the Boussinesq type.

While the Boussinesq-type models (e.g. FUNWAVE (Kirby et al., 1998), MIKE21-BOUSS (Warren and Bach, 1992) and BOSZ (Roeber and Cheung, 2012)) directly account for the dispersive properties of waves in deeper water, the NLSW models assume that waves are non-dispersive and are therefore limited to shallow-water applications (Brocchini and Dodd, 2008; Zijlema and Stelling, 2008). This limitation can be removed by taking a SS-wave averaged approach; however, at the cost of decreased accuracy (due to exclusion of SS-wave motions). The high-frequency waves are averaged, resulting in only motions at the

scale of the wave group; thus, reducing the computational demand (e.g. XBeach Surfbeat (XB-SB) model (Roelvink and Costas, 2019; Roelvink et al., 2009)).

In order to use the NLSW equations for phase-resolving simulation of SS-wave motions, Stelling and Zijlema (2003) proposed another method to account for dispersion (a result of non-hydrostatic pressure) whereby the pressure is decomposed into non-hydrostatic and hydrostatic pressure components (e.g. SWASH (Zijlema et al., 2011), NHWAVE (Ma et al., 2012) and XBeach Non-hydrostatic (XB-NH) (Smit et al., 2010) numerical models). This approach improves the dispersive properties without neglecting the higher-frequency motions; however, at the expense of more computational demand. The accuracy and range of applicability of the non-hydrostatic models may be further enhanced by coarsely dividing the model domain into a fixed number of vertical layers ($K \le 3$); thereby, improving the frequency dispersion (e.g. SWASH, NHWAVE or XB-NH in multi-layered mode (De Ridder, 2018)). By further increasing the number of vertical layers ($K \ge 10$), models like SWASH may be extended to the depth-resolving class. This approach increases the computational demand but allows processes, such as undertow and the shoreward flow near the surface, to be resolved.

Given that phase-resolving models require a grid resolution high enough to resolve the individual SS-wave components, they are generally computationally feasible only for areas of limited size. For large-scale modelling of wave motion, a phase-averaged approach is most commonly used. This type of model is constructed on the assumption that a random sea-state is composed of a superposition of linear waves whose height is a function of their frequency and direction of propagation. For an individual wave train the rate of change of wave energy (or action) flux is balanced by the wave energy transfer among different wave components in different directions and different frequencies, as well as energy input and dissipation. With the phase information filtered out, these models can use much courser computational grids and therefore be applied to large areas. However, as individual waves are not resolved, these models must be combined with empirical formulae to estimate wave run-up and overtopping (Oosterlo et al., 2018; Sierra et al., 2010). Commonly used spectral models in nearshore applications include SWAN (Booij et al., 1999) and STWAVE (Smith et al., 2001). These models are generally able to accurately reproduce higher harmonics (SS waves); however, they do not account for the interactions that force IG-wave motions (Cavaleri et al., 2007), which tend to dominate in shallow water.

With respect to previous model comparisons in shallow coastal environments, Buckley et al. (2014) assessed the performance of SWASH, SWAN and XB-SB in predicting SS wave heights ($H_{m0,SS}$), IG wave heights ($H_{m0,IG}$) and setup ($\bar{\eta}$) across a steep laboratory fringing reef profile (varying from 1:5 to 1:18.8). Results showed that each model was capable of accurately predicting $H_{m0,SS}$; however, SWAN failed to simulate the transformation of energy to lower frequencies and thus, failed to predict $H_{m0,IG}$. Likewise, SWAN showed considerably more error in its prediction of $\bar{\eta}$ compared to SWASH and XB-SB. On the other hand, XB-SB performed comparably well to its phase-resolving counterparts in the prediction of nearshore wave heights; and surprisingly the extent of wave run-up, particularly when IG-waves dominated at the shoreline (Lashley et al., 2018). From these previous studies, the points of discussion that naturally arise are:

i) Can phase-averaged models like SWAN be accurately applied under very shallow conditions, where IG waves dominate and $\bar{\eta}$ is significant?

ii) Given that IG waves dominate, are models of increasing complexity needed or is a short-wave averaged but IG-wave resolving approach all that is required? and

iii) While attempts at model comparisons for wave overtopping have been made (St-Germain et al., 2014; Vanneste et al., 2014), no study to date has the full range of model complexity (from
depth-resolving to phase-averaged) or successfully quantified the accuracy versus speed of these models under irregular wave forcing.

## 1.2. Objective

In the present study, it is our primary aim to quantify the accuracy versus speed of computation of six commonly-used nearshore wave models (Table 1) in their prediction of irregular wave overtopping at sea dikes with very shallow mildly-sloping foreshores—where IG waves and setup contribute significantly.

## 1.3. Outline

This paper is organized as follows: Section 2 provides overviews of the physical and numerical models applied, followed by descriptions of key parameters and empirical formulae used in the analysis. It ends with a description of the metrics used to quantify model accuracy. In Section 3, the results of the model-data comparisons and the overall influence of IG waves on overtopping are presented and discussed. Section 4 concludes the paper by summarising the findings, acknowledging the present study's limitations and identifying areas for future work.

# 2. Methods

This section begins with a description of the physical model tests under consideration. After which it describes the six numerical models under evaluation, including their governing equations and setup details. A description on the parameters and metrics used to assess model accuracy and computation speed is then provided. Finally, the additional numerical simulations for comparative analysis are described.

## 2.1. Description of the physical models

In the present study, we consider two specific test cases that were both performed at Flanders Hydraulics Research in a smooth, 1-m wide section of their 70-m long and 1.45-m deep wave flume (Altomare et al., 2016) with different deep water wave heights ($H_{m0.T deep}$), peak periods ($T_p$), foreshore slopes ($m$), initial water depths at the toe ($h_{toe}$), dike slopes ($\alpha$) and dike freeboards ($R_c$) (Table 2). These cases were selected to cover a wide range of deep-water wave steepness ($s_0$), from very mild ($s_0 = 0.007$, typical of swell conditions) to very steep ($s_0 = 0.047$, typical of wind-sea conditions). With relative water depths ($h_{toe}/H_{m0.T deep}$) < 1, these conditions are considered very shallow (Hofland et al., 2017). Both experiments simulated irregular spilling waves (with breaker parameter based on $m$, $\xi_{0fore} < 0.5$) with a duration approximately equal to 500 waves to obtain accurate and comparable estimates of the mean overtopping discharge ($q$) (Romano et al., 2015).

<table><caption>**Table 1**

Overview of the numerical models considered for comparative analysis.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Model Type</th><th colspan="2">Wave Propagation</th><th rowspan="2">Overtopping</th></tr><tr><th>SS<br>Waves</th><th>IG<br>Waves</th></tr></thead><tbody><tr><td>OpenFOAM<br>SWASH<sup>a</sup><br>BOSZ<br>XB-NH<br>XB-SB</td><td>Phase-resolving</td><td>Depth-resolving<br>Depth-averaged</td><td>Directly</td><td>Directly</td></tr><tr><td>SWAN</td><td>IG-wave<br>resolving;<br>SS-wave<br>averaged<br>Phase-averaged</td><td>Action-balance</td><td>Directly</td><td>Directly for<br>IG waves<sup>b</sup></td></tr><tr><td></td><td></td><td>Action-balance</td><td>Excluded</td><td>Empirically</td></tr></tbody></table>

For the mild swell-wave case, the variations of water-surface elevations were measured using 10 resistance-type gauges, all synchronously sampling at 50 Hz (Fig. 2a); while 6 gauges with a sample frequency of 20 Hz were used in the steep-wave case (Fig. 2b). In the analysis to follow, the term “offshore” is used to refer to gauges 1 to 7 and 1 to 3 of the mild swell and steep-wind wave cases, respectively; and the term “nearshore” to refer to gauges 8 to 10 and 4 to 6, respectively. In either case, the term “toe” refers to the last wave gauge (gauge 10 and gauge 6 of the mild swell and steep wind-wave cases, respectively).

In both cases, the instantaneous overtopping was measured using two Balluff “Micropulse” water sensors situated inside the overtopping box; and $q$ was then obtained by dividing the total volume of water collected at the end during the test by the total test duration.

## 2.2. Description of numerical models

In this study, six widely-used open-source numerical wave models are considered for comparative analysis. Each model is forced at its boundary with still water levels and parametric spectra (JONSWAP) to match those observed at the most offshore wave gauge during each physical experiment. Likewise, the smooth flume bottom was represented as either a Manning coefficient ($n$) of 0.01 $s/m^{1/3}$ or a Nikuradse geometrical roughness ($k_s$) of $0.3 \times 10^{-3}$ m (in the case of SWAN). A general description of each model is provided in the sections that follow. As we investigate two extremes: very mild swell and very steep wind waves, it is reasonable that some calibration was required for the depth-averaged models (BOSZ, XB-NH and XB-SB). Therefore, a description of the main calibration parameters, their optimum values and impact on model results is also provided. In general, calibration was aimed at reducing the error in $\bar{\eta}$ and $H_{m0.T}$.

### 2.2.1. OpenFOAM

The software OpenFOAM is an Open Source object-oriented library, composed by solvers and utilities (Jasak et al., 2007). The formers are designed to numerically solve continuum mechanics problems, while the latter perform tasks involving data manipulation.

For the present study, the library *waves2Foam*, a toolbox capable of generating and absorbing free surface water waves, has been adopted. Currently, the method applies the relaxation zone technique (active sponge layers) and supports a large range of wave theories (Jacobsen et al., 2012). The governing equations for the combined flow of air and water are given by the Reynolds Averaged Navier-Stokes equations (Equations (1) and (2)):

$$ \frac{\partial \rho u}{\partial t} + \nabla \cdot [\rho u u^T] = - \nabla p^* - g \cdot x \nabla \rho + \nabla \cdot [\mu \nabla u + \rho \tau] + \sigma_T \kappa_T \nabla \gamma, \quad (1) $$

coupled with the continuity equation (2) for incompressible flow:

$$ \nabla \cdot u = 0, \quad (2) $$

where $u$ is the velocity field, $p^*$ is the dynamic pressure component, $\rho$ is the density, $g$ is the acceleration due to gravity and $\mu$ is the dynamic molecular viscosity. The Reynolds stress tensor $\tau$ is defined as:

$$ \tau = \frac{2}{\rho} \mu_t S - \frac{2}{3} kl, \quad (3) $$

where $\mu_t$ is the dynamic eddy viscosity, $S$ is the strain rate tensor, $k$ is the turbulent kinetic energy per unit mass and $I$ is the identity matrix. The last term in Equation (1) is the effect of surface tension, where $\sigma_T$ is the surface tension coefficient and $\kappa_T$ is the surface curvature (Jacobsen et al., 2012). The track of the free surface is performed by using the VOF method (Hirt and Nichols, 1981), where the scalar field ($\gamma$) is 0 for air, 1 for water and any intermediate value is a mixture of the two fluids.

For the mild and the steep cases, two regular slightly different meshes have been generated, to account for the differences between the

<footnote><sup>a</sup> Does not resolve wave overturning or wave roller formation.</footnote>

<footnote><sup>b</sup> Does not include SS-wave overtopping.</footnote>
<table><caption>**Table 2**
Summary of test conditions for both the mild- and steep-wave cases.</caption><thead><tr><th>Case</th><th>H<sub>m0,T,deep</sub> (m)</th><th>T<sub>p</sub> (s)</th><th>cot m</th><th>s<sub>0</sub></th><th>kh</th><th>h<sub>toe</sub> (m)</th><th><sup>h<sub>toe</sub></sup>&frasl;<sub>H<sub>m0,deep</sub></sub></th><th>cot α</th><th>R<sub>c</sub> (m)</th></tr></thead><tbody><tr><td>Mild swell</td><td>0.06</td><td>2.29</td><td>50</td><td>0.007</td><td>0.98</td><td>0.032</td><td>0.53</td><td>2</td><td>0.06</td></tr><tr><td>Steep wind-wave</td><td>0.21</td><td>1.70</td><td>35</td><td>0.047</td><td>1.45</td><td>0.025</td><td>0.12</td><td>3</td><td>0.08</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Physical model setups showing gauge locations for both the: a) mild swell; and b) steep wind-wave cases.</figcaption></figure>

two wave conditions. The numerical domains of the mild and steep cases are respectively composed by 49021 and by 70316 cells, with a graded mesh both in the x (0.3–0.005 m for the mild, 0.1–0.01 m for the steep) and in the y (0.05–0.005 m for the mild, 0.1–0.01 m for the steep) directions. In both cases, the selected regular and constant mesh allowed for a fair compromise between the computational effort and the accuracy of the results.

### 2.2.2. SWASH

SWASH is a time domain model for simulating non-hydrostatic, free-surface and rotational flow. It solves the NLSW equations with an added non-hydrostatic pressure correction term (Smit et al., 2013):

$$ \frac{\partial \eta}{\partial t} + \frac{\partial u h}{\partial x} = 0, \qquad (4) $$

$$ \frac{\partial u}{\partial t} + \frac{\partial u u}{\partial x} + \frac{\partial w u}{\partial z} = - \frac{1}{\rho} \frac{\partial (p_h + p_{nh})}{\partial x} + \frac{\partial \tau_{xx}}{\partial x} + \frac{\partial \tau_{xz}}{\partial z}, \qquad (5) $$

$$ \frac{\partial w}{\partial t} + \frac{\partial u w}{\partial x} + \frac{\partial w w}{\partial z} = - \frac{1}{\rho} \frac{\partial p_{nh}}{\partial z} + \frac{\partial \tau_{zz}}{\partial z} + \frac{\partial \tau_{zx}}{\partial x}, \qquad (6) $$

$$ \frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0, \qquad (7) $$

where η is the free surface elevation; u(x, z, t) and w(x, z, t) are the horizontal and vertical velocities, respectively; h is the water depth; ρ is the density of water; p_h and p_nh are the hydrostatic and non-hydrostatic pressures, respectively; and τ_{xx}, τ_{xz}, τ_{zz} and τ_{zx} are the turbulent stresses.

The model exhibits good linear dispersion up to kh ≈ 8 and kh ≈ 16 with two and three equidistant (sigma) vertical layers (K), respectively; its frequency dispersion is further improved by increasing K.

Here, the model was applied with K = 20, which is sufficient for the phase velocity at the breaking wave front to be computed accurately. As such, no additional control is required to initiate or terminate wave breaking. The vertical pressure gradient was discretized by the standard central differencing scheme with the ILU pre-conditioner. The standard

k-ε turbulence model is applied to take into account vertical mixing.

A cross-shore grid spacing (Δx) of 0.025 m was specified for the mild swell cases, while a coarser resolution (Δx = 0.04 m) was needed for the steep wind-wave case to maintain model stability. This resulted in approximately 330 and 110 grid cells per deep-water wavelength (L₀/Δx) for the mild swell and steep wind-wave cases, respectively. For phase-resolving models, L₀/Δx is typically kept between 50 and 100 (by rule of thumb) to ensure that the wave components are accurately resolved; however, as waves propagate in very shallow water, the local wavelength becomes much shorter than L₀. Thus, in order to maintain a reasonable number of grid cells per local wave length, these higher-than-typical grid resolutions (L₀/Δx = 330 and 110) were specified. It should also be noted that while Hm0,T in the mild swell case showed convergence (~3% change compared to higher resolutions) at Δx = 0.04 m, the resulting q̄ increased by a factor of 6 with Δx = 0.025 m. This suggests that wave run-up and overtopping at the dike were particularly sensitive to the grid resolution.

Rijnsdorp et al. (2017) proposed a sub-grid approach to improve model efficiency, where vertical accelerations and non-hydrostatic pressures are resolved on a relative course grid while the horizontal velocities and turbulent stresses are resolved on a much finer sub-grid. This approach was attempted here, however, the simulations failed due to instabilities.

### 2.2.3. BOSZ

The BOSZ wave model—which is freely-available upon request from the developers—computes hazardous free surface flow problems ranging from near-field tsunamis to extreme swell ranges generated by hurricanes. It solves the following re-formulated, depth-integrated Boussinesq equations of Nwogu (1993), in conservative form:

$$ H_t + (HU)_x + \left[ \left( \frac{\bar{z}_\alpha}{2} - \frac{h^2}{6} \right) \cdot hU_{xx} + \left( \bar{z}_\alpha + \frac{h}{2} \right) \cdot h(hU)_{xx} \right]_x = 0 \qquad (8) $$

$$ (HU)_r + \left( \frac{HU^2}{2} + \frac{1}{2}gH^2 \right)_x - gHh_s + U \left[ \left( \frac{\bar{z}_\alpha}{2} - \frac{h^2}{6} \right) \cdot (hU)_{xx} + \left( \bar{z}_\alpha + \frac{h}{2} \right) \cdot h(hU)_{xx} \right]_x \\ + H \left[ \bar{z}_\alpha \left( \frac{\bar{z}_\alpha}{2} U_{txx} + (hU)_x \right) + \tau \right] = 0, \qquad (9) $$

where H = h + η is the flow depth, U is the horizontal flow velocity defined at a reference depth $\bar{z}_\alpha = -0.55502h$ (Simarro et al., 2013) and τ the bottom shear in terms of the Manning coefficient.

The governing equations exhibit good dispersion accuracy up to kh ≈ π. Given the difficulty of Boussinesq equations in handling flow discontinuities (such as with breaking waves), the model deactivates the dispersion terms during wave breaking and makes use of the underlying NLSW equations where the breaking wave is then approximated as a bore or hydraulic jump. Wave breaking—and the deactivation of the dispersion terms—occurs in the model based on the momentum gradient:

$$ (h + \eta) \frac{\partial u}{\partial x} > B \sqrt{g(h + \eta)}, \qquad (10) $$

where B is a calibration coefficient (by default = 0.5). Here, B = 0.8 produced the best agreement between model and observations for both cases. This suggests that under these particularly shallow conditions, the wave face becomes very steep prior to breaking. For a detailed overview
of the model's sensitivity to this parameter, the reader is pointed to (Roeber et al., 2010). All other model parameters were kept at their default values.

The grid resolution ($L_0/dx$) was set as 200 for the mild swell-wave case but was reduced to 60 for the steep wind-wave case to ensure model stability. For the steep-wave case, higher grid resolutions and lower B values led to instabilities in the form of strong oscillations in surface elevation in the breaking region. This phenomenon, explored extensively by Kazolea and Ricchiuto (2018), is due to the model's hybrid approach to handling wave breaking; that is, where the Boussinesq equations are reduced to the NLSW equations during wave breaking. It should be noted that Boussinesq wave models which take a different (eddy viscosity) approach to wave breaking reportedly show less sensitivity to the grid size (Kazolea and Ricchiuto, 2018); however, this was not evaluated here.

### 2.2.4. XBeach non-hydrostatic

Like SWASH, XB-NH solves the NLSW equations with a non-hydrostatic pressure correction term (Equations (4)–(7)). Here, XBeach version 1.235527 (also known as the “XBeachX” release) is applied in reduced (simplified) two-layer mode, where the non-hydrostatic pressure is assumed constant in the lower (first) layer (De Ridder, 2018). The water depth is divided into two layers with heights $z_1 = \alpha h$ and $z_2 = (1 - \alpha)h$, where $\alpha$ is the layer distribution. The resulting layer-averaged velocities ($u_1$ and $u_2$) are transformed to a depth-averaged velocity (U) and a velocity difference ($\Delta u$). Due to the simplified non-hydrostatic pressure in the lower layer, the vertical velocity between layers is neglected. Therefore, only the continuity relation for the upper (second) layer is required:

$$ \frac{\partial}{\partial x}[(1 + \alpha)hU + (1 - \alpha)h\alpha\Delta u] + 2w_2 - u_2\frac{\partial\eta}{\partial x} - u_1\frac{\partial z_1}{\partial x} = 0, \quad (11) $$

where $w_2$ is the vertical velocity in the upper (second) layer. To determine the water elevation, the global continuity equation is applied:

$$ \frac{\partial \eta}{\partial t} + \frac{\partial h U}{\partial x} = 0, \quad (12) $$

In order to control the computed location and magnitude of depth-limited wave breaking, a hydrostatic front approximation is applied. With this, the pressure distribution under breaking waves is considered hydrostatic when the local surface steepness exceeds a maximum prescribed value ($\lambda = 0.4$, by default):

$$ \frac{\partial \eta}{\partial t} > \lambda, \quad (13) $$

here, $\lambda = 0.9$ and $0.7$ produced the best agreement between the model and observations for the mild- and steep-wave cases, respectively. This further supports the statement that for very shallow foreshores, the waves become particularly steep before breaking. All other model parameters were kept at their default values. Additionally, the grid resolution ($L_0/\Delta x$) was set to ~200 and ~180 for the two respective cases.

### 2.2.5. XBeach surfbeat

XB-SB solves SS-wave motions using the wave-action equation with time-dependent forcing, similar to that of the HISWA model (Holthuijsen et al., 1989). The model represents the SS-wave frequency spectrum by a single frequency ($f_{rep}$) and the wave-action equation is applied at the timescale of the wave group:

$$ \frac{\partial A}{\partial t} + \frac{\partial c_{gx} A}{\partial x} = - \frac{D_w}{\sigma}, \quad (14) $$

$$ A(x, t) = \frac{S_w(x, t)}{\sigma(x, t)} \quad (15) $$

$$ \sigma = \sqrt{gk \tanh kh} \quad (16) $$

where A is the wave action, $S_w$ is the wave energy density, $\sigma$ is the intrinsic wave frequency, k is the wave number, $D_w$ is a dissipation term to account for wave breaking and $c_{gx}$ is the wave-action propagation speed in the cross-shore direction. To simulate wave breaking, XB-SB applies a dissipation model (Roelvink, 1993), by default, for use with SS-wave groups; and a roller model (Nairn et al., 1990; Svendsen, 1984) to represent momentum stored in surface rollers which results in a shoreward delay in wave forcing. The radiation stress gradients that result from these variations in wave action exert forces on the water column and drive IG waves and unsteady currents which are solved by the NLSW equations (Equations (4)–(7)). Therefore, the model directly simulates wave-driven currents and the run-up and overtopping of IG waves.

$$ \bar{D}_w = 2 \frac{\alpha}{T_{rep}} Q_b E_w \frac{H_{rms}}{h}, \quad (17) $$

$$ Q_b = 1 - \exp\left(-\left(\frac{H_{rms}}{H_{max}}\right)^{10}\right) \quad (18) $$

where $\bar{D}_w$ is the total (directionally-integrated) wave energy dissipation due to breaking, $T_{rep} = 1/f_{rep}$ is the representative wave period and $Q_b$ is the fraction of breaking waves; the root-mean-square SS-wave height, $H_{rms} = \sqrt{8E_w/\rho g}$; the maximum wave height, $H_{max} = \gamma_r h$; $E_w$ is the wave-group varying SS-wave energy; $\alpha$ is a dissipation (by default = 1) and $\gamma_r$ is the ratio of breaking waves to local water depth (by default = 0.55 but typically used for calibration).

Here, $\gamma_r = 0.45$ and $0.65$ provided the best agreement between the model and observations for the mild swell and steep wind-wave cases, respectively.

XB-SB does not directly produce the SS-wave component of the energy density spectrum, instead it computes the change in SS-wave energy as a change in the bulk $H_{rms}$ parameter, as described above. In order to produce a complete energy density ($C_{\eta\eta}$) spectrum at each gauge location, a JONSWAP distribution was assumed around the peak-frequency ($f_p$), where $\sqrt{8 \int_{f_p/2}^{2} C_{\eta\eta} df} = H_{rms}$. This SS-wave spectrum (Fig. 3b) was then combined with the IG-wave spectrum (Fig. 3a)—obtained directly from the computed surface elevation—to produce the complete spectrum (Fig. 3c).

For the mild swell-wave case, the grid resolution was varied such that it increased shoreward. This reduced computation time while ensuring that the steep dike slope was accurately captured. As such, $L_0/\Delta x$ varied from ~25 (offshore) to ~160 (at the dike) in the mild-wave case; and from ~45 to ~90 in the steep-wave case.

### 2.2.6. SWAN

SWAN is a third-generation, phase-averaged wave model used to estimate the generation (by wind), propagation and dissipation (by depth-induced breaking and bottom friction) of waves from deep water to the surf zone. This includes wave-wave interactions, in both deep and shallow water, and wave-induced setup; but neglects wave-induced currents and the generation or propagation of IG waves. Like XB-SB, SWAN computes the spectral evolution of A in space and time. This is done in a manner similar to Equation (14); however, unlike XB-SB which makes use of a single representative frequency, SWAN takes the frequency distribution of action density into account. To simulate depth-limited wave breaking, SWAN uses the following parametric dissipation model, by default (Battjes and Janssen, 1978 Battjes and Janssen):

$$ \bar{D}_w = \frac{\alpha}{4} \rho g f_{mean} Q_b H_{max}^2, \quad (19) $$

and $Q_b$ is estimated as:
<figure><img src="image_6.png" /><figcaption>Fig. 3. Example of a) the IG-wave spectrum based on the computed surface elevation; b) an assumed SS-wave spectrum (JONSWAP shape) based on the computed root-mean-square SS-wave height (H<sub>rms</sub>); and c) the total combined spectrum, for XB-SB at the dike toe (steep wind-wave case).</figcaption></figure>

$$ \frac{1 - Q_b}{\ln Q_b} = -8 \frac{E_{tot}}{H_{max}^2}, \qquad (20) $$

where $f_{mean}$ is the mean wave frequency, $H_{max} = \gamma_{bj}h$ and $E_{tot}$ is the total wave-energy variance. Here, $\gamma_{bj} = 0.73$ (default value) provided good agreement between the model and observations for both the mild swell and steep wind-wave cases. For both wave cases, a constant grid spacing of 0.25 m was applied. This corresponded to $L_0/\Delta x \approx 30$ for the mild-wave case and $L_0/\Delta x \approx 20$ for the steep-wave case.

## 2.3. Data processing and analysis

### 2.3.1. Mean water level

The mean water level ($\bar{\eta}$) was calculated by taking the average of the surface elevation, $\eta(t)$, at each gauge location, relative to the elevation of the dike toe. The wave-induced setup, $\langle \eta \rangle$, was then obtained as the difference between $\bar{\eta}$ at each gauge location and $\bar{\eta}$ at the most offshore gauge.

### 2.3.2. Separation of infragravity and sea-swell waves

The time series of $\eta(t)$ were further analysed using the Welch's average periodogram method and a Hann filter with a 50% maximum overlap. The resulting one-dimensional spectra of wave energy density, $C_{\eta\eta}(f)$—with $\sim 30^\circ$ of freedom (Percival and Walden, 1993)—were then used to determine $H_{m0,T}$, $H_{m0,SS}$ and $H_{m0,IG}$, as follows:

$$ H_{m0,T} = 4 \sqrt{\int_{0.005}^{2} C_{\eta\eta} df}, \qquad (21) $$

$$ H_{m0,SS} = 4 \sqrt{\int_{f_p/2}^{2} C_{\eta\eta} df}, \qquad (22) $$

and

$$ H_{m0,IG} = 4 \sqrt{\int_{0.005}^{f_p/2} C_{\eta\eta} df}, \qquad (23) $$

where half the peak frequency ($f_p/2 = 1/2T_p$) is taken as the cut-off to separate SS and IG motions (Roelvink and Stive, 1989). This choice of cut-off frequency is based on the tendency that, in deep water, the majority of SS-wave energy is found at frequencies $> f_p/2$ while the majority of IG-wave energy lies at frequencies $< f_p/2$.

### 2.3.3. Spectral wave period

In addition to wave heights, the spectral wave period ($T_{m-1,0}$) at each gauge location was calculated as follows:

$$ T_{m-1,0} = \frac{m-1}{m_0}, \qquad (24) $$

where,

$$ m_n = \sqrt{\int_{0.005}^{2} C_{\eta\eta} f^n df}. \qquad (25) $$

### 2.3.4. Empirical estimate of the incident infragravity waves

As SWAN neglects the contribution of IG waves to the total wave incident height at the dike toe ($H_{m0,T, toe, in}$) we apply an empirical correction, proposed by (Lashley et al., Forthcoming). Using a dataset of 672 XBeach simulations, an empirical formula for the relative magnitude of the IG waves ($\tilde{H}_{IG}$) was derived with influence factors to account for variations in offshore wave directional spreading ($\tilde{\gamma}_o$), $h_{toe}$ ($\tilde{\gamma}_h$), cot $m$ ($\tilde{\gamma}_f$), vegetation ($\tilde{\gamma}_v$) and cot $\alpha$ ($\tilde{\gamma}_d$):

$$ \tilde{H}_{IG} = 0.36 \cdot H_{m0,T,deep}^{0.5} \cdot \tilde{\gamma}_o \cdot \tilde{\gamma}_h \cdot \tilde{\gamma}_f \cdot \tilde{\gamma}_v \cdot \tilde{\gamma}_d, \qquad (26) $$

For an incident waves analysis (i.e. without the influence of the dike slope) with no directional spreading (1D flume conditions) or vegetation, $\tilde{\gamma}_o$, $\tilde{\gamma}_v$ and $\tilde{\gamma}_d = 1$; while,

$$ \tilde{\gamma}_h = 1.04 \cdot \exp(-1.4 \cdot h_{toe}) + 0.9 \cdot \exp(-0.19 \cdot h_{toe}) \qquad (27) $$

and

$$ \tilde{\gamma}_f = 1.56 - 3.09 \cdot \cot m^{-0.44}. \qquad (28) $$

As $\tilde{H}_{IG}$ represents the ratio of IG to SS waves, $H_{m0,IG, toe, in}$ can be obtained from a SWAN estimate of $H_{m0,SS, toe, in}$:

$$ H_{m0,IG, toe, in} = \tilde{H}_{IG} \cdot H_{m0,SS, toe, in}. \qquad (29) $$

Finally, a corrected estimate of $H_{m0,T, toe, in}$ was obtained as follows:

$$ H_{m0,T, toe, in} = \sqrt{H_{m0,SS, toe, in}^2 + H_{m0,IG, toe, in}^2}. \qquad (30) $$

### 2.3.5. Empirical wave overtopping

While the fully phase-averaged models like SWAN are—to some extent—able to estimate nearshore wave conditions, they cannot directly simulate wave overtopping, as this requires that the individual waves be resolved. In order to estimate wave overtopping, these models can be (and are often) combined with well-established empirical models that require wave parameters at the dike toe as input. In the present study, the EurOtop (2018) formulae based on the work of Van Gent (1999) and Altomare et al. (2016) for (very) shallow foreshores are applied in combination with SWAN. For smooth dikes under perpendicular wave attack with $h_{toe}/H_{m0,T,deep} < 1.5$:

$$ \frac{q}{\sqrt{g \cdot H_{m0,T, toe, in}^3}} = 10^{-0.79} \cdot \exp \left( - \frac{R_c}{H_{m0,T, toe, in} \cdot (0.33 + 0.022 \cdot \xi_{m-1,0})} \right), \qquad (31) $$

with
$$ \xi_{m-1,0} = \frac{\tan \alpha_{sf}}{\sqrt{H_{m0,T, toe.in}/L_{m-1,0}}} \qquad (32) $$

$$ L_{m-1,0} = \frac{g \cdot T_{m-1,0, toe.in}^2}{2\pi}, \qquad (33) $$

$$ \tan \alpha_{sf} = \frac{1.5H_{m0,T, toe.in} + R_{u2\%}}{(1.5H_{m0,T, toe.in} - h_{toe}) \cdot m + (h_{toe} + R_{u2\%}) \cdot \cot \alpha}, \qquad (34) $$

$$ \frac{R_{u2\%}}{H_{m0,T, toe.in}} = 4 - \frac{1.5}{\sqrt{\xi_{m-1,0}}}, \qquad (35) $$

where g is the gravitational constant of acceleration, α is the dike slope, α<sub>sf</sub> is an equivalent slope (to account for waves breaking on the fore-shore) and T<sub>m-1,0, toe.in</sub> is the spectral wave period at the dike toe based on the incident waves (i.e. without the influence of waves reflected at the dike). It should be noted that ξ<sub>m-1,0</sub> and R<sub>u2%</sub> are obtained iteratively (until R<sub>u2%</sub> converges), with a first estimate of R<sub>u2%</sub> = 1.5H<sub>m0,T, toe.in</sub>.

Additionally, as SWAN excludes the contribution of IG waves, corrected estimates of T<sub>m-1,0, toe</sub> are typically obtained using Equations (36) and (37) (Hofland et al., 2017), as outlined in the EurOtop (2018) manual:

$$ \frac{T_{m-1,0, toe.in}}{T_{m-1,0, deep}} = 1 - 6 \cdot \exp(-4\tilde{h}) + \exp(-\tilde{h}), \qquad (36) $$

where,

$$ \tilde{h} = \frac{h_{toe}}{H_{m0,T, deep}} \left( \frac{\cot m}{100} \right)^{0.2}. \qquad (37) $$

### 2.3.6. Error metrics

In order to compare the performance of the numerical models, we assess the mean relative accuracy in an approach similar to that of Lynett et al. (2017):

$$ \text{Mean Ratio}_{\psi} = \frac{1}{n} \sum_{i=1}^{N} \frac{\Psi_{\text{mod}}^{i}}{\Psi_{\text{obs}}^{i}}, \qquad (38) $$

where Ψ is a stand-in for the parameter under consideration (η, H<sub>m0,T</sub>, H<sub>m0,SS</sub>, H<sub>m0,IG</sub> and T<sub>m-1,0</sub>) for the N wave-gauge locations; and subscripts mod and obs refer to model predictions and observations made during the physical experiment, respectively. Here we make a distinction between gauges offshore and nearshore (Fig. 2). A mean ratio of 1 suggests perfect agreement between the model and observations, while values higher or lower than one indicate over- or under-predictions, respectively. It should be noted that all wave gauges (offshore and at the dike toe, see Fig. 2) are considered in Equation (38). While the focus of this study is primarily at the dike toe, it is important to assess the model performance offshore to ensure that: i) the boundary conditions are correctly modelled; and ii) that no (significant) numerical dissipation occurs in deep water, as a result of a coarse grid resolution for example.

Finally, the performance of each model for wave overtopping was assessed by comparing the simple ratio of modelled to observed q and the absolute relative error in the prediction of q:

$$ \text{Absolute Relative Error}_q = \left| 1 - \frac{q_{\text{mod}}}{q_{\text{obs}}} \right|. \qquad (39) $$

### 2.3.7. Computation speed

Two work stations (WS) were used to carry out this research (Table 3). Given the required computational effort, the OpenFOAM simulations were performed on WS-A, while the other models were run on WS-B. To assess computation speed, the duration of each simulation (in wall clock time) was recorded.

<table><caption>**Table 3**

Overview of work stations used to carry out the numerical simulations.</caption><thead><tr><th>Work Station (WS)</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>Operating System</td><td>Ubuntu 14.04 LTS</td><td>Windows 10</td></tr><tr><td>Memory</td><td>31.2 GB</td><td>16 GB</td></tr><tr><td>Processor</td><td>Intel Xeon® CPU ES-2690 v3 @ 2.60 GHz x 16</td><td>Intel® Core™ i7-6600 CPU @ 2.60 GHz, 2.81 GHz x 4</td></tr><tr><td>Graphics</td><td>Gallium 0.4 on NVE7</td><td>Intel® HD Graphics 520</td></tr><tr><td>Type</td><td>64-bit</td><td>64-bit</td></tr><tr><td>Disk</td><td>1.9 TB</td><td>239 GB</td></tr></tbody></table>

## 3. Results and discussions

In this section, the results of the model-data comparisons are presented and discussed. As wave overtopping is the end result of wave propagation, the performance of each model for the prediction of mean water levels, wave heights and periods is first assessed. For the models where calibration was carried out (BOSZ, XB-NH and XB-SB), both default and calibrated results are presented. Note that no parameter tuning was done for the depth-resolving models (OpenFOAM and SWASH) as wave breaking is intrinsically resolved. Likewise, SWAN with default settings showed reasonable agreement and was therefore not calibrated. Lastly, it should be noted that the BOSZ simulation of the steep wind-wave case with default settings resulted in instabilities (see Section 2.2.3) and is therefore not included in the analysis.

### 3.1. Mean water level

Each model, excluding OpenFOAM, is able to accurately (within 15% error) and consistently reproduce η̄ for both the mild swell (Fig. 4a) and steep wind-wave (Fig. 4b) cases. This includes the increase in η̄ near-shore, referred to as wave-induced setup (<η>,) highlighted in Fig. 5 with the XB-NH results representing the general behaviour of the

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Mean ratio of modelled to observed η̄ (markers) for both the a) mild swell and b) steep wind-wave cases, with error bars representing the standard deviation. Solid horizontal lines represent perfect agreement between model and observations. Dashed lines correspond to ± 15% error.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 5. Cross-shore profiles of modelled (XB-NH and OpenFOAM) and observed $<\eta>$ for both the mild swell and steep wind-wave cases.</figcaption></figure>

numerical models. Note that the vertical axes limits of Fig. 4 are set for an easy comparison with the other modelled versus observed plots that follow. While OpenFOAM agrees well with the observations for the mild swell case, it overestimates $<\eta>$ offshore (gauges 2 and 3) and underestimates $<\eta>$ nearshore (gauges 4 to 6) for the steep wind-wave case (Fig. 5). This may be indicative of premature wave breaking in OpenFOAM.

The satisfactory performance of BOSZ and XB-NH observed here (Fig. 4) is in contrast with previous studies (Lashley et al., 2018; Zhang et al., 2019), which found that depth-averaged models were unable to accurately estimate $<\eta>$ due to their lack of vertical resolution and exclusion of wave roller dynamics. However, the difference in model performance here is likely due to the spilling nature of the waves ($\xi_{0,fore} = 0.23$ and $0.13$ for the mild swell and steep wind-wave cases, respectively) compared to the plunging waves and steep fore-reef slopes assessed by Lashley et al. (2018) ($\xi_{0,fore} > 1.1$) and Zhang et al. (2019) ($\xi_{0,fore} = 1.29$). That is, while resolving the vertical structure of flow may be critical for plunging breakers, depth-averaged models are able to perform well under spilling waves.

There is also a notable difference in the observed maximum $<\eta>$ between the mild swell ($<\eta> = 0.004$ m) and steep wind-wave ($<\eta> = 0.015$ m) cases (Fig. 5). This substantial increase in $<\eta>$ as $h_{toe}/H_{m0,T,deep}$ decreases agrees with the findings of Gourlay (1996) on shallow reefs, and suggests that $<\eta>$—which contributes to wave run-up (Stockdon et al., 2006) and, by extension, overtopping—increases proportionally as foreshores become more shallow, or as deep water wave conditions become more energetic.

## 3.2. Significant wave height

SWASH, BOSZ, XB-NH and XB-SB are able to reproduce $H_{m0,T}$, both offshore and nearshore, within 15% error for two cases (Fig. 6). On the other hand, OpenFOAM and SWAN both show notable differences; with SWAN consistently and considerably underestimating $H_{m0,T}$ nearshore.

While SWAN is able to accurately simulate the propagation of high-frequency waves ($H_{m0,SS}$, Fig. 7), it does not compute the low-frequency waves ($H_{m0,IG}$, Fig. 8) and therefore underestimates $H_{m0,T}$ nearshore, where the contribution of IG waves is significant (Fig. 6). The relatively high standard deviation associated with SWAN's nearshore $H_{m0,SS}$ estimates is due to its exclusion of wave reflection. In the physical model, the superposition of the incident and reflected waves results in a nodal/ anti-nodal pattern with a maximum at the dike (outsets in Fig. 9a and Fig. 10a). As SWAN excludes the reflected component, the model underestimates $H_{m0,SS}$ immediately in front of the dike, where the incident and reflected waves add up. On the other hand, this shortcoming makes

<figure><img src="image_9.png" /><figcaption>Fig. 6. Mean ratio of modelled to observed $H_{m0,T}$ (markers) for both the a) mild swell and b) steep wind-wave cases, with error bars representing the standard deviation. Solid horizontal lines represent perfect agreement between model and observations. Dashed lines correspond to ± 15% error.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Mean ratio of modelled to observed $H_{m0,SS}$ (markers) for both the a) mild swell and b) steep wind-wave cases, with error bars representing the standard deviation. Solid horizontal lines represent perfect agreement between model and observations. Dashed lines correspond to ± 15% error.</figcaption></figure>

SWAN especially suitable for use with the empirical overtopping models that require incident-wave conditions as input.

SWAN also predicts a higher and lower maxima in $H_{m0,SS}$ (just before breaking) than XB-NH, for the mild swell (Fig. 9a) and steep wind-wave
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Mean ratio of modelled to observed *H*<sub>m0.IG</sub> (markers) for both the a) mild- and b) steep-wave cases, with error bars representing the standard deviation. Solid horizontal lines represent perfect agreement between model and observations. Dashed lines correspond to ± 15% error.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Cross-shore profiles of modelled and observed: a) *H*<sub>m0.5S</sub> and b) *H*<sub>m0.1G</sub> for the mild-swell case; with c) bed level, for reference. Outset in panel 'a' magnifies the plot area between -1 and 0 m away from the dike.</figcaption></figure>

(Fig. 10a) cases, respectively. This is likely due to the dissipation model employed by SWAN (Equation (19)). Tuning $\gamma_{bj}$—the parameter which controls the maximum wave height to water depth ratio in SWAN—would yield better agreement between the two models; however as there were no wave gauges in this region it is difficult to ascertain which model is correct here.

With respect to OpenFOAM, the model shows inconsistent results

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Cross-shore profiles of modelled and observed: a) *H*<sub>m0.5S</sub> and b) *H*<sub>m0.1G</sub> for the steep wind-wave case; with c) bed level, for reference. Outset in panel 'a' magnifies the plot area between -1 and 0 m away from the dike.</figcaption></figure>

between the two cases. Under the mild swell conditions, OpenFOAM underestimates *H*<sub>m0.T</sub> nearshore (Fig. 6a); however for the steep wind-wave case, the model overestimates *H*<sub>m0.T</sub> nearshore (Fig. 6b). In both cases, the model appears to be too dissipative, resulting in a reduction in *H*<sub>m0.5S</sub> offshore. For the mild swell case, this dissipation is minor resulting in a consistent under-prediction of *H*<sub>m0.5S</sub> (Figs. 7a and 9a) and *H*<sub>m0.1G</sub> (Figs. 8a and 9b). Under the steep wind-wave conditions, however, the dissipation is significant. This observation, combined with the overestimation of <$\eta$> offshore (Fig. 5), indicates premature wave breaking in OpenFOAM. This reduction in *H*<sub>m0.5S</sub> offshore results in unbroken SS-waves reaching the dike and the overestimation of *H*<sub>m0.5S</sub> nearshore (Figs. 7b and 10a). XB-NH also shows some numerical dissipation offshore but this is negligible compared to that of OpenFOAM (Fig. 10a). The wave damping along the flume by OpenFOAM was already found by other authors (Conde, 2019) and is popular among the user community. Previous works (Larsen et al., 2019) and this contribution suggest that a more refined mesh by increasing the number of cells per wave height may improve the wave form modelling during the wave propagation. However, a reduction in grid size for the steep wind-wave case here did not significantly improve the OpenFOAM model results.

Though SWASH is able to accurately predict *H*<sub>m0.1G</sub> it underestimated *H*<sub>m0.5S</sub> nearshore in both cases. This is possibly due to the standard k-ε turbulence model applied in multi-layered mode, which may overestimate the turbulent (vertical) viscosity. Similar to OpenFOAM, a reduction in grid size from 0.04 m to 0.025 m did not significantly improve the estimates (~3% change in *H*<sub>m0.5S</sub>). While calibration generally improved model performance for $\bar{\eta}$ (Fig. 4), *H*<sub>m0.T</sub> (Fig. 6) and *H*<sub>m0.1G</sub> (Fig. 8), it resulted in the overestimation of *H*<sub>m0.5S</sub> nearshore by XB-SB (Fig. 7). This is as a result of tuning $\gamma_r$ (Equation (18)) which affects both the maximum *H*<sub>m0.5S</sub> and *H*<sub>m0.1G</sub>. Perhaps a different approach, where $\alpha$ (Equation (17))—the parameter that controls the magnitude of dissipation—is calibrated (Lashley et al., 2018) would yield better results. However, as XB-SB predicts IG-wave overtopping only, the loss in accuracy for *H*<sub>m0.5S</sub> to improve *H*<sub>m0.1G</sub> predictions was considered acceptable.
### 3.3. Spectral wave period

SWASH, XB-NH and XB-SB show good agreement between modelled and observed $T_{m-1.0}$ predictions; while SWAN, OpenFOAM and BOSZ show notable deviations. As the accurate prediction of $T_{m-1.0}$ requires the models to correctly represent the distribution of wave energy by frequency (Equation (25)), we assess the modelled versus observed wave spectra (Fig. 12 and Fig. 13). SWASH, BOSZ, XB-NH and XB-SB correctly capture the shift in peak energy density ($C_{\eta\eta}$) from the SS-wave (Figs. 12a and 13a) to the IG-wave band (Figs. 12b and 13b); however BOSZ overestimates the magnitude of the IG peak and shows it at slightly lower frequencies than observed. This, coupled with a minor underestimation of the SS-wave energy—most evident for the mild swell case (Fig. 12b)—results in an overestimation of $T_{m-1.0}$.

The consistent underestimation of $T_{m-1.0}$ by SWAN is expected due to its exclusion of $C_{\eta\eta}$ at IG frequencies (Figs. 12b and 13b). OpenFOAM, on the other hand, does show a shift in energy from offshore to the dike toe; however, it shows two distinct IG peaks (Figs. 12b and 13b), not present in the observations. In the mild swell case, this misrepresentation of $C_{\eta\eta}$ at IG frequencies couple with the underestimation of $C_{\eta\eta}$ in the SS-wave band resulted in the significant overestimation of $T_{m-1.0}$ nearshore (Fig. 11a). Under the steep wind-wave conditions, OpenFOAM also shows considerable $C_{\eta\eta}$ in the SS-wave band (nearshore) while the observed spectra shows very little (Fig. 13b). This further supports the argument that due to premature wave decay in the model, some unbroken SS waves are able to reach the dike.

### 3.4. Wave overtopping

#### 3.4.1. General

An important remark is the difference in the observed $q$ between the mild swell (0.094 l/s per m) and steep wind-wave (0.205 l/s per m) cases. Despite having similar $H_{m0,SS, toe}$ values—0.039 m (Figs. 9) and 0.038 m (Fig. 10a) for the mild swell and steep wind-wave cases, respectively—the steep-wave case with a higher $R_c$ and lower $h_{toe}$

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Mean ratio of modelled to observed $T_{m-1.0}$ (markers) for both the a) mild swell and b) steep wind-wave cases, with error bars representing the standard deviation. Solid horizontal lines represent perfect agreement between model and observations. Dashed lines correspond to ± 15% error.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Model-data comparison of wave spectra: a) offshore (at gauge 1) and b) nearshore (at gauge 10) for the mild swell case. Dashed vertical lines indicate the frequency separating SS- and IG-wave motions.</figcaption></figure>

(Table 2), produces double the $q$. This observation suggests that the nonlinear effects of wave breaking—that is, the generation of IG-waves and wave-induced setup—contribute significantly to the resulting overtopping discharge. While the effects of vegetation are not considered here, this observation highlights a potential limitation in studies

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Model-data comparison of wave spectra: a) offshore (at gauge 1) and b) nearshore (at gauge 6) for the steep wind-wave case. Dashed vertical lines indicate the frequency separating SS- and IG-wave motions.</figcaption></figure>
that assess the effectiveness of shallow foreshores but focus only on the attenuation of $H_{m0,SS}$ and neglect the contribution of $H_{m0,IG}$ (Vuik et al., 2016; Yang et al., 2012).

To further investigate the influence of the IG-waves, we compare the overtopping estimated using SWAN and EurOtop (2018) with and without the corrections to $H_{m0,IG, toe, in}$ and $T_{m-1,0, toe, in}$ obtained through Equations (29) and (36), respectively:

Under the mild swell-wave conditions, the ratio $\tilde{H}_{IG} = 0.5$ (Equation (29)) and the contribution of $H_{m0,IG, toe, in}$ to $H_{m0,T, toe, in}$ is minor (Table 4). On the other hand, including the IG waves resulted in a 2.5-fold increase in $T_{m-1,0, toe, in}$ and magnitude 4-fold increase in the predicted $q$, compared to the original SWAN estimates. The difference is more striking for the steep-wave case where the IG waves dominate at the dike toe ($\tilde{H}_{IG} = 0.92$). The inclusion of the IG waves resulted in 36% increase in $H_{m0,T, toe, in}$, a 5-fold increase in $T_{m-1,0, toe, in}$ and an order of magnitude increase in the predicted $q$. Furthermore, the original SWAN estimates—without any corrections to $T_{m-1,0, toe, in}$ and $H_{m0,T, toe, in}$—erroneously show a decrease in $q$ between the mild swell and steep wind-wave cases, while the observations show a notable increase. These results further emphasize the danger of neglecting the IG-wave contribution—demonstrated here by the correction of input to the empirical formulae—in the design and assessment of coastal structures with very shallow foreshores.

Considering the wider model comparison, each model—with the exception of BOSZ—fails to reproduce the overtopping for the mild swell case (Fig. 14). This is particularly evident for XB-SB which significantly underestimates $q$ for both wave cases. This suggests that while XB-SB may estimate wave run-up accurately in IG-wave dominant environments (Lashley et al., 2018), it's exclusion of the SS-wave component considerably limits its performance for wave overtopping.

Compared to the depth-resolving models (XB-NH and BOSZ), both SWASH and OpenFOAM show larger errors in wave overtopping predictions (Fig. 14). While the performance of OpenFOAM may be attributed to its misrepresentation of nearshore waves and periods (Sections 3.2 and 3.3), the performance of SWASH—particularly under mild-swell conditions—is somewhat surprising since it performed reasonably well in the prediction of $\bar{\eta}_i$, $H_{m0}$ and $T_{m-1,0}$ in both cases here and has been previously successful in one-layered mode (Suzuki et al., 2017). However, Suzuki et al. (2017) focused on obtaining good agreement at the toe (the last wave gauge only) and the resulting $\bar{q}$ in their tuning of SWASH; therefore $q$ was not assessed unless wave heights and periods at the toe were within a certain accuracy range, regardless of the input (offshore) conditions. Whereas here, we assess the model's general performance for wave propagation (both offshore and near-shore), in addition to $q$. The models do, however, perform considerably better for the steep-wave case. This is consistent with the findings of Roelvink et al. (2018) and Suzuki et al. (2017) who showed that XB-NH and SWASH, respectively, were more accurate for higher overtopping rates, but suffered for rates below 0.08–0.16 l/s per m (in model scale).

The improvement in SWAN with the corrections is most evident for the steep-wave case, with the estimated $q$ now on par with that of BOSZ and outperforming the other more physically-complex models. Fig. 15 shows the modelled relative overtopping discharge ($q/\sqrt{gH_{m0,T, toe, in}^3}$) versus the relative freeboard ($R_c/(H_{m0,T, toe, in}(0.33 + 0.022\xi_{m-1,0,sf}))$), where $H_{m0,T, toe, in}$ and $T_{m-1,0}$ (to compute $\xi_{m-1,0,sf}$) are taken from Table 4.

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Ratio of modelled to observed $q$ (markers) for both the a) mild swell and b) steep wind-wave cases. Solid horizontal lines represent perfect agreement between model and observations. Dashed lines correspond to a factor of 4 larger and lower than the observations.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Relative overtopping discharge versus relative freeboard with triangles and circles representing the mild swell and steep wind-wave cases, respectively. Solid line corresponds to Equation (31) with dashed lines representing ± 5% exceedance. Note that where applicable $q$ corresponds to the calibrated model results.</figcaption></figure>

<caption>**Table 4**</caption>

SWAN results with and without the empirical corrections for $H_{m0,IG, toe, in}$ (Equation (29)) and $T_{m-1,0, toe, in}$ (Equation (36)).

<table><thead><tr><th rowspan="2">Case</th><th rowspan="2">SWAN</th><th rowspan="2">H<sub>m0,SS, toe, in</sub> (m)</th><th rowspan="2">H<sub>m0,IG, toe, in</sub> (m)</th><th rowspan="2">H<sub>m0,T, toe, in</sub> (m)</th><th rowspan="2">T<sub>m-1,0, toe, in</sub> (s)</th><th colspan="2">q (l/s per m)</th></tr><tr><th>Modelled</th><th>Observed</th></tr></thead><tbody><tr><td rowspan="2">Mild swell</td><td>Original</td><td>0.029</td><td>0</td><td>0.029</td><td>2.17</td><td>0.014</td><td>0.094</td></tr><tr><td>Corrected</td><td>0.029</td><td>0.014</td><td>0.032</td><td>5.38</td><td>0.053</td><td></td></tr><tr><td rowspan="2">Steep wind wave</td><td>Original</td><td>0.033</td><td>0</td><td>0.033</td><td>1.67</td><td>0.003</td><td>0.205</td></tr><tr><td>Corrected</td><td>0.033</td><td>0.03</td><td>0.045</td><td>8.37</td><td>0.089</td><td></td></tr></tbody></table>
If we take the +/-5% exceedance lines of the EurOtop (2018) formula (Equation (31), Fig. 15) as the general range of acceptable overtopping predictions, OpenFOAM (mild swell case), SWASH (steep wind-wave case) BOSZ, XB-NH and SWAN (with corrections) are all reasonable.

As most of the models performed reasonably well for wave propa-
gation, the excellent agreement between BOSZ and the observed *q* is
likely not dependent on underlying governing equations (Boussinesq
versus NLSW) but more to do with how the shoreline and wave run-up
are treated numerically. However, an in-depth analysis of the various
numerical schemes implemented in each numerical model was beyond
the scope of this study.

3.4.2. Accuracy versus speed

In contrast with the general assumption that models of increasing
physical complexity produce more accurate results, Fig. 16 shows no
clear relationship between computational demand (simulation time)
and the absolute relative error in overtopping. Furthermore, the depth-
resolving models (SWASH and OpenFOAM), which have significantly
higher simulation times show larger errors than the depth-averaged
models (XB-NH and BOSZ). The phase-averaged models (XB-SB and
SWAN (original)), despite their considerable speed advantage, signifi-
cantly underestimated the overtopping discharge due to their exclusion
of higher- and lower-frequency wave components, respectively. How-
ever, by including the IG-waves empirically, SWAN's performance
improved significantly; now within acceptable limits and on par with
those of XB-NH and BOSZ but at little to no computational cost (Fig. 16).
It should be noted that the use of SWAN with Equation (36) is already
the recommended approach in EurOtop (2018); the novelty here is the
further improvement in results offered by Equation (29).

4. Conclusion

In the present study we assess the ability of 6 widely used numerical
models to simulate waves overtopping steep dikes with mildly-sloping
shallow foreshores. However, with the exception of OpenFOAM and to
some extent SWASH (K ≥ 10), the above (phase-resolving) models were
originally developed to simulate wave evolution over mildly-sloping
foreshores; and not specifically for wave run-up and overtopping of
steep structure slopes. Since their development, the phase-resolving
models have each been successfully applied to simulate wave propaga-
tion over steep reefs and run-up of relatively steep beaches. Likewise,
depth-resolving models like OpenFOAM and SWASH (K ≥ 10) were
originally developed to simulate wave-structure interaction and not
specifically for wave propagation. In the present study we tested the
ability of these models in both applications: i) wave evolution over a
shallow mildly-sloping foreshore; and ii) the resulting overtopping
discharge.

Overall, BOSZ and XB-NH (under steep wind-waves) showed high
skill in both applications with a reasonable computational demand;
while OpenFOAM, with a much higher computational demand—showed
difficulty in performing both functions. The broad implication of the
present work is that higher-resolution, more computationally-
demanding wave models may simply not be needed; specifically
where the analysis is focused on bulk, time-averaged physical quantities
(Hm0, Tm-1.0 and q), as shown here. Should more detail be required—for
example, estimates of the vertical velocity profile or turbulence—then a
depth-resolving model such as SWASH (K ≥ 10) or OpenFOAM should
be applied. Moreover, SWASH (K ≥ 10) and OpenFOAM are likely to
perform well if the computational domain begins at the dike toe and
ends at the overtopping box; i.e., where simulating wave propagation
over a large distance is not required.

In addition, our results showed that with simple empirical correc-
tions, phase-averaged models like SWAN can perform on par—if not
better than—phase-resolving models, with much less computational
effort. Importantly, our work emphasizes the importance of including IG
waves in the design and assessment of coastal dikes; as neglecting their

<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Accuracy versus speed of the numerical wave models for wave over-
topping with triangles and circles representing the mild swell and steep wind-
wave cases, respectively.</figcaption></figure>

contribution to Hm0,T.toe and Tm-1.0.toe can lead to under-predictions in q
of up to two orders of magnitude.

Given the scope of the model comparison, including both phase-
resolving and phase-averaged, a detailed wave-by-wave comparison of
the higher-resolution models was not carried out. Future work should
address this and investigate the influence of the various numerical
schemes implemented in the respective numerical models, as this was
not within the scope of the present work. Additionally, Equations (31)
and (36) were developed (in part) using the wider dataset from which
these cases were taken; therefore their performance under different
conditions is still to be confirmed. Despite these limitations, the findings
here can aid practitioners in their decision making; specifically in
deciding which numerical model should be applied based on the level of
accuracy required.

Software availability

* OpenFOAM - developed by OpenCFD Ltd, the software package is freely available from: https://www.openfoam.com/

* SWASH - developed at Delft University of Technology, the model is available freely from: http://swash.sourceforge.net/download/download.htm

* BOSZ - developed at the Université de Pau et des Pays de l'Adour, the model is freely available under request from: Volker.roeber@univ-pau.fr

* XBeach Non-hydrostatic and XBeach Surfbeat - developed by IHE Delft, Deltares, Delft University of Technology and the University of Miami, both models are freely available from: https://oss.deltares.nl/web/xbeach/download

* SWAN - developed at Delft University of Technology, the third-generation wave model is freely available from: http://swanmodel.sourceforge.net/download/download.htm

Declaration of competing interest

□ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work is part of the *Perspectief* research programme *All-Risk* with
project number P15-21 which is (partly) financed by NWO Domain Applied and Engineering Sciences, in collaboration with the following private and public partners: the Dutch Ministry of Infrastructure and Water Management (RWS); Deltares; STOWA; the regional water authority, Noorderzijlvest; the regional water authority, Vechtmorgan; It Fryske Gea; HKV consultants; Natuurmonumenten; and waterboard HHNK. Dr. Corrado Altomare acknowledges funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No.: 792370. Volker Roeber acknowledges financial support from the Isite program Energy Environment Solutions (E2S), the Communauté d'Agglomération Pays Basque (CAPB) and the Communauté Région Nouvelle Aquitaine (CRNA) for the chair position HPC-Waves and the European Union's Horizon 2020 Research and Innovation Programme. The authors also acknowledge Akshay Patil for his incite and efforts towards the success of the research.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.envsoft.2020.104740.

## References

Akbar, M., Aliabadi, S., 2013. Hybrid numerical methods to solve shallow water equations for hurricane induced storm surge modeling. Environ. Model. Software 46, 118-128.

Altomare, C., Domínguez, J.M., Crespo, A.J.C., Suzuki, T., Caceres, I., Gómez-Gesteira, M., 2015. Hybridization of the wave propagation model SWASH and the meshfree particle method SPH for real coastal applications. Coast Eng. J. 57 (4), 1550024.

Altomare, C., Suzuki, T., Chen, X., Verwaest, T., Kortenhaus, A., 2016. Wave overtopping of sea dikes with very shallow foreshores. Coast Eng. 116, 236-257.

Altomare, C., Tagliafierro, B., Domínguez, J.M., Suzuki, T., Viccione, G., 2018. Improved relaxation zone method in SPH-based model for coastal engineering applications. Appl. Ocean Res. 81, 15-33.

Battjes, J. A., and Janssen, J. "Energy loss and set-up due to breaking of random waves." Proc., 16th Conference on Coastal Engineering, ASCE, 569-587.

Booj, N., Ris, R.C., Holthuijsen, L.H., 1999. A third-generation wave model for coastal regions - 1. Model description and validation. J Geophys Res-Oceans 104 (C4), 7649-7666.

Brocchini, M., Dodd, N., 2008. Nonlinear shallow water equation modeling for coastal engineering. J Waterw Port C-Asce 134 (2), 104-120.

Buckley, M., Lowe, R., Hansen, J., 2014. Evaluation of nearshore wave models in steep reef environments. Ocean Dyn. 64 (6), 847-862.

Cavaleri, L., Alves, J.H.G.M., Ardhuin, F., Babanin, A., Banner, M., Belibassakis, K., Benoit, M., Donelan, M., Groeneweg, J., Herbers, T.H.C., Hwang, P., Janssen, P.A.E., Janssen, T., Lavrenov, I.V., Magne, R., Monbaliu, J., Onorato, M., Polnikov, V., Resio, D., Rogers, W.E., Sheremet, A., McKee Smith, J., Tolman, H.L., van Vledder, G., Wolf, J., Young, I., 2007. Wave modelling – the state of the art. Prog. Oceanogr. 75 (4), 603-674.

Conde, J.M.P., 2019. Comparison of different methods for generation and absorption of water waves. Revista de Engenharia Térmica 18 (1), 71-77.

Crespo, A.J.C., Domínguez, J.M., Rogers, B.D., Gómez-Gesteira, M., Longshaw, S., Canelas, R., Vacondio, R., Barreiro, A., García-Feal, O., 2015. DualSPHysics: open-source parallel CFD solver based on smoothed particle hydrodynamics (SPH). Comput. Phys. Commun. 187, 204-216.

De Ridder, M., 2018. Non-hydrostatic Wave Modelling of Coral Reefs with the Addition of a Porous In-canopy Model. MSc, Delft University of Technology. Retrieved from. http://resolver.tudelft.nl/uuid:4f738968-be52-4484-af1d-60a0e74.

EurOtop, 2018. In: Van der Meer, J., Allsop, N., Bruce, T., De Rouck, J., Kortenhaus, A., Pullen, T., Schuttrumpf, H., Troch, P., Zanuttigh, B. (Eds.), Manual on Wave Overtopping of Sea Defences and Related Structures. An Overtopping Manual Largely Based on European Research, but for Worldwide Application. Retrieved from. www.overtopping-manual.com.

Gourlay, M.R., 1996. Wave set-up on coral reefs. 1. Set-up and wave-generated flow on an idealised two dimensional horizontal reef. Coast Eng. 27 (3), 161-193.

Higuera, P., Lara, J.L., Losada, I.J., 2013. Simulating coastal engineering processes with OpenFOAM®. Coast Eng. 71, 119-134.

Hirt, C.W., Nichols, B.D., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39 (1), 201-225.

Hofland, B., Chen, X., Altomare, C., Oosterlo, P., 2017. Prediction formula for the spectral wave period Tm-1,0 on mildly sloping shallow foreshores. Coast Eng. 123 (Suppl. C), 21-28.

Holthuijsen, L.H., Booj, N., Herbers, T.H.C., 1989. A prediction model for stationary, short-crested waves in shallow water with ambient currents. Coast Eng. 13 (1), 23-54.

Jacobsen, N.G., Fuhrman, D.R., Fredsøe, J., 2012. A wave generation toolbox for the open-source CFD library: OpenFoam®. Int. J. Numer. Methods Fluid. 70 (9), 1073-1088.

Jasak, H., Jemcov, A., Tukovic, Z., 2007. OpenFOAM: a C++ library for complex physics simulations. In: Proc., International Workshop on Coupled Methods in Numerical Dynamics. IUC Dubrovnik, Croatia, pp. 1-20.

Kazolea, M., Ricciutto, M., 2018. On wave breaking for Boussinesq-type models. Ocean Model. 123, 16-39.

Kirby, J.T., Wei, G., Chen, Q., Kennedy, A.B., Dalrymple, R.A., 1998. Funwave 1.0: fully nonlinear Boussinesq wave model-Documentation and user's manual research report NO. CACR-98-06.

Larsen, B.E., Fuhrman, D.R., Roenby, J., 2019. Performance of interFoam on the simulation of progressive waves. Coast Eng. J. 61 (3), 380-400.

Lashley, C. H., Bricker, J. D., Van der Meer, J., Altomare, C., and Suzuki, T. (Forthcoming). "Relative magnitude of infragraft waves at coastal dikes with shallow foreshores: a prediction tool." J. Waterw. Port, Coast. Ocean Eng.

Lashley, C.H., Roelvink, D., van Dongeren, A., Buckley, M.L., Lowe, R.J., 2018. Nonhydrostatic and surfbeat model predictions of extreme wave run-up in fringing reef environments. Coast Eng. 137, 11-27.

Lowe, R.J., Buckley, M.L., Altomare, C., Rijnsdorp, D.P., Yao, Y., Suzuki, T., Bricker, J.D., 2019. Numerical simulations of surf zone wave dynamics using Smoothed Particle Hydrodynamics. Ocean Model. 144, 101481.

Lynett, P.J., Gately, K., Wilson, R., Montoya, L., Arcas, D., Aytore, B., Bai, Y., Bricker, J. D., Castro, M.J., Cheung, K.F., David, C.G., Dogan, G.G., Escalante, C., González-Vida, J.M., Grilli, S.T., Heitmann, T.W., Horrillo, J., Kánoğlu, U., Kian, R., Kirby, J. T., Li, W., Macías, J., Nicolsky, D.J., Ortega, S., Pampell-Manis, A., Park, Y.S., Roeber, V., Sharghivand, N., Shelby, M., Shi, F., Tehranirad, B., Tolkova, E., Thio, H. K., Velioğlu, D., Yalcqner, A.C., Yamazaki, Y., Zaytsev, A., Zhang, Y.J., 2017. Inter-model analysis of tsunami-induced coastal currents. Ocean Model. 114, 14-32.

Ma, G., Shi, F., Kirby, J.T., 2012. Shock-capturing non-hydrostatic model for fully dispersive surface wave processes. Ocean Model. 43-44, 22-35.

Mase, H., Tamada, T., Yasuda, T., Hedges, T.S., Reis, M.T., 2013. Wave runup and overtopping at seawalls built on land and in very shallow water. J. Waterw. Port, Coast. Ocean Eng. 139 (5), 346-357.

Nairn, R.B., Roelvink, J., Southgate, H.N., 1990. Transition zone width and implications for modelling surfzone hydrodynamics. In: Proc., 22nd International Conference on Coastal Engineering. ASCE, pp. 68-81.

Nwogu, O., 1993. Alternative form of Boussinesq equations for nearshore wave-propagation. J Waterw Port C-Asce 119 (6), 618-638.

Oosterlo, P., McCall, R., Vuik, V., Hofland, B., van der Meer, J., Jonkman, S., 2018. Probabilistic assessment of overtopping of sea dikes with foreshores including infragraft waves and morphological changes: westkapelle case study. J. Mar. Sci. Eng. 6 (2), 48.

Percival, Donald, Walden, Andrew, 1993. Spectral analysis for physical applications. Cambridge University Press.

Rijnsdorp, D.P., Smit, P.B., Zijlema, M., Reniers, A.J.H.M., 2017. Efficient non-hydrostatic modelling of 3D wave-induced currents using a subgrid approach. Ocean Model. 116, 118-133.

Roeber, V., Cheung, K., 2012. BOSZ (Boussinesq ocean and surf zone model). In: Proc., 2011 NTHMP Model Benchmarking Workshop. NOAA.

Roeber, V., Cheung, K.F., Kobayashi, M.H., 2010. Shock-capturing Boussinesq-type model for nearshore wave processes. Coast Eng. 57 (4), 407-423.

Roelvink, D., Costas, S., 2019. Coupling nearshore and aeolian processes: XBeach and duna process-based models. Environ. Model. Software 115, 98-112.

Roelvink, D., McCall, R., Mehvar, S., Nederhoff, K., Dastgheib, A., 2018. Improving predictions of swash dynamics in XBeach: the role of groupiness and incident-band runup. Coast Eng. 134, 103-123.

Roelvink, D., Reniers, A., van Dongeren, A., van Thiel de Vries, J., McCall, R., Lescinski, J., 2009. Modelling storm impacts on beaches, dunes and barrier islands. Coast Eng. 56 (11-12), 1133-1152.

Roelvink, J.A., 1993. Dissipation in random wave groups incident on a beach. Coast Eng. 19 (1), 127-150.

Roelvink, J.A., Stive, M.J.F., 1989. Bar-generating cross-shore flow mechanisms on a beach. J. Geophys. Res.: Oceans 94 (C4), 4785-4800, 1978-1982.

Romano, A., Bellotti, G., Briganti, R., Franco, L., 2015. Uncertainties in the physical modelling of the wave overtopping over a rubble mound breakwater: the role of the seeding number and of the test duration. Coast Eng. 103, 15-21.

Sierra, J.P., González-Marco, D., Mestres, M., Gironella, X., Oliveira, T.C.A., Cáceres, I., Mósco, C., 2010. Numerical model for wave overtopping and transmission through permeable coastal structures. Environ. Model. Software 25 (12), 1897-1904.

Simarro, G., Orfila, A., Galan, A., 2013. Linear shoaling in Boussinesq-type wave propagation models. Coast Eng. 80, 100-106.

Smit, P., Stelling, G., Roelvink, J., Van Thiel de Vries, J., McCall, R., Van Dongeren, A., Zwinkels, C., Jacobs, R., 2010. XBeach: Non-hydrostatic Model: Validation, Verification and Model Description. Delft Univ. Technol.

Smit, P., Zijlema, M., Stelling, G., 2013. Depth-induced wave breaking in a non-hydrostatic, near-shore wave model. Coast Eng. 76, 1-16.

Smith, J.M., Sherlock, A.R., Resio, D.T., 2001. STWAVE: Steady-State Spectral Wave Model User's Manual for STWAVE, Version 3.0. ENGINEER RESEARCH AND DEVELOPMENT CENTER VICKSBURG MS COASTAL AND HYDRAULICSLAB.

Smith, R.A.E., Bates, P.D., Hayes, C., 2012. Evaluation of a coastal flood inundation model using hard and soft data. Environ. Model. Software 30, 35-46.

St-Germain, P., Nistor, I., Readshaw, J., Lamont, G., 2014. Numerical modeling of coastal dike overtopping using sph and non-hydrostatic nlsw equations. Coast Eng. Proc. 1 (34).
Stelling, G., Zijlema, M., 2003. An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation. Int. J. Numer. Methods Fluid. 43 (1), 1–23.

Stockdon, H.F., Holman, R.A., Howd, P.A., Sallenger, A.H., 2006. Empirical parameterization of setup, swash, and runup. Coast Eng. 53 (7), 573–588.

Suzuki, T., Altomare, C., Veale, W., Verwaest, T., Trouw, K., Troch, P., Zijlema, M., 2017. Efficient and robust wave overtopping estimation for impermeable coastal structures in shallow foreshores using SWASH. Coast Eng. 122 (Suppl. C), 108–123.

Svendsen, I.A., 1984. Wave heights and set-up in a surf zone. Coast Eng. 8 (4), 303–329.

Van Gent, M., 1999. Physical Model Investigations on Coastal Structures with Shallow Foreshores: 2D Model Tests with Single and Double-Peaked Wave Energy Spectra. Hydraulic Engineering Reports, Delft.

Vanneste, D.F.A., Altomare, C., Suzuki, T., Troch, P., Verwaest, T., 2014. Comparison of numerical models for wave overtopping and impact on a sea wall. Coast Eng. Proc. 1 (34).

Verbrugghe, T., Domínguez, J.M., Crespo, A.J.C., Altomare, C., Stratigaki, V., Troch, P., Kortenhaus, A., 2018. Coupling methodology for smoothed particle hydrodynamics modelling of non-linear wave-structure interactions. Coast Eng. 138, 184–198.

Vuik, V., Jonkman, S.N., Borsje, B.W., Suzuki, T., 2016. Nature-based flood protection: the efficiency of vegetated foreshores for reducing wave loads on coastal dikes. Coast Eng. 116 (Suppl. C), 42–56.

Vyzikas, T., Greaves, D., 2018. Numerical Modelling. Wave and Tidal Energy, John Wiley & Sons, Ltd, pp. 289–363.

Warren, I.R., Bach, H.K., 1992. Mike 21: a modelling system for estuaries, coastal waters and seas. Environ. Software 7 (4), 229–240.

Yang, S.L., Shi, B.W., Bouma, T.J., Ysebaert, T., Luo, X.X., 2012. Wave attenuation at a salt marsh margin: a case study of an exposed coast on the yangtze estuary. Estuar. Coast 35 (1), 169–182.

Zhang, S.J., Zhu, L.S., Zou, K., 2019. A comparative study of numerical models for wave propagation and setup on steep coral reefs. China Ocean Eng. 33 (4), 424–435.

Zijlema, M., Stelling, G., Smit, P., 2011. SWASH: an operational public domain code for simulating wave fields and rapidly varied flows in coastal waters. Coast Eng. 58 (10), 992–1012.

Zijlema, M., Stelling, G.S., 2008. Efficient computation of surf zone waves using the nonlinear shallow water equations with non-hydrostatic pressure. Coast Eng. 55 (10), 780–790.