# Numerical modelling of nearshore wave transformation, breaking and overtopping of coastal protections with the enhanced Serre-Green–Naghdi equations

Guillaume Coulaud, Maria Teles, Michel Benoit

► To cite this version:

Guillaume Coulaud, Maria Teles, Michel Benoit. Numerical modelling of nearshore wave transformation, breaking and overtopping of coastal protections with the enhanced Serre-Green-Naghdi equations. Coastal Engineering, 2025, 10.1016/j.coastaleng.2025.104857. http://dx.doi.org/10.1016/j.coastaleng.2025.104857

HAL Id: hal-05238645

https://edf.hal.science/hal-05238645v1

Submitted on 3 Sep 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>

# Numerical modelling of nearshore wave transformation, breaking and overtopping of coastal protections with the enhanced Serre-Green-Naghdi equations

Guillaume Coulaud<sup>a,b</sup>,<sup>DID</sup>*, Maria Teles<sup>a</sup>,<sup>DID</sup>, Michel Benoit<sup>a,b</sup><sup>DID</sup>

## ARTICLE INFO

**Keywords:**

Wave overtopping
Serre-Green-Naghdi
Wave breaking
Nonlinear waves

## ABSTRACT

Admissible average overtopping discharges in given storm conditions are typically used to design coastal protections, in particular dykes or breakwaters. These discharges are usually estimated using semi-empirical formulas relying on wave conditions at the toe of the structure. These formulas, unfortunately, only work for simple configurations, invariant alongshore, and can be insufficient for complex sea states. Therefore, numerical modelling could be a more flexible alternative for estimating these discharges. This work presents the development and validation of a Boussinesq-type numerical model solving the fully-nonlinear weakly-dispersive enhanced Serre-Green-Naghdi equations for the simulation of random wave overtopping over impermeable structures in one horizontal dimension. Wave breaking is modelled with an eddy viscosity approach based on the turbulent kinetic energy, which is robust and accurate at describing energy dissipation in the surf zone. Two distinct experimental datasets, with 184 trials in total and very dissimilar wave conditions and foreshore seabed profiles, are used to validate the model regarding wave propagation, shoaling, breaking and overtopping. Both unimodal and bimodal sea states are considered. Average overtopping discharges in configurations with deep and very shallow foreshores, as well as for breaking and non-breaking waves, are well reproduced by the model. For instance, typical mean relative errors on the simulated mean overtopping rates are found to lie within ±20% compared with the measurements, at least for the largest discharges of the considered campaigns. The scatter of simulated discharges is somewhat higher for lower discharges, but the results remain in an acceptable range.

## 1. Introduction

The design and safety assessment of coastal protections is typically done by specifying a maximum admissible average overtopping rate for given storm conditions. The discharges are usually estimated from semi-empirical relations (TAW, 2002; EurOtop, 2007; Goda, 2009; Mase et al., 2013; EurOtop, 2018; Lashley et al., 2021), determined using a large number of laboratory experiments. They relate the over- topping rate over a given structure to the incident wave properties (typically the significant wave height $H_{1/3}$ or $H_{m0}$, the peak or mean wave period and the mean incoming direction) in deep water or at its toe, and some characteristics of the structure, including crest height, cross-sectional profile, porosity of the armour, filter layers and core, etc. Among the most widely used formulas are those from the EurOtop (2018) guide, which are based on wave conditions at the toe of the structures. However, the applicability of all formulas is limited to sim- ple configurations, invariant along the shore. Also, if no measurements

are available, then the wave conditions at the toe of the structure must be estimated somehow. When propagating in varying water depth, incident waves undergo significant transformations from deep waters to the structure's toe, due to shoaling, refraction, diffraction, nonlinear interactions and wave breaking. These processes must be carefully accounted for. Physical (scale) models allow for representing most of these processes accurately, although considering large-scale domains may sometimes not be feasible or alternatively require a large scaling factor (then introducing some spurious effects in the laboratory model). Another complementary, flexible and accurate way to do so is to use numerical models. If a phase-resolving wave model is used, pro- vided it can model moving shorelines and run-up/run-down, average overtopping rates can also be directly estimated.

Among phase-resolving models, those based on the Nonlinear Shal- low Water Equations (NSWE) are the most commonly used due to their

<footnote><sup>a</sup> EDF R&D, Laboratoire National d'Hydraulique et Environnement (LNHE), 6 quai Watier, Chatou, 78400, France</footnote>

<footnote><sup>b</sup> LHSV, Saint-Venant Hydraulics Laboratory (ENPC, Institut Polytechnique de Paris, EDF R&D), 6 quai Watier, Chatou, 78400, France</footnote>

<footnote>* Corresponding author at: EDF R&D, Laboratoire National d'Hydraulique et Environnement (LNHE), 6 quai Watier, Chatou, 78400, France. E-mail addresses: guillaume.coulaud@edf.fr (G. Coulaud), maria.teles@edf.fr (M. Teles), michel.benoit@edf.fr (M. Benoit).</footnote>
simplicity and the availability of accurate and robust numerical methods suited for modelling moving shorelines. They are depth-averaged equations, appropriate to study long wave run-up, overtopping and flooding, see e.g. Kobayashi and Wurjanto (1989), Van Gent (1994), Dodd (1998), Hu et al. (2000), Hubbard and Dodd (2002), Shiach et al. (2004), Briganti and Dodd (2009a), Tuan and Oumeraci (2010). The equations, however, do not model frequency dispersion and are thus only valid in shallow waters and for very long waves (wavelength greater than approximately 5–6 times the water depth Pullen and Allsop, 2003). More complex wave models are required to study wave transformation up to the toe of the structure.

Weakly dispersive long wave models are extensions of the NSWE, retaining the dispersive effects of lowest order(s) by including terms with high-order derivatives of the flow variables, while still being depth-averaged. These models, denoted here as Boussinesq-type (BT) models, include various forms of the Boussinesq (Peregrine, 1967; Madsen and Sørensen, 1992; Nwogu, 1993) or Serre-Green-Naghdi (SGN) (Serre, 1953; Green and Naghdi, 1976) equations. They can be regarded as weakly to mildly dispersive, and weakly to fully nonlinear depending upon the terms neglected in the course of their mathematical derivation. Boussinesq equations are usually weakly nonlinear, i.e. they include only linear dispersive terms, while the SGN equations are fully nonlinear, i.e. they include all nonlinear dispersive terms at the considered order in dispersion. They can model wave propagation in intermediate and shallow waters, and can be applied to study wave run-up and overtopping (Stansby, 2003; Lynett et al., 2010; Tonelli and Petti, 2010; Tsung et al., 2012; McCabe et al., 2013; Whittaker et al., 2018; Pinault et al., 2020; Kalisch et al., 2024). Comprehensive reviews of BT models are available in Kirby (2003), Brocchini (2013), Kazolea and Ricchiuto (2024) and Coulaud et al. (2025).

A possible alternative to BT models is non-hydrostatic wave models, which compute the non-hydrostatic component of the pressure to include dispersion. SWASH (Zijlema et al., 2011) and XBeach Non-Hydrostatic (Smit et al., 2010) are well-known examples of this class. The former is increasingly applied to model wave overtopping (Suzuki et al., 2014, 2017, 2020; Zhang et al., 2020; Henderson et al., 2022; Tuozzo et al., 2024).

Finally, the most complex phase-resolving models solve the full 3D Euler or Navier-Stokes equations. They can theoretically provide the most detailed and accurate information on wave propagation and overtopping characteristics, at a high computational cost, though. Performing, in a reasonable time, computations over sufficient durations to obtain reliable averaged quantities can thus only be done in 2DV cross-shore sections. Examples of wave overtopping computations with such models can be found in Losada et al. (2008), Higuera et al. (2014), Akbari (2017), Altomare et al. (2021), Chen et al. (2021) and Irías Mata and Van Gent (2023).

Lashley et al. (2020) compared the performance of six different numerical models for simulating random wave propagation and overtopping. They considered two cases of long-crested random waves overtopping a dyke with smooth slope, preceded by a foreshore with a mild slope, inducing significant wave transformation and wave breaking. The models ranged from the phase-averaged model SWAN (Booij et al., 1999) to the Navier-Stokes solver OpenFOAM (Jasak et al., 2007). Their results show that the BT model BOSZ (Roeber and Cheung, 2012) provides the best results and an almost perfect agreement on the average discharge, while having a moderate computational cost. The non-hydrostatic models XBeach Non-Hydrostatic and SWASH, which have similar computational costs, provide accurate results as well. This indicates that BT models are a good compromise between accuracy and efficiency for estimating overtopping rates.

As already mentioned, various such BT models exist, with different dispersive and nonlinear properties and showing different results depending on the situation. In particular, fully nonlinear models include nonlinear dispersive effects and theoretically model wave shoaling more accurately up to the breaking point. This paper presents the

development of a one-dimensional numerical model, based on the enhanced SGN (eSGN) equations of Bonneton et al. (2011), and assesses its performance for the simulation of random waves overtopping smooth structures. The selection of this particular model is based on the extensive comparison done by Coulaud et al. (2025) who studied 8 variants of BT models, both from a theoretical perspective (e.g. linear dispersion relation, second-order coupling coefficient, linear shoaling gradient, etc.) and through numerical simulations of three challenging test cases involving strong depth variations. They concluded that this weakly dispersive, fully nonlinear model accurately describes coastal wave propagation while being robust and stable. One crucial aspect of wave modelling with BT models is the treatment of wave breaking, as they do not model turbulence or energy dissipation. For our application, the robustness of the breaking model is paramount, as it needs to be able to deal with random wave breaking, violent uprush and backwash, and overtopping. Therefore, it was chosen to work with an eddy viscosity approach based on the turbulence kinetic energy (TKE), which has been proven accurate and robust in different situations, including similar applications of irregular wave run-up (Kalisch et al., 2024).

The remainder of this article is divided as follows. Section 2 introduces the chosen mathematical and wave breaking models. In Section 3, the numerical discretisation of the model is presented in great detail. The results obtained when simulating 184 cases belonging to two experimental datasets of long-crested irregular wave overtopping are reported in Section 4 and Section 5. Finally, Section 6 discusses the results and provides conclusions and perspectives.

## 2. Mathematical model

### 2.1. The enhanced Serre-Green-Naghdi equations

The numerical model developed and used to simulate wave propagation and overtopping is based on the one-dimensional enhanced Serre-Green-Naghdi (eSGN) equations, derived by Bonneton et al. (2011). This set of equations is equivalent at $O((kd)^2)$ to the original SGN equations from Serre (1953) and Green and Naghdi (1976), with $k$ and $d$ the representative wavenumber and water depth, respectively. The eSGN model has better dispersive properties than the original equations, making it applicable up to the deep water limit, i.e. $kd \approx \pi$, while keeping the relative error on the linear phase speed below 1%, with respect to that predicted by the linear wave theory. It has been widely used to study coastal wave propagation, showing satisfactory results in a wide range of configurations, both in one and two horizontal dimensions, see e.g. Tissier et al. (2012), Popinet (2015), Filippini et al. (2016), Kazolea and Ricchiuto (2018), Chassagne et al. (2019), Li et al. (2019), Kazolea and Ricchiuto (2024) and Coulaud et al. (2025).

The equations are solved under the following form, proposed by Kazolea et al. (2023), separating the hyperbolic part of the model (i.e. the NSWE) from the dispersive part:

$$h_t + (hu)_x = 0 \qquad (1)$$

$$q_t + (hu^2 + \frac{1}{2}gh^2)_x = ghd_x + S_f + \Phi \qquad (2)$$

$$\Phi = \Psi + \frac{gh}{\alpha}\eta_x \qquad (3)$$

$$(I + \alpha T)\Psi = -\frac{gh}{\alpha}\eta_x - hQ + D_{wb} \qquad (4)$$

Eqs. (1) and (2) express conservation of mass and momentum, respectively. $d(x)$ is the still water depth, $\eta(x, t)$ the free surface elevation, $h(x, t) = d + \eta$ the total water depth, $u(x, t)$ the depth-averaged horizontal velocity, $g$ the acceleration of gravity and $q(x, t) = hu$ is the momentum and horizontal volumic discharge. Subscripts $t$ and $x$ denote partial derivatives with respect to time and space, respectively. $\Phi(x, t)$, defined by Eq. (3), is a source term accounting for all dispersive effects, obtained from the resolution of the elliptic Eq. (4). Setting $\Phi = 0$ in
Eqs. (1) and (2), one recovers the hyperbolic NSWE. $S_f$ is a source term modelling bottom friction. Manning's law is used here:

$$S_f = -\frac{gn^2}{h^{1/3}}|u|u| \quad (5)$$

The expressions for the linear operator $T(x, t)$ and the term $Q(x, t)$ in Eq. (4) are:

$$T(\Psi) = \left[ \frac{1}{3} (h_x^2 + h h_{xx}) + d_x^2 - h_x d_x - \frac{1}{2} h d_{xx} \right] \Psi - \frac{1}{3} h h_x \Psi_x - \frac{1}{3} h^2 \Psi_{xx} \quad (6)$$

$$Q = 2hh_x u_x^2 + \frac{4}{3}h^2 u_x u_{xx} - hd_x u_x^2 - hd_{xx} u_x u_x - \left[ h_x d_{xx} + \frac{1}{2}hd_{xxx} - d_x d_{xx} \right] u^2 \quad (7)$$

$D_{wb}(x, t)$ is an additional term to model wave breaking and is described in the next section. The parameter $\alpha$ in Eqs. (3) and (4) is used to tune the dispersive properties of the model. $\alpha = 1.159$ was suggested by Bonneton et al. (2011) to optimise the linear phase and group speeds of the model with respect to linear wave theory, and is used in this work. Setting $\alpha = 1$, one recovers the original SGN equations.

## 2.2. Wave breaking modelling

Similar to all depth-integrated models, the eSGN equations cannot model the overturning of the free surface and do not include any energy dissipation term, which means additional techniques are needed to deal with wave-breaking dissipation. Two methods are usually employed: (i) locally neglecting dispersion to solve the NSWE and model the breaking waves as shocks, which dissipate energy (the so-called hybrid approach); (ii) adding a diffusive-like term based on an eddy viscosity to dissipate energy (the eddy viscosity approach). Several definitions for this eddy viscosity have been proposed, but the one based on a one-equation turbulence model on the TKE, initially proposed by Nwogu (1996), has recently received a lot of interest (Zhang et al., 2014; Kazolea and Ricchiuto, 2018; Joshi et al., 2022; Kalisch et al., 2024).

In their thorough comparison between the hybrid and TKE approaches, Kazolea and Ricchiuto (2018) showed that both can correctly describe the spatial evolution of integral wave properties due to wave breaking, such as the wave height, for instance. Although the former is seemingly simpler and requires less parameter calibration, the switching between the dispersive equations and the NSWE introduces spurious oscillations, which can compromise the stability of the solution when a fine grid size is used. Galaz et al. (2025) evidenced that these oscillations are wave reflections at the switching interface. The TKE approach, however, allows for the computation of mesh-independent solutions. This approach seems more appropriate because the numerical modelling of wave run-up and overtopping requires fine grid sizes. Furthermore, a sensitivity analysis performed by Joshi et al. (2022) highlighted that the maximum run-up is highly influenced by the predicted breaking point when using the hybrid approach. In contrast, it is almost insensitive to it with the TKE model. Wave breaking in BT models is detected with wave breaking criteria, with onset and termination thresholds that need to be calibrated and may vary from one case to another. This considerable sensitivity of the run-up to the breaking point is thus a significant drawback of the hybrid approach. For those reasons, the eddy viscosity approach based on the TKE was selected for this work.

The following diffusive-like term is added to the right-hand side (RHS) of the elliptic Eq. (4):

$$D_{wb} = (2v_T h u_x)_x \quad (8)$$

To specify the eddy viscosity $v_T$, a mixing length and a TKE model are used (Pope, 2000):

$$v_T = C_v \sqrt{k} l_T \quad (9)$$

$$l_T = \kappa h \quad (10)$$

$C_v = 0.55$ is a constant (Pope, 2000), $k(x, t)$ [m$^2$/s$^2$] is the TKE, and $l_T(x, t)$ [m] is the mixing length, which is assumed proportional to

the local total water depth $h$ as in Kazolea and Ricchiuto (2018). $\kappa$ is a constant parameter with typical values in the range [0.1,1], and is used to tune the intensity of the energy dissipation. The TKE is assumed constant over the water column, and is computed from the following depth-integrated evolution equation, similar to that proposed by Kazolea and Ricchiuto (2018):

$$(hk)_t + (huk)_x = h v_T k_{xx} + P - \varepsilon \quad (11)$$

The three terms on the RHS of Eq. (11) correspond to diffusion, production and destruction of TKE, respectively. The destruction term is (Pope, 2000):

$$\varepsilon = h C_v^3 \frac{k^{3/2}}{l_T} \quad (12)$$

Finally, the expression for the production term suggested by Nwogu (1996) is:

$$P = h B \frac{l_T^2}{C_v^{3/2}} |U_z(z = \eta)|^3 \quad (13)$$

When waves break, TKE is generated by the vertical gradient of the horizontal velocity $U(x, z, t)$ at the free surface. Only the depth-averaged horizontal velocity $u$ is directly available from the eSGN model. However, a vertical expansion of the horizontal velocity $U(x, z, t)$ accurate at $O((kd)^2)$ can be derived from the irrotationality condition, and is given by Khakimzyanov et al. (2020) and Lannes (2013):

$$U(x, z, t) = u(x, t) + \left[ \frac{1}{2}h - (z + d) \right] (2u_x d_x + u d_{xx}) + \left[ \frac{1}{6}h^2 - \frac{1}{2}(z + d)^2 \right] u_{xx} \quad (14)$$

The vertical gradient of the horizontal velocity at the free surface is then:

$$U_z(x, z = \eta, t) = -h u_{xx} - 2u_x d_x - u d_{xx} \quad (15)$$

The coefficient $B(x, t)$ in Eq. (13) equals 1 or 0, whether wave breaking is activated or not. TKE is thus produced locally where and when wave breaking happens. Wave breaking is detected with the following criterion, ratio of the time derivative of the free surface to the shallow water phase speed, used by Kennedy et al. (2000):

$$\frac{\eta_t}{\sqrt{gh}} > \gamma \quad (16)$$

This criterion is commonly used to detect depth-induced wave breaking in phase-resolving wave models, see e.g. Borthwick et al. (2006), Lynett (2006), McCabe et al. (2013), Smit et al. (2013), Simon et al. (2019) and Kazolea et al. (2023). The onset of wave breaking is detected when this ratio exceeds a threshold $\gamma$, with values in [0.35,0.65] usually allowing to predict the correct breaking location (Kazolea et al., 2023). Breaking is then deactivated when the ratio falls below 0.15 (Kennedy et al., 2000).

## 3. Numerical methods

The eSGN Eqs. (1)–(4) and the TKE Eq. (11) are solved with a hybrid finite volume (FV)/finite difference (FD) scheme, which has become a standard approach for discretising BT models in the last two decades (Erduran et al., 2005; Tonelli and Petti, 2009; Bonneton et al., 2011; Roeber and Cheung, 2012; Shi et al., 2012; Kazolea and Delis, 2013; Lannes and Marche, 2015; Lteif, 2024). The hyperbolic NSWE part of the system is discretised with a Godunov-type FV scheme, while spatial derivatives of order up to three in the additional dispersive terms are computed with standard second-order FD formulas. The popularity of this approach comes from the relative simplicity of the numerical scheme and the well-established techniques used in FV schemes to model wetting and drying and obtain a scheme which preserves rest steady states (i.e. well-balanced scheme, or preserving the C-property).
A one-dimensional domain $[x_1, x_N]$ is considered, discretised with N grid points $(x_i)_{i=1,N}$ with a regular grid size $\Delta x = x_{i+1} - x_i$. A FV discretisation is applied to Eqs. (1) and (2) on dual cells $[x_{i-1/2}, x_{i+1/2}]$, with $x_{i\pm 1/2} = x_i \pm \frac{1}{2}\Delta x$.

### 3.1. FV discretisation of the hyperbolic part

The Eqs. (1) and (2) are written under the following vector conservative form:

$$U_t + F(U)_x = S_b + S_f + \Phi \quad (17)$$

with $U$ the vector of conserved variables, $F$ the flux, $S_b$ the shallow water bathymetry source term, $S_f$ the bottom friction term and $\Phi$ the dispersive source term:

$$ 
\begin{aligned}
U &= \begin{pmatrix} h \\ hu \end{pmatrix}, \quad F = \begin{pmatrix} hu \\ hu^2 + \frac{1}{2}gh^2 \end{pmatrix}, \\
S_b &= \begin{pmatrix} 0 \\ ghd_x \end{pmatrix}, \quad S_f = \begin{pmatrix} 0 \\ S_f \end{pmatrix}, \quad \Phi = \begin{pmatrix} 0 \\ \Phi \end{pmatrix}
\end{aligned} 
\quad (18) $$

The FV method will be used to evolve in time averaged values over the dual cells of the variables:

$$U_i = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} U \, dx \quad (19)$$

Integration of Eq. (17) over the dual cells gives the following semi-discrete scheme:

$$\frac{dU_i}{dt} = -\frac{1}{\Delta x} (F_{i+1/2} - F_{i-1/2}) + S_{bi} + S_{fi} + \bar{\Phi}_i \quad (20)$$

with $F_{i\pm 1/2}$ the fluxes at the cell interfaces, $S_{bi}$ the cell-averaged bathymetry source term, $S_{fi}$ the friction term, the treatment of which is described in Section 3.3, and $\bar{\Phi}_i$ the cell-averaged dispersive source term, approximated by its nodal value in $x_i$:

$$\bar{\Phi}_i \approx \begin{pmatrix} 0 \\ \Phi_i \end{pmatrix} \quad (21)$$

The fluxes at an interface with left and right states $U_{i\pm 1/2}^L$ and $U_{i\pm 1/2}^R$ are computed with the HLL Riemann solver (Harten et al., 1983):

$$F_{i\pm 1/2} = F \begin{pmatrix} U_{i\pm 1/2}^L & U_{i\pm 1/2}^R \end{pmatrix} \quad (22)$$

$$F(U^L, U^R) = \begin{cases} F(U^L) & \text{if } S_L \ge 0 \\ F^{\text{HLL}} & \text{if } S_L \le 0 \le S_R \\ F(U^R) & \text{if } S_R \le 0 \end{cases} \quad (23)$$

$$F^{\text{HLL}} = \frac{S_R F(U^L) - S_L F(U^R) + S_R S_L (U^R - U^L)}{S_R - S_L} \quad (24)$$

The wave speeds $S_{L/R}$ are computed with the following estimates, where $u_{L/R}$ and $h_{L/R}$ are the horizontal velocity and water depth on each side of the interface:

$$S_L = \min \left( u_L - \sqrt{gh_L}, u_* - \sqrt{gh_*} \right), \quad (25)$$

$$S_R = \max \left( u_R + \sqrt{gh_R}, u_* + \sqrt{gh_*} \right), \quad (26)$$

$$u_* = \frac{1}{2} (u_L + u_R) + \sqrt{gh_L} - \sqrt{gh_R}, \quad (27)$$

$$h_* = \frac{1}{g} \left[ \frac{1}{2} \left( \sqrt{gh_L} + \sqrt{gh_R} \right) + \frac{1}{4} (u_L - u_R) \right]^2. \quad (28)$$

If at least one of the two cells is dry, the wave speeds are instead computed as (Toro et al., 2020):

$$S_L = u_L - \sqrt{gh_L}, \quad S_R = u_L + 2\sqrt{gh_L} \quad \text{if } h_R \le \varepsilon_h \quad (29)$$

$$S_L = u_R - 2\sqrt{gh_R}, \quad S_R = u_R + \sqrt{gh_R} \quad \text{if } h_L \le \varepsilon_h$$

$\varepsilon_h$ is a small threshold used to identify dry cells. The appropriate values for this parameter are discussed later.

A third-order MUSCL reconstruction (Van Leer and Nishikawa, 2021) of $h$, $u$ and $\eta$ on each side of the cell interfaces is performed

to achieve sufficient accuracy. In a given cell and for variable $v$, the reconstruction on the interfaces is:

$$v_{i-1/2}^R = v_i - \frac{1}{6} [2(v_i - v_{i-1}) + v_{i+1} - v_i] \quad (30)$$

$$v_{i+1/2}^L = v_i + \frac{1}{6} [v_i - v_{i-1} + 2(v_{i+1} - v_i)] \quad (31)$$

Since dispersive equations are solved, and thanks to the addition of the energy dissipation term for wave breaking, no discontinuities or shocks are expected to arise in the solution. Thus, no use is made of any slope limiter, which would introduce unnecessary numerical diffusion (Jouy et al., 2024). The reconstruction is performed only if all cells in the stencil are wet. If at least one cell is dry, the reconstructed variables are the cell-averaged values $v_i$, and the scheme is locally of first order.

To guarantee the positivity of the water depth and that the scheme is well-balanced, the hydrostatic reconstruction method of Chen and Noelle (2017) and Dong and Li (2020) is used. From the reconstructed values of water depth and surface elevation on cell interfaces, the reconstruction of the still water depth $d$ is computed as:

$$d_{i+1/2}^{L/R} = h_{i+1/2}^{L/R} - \eta_{i+1/2}^{L/R} \quad (32)$$

A single still water depth is then defined at the interface, and new reconstructed water depths are computed as:

$$d_{i+1/2} = \max \left[ \min \left( d_{i+1/2}^L, d_{i+1/2}^R \right), -\min \left( \eta_{i+1/2}^L, \eta_{i+1/2}^R \right) \right], \quad (33)$$

$$h_{i+1/2}^{-} = \min \left( h_{i+1/2}^{L}, \eta_{i+1/2}^{L} + d_{i+1/2} \right), \quad (34)$$

$$h_{i+1/2}^{+} = \min \left( h_{i+1/2}^{R}, \eta_{i+1/2}^{R} + d_{i+1/2} \right). \quad (35)$$

The HLL Riemann solver is applied with left and right states:

$$U_{i+1/2}^{L} = \begin{pmatrix} h_{i+1/2}^{-} \\ h_{i+1/2}^{-} u_{i+1/2}^{L} \end{pmatrix}, \quad U_{i+1/2}^{R} = \begin{pmatrix} h_{i+1/2}^{+} \\ h_{i+1/2}^{+} u_{i+1/2}^{R} \end{pmatrix} \quad (36)$$

A second-order well-balanced discretisation of the topography source term is (Dong and Li, 2020):

$$S_{bi} = \begin{pmatrix} 0 \\ S_{b,i-1/2}^{+} + S_{bc,i} + S_{b,i+1/2}^{-} \end{pmatrix} \quad (37)$$

$$S_{bc,i} = \frac{g}{2\Delta x} \left( h_{i-1/2}^{R} + h_{i+1/2}^{L} \right) \left( d_{i+1/2}^{L} - d_{i-1/2}^{R} \right) \quad (38)$$

$$S_{b,i-1/2}^{+} = \frac{g}{2\Delta x} \left( h_{i-1/2}^{R} + h_{i-1/2}^{+} \right) \left( d_{i-1/2}^{R} - d_{i-1/2}^{L} \right) \quad (39)$$

$$S_{b,i+1/2}^{-} = \frac{g}{2\Delta x} \left( h_{i+1/2}^{L} + h_{i+1/2}^{-} \right) \left( d_{i+1/2}^{L} - d_{i+1/2}^{R} \right) \quad (40)$$

### 3.2. FD discretisation of the elliptic equation

In order to compute the dispersive source term $\Phi_i$, the spatial derivatives in Eqs. (3) and (4) are discretised with the following second-order centred FD formulas:

$$v_{x,i} = \frac{v_{i+1} - v_{i-1}}{2\Delta x} \quad (41)$$

$$v_{xx,i} = \frac{v_{i+1} - 2v_i + v_{i-1}}{\Delta x^2} \quad (42)$$

$$v_{xxx,i} = \frac{v_{i+2} - 2v_{i+1} + 2v_{i-1} - v_{i-2}}{2\Delta x^3} \quad (43)$$

Spatial derivatives appearing in the advection, diffusion and production terms of the TKE Eq. (11) are discretised with the same FD formulas. The discretisation of the coefficients in front of $\Psi$, $\Psi_x$ and $\Psi_{xx}$ in the operator T (6) leads to:

$$T(\Psi)_i = C_{\psi} \Psi_i + C_{\psi_x} \Psi_{x,i} + C_{\psi_{xx}} \Psi_{xx,i} \quad (44)$$

$$C_{\psi} = \frac{1}{3} (h_{x,i}^2 + h_i h_{xx,i}) + d_{x,i}^2 - h_{x,i} d_{x,i} - \frac{1}{2} h_i d_{xx,i} \quad (45)$$

$$C_{\psi_x} = -\frac{1}{3} h_i h_{x,i} \quad (46)$$
$$C_{\psi_{xx}} = -\frac{1}{3}h_l^2 \qquad (47)$$

By discretising the derivatives of $\Psi'$, the complete discretisation of the elliptic Eq. (4) eventually gives the following linear system:

$$a_i \Psi_{i-1} + b_i \Psi'_{i} + c_i \Psi_{i+1} = R H S_i \qquad (48)$$

$$\text{with } b_i = 1 + \alpha \left( C_{\psi} - \frac{2}{\Delta x^2} C_{\psi_{xx}} \right) \qquad (49)$$

$$a_i = \alpha \left( -\frac{1}{2\Delta x} C_{\psi_x} + \frac{1}{\Delta x^2} C_{\psi_{xx}} \right) \qquad (50)$$

$$c_i = \alpha \left( \frac{1}{2\Delta x} C_{\psi_x} + \frac{1}{\Delta x^2} C_{\psi_{xx}} \right) \qquad (51)$$

$$RHS_i = -\frac{gh_i}{\alpha} \eta_{x,i} - h_i Q_i + D_{wb,i} \qquad (52)$$

This tridiagonal system is solved with the Thomas algorithm to get the values of $\Psi_i$ in all mesh nodes. The final dispersive source term $\Phi_i$ is readily obtained from Eq. (3):

$$\Phi_i = \Psi_i + \frac{gh_i}{\alpha} \eta_{x,i} \qquad (53)$$

### 3.3. Time integration, bottom friction

From the spatial discretisation of Eqs. (11) and (20) (omitting bottom friction so far), the following equation must be integrated in time:

$$\frac{dV}{dt} = L(V), \quad \text{with } V = \begin{pmatrix} h \\ hu \\ hk \end{pmatrix} \qquad (54)$$

This is done with an explicit third-order strong-stability-preserving Runge-Kutta scheme (Gottlieb et al., 2001), with adaptive time-step $\Delta t$. Starting from the solution $V^n$ at time $t^n$, the three following steps are applied to compute $V^{n+1}$ at time $t^n + \Delta t$:

$$V^{(1)} = V^n + \Delta t L(V^n), \qquad (55)$$

$$V^{(2)} = \frac{3}{4}V^n + \frac{1}{4}[V^{(1)} + \Delta t L(V^{(1)})], \qquad (56)$$

$$V^{n+1} = \frac{1}{3}V^n + \frac{2}{3}[V^{(2)} + \Delta t L(V^{(2)})]. \qquad (57)$$

The variable time step $\Delta t$ is computed from the following Courant-Friedrichs-Lewy (CFL) condition as:

$$\Delta t = \text{CFL min}_{i=1:N} \left( \frac{\Delta x}{|u_i| + \sqrt{gh_i}}, \frac{\Delta x^2}{v_{T,i}} \right), \qquad (58)$$

where CFL is fixed by the user. The first term in the brackets of Eq. (58) is the classical stability condition inherited from NSWE-like models, while the second one is a supplementary parabolic stability condition due to the presence of the diffusive-like wave-breaking term in Eq. (4).

Once $V^{n+1}$ is known, a first-order implicit scheme is used to discretise the bottom friction term (5). Following Varra et al. (2024), the new value of the momentum $q_{\text{new}}^{n+1}$ is computed as:

$$q_{\text{new}}^{n+1} = \frac{2q^{n+1}}{1 + \sqrt{1 + \frac{4g n^2 \Delta t |q^{n+1}|}{(h^{n+1})^{7/3}}}}. \qquad (59)$$

### 3.4. Treatment of wetting and drying

Because conservative variables are used, $u$ and $k$ must be computed from $q = hu$ and $hk$ by dividing them by the water depth $h$. Unphysically high values of $u$ and $k$ might be computed near the shoreline because of the division by small water depths. To avoid this, similarly to Filippini et al. (2016) a threshold $\epsilon_u$ is introduced, such that when in one node $h \le \epsilon_u$, all variables $u$, $q$, $k$ and $hk$ are set to 0. Furthermore, when $h \le \epsilon_h$, $h$ is also set to 0. $\epsilon_h$ is a second threshold used to identify dry cells, already mentioned in Section 3.1, and is smaller than

$\epsilon_u$. Both thresholds should be small, so that they do not perturb the advancement of the wet/dry front, but not too small or else spurious velocities and TKE values near the shoreline can appear. As pointed out by Kazolea et al. (2012) and Tonelli and Petti (2013), the choice of these thresholds is thus not obvious and depends considerably on the problem and mesh size. Working with equal thresholds $\epsilon_h = \epsilon_u$, Briganti and Dodd (2009b) and Roeber and Cheung (2012) suggested that a reasonable value for the threshold is:

$$\epsilon_u \approx \Delta x/1000. \qquad (60)$$

For our application, i.e. random wave overtopping, special care must be taken because the accurate modelling of the movements of the wet/dry fronts is necessary, as well as robustness of the numerical method. Indeed, the aim is to study the potentially violent and frequent run-up/run-down of thousands of waves. Because of the presence of first- and second-order spatial derivatives of the velocity in the equations, especially in the TKE production term, if spurious velocities appear near the shoreline, numerical difficulties and strong perturbations of the solution can occur. For those reasons, values of $\epsilon_u$ higher than the one given by the estimate (60) might be needed. In that case, $\epsilon_h$ should be chosen at least an order of magnitude smaller than $\epsilon_u$ not to affect the run-up. The values of $\epsilon_h$ and $\epsilon_u$ used in practical applications are given in Sections 4 and 5.

### 3.5. Boundary conditions and wave generation

Two boundary conditions are considered: fully reflective walls or wave generation/absorption with relaxation zones. For the fully reflective condition, a homogeneous Neumann condition on $h$ and a homogeneous Dirichlet condition on $u$ must be imposed. Ghost cells are introduced on each side of the domain to take these conditions into account in the FV scheme. The boundary condition for the dispersive source term $\Phi$ is (Kazolea et al., 2023):

$$\Phi = 0 \Leftrightarrow \Psi = -\frac{gh}{\alpha} \eta_x. \qquad (61)$$

This condition is enforced by replacing the boundary equations of the linear system (48) by a discrete version of Eq. (61).

Relaxation zones are used to generate waves and absorb reflected ones. Over a length $L$ next to the boundary, the target free surface and velocity signals are progressively imposed, while the computed ones are attenuated:

$$\eta_{\text{new}}(x, t) = C_r(x) \eta_{\text{imp}}(x, t) + (1 - C_r(x)) \eta(x, t) \qquad (62)$$ $$u_{\text{new}}(x, t) = C_r(x) u_{\text{imp}}(x, t) + (1 - C_r(x)) u(x, t)$$

The coefficient $C_r$ goes from 1 at the domain's boundary to 0 at the other end of the relaxation zone. The ramp function proposed by Engsig-Karup et al. (2006) is used (with $\xi$ the distance from the boundary of the domain):

$$C_r(\xi) = \left(1 - \frac{\xi}{L}\right)^5. \qquad (63)$$

Incident irregular sea states are generated by the linear superposition of a high number of wave components:

$$\eta_{\text{imp}}(x, t) = \sum_j a_j \cos(\omega_j t - k_j x + \phi_j), \qquad (64)$$

$$u_{\text{imp}}(x, t) = \sum_j \frac{\omega_j}{k_j d} a_j \cos(\omega_j t - k_j x + \phi_j). \qquad (65)$$

$a_j$, $\omega_j$ and $\phi_j$ are the input amplitudes, angular frequencies and phases of the wave components, and $d$ is the still water depth in the generation zone. The wave-numbers $k_j$ are computed from the dispersion relation of the linearised eSGN equations:

$$\frac{\omega_j^2}{k_j^2 g d} = \frac{1 + \frac{\alpha-1}{3}(k_j d)^2}{1 + \frac{\alpha}{3}(k_j d)^2}. \qquad (66)$$
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Exp. 1 – Characteristic frequency spectra for the eight sea states generated in the experiments from Villefer et al. (2024).</figcaption></figure>

The length *L* of the relaxation zone should be one to two times the typical wavelength of the dominant waves to absorb or generate.

The validation of our approach for the simulation of random waves overtopping impermeable dykes with smooth slopes will now be presented. Two experimental datasets, totalling 184 cases, with both non-breaking and breaking waves, and unimodal or bimodal sea states, are considered.

# 4. Experiments 1: Wave overtopping for non-breaking waves and complex sea-states

## 4.1. Overview of the experimental set-up

The first experiments considered were conducted in the OSU Institut Pytheas wind wave tank in Marseilles, France, and are described in Villefer (2022) and Villefer et al. (2024). The goal of the experiments was to study the average overtopping rate over smooth and rubble-mound breakwaters with a slope of 2:3 (V:H), for unimodal and bimodal sea states. Only the tests performed with the smooth slope, made of aluminium, are considered here. There are 73 such tests. The wave tank is 40 m long, 2.6 m wide, and the still water depth can vary between 72 and 74 cm. The bottom profile before the breakwater is horizontal, and so is the crest of the breakwater. The waves were not breaking before reaching the breakwater and surged up its slope.

Eight different sea states, termed A to H, were generated linearly with a piston-type wavemaker and without active reflection compensation. Sea states A and B are unimodal swell and wind wave sea states, respectively, while sea states C to H are bimodal cases, obtained from

<caption>**Table 1**</caption>

Exp. 1 – Target wave conditions at model's scale in the experiments from Villefer et al. (2024). Corresponding representative measured spectra are plotted in Fig. 1.

<table><thead><tr><th rowspan="2">Cases</th><th colspan="2">Swell</th><th colspan="2">Wind waves</th></tr><tr><th>H<sub>m0.5</sub> (cm)</th><th>T<sub>p.5</sub> (s)</th><th>H<sub>m0.5,ww</sub> (cm)</th><th>T<sub>p,ww</sub> (s)</th></tr></thead><tbody><tr><td>A</td><td>5.6</td><td>1.67</td><td>-</td><td>-</td></tr><tr><td>B</td><td>-</td><td>-</td><td>5.2</td><td>1</td></tr><tr><td>C</td><td>4.8</td><td>1.67</td><td>4.8</td><td>1</td></tr><tr><td>D</td><td>3.6</td><td>1.67</td><td>4.4</td><td>1</td></tr><tr><td>E</td><td>3.6</td><td>1.67</td><td>4.0</td><td>0.75</td></tr><tr><td>F</td><td>3.6</td><td>1.67</td><td>2.4</td><td>0.66</td></tr><tr><td>G</td><td>3.6</td><td>1.67</td><td>2.0</td><td>0.56</td></tr><tr><td>H</td><td>3.6</td><td>1.67</td><td>5.4</td><td>1</td></tr></tbody></table>

the superposition of two JONSWAP spectra (Hasselmann et al., 1973). The characteristics of the sea states (target significant wave heights and peak periods of the swell and wind wave components) are reported in Table 1, while one typical frequency spectrum from each sea state is plotted in Fig. 1.

The spectra are coloured by their representative wave steepness $s_{m-1,0}$ computed with the incident significant wave height $H_{m0}$ and the deep-water wavelength $L_{m-1,0}$ associated with the energy period $T_{m-1,0} = m_{-1}/m_0$:

$$s_{m-1,0} = \frac{H_{m0}}{L_{m-1,0}}, \qquad (67)$$

$$L_{m-1,0} = \frac{gT_{m-1,0}^2}{2\pi}, \qquad (68)$$

$$m_n = \int_0^\infty f^n E(f) df. \qquad (69)$$

The total incident significant wave height being roughly constant among the various tests, around 5 – 6 cm, the wave steepness allows for distinguishing between the different sea states, since a sea state with a higher proportion of wind waves will have a lower energy period and a higher steepness. Five irregularly spaced wave gauges (WG) were placed before the breakwater to compute the incident wave characteristics, using the linear reflection analysis method from Zelt and Skjelbreia (1992).

Configurations with various freeboard heights $R_c$ were studied, with dimensionless freeboards $R_c^* = R_c/H_{m0}$ ranging between 0.4 and 1.5. Overtopping volumes were collected in a container placed behind the crest of the dyke, and the total volume was measured at the end of each run to deduce the average overtopping rate. Each trial lasted around $D = 30$ min, so at least 1000 waves based on the peak period interacted with the breakwater ($D/T_{p,s} \approx 1080$). The measured overtopping rates ranged between 0.13 L/s/m and 2.66 L/s/m.

## 4.2. Numerical set-up and parameters

In the numerical model, a relaxation zone of 4 m is used near the left boundary to generate incident waves and absorb those reflected by the breakwater (the wavelength associated with the swell peak period $T_{p,s} = 1.67$ s is approximately 3.7 m). First-order (linear) wave generation is used, as in the experiments. As pointed out by Williams et al. (2014) and Romano et al. (2015), different incident wave time series characterised by the same wave spectrum can lead to different overtopping statistics. To compare the computed average overtopping discharges to the measured ones, the exact amplitudes and phases of the incident wave components in the experiments are used to build the incident wave signal. Between 4600 and 6000 frequency components are prescribed, in the range [0.2, 2.5] Hz. Five numerical WG are set at the same location as in the experiments to compute and compare the incident wave fields. A numerical gauge at the end of the breakwater crest records the flow rate time series to calculate the
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Exp. 1 – Numerical setup for simulating one case from Villefer et al. (2024) with a bimodal sea state D ($H_{m0}$ = 6.31 cm, $R_c$ = 3.2 cm, see Table 1 for detailed wave parameters). Wave gauges (WG) 1 to 5, used for the reflection analysis, are located at $x = -3.22$ m, $-3.06$ m, $-2.85$ m, $-2.3$ m and $-2$ m, respectively. The breakwater toe is located at $x = 0$ m.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Exp. 1 – Relative error on the overtopping rates for all cases with sea state F and five decreasing grid sizes.</figcaption></figure>

average discharge. A relaxation zone is used behind the breakwater to absorb overtopping waves. Fig. 2 shows the setup for one run with a bimodal spectrum of case D, with a low freeboard height $R_c^* = 0.51$.

All 73 cases were simulated with the same duration as the experiments, i.e. $D = 30$ min. A Manning coefficient of $n = 0.01$ s/m$^{1/3}$ is used to model bottom friction on the smooth aluminium slope. Computations with $n = 0.012$ s/m$^{1/3}$ and $n = 0.014$ s/m$^{1/3}$ were also performed, leading to average overtopping discharges 3.3% and 7.1% lower, respectively. The grid size $\Delta x$ and CFL number were selected after dedicated convergence studies were done on the mean overtopping discharge $q$ computed over the test duration $D$. As an example, Fig. 3 shows the convergence with 5 decreasing $\Delta x$ (from 8 cm to 0.5 cm) of the dimensional $q$ for all cases of sea state F. The relative error on the overtopping rate is represented as a function of the ratio of the deep-water wavelength $L_{m-1,0}$ and the grid size $\Delta x$. The stabilisation of the error with the smallest grid sizes indicates that convergence has been reached. Note that, in this figure, the points are coloured by the error on the representative steepness, namely $s_{m-1,0,comp} - s_{m-1,0,exp}$ (in %); this provides an additional indication of the convergence of the incident wave conditions and spectrum, through the integral measure of wave steepness. The incident wave spectra are already well resolved with $\Delta x = 2$ cm since the error on $s_{m-1,0}$ does not change significantly with lower $\Delta x$. However, smaller grid sizes are necessary to model the run-up and overtopping accurately. Similar or faster convergence of $q$

is observed for the other sea states. A similar convergence study was done regarding the time discretisation, allowing us to finally select a grid size of $\Delta x = 5$ mm, and a Courant number CFL = 0.3, subsequently used for all simulations. Performing computations with such fine $\Delta x$ is permitted by the TKE breaking model, and its robustness is thereby demonstrated.

The wave breaking model is needed to deal with the bores generated behind the crest when the overtopping volumes run down the back slope of the breakwater. A small amount of TKE is sometimes also produced near the shoreline in front of the breakwater. The simulations are run with the breaking detection and strength parameters $\gamma = 0.5$, taken in the middle of the typical range of values [0.35,0.65], and $\kappa = 0.1$. A sensitivity study (not shown here) has demonstrated that the overtopping discharge in these experiments is insensitive to the value of $\gamma$, leading to the above choice. The small value $\kappa = 0.1$ was found to be necessary, otherwise the overtopping discharge is significantly underestimated. For example, when using $\kappa = 0.5$, which is on the lower bound of values we typically use to simulate the propagation of breaking waves, the computed overtopping rate is lower by 29%. It can be concluded that when wave breaking is not expected to happen on the breakwater slope, a small $\kappa$ of the order of 0.1 should be used.

The thresholds for wetting and drying were set as $\epsilon_u = \Delta x/100$ and $\epsilon_h = \epsilon_u/10$. Using values of magnitude higher slightly decreased the accuracy of the overtopping rate estimations, while using values one or several orders of magnitude lower did not improve the results, and even led to some computations blowing up. As explained in Section 3.4, solution blow up might occasionally occur when too low values of $\epsilon_u/h$ are used, because spurious values of the velocity $u$ could be computed near the shoreline, causing large values of the TKE production term which disturb the solution.

### 4.3. Wave propagation up to the breakwater

First, the ability of the model to generate and propagate the wave spectra correctly is checked. Fig. 4 compares the incident spectra of one case of each sea state, computed from a linear reflection analysis performed using the wave signals recorded at WG 1 to 5 with the least squares method from Zelt and Skjelbreia (1992). The agreement is very satisfactory, with only local differences. Note that for frequencies above 1.1 Hz approximately, the associated wavenumbers $k$ are such that $kd \ge 3.5$, therefore, these wave components are deep water waves. The model theoretically shows increasing dispersion errors for such waves, which explains the discrepancies at high frequencies. Table 2 shows the mean ratio of the incident $H_{m0}$ in the simulations and experiments
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Exp. 1 – Incident wave spectra in the deep region, in the experiments (blue lines) and the simulations (red lines) for one case of each sea state. The spectra are obtained from a linear reflection analysis performed with gauges 1 to 5.</figcaption></figure>

**Table 2**

<table><caption>Exp. 1 – Mean ratio of incident significant wave heights in the computations and experiments, mean index of agreement *d*<sub>r</sub> (defined by eq. (75)) of the total and incident free surface elevation time series at WG 5, for all the sea states.</caption><thead><tr><th>Cases</th><th>Mean<br>H<sub>m0,comp</sub>/H<sub>m0,meas</sub></th><th>Mean d<sub>r</sub> (total wave<br>signal)</th><th>Mean d<sub>r</sub> (incident<br>wave signal)</th></tr></thead><tbody><tr><td>A</td><td>0.97</td><td>0.53</td><td>0.90</td></tr><tr><td>B</td><td>0.98</td><td>0.45</td><td>0.87</td></tr><tr><td>C</td><td>0.97</td><td>0.51</td><td>0.87</td></tr><tr><td>D</td><td>0.97</td><td>0.49</td><td>0.87</td></tr><tr><td>E</td><td>0.96</td><td>0.51</td><td>0.85</td></tr><tr><td>F</td><td>0.95</td><td>0.49</td><td>0.84</td></tr><tr><td>G</td><td>0.95</td><td>0.49</td><td>0.84</td></tr><tr><td>H</td><td>0.96</td><td>0.46</td><td>0.86</td></tr></tbody></table>

for each sea state. The significant wave height is systematically slightly underestimated by 2 to 5% depending on the sea state.

**Fig. 5** shows samples of free surface elevation time series of the same cases at the WG closest to the dyke (WG 5), to check that the phases of the wave components are correct. Overall, the phases of the long waves are well correlated in the experiments and simulations, although significant differences in the wave amplitudes are observed. The limited dispersive properties of the model partly cause this, but the main reason behind these discrepancies is the differences in wave components reflected by the breakwater. This was also found to be the source of discrepancies in free surface elevation time series in Vanneste et al. (2014). Here, the waves being essentially non-breaking, the incident wave energy not transmitted behind the breakwater due to overtopping is mostly reflected. The overtopping not being perfectly reproduced by the model, as will be seen later, reflected waves are
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Exp. 1 – Examples of free surface elevation time series at WG 5 in the experiments (solid blue lines) and the simulations (dashed red lines) for one case of each sea state.</figcaption></figure>

partially different in the simulations and the experiments, explaining the differences in Fig. 5. To confirm this, the incident free surface time series at WG 5 were reconstructed, using the wave amplitudes and phases computed in the reflection analysis. The results are presented in Fig. 6, for the same cases as in Fig. 5. One can see that the experimental and simulated time series now match very well, confirming that the differences noticed previously are primarily due to wave reflections.

To perform a quantitative analysis of the time series and as in Al-
tomare et al. (2021), the Willmott index $d_r$ was computed (Willmott
et al., 2012). It measures the general accuracy of the simulations and
is defined in Appendix (Eq. (75)). $d_r$ is bounded by -1 and 1, and
the closer it is to 1, the more accurate the results are. Table 2 shows
the mean index $d_r$ for the total wave signal (incident + reflected) at
WG 5, and the reconstructed signal at the same WG, and each sea
state. Gruwez et al. (2020) proposed a classification of the results
depending on the values of $d_r$. According to it, our results on the
total time series are rated as “reasonable” ($d_r \in [0.5, 0.7]$) or “poor”

$(d_r \in [0.3, 0.5])$ depending on the sea state. When it comes to the
reconstructed incident time series, however, the results are rated as
"very good" $(d_r \in [0.8, 0.9])$ for all sea states.

4.4. Mean wave overtopping discharges

The results concerning the estimation of the average overtopping
discharge are discussed now. An issue for the study of wave overtopping
in an operational context is the long physical time that must be simu-
lated to get a reasonable estimation of mean statistics. Before validating
the computed average overtopping discharges with the experimental
data, an analysis of the time evolution of these discharges in the simula-
tions is presented, to determine how many waves must interact with the
structure for them to converge. The objective is to estimate the length
of the time series that will be required to obtain a reliable estimation
of the average discharge. Fig. 7 shows for all cases the minimum
duration of the time series (expressed in number of waves based on
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Exp. 1 – Examples of incident free surface elevation time series reconstructed at WG 5 in the experiments (solid blue lines) and the simulations (dashed red lines) for one case of each sea state.</figcaption></figure>

the peak period of the spectrum) so that the overtopping rate stays within 5% of the one computed with the whole time series. Significant differences between the different sea states can be observed. For the lowest freeboard heights, the required number of waves decreases when the proportion of swell increases. Only 300 waves are needed for sea state A (pure swell), compared to around 1000 for case B (wind waves only). The required number increases by a few hundred waves for high freeboards, as observed experimentally by Romano et al. (2015). The recommendation from EurOtop (2018) is to include 500 to 1000 waves in the time series to have converged statistics. From our results, using 800 to 1000 waves indeed appears necessary to estimate correctly the discharge in most cases.

The relative errors on the average overtopping discharge for all cases are represented in Fig. 8, considering both the dimensional results in Fig. 8(a), and the non-dimensional discharges $q^* = q/\sqrt{gH_{m0}^3}$ in Fig. 8(b). The incident $H_{m0}$ from the reflection analysis of the experiments and the simulations are used to scale the measured and

computed discharges, respectively. The dimensional overtopping rates are underestimated by -22% on average. Slight underestimations of the incident $H_{m0}$ in the numerical simulations (as discussed in Section 4.3 and which can also be noticed in the plots of Fig. 4) account partially for this underestimation. This is confirmed by the reduced value of the average error on the dimensionless discharge $q^*$ of -18% in Fig. 8(b). The agreement is satisfactory over the whole range of measured overtopping rates, although the scatter of the results is slightly more important for low values. In such cases, bottom friction highly influences overtopping and is only approximately modelled with depth-averaged equations. Suzuki et al. (2017) suggested that surface tension, not modelled here, may also significantly influence overtopping at the physical model's scale.

Finally, Fig. 9 shows, for the experiments (grey points) and the computations (coloured points), the dimensionless overtopping rate as a function of the dimensionless freeboard height. The reference formula from the EurOtop (2018) manual to predict the overtopping discharge
<figure><img src="image_9.png" /><figcaption>Fig. 7. Exp. 1 – Number of waves necessary so that the estimated average discharge *q* stays within 5% of the discharge computed with the whole time series.</figcaption></figure>

in this configuration is:

$$ \frac{q}{\sqrt{gH_{m0}^3}} = 0.09 \exp \left( - \left( 1.5 \frac{R_c}{H_{m0}} \right)^{1.3} \right). \quad (70) $$

The formula is plotted in solid black, while the upper and lower bounds of its 90% confidence interval are plotted in dashed black. The experimental data agree very well with the formula, with a mean normalised bias (MNB) of -0.3% and a root mean squared percentage error (RMSPE) of 10.7% (both defined in Appendix). The MNB and RMSPE computed from the simulation results are -11.3% and 15.7% respectively, a bit larger as expected from the previous comparisons, but still low considering the difficulty of predicting mean overtopping rates numerically. It was observed in Villefer et al. (2024) that the overtopping rate increases with the proportion of swell in the sea state when the relative freeboard height is low, and that the points with a higher proportion of wind waves (i.e. characterised by a larger value of the representative wave steepness) follow a less steep slope as a function of $R_c^*$. This behaviour is also clearly observable in our simulation results.

### 4.5. Summary of model's accuracy and efficiency

These 73 cases showed that the proposed eSGN model can simulate the overtopping of a smooth dyke for complex sea states (i.e. unimodal and bimodal spectra) and essentially non-breaking waves, for mean overtopping discharges of the order 0.1–3 L/s/m at model's scale. The accuracy of the numerical predictions is demonstrated by a low underestimation of the mean overtopping discharges, by about -20% on average in comparison with the measured discharges, and by -11.3% only on average in comparison with the reference formula (70) from EurOtop (2018).

To illustrate the efficiency of the current version of the code (without any parallelisation at this stage), the simulation of 30 min of physical time at model's scale (i.e. ≈1100 waves with period $T_p$) took approximately 5 h with $\Delta x = 1$ cm and CFL = 0.3 on one core of an Intel Xeon 6140 processor, and around 16 h with $\Delta x = 5$ mm.

## 5. Experiments 2: Wave overtopping of breaking waves over dykes with very shallow foreshores

### 5.1. Overview of the experimental set-up

The second set of experiments was conducted at Flanders Hydraulics and is described in Altomare et al. (2016) and Suzuki et al. (2017). The simulated cases belong to the experimental campaigns referred

<caption>Table 3</caption>

Exp. 2 – Range of water depths, wave conditions, freeboard heights and measured mean overtopping rates in the selected cases from the overtopping experiments of Suzuki et al. (2017).

<table><thead><tr><th>Dataset</th><th><i>h</i> (m)</th><th><i>T</i><sub>p</sub> (s)</th><th><i>H</i><sub>m0</sub> (m)</th><th><i>R</i><sub>c</sub> (cm)</th><th><i>q</i> (L/s/m)</th></tr></thead><tbody><tr><td>00–025</td><td>0.935–0.989</td><td>2.16–2.56</td><td>0.18–0.19</td><td>4.48–12.6</td><td>0.1212–1.824</td></tr><tr><td>13–116</td><td>0.9–0.935</td><td>1.7–2.5</td><td>0.1–0.25</td><td>7.6–11.5</td><td>0.0013–0.25</td></tr></tbody></table>

to as datasets ‘00–025’ and ‘13–116’ from Suzuki et al. (2017). The experiments studied the overtopping of irregular waves over smooth impermeable dykes with very shallow foreshores. The slope of the dyke is either 1:2 for dataset 00–025, or 1:3 or 1:6 for dataset 13–116. The dyke is preceded by a mildly sloping foreshore of slope 1:35, which induces wave shoaling and breaking. For some of the cases of dataset 00–025, a promenade of slope 1:100 is placed at the crest of the breakwater. The crest is horizontal in the rest of the cases. These experiments are used to validate the model for irregular wave propagation over varying bottoms, wave breaking and overtopping. 21 cases from dataset 00–025 and 90 from dataset 13–116 were simulated, with varying water depths and wave conditions. In the experiments, irregular waves were generated linearly with a piston-type wavemaker. The main characteristics of these cases are reported in Table 3.

### 5.2. Numerical set-up and parameters

Fig. 10 shows examples of the numerical setups for one case of each dataset, with numerical WG located at the same position as the resistance-type WG used in the experiments. The three most offshore WG for dataset 00–025, and the two most offshore for dataset 13–116, are used to compute the experimental incident wave field provided as input to the numerical model. Between 5000 and 5300 wave components are prescribed in the range [0.075, 1.7] Hz, and generated with a linear superposition method.

All the trials lasted around $D = 40$ min at the model's scale, and the same duration is simulated. Depending upon the peak wave period, this corresponds to a range of 960 waves ($T_p = 2.5$ s) to 1410 waves ($T_p = 1.7$ s). A Manning coefficient of $n = 0.012 \text{ s/m}^{1/3}$ is used, as in Suzuki et al. (2017). This value is in the recommended range ([0.01, 0.014] s/m$^{1/3}$) to model bottom friction on smooth wood (Chow, 1959), which is the material the dyke was made of. Computations with $n = 0.01 \text{ s/m}^{1/3}$ and $n = 0.014 \text{ s/m}^{1/3}$ were also performed, leading to 22% higher and 19% lower overtopping discharges, respectively.

Convergence studies similar to the ones of the previous section were repeated, leading to a selected grid size of $\Delta x = 1$ cm, and a Courant number CFL = 0.3. As an illustrative example, Fig. 11 shows the grid convergence of the average discharge $q$ for all 21 cases of dataset 00–025. The representative wavelength at the toe is here computed as in Suzuki et al. (2017) by $L_{m-1,0} = T_{m-1,0} \sqrt{g(d + H_{m0})}$ based on incident wave characteristics at the toe. It is observed that the errors do not evolve much with the smallest grid sizes. The mean ratio of incident wave height at the toe of the dyke in the simulation and experiments is 0.93 with $\Delta x = 8$ cm, 1.02 with $\Delta x = 4$ cm, and 1.04 for lower values. As for the previous case, the smallest grid sizes are not necessary to have converged results regarding wave properties, but they are for the run-up and overtopping estimations. Similar results were observed for dataset 13–116.

The value of the breaking detection parameter is again $\gamma = 0.5$. The breaking strength parameter $\kappa$ was set to 0.75, in the middle of the range [0.5, 1] that, in our experience, provided good results for the simulation of various types of breaking waves. Both parameters have a significant influence on the average overtopping rate. For instance, when using $\gamma = 0.4$ and 0.6, the computed average discharges are 30% lower and 31% higher on average. When using $\kappa = 0.5$ and 1, the computed average discharges are 32% higher and 22% lower on average. The wetting and drying thresholds are set to $\epsilon_u = \Delta x/100$
<figure><img src="image_10.png" /><figcaption>Fig. 8. Exp. 1 – Relative error on the mean overtopping discharges for all cases as a function of the measured mean overtopping discharges.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Exp. 1 – Dimensionless average discharge as a function of the dimensionless freeboard, in the simulations (symbols coloured by the representative wave steepness) and the experiments from Villefer et al. (2024) (grey symbols).</figcaption></figure>

(same value as in Section 4) and $\epsilon_h = \epsilon_u/100$. Waves are generated in relaxation zones of length equal to one peak wavelength of the incident spectrum.

### 5.3. Wave propagation up to the breakwater

Fig. 12 shows the evolution of the total wave spectrum at several WG for the case WEN_004, chosen as a representative case ($h = 0.942$ m, $T_p = 2.13$ s, $H_{m0} = 18.9$ cm, $R_c = 6.88$ cm). At the first WG, the amplitude of the peak is slightly underestimated. The model captures the energy transfers to lower and higher frequencies during shoaling. The same can be said about the energy dissipation and transfer to the infragravity (IG) range due to wave breaking. At the toe of the dyke, all the energy from the offshore spectrum as been dissipated or transferred to IG waves. Only long waves thus overtop the structure. Fig. 13 shows the evolution of the mean water level $\langle \eta \rangle$ (i.e. wave set-down and setup), the energy period $T_{m-1,0}$ and the total significant wave height

of the IG waves and the “short” waves, defined as (with $f_p$ the offshore peak frequency):

$$H_{m0,IG} = 4\sqrt{\int_{0}^{f_p/2} E(f) \mathrm{d}f}, \qquad (71)$$

$$H_{m0,SW} = 4\sqrt{\int_{f_p/2}^{\infty} E(f) \mathrm{d}f}. \qquad (72)$$

The skewness $Sk$ of the free surface elevation is also shown, defined as:

$$Sk = \frac{\langle (\eta - \langle \eta \rangle)^3 \rangle}{\sigma^3}, \qquad (73)$$

$$\sigma = \sqrt{\langle (\eta - \langle \eta \rangle)^2 \rangle}. \qquad (74)$$

The agreement between simulations and experiments is satisfactory, although the wave setup at the toe of the structure is slightly underestimated. The skewness is slightly overestimated at the end of
<figure><img src="image_12.png" /><figcaption>Fig. 10. Exp. 2 – Numerical setup for one case of each dataset of the experiments from Suzuki et al. (2017), with a visualisation of the TKE profile.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Exp. 2 – Relative error on the overtopping rates for all cases of dataset 00–025 and four decreasing grid sizes.</figcaption></figure>

the foreshore slope. Table 4 shows a summary of the estimated wave properties at the toe of the dykes for all cases. For dataset 00–025, wave properties at the toe are quite well estimated. For dataset 13–116, the wave heights are underestimated by 9%.

The free surface profiles at WG 03 and 14 for case WEN_004 are compared in Fig. 14. The agreement is very satisfactory, both regarding wave amplitude and phases, even at the toe of the dyke. The index of agreement $d_r$ of the time series at each WG (except the one at the edge

of the crest used to estimate the average discharge) was computed and averaged over all the cases. Table 4 shows the range of these averaged $d_r$. The rating of the free surface time series based on these values range from “very good” offshore, to “reasonable” at the toe of the dyke for both datasets.

### 5.4. Mean wave overtopping discharges

An analysis of the time evolution of the average discharge, similar to that presented before, is performed. The results are shown in Fig. 15. An estimation of the incident significant wave height at the toe of the dyke, obtained from a reflection analysis with the method of Zelt and Skjelbreia (1992) with three WG close to the toe, is used to scale the freeboard. The number of waves required (based on the offshore peak period) increases more notably with the freeboard height than in the previous case, from approximately 900 when $R_c/H_{m0} \le 1$ up to 1300 when $R_c/H_{m0} > 3$. This number also tends to be lower for longer incident waves, as was observed in the previous case.

Finally, Fig. 16 shows the computed overtopping rates as a function of the measured ones, for all 111 cases in Fig. 16(a), and only the 21 cases from dataset 00–025 in Fig. 16(b). In the first figure, the dashed lines show overestimation or underestimation by a factor of two. All but 21 points, all from dataset 13–116, lie between these lines. Overall, the overtopping rates are well reproduced over the three orders of magnitude of measured $q$, with a tendency to overestimate $q$ for cases with slope 1:3, and underestimate it for cases with slope 1:2. Following Suzuki et al. (2017), the geometric mean (Geo) and standard deviation (GSD) of the simulation results have been computed. Their definitions are given in Appendix. Table 5 summarises these statistics of the discharge estimation. The mean absolute relative error on $q$ for all 111 cases is 57%. For cases of dataset 00–025 in Fig. 16(b), for
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Exp. 2 – Total wave spectra in the experiments (blue lines) and the simulations (red lines) at 6 WG for case WEN_004 from dataset 00–025. The spectrum at WG 03 is duplicated with a grey dotted line in all other panels to better visualise the magnitude of spectral changes as waves approach the dyke’s toe.</figcaption></figure>

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

Exp. 2 – Mean ratio of total significant wave heights and energy period, mean difference on the wave setup at the toe of the dyke, and range of mean $d_r$ at the different WG, for datasets 00-025 and 13–116.

<table><thead><tr><th>Datasets</th><th>H<sub>m0,comp</sub>/H<sub>m0,meas</sub></th><th>T<sub>m-1,0,comp</sub>/T<sub>m-1,0,meas</sub></th><th>⟨η⟩<sub>comp</sub> - ⟨η⟩<sub>meas</sub> (m)</th><th>Range of mean d<sub>r</sub> at all WG</th></tr></thead><tbody><tr><td>00-025</td><td>1.02</td><td>0.90</td><td>-0.0033</td><td>0.54-0.87</td></tr><tr><td>13-116</td><td>0.91</td><td>1.14</td><td>-0.0028</td><td>0.53-0.86</td></tr></tbody></table>

which *q* is the largest, the mean error is only 13%, and exceeds 20% for seven cases. The estimation of large overtopping rates, above 0.3 L/s/m here, is thus quite accurate. As in the previous case, the scatter of the results is higher for low discharges. Geo = 1.18 is close to unity, which confirms that overall overtopping discharges are well estimated. GSD measures the scatter of the data, and is equal to 1.67 for the ensemble of 111 cases. For a normally distributed quantity, 90% of the data lies between Geo multiplied and divided by 1.64 × GSD.

A visual comparison of our Fig. 16(a) with the Fig. 7 of Suzuki et al. (2017), obtained when simulating cases of the same datasets with SWASH, indicates a comparable level of model performance, and an apparent weaker scatter in our case (noting that the dashed lines correspond to a factor 3 in their case). In their table 8, they report values of GSD of 1.4 for dataset 00–025, and of 2.12 and 2.04 for cases of dataset 13–116 with dyke slope 1:3 and 1:6 respectively, confirming that the scatter of our results is weaker. Their results for dataset 13–116 (Geo = 1.11 and 0.95 for dyke slopes 1:3 and 1:6, respectively) are more accurate than ours. For dataset 00–025, it is the opposite, as they report Geo = 0.7.

<caption>**Table 5**</caption>

Exp. 2 – Number of simulated cases *N* and statistics of the overtopping discharge estimations for the experiments from Suzuki et al. (2017). Geo and GSD are defined in Appendix.

<table><thead><tr><th>Dataset</th><th>N</th><th>Relative error (%)</th><th>Geo</th><th>GSD</th></tr></thead><tbody><tr><td>00-025</td><td>21</td><td>13</td><td>0.87</td><td>1.12</td></tr><tr><td>13-116</td><td>90</td><td>67</td><td>1.27</td><td>1.72</td></tr><tr><td>All</td><td>111</td><td>57</td><td>1.18</td><td>1.67</td></tr></tbody></table>

### 5.5. Summary of model’s accuracy and efficiency

This second set of experiments has proven the ability of the model to simulate irregular wave transformations over mildly sloping foreshores, with very intense wave breaking dissipation and energy transfers to the IG range. As illustrated in Figs. 12 (frequency spectra) and 13 (integral wave parameters), the incident sea states undergo tremendous changes before reaching the toe of the breakwater, and the numerical model succeeds in modelling the various physical processes at play. The
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Exp. 2 – Evolution of the total infragravity and short wave heights ($H_{m0,IG}$, $H_{m0,SW}$), mean water level $\langle\eta\rangle$ (wave set-down/setup), energy period $T_{m-1,0}$ and skewness $Sk$ in the experiments and the simulation for case WEN_004 of dataset 00–025.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Exp. 2 – Examples of free surface elevation time series at the most offshore WG and the one closest to the toe of the dyke for case WEN_004.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Exp. 2 – Number of waves necessary so that the average discharge $q$ stays within 5% of the discharge computed with the whole time series. Slope 1:2 corresponds to dataset 00–025, and the other points correspond to cases from dataset 13–116 with dyke slope 1:3 or 1:6.</figcaption></figure>

overtopping rates over steep dykes with shallow foreshores are well estimated over three orders of magnitude, between 0.01 and 1.8 L/s/m at model's scale as summarised in Fig. 16. In particular, the numerical predictions of the largest overtopping discharges (dataset 00–025) show an average relative error of 13% in comparison with the measurements, and this relative error barely exceeds ±20%.

A single set of parameter values ($\kappa$ = 0.75, $\gamma$ = 0.5) was applied to simulate all 111 cases from both datasets, 00–025 and 13–116. While using distinct parameter values for each dataset would certainly enhance the accuracy of the results, a unified parameterisation was deliberately chosen to ensure consistency across simulations. Using the experimental data, $\gamma$ could be further calibrated by comparing the location of the breaking point in the simulations and the experiments, while $\kappa$ could be calibrated by comparing the rate at which the significant wave height decreases in the breaking zone along the slope, or by comparing the overtopping rates directly. However, and as can be seen from the literature (e.g. Kazolea and Ricchiuto, 2018; Joshi et al., 2022), the appropriate values for these breaking parameters are somewhat case-dependent; in particular, they are expected to depend on the geometry of the problem (seabed shape and slope) and the incident wave conditions. Empirical laws to compute these values
<figure><img src="image_18.png" /><figcaption>(a) Datasets 00-025 and 13-116 (111 tests in total), plotted using a logarithmic scale. The dashed lines indicate overestimation or underestimation by a factor of 2.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Dataset 00-025 only (21 tests), plotted using a linear scale. The dashed lines indicate overestimation or underestimation by 20%.</figcaption></figure>

<caption>**Fig. 16.** Exp. 2 – Computed vs. measured averaged dimensional overtopping rates for the experiments of Suzuki et al. (2017).</caption>

are not available currently, and efforts should be made to provide more quantitative guidelines or parameterisations to estimate them by considering a more comprehensive and diverse set of experimental conditions.

In terms of computational resources, the simulations of $D = 40$ min of physical time at the model's scale (roughly 1000–1400 waves depending upon $T_p$) took approximately 13 h with $\Delta x = 1$ cm and CFL = 0.3 on one core of an Intel Xeon 6140 processor.

## 6. Conclusions and outlook

A one-dimensional numerical model solving the enhanced Serre-Green-Naghdi (eSGN) equations of Bonneton et al. (2011) with a hybrid FV/FD scheme was developed. Some salient features of this model are: (i) the use of a third-order MUSCL reconstruction of the flow variables in the FV phase to achieve low discretisation errors and of second-order FD schemes to compute the dispersive terms, (ii) a robust treatment of wetting and drying using the hydrostatic reconstruction of Dong and Li (2020), and (iii) the modelling of wave breaking with an accurate and robust eddy viscosity approach, based on the solution of an evolution equation for the TKE. Dedicated convergence studies have allowed us to study the effect of spatial and temporal discretisation parameters, and to select optimal values for subsequent simulations. The model was then applied to reproduce an ensemble of 184 tests from two distinct (and very dissimilar) experimental campaigns of irregular long-crested wave propagation and overtopping, with mean overtopping discharges varying over three orders of magnitude. While only the average overtopping discharge was studied here, this deterministic model can theoretically reproduce the characteristics of individual overtopping events, such as the thickness of the overtopping layer, the overtopping velocity or volume. As no data on instantaneous overtopping were available in the selected experimental campaigns, results on individual overtopping events were not discussed. The validation of the model for that purpose is then left for future work.

The 73 cases of the experiments from Villefer et al. (2024) demonstrate the model's ability to estimate overtopping rates of essentially non-breaking waves characterised by complex sea states (i.e. unimodal and bimodal spectra, with different peak periods and representative wave steepness). The results are very satisfactory, as summarised in Section 4.5, even though a significant part of the energy was carried by deep-water waves, which are theoretically too short to be propagated correctly by a weakly dispersive BT model. Overall, the mean relative errors of the simulated overtopping discharges are about -20% compared with the measurements, and about -11.3% compared with the relevant formula from EurOtop (2018).

The second set of experiments from Suzuki et al. (2017) shows that the model can simulate complex wave transformation over varying bathymetries, including the transfer of energy towards the IG range. The TKE breaking model describes well the observed strong energy dissipation due to wave-breaking. Variations of the mean water level (i.e wave set-down and setup) are also fairly well simulated. These cases are very demonstrative of the need to use an accurate wave model to transfer offshore wave conditions to the toe of the structure, as summarised in Section 5.5. Again, mean overtopping rates are very well estimated, with a 13% relative error on average for cases with the highest discharges (dataset 00-025), and a slightly higher scatter for lower discharges.

The model's robustness is highlighted by the ability to use very fine grid sizes when needed to describe the run-up and overtopping processes accurately. It was shown that wave properties converge faster than overtopping rates with decreasing grid sizes. This suggests that a coarser resolution in the wave propagation regions could boost the model's efficiency. Irregular meshes, or techniques such as grid nesting, should thus be used. In addition, the efficiency of the model could be improved by parallelising it or using GPU computations, for instance. Nevertheless, the current sequential version can already be used in one-dimensional engineering applications with typical CPU times in the order of 5–10 h to simulate sea states of more than 1000 waves on standard scientific computers.

This eSGN model is currently being extended to two horizontal dimensions and parallelised in order to test its applicability to study wave transformation over complex nearshore seabeds and wave overtopping over real coastal structures. Future work will also address the modelling of porous coastal structures to better model rubble-mound breakwaters and overtopping of coastal structures equipped with a vertical crown wall.
## CRediT authorship contribution statement

**Guillaume Coulaud:** Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization.
**Maria Teles:** Writing – review & editing, Supervision.
**Michel Benoit:** Writing – review & editing, Supervision, Software, Conceptualization.

## 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 was carried out as part of EDF R&D's MOISE-2 research project on coastal and river flood hazard assessment, whose support the authors gratefully acknowledge. Financial support for Guillaume Coulaud's Ph.D. research program was provided by the French ANRT (Association Nationale de la Recherche et de la Technologie) through the CIFRE grant 2022/1017. The authors are grateful to Flanders Hydraulics (Belgium), and especially Toon Verwaest and Tomohiro Suzuki (now at KU Leuven, Belgium), for providing the data and additional explanations about the wave overtopping laboratory experiments described in Suzuki et al. (2017). Finally, the authors thank the two anonymous reviewers for providing insightful comments and suggestions on an earlier version of this article.

## Appendix. Statistical quantities used to measure model performance

To assess the model performance, the following statistical indices are defined, for two series of N data points, $P_i$ denoting the Prediction (here simulation results) and $O_i$ denoting the Observation (either measurements or other reference results):

* Willmott index $d_r$ (Willmott et al., 2012):

$$ 
\begin{aligned}
d_r &= \begin{cases} 1 - \frac{MAE}{2 \times MAD} & \text{if } MAE \le 2 \times MAD \\ \frac{2 \times MAD}{MAE} - 1 & \text{if } MAE > 2 \times MAD \end{cases} \\ 
MAE &= \frac{1}{N} \sum_{i=1}^{N} |P_i - O_i| \\ 
MAD &= \frac{1}{N} \sum_{i=1}^{N} |O_i - \bar{O}|
\end{aligned} 
\quad (75) 
$$

MAE is the mean absolute error and MAD is the mean absolute deviation of the observation. $d_r$ is bounded by -1 and 1, and the closer to 1 is its value, the better is the accuracy of the results.

* Mean normalised bias (MNB):

$$ MNB = \frac{1}{N} \sum_{i=1}^{N} \frac{P_i - O_i}{O_i} \quad (76) $$

* Root mean squared percentage error (RMSPE):

$$ RMSPE = \sqrt{\frac{1}{N} \sum_{i=1}^{N} \left( \frac{P_i - O_i}{O_i} \right)^2} \quad (77) $$

* Geometric mean (Geo):

$$ Geo = \exp \left( \frac{1}{N} \sum_{i=1}^{N} \ln \left( \frac{P_i}{O_i} \right) \right) \quad (78) $$

* Geometric standard deviation (GSD):

$$ GSD = \exp \left( \sqrt{\frac{1}{N} \sum_{i=1}^{N} \left( \ln \left( \frac{P_i}{O_i} \right)^2 - \ln(Geo)^2 \right)} \right) \quad (79) $$

## Data availability

The authors do not have permission to share the experimental data of the two campaigns used to validate the numerical model. The results of the numerical model can be partially made available by the authors upon reasonable request.

## References

Akbari, H., 2017. Simulation of wave overtopping using an improved SPH method. Coast. Eng. 126, 51–68. http://dx.doi.org/10.1016/j.coastaleng.2017.04.010.

Altomare, C., Gironella, X., Crespo, A.J.C., 2021. Simulation of random wave overtopping by a WCSPH model. Appl. Ocean Res. 116, 102888. http://dx.doi.org/10.1016/j.apor.2021.102888.

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. http://dx.doi.org/10.1016/j.coastaleng.2016.07.002.

Bonneton, P., Chazel, F., Lannes, D., Marche, F., Tissier, M., 2011. A splitting approach for the fully nonlinear and weakly dispersive Green–Naghdi model. J. Comput. Phys. 230 (4), 1479–1498. http://dx.doi.org/10.1016/j.jcp.2010.11.015.

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.: Ocean. 104 (C4), 7649–7666. http://dx.doi.org/10.1029/98JC02622.

Borthwick, A.G.L., Ford, M., Weston, B.P., Taylor, P.H., Stansby, P.K., 2006. Solitary wave transformation, breaking and run-up at a beach. Proc. Inst. Civ. Eng. - Marit. Eng. 159 (3), 97–105. http://dx.doi.org/10.1680/maen.2006.159.3.97.

Briganti, R., Dodd, N., 2009a. On the role of shoreline boundary conditions in wave overtopping modelling with non-linear shallow water equations. Coast. Eng. 56 (10), 1061–1067. http://dx.doi.org/10.1016/j.coastaleng.2009.06.011.

Briganti, R., Dodd, N., 2009b. Shoreline motion in nonlinear shallow Water Coastal models. Coast. Eng. 56 (5), 495–505. http://dx.doi.org/10.1016/j.coastaleng.2008.10.008.

Brocchini, M., 2013. A reasoned overview on Boussinesq-type models: The interplay between physics, mathematics and numerics. Proc. R. Soc. A: Math. Phys. Eng. Sci. 469 (2160), 20130496. http://dx.doi.org/10.1098/rspa.2013.0496.

Chassagne, R., Filippini, A.G., Ricchiuto, M., Bonneton, P., 2019. Dispersive and dispersive-like bores in channels with sloping banks. J. Fluid Mech. 870, 595–616. http://dx.doi.org/10.1017/jfm.2019.287.

Chen, G., Noelle, S., 2017. A new hydrostatic reconstruction scheme based on subcell reconstructions. SIAM J. Numer. Anal. 55 (2), 758–784. http://dx.doi.org/10.1137/15M1053074.

Chen, W., Warmink, J.J., van Gent, M.R.A., Hulscher, S.J.M.H., 2021. Numerical modelling of wave overtopping at dikes using OpenFOAM®. Coast. Eng. 166, 103890. http://dx.doi.org/10.1016/j.coastaleng.2021.103890.

Chow, V.T., 1959. Open-Channel Hydraulics. McGraw-Hill, New-York.

Coulaud, G., Teles, M., Benoit, M., 2025. A comparison of eight weakly dispersive Boussinesq-type models for non-breaking long-wave propagation in variable water depth. Coast. Eng. 195, 104645. http://dx.doi.org/10.1016/j.coastaleng.2024.104645.

Dodd, N., 1998. Numerical model of wave run-up, overtopping, and regeneration. J. Waterw. Port Coast. Ocean. Eng. 124 (2), 73–81. http://dx.doi.org/10.1061/(ASCE)0733-950X(1998)124:2(73).

Dong, J., Li, D.F., 2020. A reliable second-order hydrostatic reconstruction for shallow water flows with the friction term and the bed source term. J. Comput. Appl. Math. 376, 112871. http://dx.doi.org/10.1016/j.cam.2020.112871.

Engsig-Karup, A.P., Hesthaven, J.S., Bingham, H.B., Madsen, P.A., 2006. Nodal DG-FEM solution of high-order Boussinesq-type equations. J. Engrg. Math. 56 (3), 351–370. http://dx.doi.org/10.1007/s10665-006-9064-z.

Erduran, K.S., Ilic, S., Kutija, V., 2005. Hybrid finite-volume finite-difference scheme for the solution of Boussinesq equations. Internat. J. Numer. Methods Fluids 49 (11), 1213–1232. http://dx.doi.org/10.1002/fld.1021.

EurOtop, 2007. Wave overtopping of sea defences and related structures: Assessment manual. Pullen, T., Allsop, N.W.H., Bruce, T., Kortenhaus, A., Schüttrumpf, H., van der Meer, J.W., www.overtopping-manual.com.

EurOtop, 2018. Manual on wave overtopping of sea defences and related structures. An overtopping manual largely based on European research, but for worldwide application. van der Meer, J.W., Allsop, N.W.H., Bruce, T., De Rouck, J., Kortenhaus, A., Pullen, T., Schüttrumpf, H., Troch, P., Zanuttigh, B., www.overtopping-manual.com.

Filippini, A.G., Kazolea, M., Ricchiuto, M., 2016. A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up. J. Comput. Phys. 310, 381–417. http://dx.doi.org/10.1016/j.jcp.2016.01.027.

Galaz, J., Kazolea, M., Rousseau, A., 2025. Analysis of linear Boussinesq-type models coupled with static interfaces. http://dx.doi.org/10.48550/arXiv.2503.10300, ArXiv Preprint. arXiv:2503.10300.

Goda, Y., 2009. Derivation of unified wave overtopping formulas for seawalls with smooth, impermeable surfaces based on selected CLASH datasets. Coast. Eng. 56 (4), 385–399. http://dx.doi.org/10.1016/j.coastaleng.2008.09.007.
Gottlieb, S., Shu, C.-W., Tadmor, E., 2001. Strong stability-preserving high-order time discretization methods. SIAM Rev. 43 (1), 89–112. http://dx.doi.org/10.1137/S003614450036757X.

Green, A.E., Naghdi, P.M., 1976. A derivation of equations for wave propagation in water of variable depth. J. Fluid Mech. 78 (2), 237–246. http://dx.doi.org/10.1017/S0022112076002425.

Gruwez, V., Altomare, C., Suzuki, T., Streicher, M., Cappietti, L., Kortenhaus, A., Troch, P., 2020. Validation of RANS modelling for wave interactions with sea dikes on shallow foreshores using a large-scale experimental dataset. J. Mar. Sci. Eng. 8 (9), 650. http://dx.doi.org/10.3390/jmse8090650.

Harten, A., Lax, P.D., van Leer, B., 1983. On upstream differencing and Godunov-Type schemes for hyperbolic conservation laws. SIAM Rev. 25 (1), 35–61. http://dx.doi.org/10.1137/1025002.

Hasselmann, K., Barnett, T.P., Bouws, E., Carlson, H., Cartwright, D.E., Enke, K., Ewing, J.A., Gienapp, A., 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). Ergänzungsheft Zur Dtsch. Hydrogr. Z. Reihe A Nr. 12.

Henderson, C.S., Fiedler, J.W., Merrifield, M.A., Guza, R.T., Young, A.P., 2022. Phase resolving runup and overtopping field validation of SWASH. Coast. Eng. 175, 104128. http://dx.doi.org/10.1016/j.coastaleng.2022.104128.

Higuera, P., Lara, J.L., Losada, I.J., 2014. Three-dimensional interaction of waves and Porous Coastal structures using OpenFoAM®. Part II: Application. Coast. Eng. 83, 259–270. http://dx.doi.org/10.1016/j.coastaleng.2013.09.002.

Hu, K., Mingham, C.G., Causon, D.M., 2000. Numerical simulation of wave overtopping of coastal structures using the non-linear shallow water equations. Coast. Eng. 41 (4), 433–465. http://dx.doi.org/10.1016/S0378-3839(00)00040-5.

Hubbard, M.E., Dodd, N., 2002. A 2D numerical model of wave run-up and overtopping. Coast. Eng. 47 (1), 1–26. http://dx.doi.org/10.1016/S0378-3839(02)00094-7.

Irías Mata, M., Van Gent, M.R., 2023. Numerical modelling of wave overtopping discharges at rubble mound breakwaters using OpenFoAM®. Coast. Eng. 181, 104274. http://dx.doi.org/10.1016/j.coastaleng.2022.104274.

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

Joshi, S.M., Kazolea, M., Ricchiuto, M., 2022. Parameter sensitivity for wave-breaking closures in Boussinesq-Type Models. Water Waves 4 (3), 491–515. http://dx.doi.org/10.1007/s42286-022-00068-2.

Jouy, B., Violeau, D., Ricchiuto, M., Le, M., 2024. One dimensional modelling of Favre waves in channels. Appl. Math. Model. 133, 170–194. http://dx.doi.org/10.1016/j.apm.2024.05.020.

Kalisch, H., Lagona, F., Roeber, V., 2024. Sudden wave flooding on steep rock shores: A clear but hidden danger. Nat. Hazards 120 (3), 3105–3125. http://dx.doi.org/10.1007/s11069-023-06319-w.

Kazolea, M., Delis, A.I., 2013. A well-balanced shock-capturing hybrid finite volume–finite difference numerical scheme for extended 1D Boussinesq models. Appl. Numer. Math. 67, 167–186. http://dx.doi.org/10.1016/j.apnum.2011.07.003.

Kazolea, M., Delis, A.I., Nikolos, I.K., Synolakis, C.E., 2012. An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations. Coast. Eng. 69, 42–66. http://dx.doi.org/10.1016/j.coastaleng.2012.05.008.

Kazolea, M., Filippini, A.G., Ricchiuto, M., 2023. Low dispersion finite volume/element discretization of the enhanced Green–Naghdi equations for wave propagation, breaking and runup on unstructured meshes. Ocean. Model. 182, 102157. http://dx.doi.org/10.1016/j.ocemod.2022.102157.

Kazolea, M., Ricchiuto, M., 2018. On wave breaking for Boussinesq-type models. Ocean. Model. 123, 16–39. http://dx.doi.org/10.1016/j.ocemod.2018.01.003.

Kazolea, M., Ricchiuto, M., 2024. Full nonlinearity in weakly dispersive Boussinesq Models: Luxury or necessity. J. Hydraul. Eng. 150 (1), 04023061. http://dx.doi.org/10.1061/JHEND8.HYENG-13718.

Kennedy, A.B., Chen, Q., Kirby, J.T., Dalrymple, R.A., 2000. Boussinesq modeling of wave transformation, breaking, and runup. I: 1D. J. Waterw. Port Coast. Ocean. Eng. 126 (1), 39–47. http://dx.doi.org/10.1061/(ASCE)0733-950X(2000)126:1(39).

Khakimzyanov, G., Dutykh, D., Fedotova, Z., Gusev, O., 2020. Dispersive Shallow Water Waves: Theory, Modeling, and Numerical Methods. In: Lecture Notes in Geosystems Mathematics and Computing, Springer International Publishing, http://dx.doi.org/10.1007/978-3-030-46267-3.

Kirby, J.T., 2003. Boussinesq models and applications to nearshore wave propagation, surf zone processes and wave-induced currents. In: Lakhan, V.C. (Ed.), Advances in Coastal Modeling. In: Elsevier Oceanography Series, vol. 67, Elsevier, pp. 1–41.

Kobayashi, N., Wurjanto, A., 1989. Wave overtopping on coastal structures. J. Waterw. Port Coast. Ocean. Eng. 115 (2), 235–251. http://dx.doi.org/10.1061/(ASCE)0733-950X(1989)115:2(235).

Lannes, D., 2013. The Water Waves Problem: Mathematical Analysis and Asymptotics. In: Mathematical Surveys and Monographs, Vol. 188, American Mathematical Society, http://dx.doi.org/10.1090/surv/1880.

Lannes, D., Marche, F., 2015. A new class of fully nonlinear and weakly dispersive Green–Naghdi models for efficient 2D simulations. J. Comput. Phys. 282, 238–268. http://dx.doi.org/10.1016/j.jcp.2014.11.016.

Lashley, C.H., van der Meer, J., Bricker, J.D., Altomare, C., Suzuki, T., Hirayama, K., 2021. Formulating wave overtopping at vertical and sloping structures with shallow foreshores using deep-water wave characteristics. J. Waterw. Port Coast. Ocean. Eng. 147 (6), 04021036. http://dx.doi.org/10.1061/(ASCE)WW.1943-5460.0000675.

Lashley, C.H., Zanuttigh, B., Bricker, J.D., van der Meer, J., Altomare, C., Suzuki, T., Roeber, V., Oosterlo, P., 2020. Benchmarking of numerical models for wave overtopping at dikes with shallow mildly sloping foreshores: Accuracy versus speed. Environ. Model. Softw. 130, 104740. http://dx.doi.org/10.1016/j.envsoft.2020.104740.

Li, M., Xu, L., Cheng, Y., 2019. A CDG-FE method for the two-dimensional Green–Naghdi model with the enhanced dispersive property. J. Comput. Phys. 399, 108953. http://dx.doi.org/10.1016/j.jcp.2019.108953.

Losada, I.J., Lara, J.L., Guanche, R., Gonzalez-Ondina, J.M., 2008. Numerical analysis of wave overtopping of rubble mound breakwaters. Coast. Eng. 55 (1), 47–62. http://dx.doi.org/10.1016/j.coastaleng.2007.06.003.

Lteif, R., 2024. An operator-splitting approach with a hybrid finite volume/finite difference scheme for extended Boussinesq models. Appl. Numer. Math. 196, 159–182. http://dx.doi.org/10.1016/j.apnum.2023.10.009.

Lynett, P.J., 2006. Nearshore wave modeling with high-order Boussinesq-Type equations. J. Waterw. Port Coast. Ocean. Eng. 132 (5), 348–357. http://dx.doi.org/10.1061/(ASCE)0733-950X(2006)132:5(348).

Lynett, P.J., Melby, J.A., Kim, D.-H., 2010. An application of Boussinesq modeling to Hurricane wave overtopping and inundation. Ocean Eng. 37 (1), 135–153. http://dx.doi.org/10.1016/j.oceaneng.2009.08.021.

Madsen, P.A., Sørensen, O.R., 1992. A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry. Coast. Eng. 18 (3–4), 183–204. http://dx.doi.org/10.1016/0378-3839(92)90019-Q.

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. http://dx.doi.org/10.1061/(ASCE)WW.1943-5460.0000199.

McCabe, M.V., Stansby, P.K., Apsley, D.D., 2013. Random wave runup and overtopping of a steep sea wall: Shallow-water and Boussinesq modelling with generalised breaking and wall impact algorithms validated against laboratory and field measurements. Coast. Eng. 74, 33–49. http://dx.doi.org/10.1016/j.coastaleng.2012.11.010.

Nwogu, O., 1993. Alternative form of Boussinesq equations for nearshore wave propagation. J. Waterw. Port Coast. Ocean. Eng. 119 (6), 618–638. http://dx.doi.org/10.1061/(ASCE)0733-950X(1993)119:6(618).

Nwogu, O.G., 1996. Numerical prediction of breaking waves and currents with a Boussinesq model. In: Proceedings of 25th International Conference on Coastal Engineering 1996. American Society of Civil Engineers, Orlando, Florida, United States, pp. 4807–4820. http://dx.doi.org/10.1061/9780784402429.374.

Peregrine, D.H., 1967. Long waves on a beach. J. Fluid Mech. 27 (4), 815–827. http://dx.doi.org/10.1017/S0022112067002605.

Pinault, J., Morichon, D., Roeber, V., 2020. Estimation of irregular wave runup on intermediate and reflective beaches using a phase-resolving numerical model. J. Mar. Sci. Eng. 8 (12), 993. http://dx.doi.org/10.3390/jmse8120993.

Pope, S.B., 2000. Turbulent Flows, first ed. Cambridge University Press, http://dx.doi.org/10.1017/CBO9780511840531.

Popinet, S., 2015. A quadtree-adaptive multigrid solver for the Serre–Green–Naghdi equations. J. Comput. Phys. 302, 336–358. http://dx.doi.org/10.1016/j.jcp.2015.09.009.

Pullen, T., Allsop, N., 2003. Use of Numerical Models of Wave Overtopping: A Summary of Current Understanding. Technical Report R&D Interim Guidance Note FD2410/GN1, Department for Environment Food and Rural Affairs / Environment Agency.

Roeber, V., Cheung, K.F., 2012. Boussinesq-type model for energetic breaking waves in fringing reef environments. Coast. Eng. 70, 1–20. http://dx.doi.org/10.1016/j.coastaleng.2012.06.001.

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. http://dx.doi.org/10.1016/j.coastaleng.2015.05.005.

Serre, F., 1953. Contribution à l'étude des Écoulements permanents et variables dans les canaux. Houille Blanche 39 (6), 830–872. http://dx.doi.org/10.1051/lhb/1953058.

Shi, F., Kirby, J.T., Harris, J.C., Geiman, J.D., Grilli, S.T., 2012. A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation. Ocean. Model. 43–44, 36–51. http://dx.doi.org/10.1016/j.ocemod.2011.12.004.

Shiach, J.B., Mingham, C.G., Ingram, D.M., Bruce, T., 2004. The applicability of the shallow water equations for modelling violent wave overtopping. Coast. Eng. 51 (1), 1–15. http://dx.doi.org/10.1016/j.coastaleng.2003.11.001.

Simon, B., Papoutsellis, C.E., Benoit, M., Yates, M.L., 2019. Comparing methods of modeling depth-induced breaking of irregular waves with a fully nonlinear potential flow approach. J. Ocean. Eng. Mar. Energy 5 (4), 365–383. http://dx.doi.org/10.1007/s40722-019-00154-7.

Smit, P., Stelling, G., Roelvink, D.J., Thiel de Vries, J., McCall, R., van Dongeren, A., Zwinkels, C., Jacobs, R., 2010. XBeach: Non-hydrostatic Model. Validation, Verification and Model Description. Technical Report, Delft University of Technology and Deltares, The Netherlands.
Smit, P., Zijlema, M., Stelling, G., 2013. Depth-induced wave breaking in a non-hydrostatic, near-shore wave model. Coast. Eng. 76, 1–16. http://dx.doi.org/10.1016/j.coastaleng.2013.01.008.

Stansby, P.K., 2003. Solitary wave run up and overtopping by a semi-implicit finite-volume shallow-water Boussinesq model. J. Hydraul. Res. 41 (6), 639–647. http://dx.doi.org/10.1080/00221680309506896.

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, 108–123. http://dx.doi.org/10.1016/j.coastaleng.2017.01.009.

Suzuki, T., Altomare, C., Verwaest, T., Trouw, K., Zijlema, M., 2014. Two-dimensional wave overtopping calculation over a dike in shallow foreshore by SWASH. Coast. Eng. Proc. (34), 3. http://dx.doi.org/10.9753/icce.v34.structures.3.

Suzuki, T., Altomare, C., Yasuda, T., Verwaest, T., 2020. Characterization of overtopping waves on sea dikes with gentle and shallow foreshores. J. Mar. Sci. Eng. 8 (10), 752. http://dx.doi.org/10.3390/jmse8100752.

TAW, 2002. Wave Run-up and Wave Overtopping at Dikes. Technical Report, Technical Advisory Committee on Flood Defence, The Netherlands.

Tissier, M., Bonneton, P., Marche, F., Chazel, F., Lannes, D., 2012. A new approach to handle wave breaking in fully non-linear Boussinesq models. Coast. Eng. 67, 54–66. http://dx.doi.org/10.1016/j.coastaleng.2012.04.004.

Tonelli, M., Petti, M., 2009. Hybrid finite volume – finite difference scheme for 2DH improved Boussinesq equations. Coast. Eng. 56 (5), 609–620. http://dx.doi.org/10.1016/j.coastaleng.2009.01.001.

Tonelli, M., Petti, M., 2010. Finite volume scheme for the solution of 2D extended Boussinesq equations in the surf zone. Ocean Eng. 37 (7), 567–582. http://dx.doi.org/10.1016/j.oceaneng.2010.02.004.

Tonelli, M., Petti, M., 2013. Numerical simulation of wave overtopping at coastal dikes and low-crested structures by means of a shock-capturing Boussinesq model. Coast. Eng. 79, 75–88. http://dx.doi.org/10.1016/j.coastaleng.2013.04.007.

Toro, E.F., Müller, L.O., Siviglia, A., 2020. Bounds for wave speeds in the Riemann Problem: Direct theoretical estimates. Comput. & Fluids 209, 104640. http://dx.doi.org/10.1016/j.compfluid.2020.104640.

Tsung, W.-S., Hsiao, S.-C., Lin, T.-C., 2012. Numerical simulation of solitary wave run-up and overtopping using Boussinesq-Type model. J. Hydrodyn. 24 (6), 899–913. http://dx.doi.org/10.1016/S1001-6058(11)60318-1.

Tuan, T.Q., Oumeraci, H., 2010. A numerical model of wave overtopping on seadikes. Coast. Eng. 57 (8), 757–772. http://dx.doi.org/10.1016/j.coastaleng.2010.04.007.

Tuozzo, S., Calabrese, M., Buccino, M., 2024. An overtopping formula for shallow water vertical seawalls by SWASH. Appl. Ocean Res. 148, 104009. http://dx.doi.org/10.1016/j.apor.2024.104009.

Van Gent, M.R.A., 1994. The modelling of wave action on and in coastal structures. Coast. Eng. 22 (3), 311–339. http://dx.doi.org/10.1016/0378-3839(94)90041-8.

Van Leer, B., Nishikawa, H., 2021. Towards the ultimate understanding of MUSCL: Pitfalls in achieving third-order accuracy. J. Comput. Phys. 446, 110640. http://dx.doi.org/10.1016/j.jcp.2021.110640.

Vanneste, D.F., Altomare, C., Suzuki, T., Troch, P., Verwaest, T., 2014. Comparison of numerical models for wave overtopping and impact on a sea wall. In: Proceedings of 33rd International Conference on Coastal Engineering 2014. Vol. 1, American Society of Civil Engineers, Seoul, Korea, p. 5. http://dx.doi.org/10.9753/icce.v34.structures.5.

Varra, G., Pepe, V., Della Morte, R., Cozzolino, L., 2024. A novel efficient and robust treatment of the friction source term in 2D shallow water inundation models. J. Hydrol. 634, 131045. http://dx.doi.org/10.1016/j.jhydrol.2024.131045.

Villefer, A., 2022. Study of wave overtopping processes of coastal protections in complex sea-states. PhD thesis of École Nationale des Ponts et Chaussées (ENPC). Champs-sur-Marne, France, 184 p. https://pastel.hal.science/tel-04076455.

Villefer, A., Benoit, M., Violeau, D., Teles, M., Luneau, C., 2024. Experimental study of wave overtopping of coastal protections for sea state combining swell and wind waves. In: XVIIèmes Journées Nationales Génie Côtier - Génie Civil, Chatou. Vol. 14, Revue Paralia, pp. s02.13-s02.24. http://dx.doi.org/10.5150/revue-paralia.2024.s02.

Whittaker, C.N., Fitzgerald, C.J., Raby, A.C., Taylor, P.H., Borthwick, A.G.L., 2018. Extreme Coastal responses using focused wave groups: Overtopping and horizontal forces exerted on an inclined seawall. Coast. Eng. 140, 292–305. http://dx.doi.org/10.1016/j.coastaleng.2018.08.004.

Williams, H.E., Briganti, R., Pullen, T., 2014. The role of offshore boundary conditions in the uncertainty of numerical prediction of wave overtopping using non-linear shallow water equations. Coast. Eng. 89, 30–44. http://dx.doi.org/10.1016/j.coastaleng.2014.03.003.

Willmott, C.J., Robeson, S.M., Matsuura, K., 2012. A refined index of model performance. Int. J. Climatol. 32 (13), 2088–2094. http://dx.doi.org/10.1002/joc.2419.

Zelt, J.A., Skjelbreia, J.E., 1992. Estimating incident and reflected wave fields using an arbitrary number of wave gauges. In: Proceedings of 23rd International Conference on Coastal Engineering 1992. American Society of Civil Engineers, Venice, Italy, pp. 777–789. http://dx.doi.org/10.1061/9780872629332.058.

Zhang, Y., Kennedy, A.B., Donahue, A.S., Westerink, J.J., Panda, N., Dawson, C., 2014. Rotational surf zone modeling for $O(\mu^4)$ Boussinesq-Green-Naghdi systems. Ocean. Model. 79, 43–53. http://dx.doi.org/10.1016/j.ocemod.2014.04.001.

Zhang, N., Zhang, Q., Wang, K.-H., Zou, G., Jiang, X., Yang, A., Li, Y., 2020. Numerical simulation of wave overtopping on breakwater with an armor layer of accropode using SWASH model. Water 12 (2), 386. http://dx.doi.org/10.3390/w12020386.

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. http://dx.doi.org/10.1016/j.coastaleng.2011.05.015.