# Full nonlinearity in weakly dispersive Boussinesq models: luxury or necessity ?

Maria Kazolea, Mario Ricchiuto

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

## ► To cite this version:

Maria Kazolea, Mario Ricchiuto. Full nonlinearity in weakly dispersive Boussinesq models: luxury or necessity?. Journal of Hydraulic Engineering, 2023, 150 (1), 10.1061/JHEND8.HYENG-13718 . hal-04334853v2

HAL Id: hal-04334853

https://inria.hal.science/hal-04334853v2

Submitted on 11 Dec 2023

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_2.png" /></figure>
# Full nonlinearity in weakly dispersive Boussinesq models: luxury or necessity ?

Maria Kazolea<fnref n="1" />, and Mario Ricchiuto<fnref n="2" />

## Abstract

Boussinesq type (BT) models play a major role in coastal engineering applications. Due to their asymptotic nature, these models have lent themselves to constant study, and improvement based on the inclusion/removal of asymptotically small (or not) terms. In this work we go back to the fundamental question of whether full nonlinear in weakly dispersive BT models is a necessity. We re-consider the tests by [67] to address this issue, as well as a number of more demanding ones. We also consider different families of weakly nonlinear BT models, with different shoaling characteristics, especially when nonlinear waves are involved. Our study allows to point out that for many cases it is quite hard to conclude in favor of one class of model or another. There are a few discriminant cases, which are unfortunately not those mostly used in the literature proposing new models or new numerical methods.

## 1 Introduction

No one can dispute the importance of the Boussinesq type (BT) models in coastal engineering. In the last 50 years significant developments have been achieved both at the level of modeling, as well as in the development of related numerical methods. Boussinesq type models are depth-integrated, phase resolving models for surface wave propagation. Two main parameters characterize these type of models. The first one is the dispersion parameter $\mu$, which is the ratio of the reference water depth $d$ over the characteristic wavelength $L$. The second is nonentity parameter $\epsilon$, defined as the ratio of the characteristic wave amplitude $\alpha_0$ over $d$. Very long waves are characterized by negligible values of $\mu^2$. In this case, the very well known system of shallow water, or Saint-Venant, equations can be recovered. In near shore applications, involving wave transformations and their impact in coastal structures, the most widely used models are the so-called weakly-dispersive ones, which are obtained in the hypothesis that $\mu^2 \ll 1$, and neglect all effects of order $O(\mu^4)$. This paper focuses on this class of Boussinesq models. Note that all these PDE systems involve weak dispersive perturbations of the shallow water equations, and usually involve the inversion of an elliptic operator in the time advancement of the depth averaged velocity.

Within weakly dispersive BT equations, we can further distinguish two important sub-families. Weakly nonlinear Boussinesq models are obtained assuming that $\epsilon \approx O(\mu^2)$, which allows to neglect all effects of order $O(\mu^4, \epsilon\mu^2)$. This ansatz allows to obtain PDE systems which can be written as the nonlinear shallow water equations with a linear dispersive perturbation. The linearity of the dispersive operator translates, in the time independent nature of the elliptic equation. A large number of the most well known models in the coastal engineering community belong to this class: from the simple model by Boussinesq [7] to Peregrine's equations [50], to the dispersion enhanced models by Nwogu [45], Madsen and Sørensen [38] or Beji and Nadaoka [4]. The linear character of the dispersive operator, translates in its time independent nature which allows efficient implementation strategies for its inversion [54]. There exist a number of known codes implementing these models, e.g. BOUSS-2D

<footnote n="1">¹Junior researcher, INRIA, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, 200 Avenue de la Vieille Tour, 33405 Talence cedex, France, maria.kazolea@inria.fr</footnote>

<footnote n="2">²Senior researcher, INRIA, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, 200 Avenue de la Vieille Tour, 33405 Talence cedex, France, mario.ricchiuto@inria.fr</footnote>
45

46

by the US Army [44], DHI's commercial platform MIKE21 [12], the code BOSZ initially developed at U. of Hawaii [57], or the more recent open source platform CELERIS [61, 62].

47

48

49

50

51

52

53

54

The BT models obtained assuming that $\epsilon = O(1)$ are referred to as fully nonlinear. In this case the depth averaged equations involve a strongly nonlinear dispersive operator, making the associated numerical resolution inherently more expensive. Despite of this fact, this choice is without doubt physically more relevant, and provides the well known models of Serre [58], Green Naghdi [21], and enhanced variants discussed for example in [67], [32, 10], or [11]. Model of this type is implemented in many operational codes among which the most well known ones are probably FUNWAVE [59, 14], and Coulwave [37, 15], but similar models have been implemented in other platforms as e.g. the open source projects Basilisk [16], and other research platforms less visible than the above.

55

56

57

58

59

Note that, other parameters may be introduced to characterize depth averaged equations in terms of mild (or not) slope, presence (or not) of surface tension effects etc. For an extensive review on the BT models we refer to [34, 31, 8] and references therein. In this work, we will limit our study to the impact of full/weak nonlinearity.

60

61

62

63

64

65

66

67

68

As already mentioned, fully nonlinear models are more complex than weakly nonlinear ones, and require the inversion of a nonlinear, and thus time dependent, elliptic operator at every time step. The natural question arising in this case is whether this additional complexity and computational cost is necessary. Many think that this issue is fully answered in the landmark paper by Wei and collaborators [67] which concludes that full nonlinearity is absolutely necessary to correctly predict water heights. Yet, the literature is still filled with simulations using weakly nonlinear models, showing accurate predictions of wave heights and phase, for now almost standard cases, in both small amplitude and configurations close to/or in incipient breaking.

69

70

71

72

73

74

75

So in this paper we go back to the question of whether full nonlinearity in weakly dispersive BT models is a necessity. We re-consider the tests by [67], as well as a number of more demanding ones. We also consider different families of weakly nonlinear BT models, with different shoaling characteristics. As we will demonstrate, for many cases it is quite hard to conclude in favor of one class of model or another. There are however, a few discriminant cases in which the height predictions are very negatively affected by the lack of full nonlinearity. These cases are however seldom used in the literature proposing new numerical methods or, more importantly, dispersive PDE models.

76

77

78

79

80

The structure of the paper is the following. We first discuss the BT models considered, their main characteristics, and briefly recall their numerical treatment. Then we review the results on benchmarks of increasing complexity in one and two space dimensions. We end the paper with a conclusive section that tries to provide a few recommendation concerning the choice and benchmarking of the models.

81

## 2 Weakly dispersive Boussinesq PDE models

82

83

84

In this section we briefly describe the Boussinesq systems used in this work, and recall some of their significant properties. As mentioned in the introduction the Boussinesq-type models derived since the 90's are numerous. In this work we use some of the most well known among these.
## 2.1 Weakly nonlinear BT equations

Among the weakly nonlinear models the ones of Madsen and Sørensen [40, 39] and Nwogu's [45] are two typical and widely used models. These two models have quite comparable linear dispersion properties, but different shoaling characteristics. For this reason, we also consider the model Nwogu-Abbot introduced in [18], which has the exact same linear properties of the Nwogu equations but is only equal to the latter within an $O(\epsilon\mu^2, \mu^4)$ correction allowing to write all the dispersive terms using derivatives of the volume flux rather than of the velocity. As shown in [18], this leads to a system which behaves like the Nwogu one for small amplitude waves, but much closer to the one of Madsen and Sørensen in more nonlinear regimes.

All these models can be written in a pseudo balanced law form:

$$ 
\begin{gathered}
h_t + \nabla \cdot \mathbf{q} = \mathcal{D}_h \\
\mathbf{U}_t + \nabla \cdot (\mathbf{u} \otimes \mathbf{q}) + gh\nabla\eta + c_f(h, \mathbf{q})\mathbf{q} = \mathcal{D}_q \\
\text{on } \Omega \times [0, t] \subset \mathbb{R}^2 \times \mathbb{R}^+
\end{gathered}
\quad (1) $$

where $h$ is the total water depth and $\mathbf{q}$ the volume flux $hu$, with $\mathbf{u}$ the depth averaged velocity. The free surface is denoted here by $\eta = h + b$ with $b$ the bathymetry. The friction term $c_f \mathbf{q}$ in this work is approximated with the usual Manning law $c_f = gn_M^2 \|\mathbf{u}\|/h^{4/3}$ with $n_M$ the Manning coefficient. Finally, all the dispersive terms are on the right hand sides of (1). The definition of $\mathcal{D}_h$, of $\mathcal{D}_q$, and of the pseud-momentum variables $\mathbf{U}$ is what allows to distinguish between different models.

**The equations of Madsen and Sørensen** This system can be seen as a dispersion-enhanced variant of the Abbot equations, which are an equivalent of Peregrine's [50] within an $O(\epsilon\mu^2, \mu^4)$ correction allowing to write all the dispersive terms using derivatives of the volume flux rather than of the velocity [18]. In two space dimensions, the equations are obtained from (1) setting $\mathcal{D}_h = 0$, and

$$ \mathcal{D}_{\mathbf{q}} = Bgd^2\nabla\nabla \cdot (d\nabla\eta) \quad (2) $$

where $d = d(x, y) = h_0 - b(x, y)$ is the still water depth. The momentum variable is defined by

$$ \mathbf{U} = \mathbf{q} - \left( B + \frac{1}{3} \right) d^2 \nabla (\nabla \cdot \mathbf{q}) - d \left[ \frac{1}{3} (\nabla d \otimes \nabla \mathbf{q}) + \frac{1}{6} (\nabla d \otimes (\nabla \mathbf{q})^{\perp}) - \frac{1}{6} ((\nabla d)^{\perp} \otimes \nabla \mathbf{q}) \right] \quad (3) $$

The value of $B = 1/15$ is the one for which the phase relation of the system matches a second-order Padé approximation of the exact phase of the Euler equations [40].

**The equations of Nwogu** Introduced in [45] the BT model due to Nwogu can be obtained from Peregrine's by replacing the depth averaged velocity with a velocity at an arbitrary distance from the still water level $z_\theta$. The value $z_\theta = -0.531d$ is the one that optimizes the dispersion relation of the system wrt the one of the Euler equations. The model is defined as in (1) with

$$ \mathcal{D}_h = -\nabla \cdot \left[ \left( \frac{z_\theta^2}{2} - \frac{d^2}{6} \right) d\nabla(\nabla \cdot \mathbf{u}) + \left( z_\theta + \frac{d}{2} \right) d\nabla[\nabla \cdot (d\mathbf{u})] \right], \quad (4) $$

$$ \mathcal{D}_{\mathbf{q}} = -u\mathcal{D}_h + h_t \frac{z_{\theta}^2}{2} \nabla(\nabla \cdot \mathbf{u}) + h_t z_{\theta} \nabla(\nabla \cdot d\mathbf{u}) \quad \text{and} \quad (5) $$

$$ \mathbf{U} = \mathbf{q} + h \left[ \frac{z_{\theta}^2}{2} \nabla(\nabla \cdot \mathbf{u}) + z_{\theta} \nabla(\nabla \cdot d\mathbf{u}) \right]. \quad (6) $$
113 **The equations of Nwogu-Abbott** These equations can be derived starting from the ones of Nwogu
114 and in the aim of find an asymptotically equivalent system in an amplitude-volume flux form which
degenerates to the same linearized equations of [45]. This derivation, detailed in [18], boils down to
115 add order O(εμ², μ⁴) correction terms to pass from derivatives of du to derivatives of q. In this case,
116 the equations are written as in (1) but now we have Dq = 0 and:
117

$$
\mathcal{D}_h = -\nabla \cdot \left[ \left( \frac{z_\theta^2}{2} - \frac{d^2}{6} \right) d\nabla \cdot \left( \frac{\mathbf{q}}{d} \right) + \left( z_\theta + \frac{d}{2} \right) d\nabla \cdot \mathbf{q} \right], \quad (7)
$$

$$
\mathbf{U} = \mathbf{q} + \frac{dz_{\theta}^{2}}{2} \nabla \left( \nabla \cdot \left( \frac{\mathbf{q}}{d} \right) \right) + dz_{\theta} \nabla (\nabla \cdot \mathbf{q}) \quad (8)
$$

118 **2.2 Fully nonlinear BT modeling: the enhanced Green-Naghdi equations**

Maybe one of the most well known fully nonlinear Boussineq models are the equations of Green-Naghdi
[21]. The study of improved formulations of these equations and their numerical approximation has
been the object of intense efforts in the last 15 years, one can refer for example to [6, 20, 42, 46, 35,
52, 11, 13, 36, 47, 34, 51, 63] (and many others, we are certainly forgetting many) and references
therein. Here we consider the enhanced equations in the form used in [19, 28] which reads as in (1)
with $\mathcal{D}_h = 0$, $\mathbf{U} = [h, \mathbf{q}]^T$ and

$$
\mathcal{D}_{\mathbf{q}} = \frac{gh\nabla\eta}{\alpha} + \mathbf{\Psi}
$$

119 where $\mathbf{\Psi}$ is the solution of the elliptic problem

$$
(I + \alpha \mathcal{T}) \mathbf{\Psi} = \mathcal{W} - \mathcal{R}. \tag{9}
$$

We define

$$
\mathcal{W} = -\frac{gh}{\alpha}\nabla\eta \quad \text{and} \quad \mathcal{R} = hQ\left(\frac{\mathbf{q}}{h}\right)
$$

120 where the operators T, Q and R are

$$
T(\cdot) = -\frac{1}{3}\nabla\left(h^3\nabla\cdot\left(\frac{\cdot}{h}\right)\right) - \frac{h^2}{2}\left(\nabla\cdot\left(\frac{\cdot}{h}\right)\right)\nabla b + \frac{1}{2}\nabla\left(h^2\nabla b\cdot\left(\frac{\cdot}{h}\right)\right) + h\left(\nabla b\cdot\left(\frac{\cdot}{h}\right)\right)\nabla b, (10)
$$

$$
\begin{align}
Q(\cdot) &= \frac{2}{3h} \nabla \left( h^3 \left( \nabla(\cdot)_1 \cdot \nabla^{\perp}(\cdot)_2 + (\nabla(\cdot))^2 \right) \right) + h^2 \left( \nabla(\cdot)_1 \cdot \nabla^{\perp}(\cdot)_2 + (\nabla(\cdot))^2 \right) \nabla b \nonumber \\
&\quad + \frac{1}{2h} \nabla \left( h^2 \left( (\cdot) \cdot (\nabla(\cdot)) \nabla b \right) \right) + \left( (\cdot) \cdot (\nabla(\cdot)) \nabla b \right) \nabla b. \tag{11}
\end{align}
$$

121 α is a parameter used to improve the dispersion properties of the model w.r.t those of the full Euler
122 equations [32]. By setting α = 1 we recover the classical Green-Naghdi model. We can write the
123 operator T in such a way as to guarantee its inversion properties and the coersivity of the operator
124 (I + αT), see [18],[30] and references therein for further details.

125 **2.3 A short summary on the properties of the models**

126 We provide here a small recap of the main properties of the models considered. In terms of dispersive
127 propagation, the linear dispersion relation for the nominal values of the models' constants proposed in
128 literature is reported on the leftmost picture on figure 1. In the classical near-shore window of values
129 kd ∈ [0, π], all models provide errors below 2% wrt the Airy theory [60]. We highlight that Nwogu
and Nwogu-Abbot models have the same linearized form so the same phase relation. The usual choice of parameters provides a monotone relation only for the MS and GN model. The Nwogu and eGN equations pass from phase lag to phase advance when going from long to short waves, the crossing taking place within the frequency window considered for the Nwogu model.

An other relevant parameter to consider is the shoaling coefficient ($S$) which relates the rate of change in wave amplitude to the rate of change in water depth. The procedure to produce explicit expressions of the coefficient can be found in [25, 39] among others. The middle picture in figure 1 presents graphically the value of $S$ for each model. The standard values for the MS equations are optimized in terms of linear shoaling, and this is clearly visible in the picture. In the linear case, the eGN shows an over shoaling behavior for short waves, while the Nwogu and GN models under-predict the linear increase in height.

The situation is actually much different in the nonlinear case, as the tests in this paper show. To advance some of the results, we also report in the rightmost picture amplification of second order Stoke-Waves obtained with a second order analysis accounting for $O(\epsilon)$ effects (see e.g. [18] and references therein). The plot shows that amplitudes are systematically under-predicted for the models in amplitude-flux form, and in particular for the MS and NA models considered here. When some nonlinearity is accounted for the Nwogu model, an over-prediction in the amplitudes, up to values $kd \approx 1.7$, is observed. The same holds for all the well known amplitude-velocity models as shown in [18] and in contrast with what the linear shoaling shows. Clearly this behavior is not universal for BT models, which is also a part of the motivation of this work.

<figure><img src="image_3.png" /><figcaption>Figure 1: Linear dispersion with respect to the one of the Airy theory and linear shoaling properties for the models considered</figcaption></figure>

## 2.4 Numerical approximation, wet/dry treatment, and breaking closure

**Numerical framework** The numerical discretization of system (1) is performed on unstructured mesh discretizations of the physical domain, and brings together work from the present authors in the last 10 years. It relies on two main building blocks: a high order node centered finite volume solver on dual cells to evolve in time vector ($h$ U) at all mesh nodes; a nodal elliptic solver used to recover the velocities/fluxes at all mesh nodes starting from the values of U, or evaluate the dispersive forcing solving (9) for the GN equations. The nodal finite volume solver, described in detail in [26, 19, 28] uses classical Roe fluxes, with appropriate modifications to deal with wet/dry fronts and to guarantee the well balanced approximation of the bathymetric source. Third order polynomial reconstructions are used to evaluate the fluxes, based either on the MUSCL approach, or on more recent compact iterative Green-Gauss polynomial recovery [28]. Where necessary, the high order polynomial reconstructions
also allow to evaluate high order derivatives in the dispersive terms on the right hand side of (1). As sections §2.1 and §2.2 show these terms are considerably more cumbersome for the fully nonlinear equations than for the weakly nonlinear models, representing one of the computational overheads of accounting for full nonlinearity. Where necessary, as e.g. for some wave breaking closures, discontinuity capturing is guaranteed by decreasing the degree of the polynomial approximation using either the van Albada limiter, of the high order limiter proposed in [41]. Time integration is achieved through a third order explicit Strong Stability Preserving Runge-Kutta method.

The elliptic systems are discretized with linear finite elements. The coupling with the FV method to express the necessary quantities at nodes is performed via the cell-wise high order polynomial (see [28] for details). Note that for fully nonlinear models these systems involve the depth $h$ and are thus time dependent. They require a new assembly step and inversion at each time step, which is another computational overhead of accounting for full nonlinearity. The appropriate formulation of boundary condition for this step is an open research area, especially on non-straight boundaries. Usually ad-hoc techniques involving sponge layers, and internal wave generators are used in combination with simple treatments of the elliptic component. The interested reader can refer to [33, 24, 43] and references therein for more on this aspect. Note also that we have not considered here less classical variants of the GN model as those proposed in [35] which involve two elliptic solves at each time step, albeit of constant (in time) elliptic operators. The general boundary treatment for the extra auxiliary variable is even less clear. We refer however to the reference for details on this promising approach.

As a final note that we would like to stress that *unless otherwise specified, all results presented in this paper are grid independent*.

**Wet/dry transitions** In wet/dry fronts special care needs to be taken in to account in order to keep the well balancing of the scheme, and avoid the blow up of the velocities. This involves modifications to the states used to evaluate the numerical fluxes, as well as the bathymetry, which are somewhat classical in the approximation of the shallow water equations. We refer the interested reader to [55, 26, 27, 53, 1] and references therein for details.

For dispersive models, an important issue is related to the so-called “dispersive run-up/rundown”. Weakly nonlinear models include in their dispersion operator the still water level $d(x, y)$, which is a somewhat problem dependent variable. Above the level $d = 0$, usually corresponding to the initial position of the shorelines, these terms do not make sense, as $d$ itself changes sign. For this reason, in all weakly nonlinear models dispersion is (smoothly) turned off in initially dry areas. For some cases involving moving waves interacting with shorelines this transition can be seen to affect the results. In contrast, the fully nonlinear GN model involves the actual depth $h$ in the dispersive terms. These terms are in practice numerically limited only close to the current location of the wet/dry front, thus also computed in initially dry regions.

**Wave breaking** Last but not least, an appropriate closure is required to account for the dissipation in breaking waves. There are many approaches to obtain such closure. The interested reader is referred to [27, 30] and references therein for an overview. A wave breaking closure consists of two parts: a detection step; a dissipation mechanism. In this work we use the simple, local, and robust detection criteria proposed in [27, 19]. A cell is flagged as breaking when:

$$ |\eta_t| \ge \gamma \sqrt{gH} \quad \begin{array}{l} \text{surface variation criterion} \\ \text{or/and} \end{array} $$

$$ \|\nabla \eta\| \ge \tan(\phi_c) \quad \text{local slope criterion} $$
Usually, the flagging is extended according to some physical and numerical criteria discussed e.g. in [27, 19, 2]. The parameters $\gamma$ and $\phi_c$ depend on the type of breaking regime. Some values are provided in the last references. In this work we use $\gamma = 0.6$ and $\phi_c = 30^\circ$ unless otherwise stated.

Concerning the dissipation mechanism the results presented in this paper have been obtained with two approaches discussed in [30]: the hybrid closure in which dispersion is turned off in breakers which eventually become shallow water (hydrostatic) bores; an eddy viscosity closure involving the solution of a balance equation for the turbulent kinetic energy (TKE). As shown in [30] both approaches allow to describe correctly wave transformation and breaking in large scales. But when we want to obtain grid converged solutions, as in most of the 1d results presented here, only the TKE closure can be used, the hybrid closure introducing unbounded physical oscillations at the transition region when the mesh is refined (see also [2]).

# 3 One dimensional benchmarks

## 3.1 Shoaling of a solitary wave

This is a benchmark used in two landmark papers [22, 67]. In the first the authors discuss comparisons between experimental data and a fully nonlinear/dispersive potential model. The second paper uses this benchmark to compare fully and weakly nonlinear variants of a Boussinesq model, based on Nwogu's equations [45], with the experimental and potential results. The benchmark setup is quite simple. The computational domain consists of a ramp of constant slope, whose toe is located at $x = 0m$. The depth on the left of the ramp is of $d = 0.44m$. A solitary wave of amplitude $A = \epsilon d$ propagates and shoals on the ramp. Here, differently from [67], we compare weakly Boussinesq models with different linear and non-linear shoaling characteristics, as well as the fully nonlinear Green-Naghdi one. For all the models, the usual analytical $\sec h^2$ solitary wave shape is used to initialize the computation.

<figure><img src="image_4.png" /><figcaption>Figure 2: Free surface elevation on the wave gauges 1, 3, 5, 7, 9 for GN and Nwogu (left), and for MS, NA, and GN (right).</figcaption></figure>

We start from the case $\epsilon = 0.2$, with a slope $s = 1:35$, and report the wave height measured at different wave gauges along the slope in figure 2. On the left we compare the experimental data to the fully nonlinear Green-Naghdi results as well as to those obtained with the weakly nonlinear model of Nwogu. On the right we compare the data with the results of the fully nonlinear equations and with
those obtained with the equations of Madsen and Sørensen, and Nwogu-Abbot. The result on the left is essentially the one discussed in [67]. We can clearly see that the weakly nonlinear model provides a large overestimation of the wave heights. The right figure however is much different, the weakly nonlinear MS and NA models providing a small underestimation of the elevations. Similar results have been obtained in [18] using as initial data for each model its own semi-analytical soliton shape.

227

<figure><img src="image_5.png" /><figcaption>Figure 3: Relative wave height computed with all the models on a 1:15 slope.</figcaption></figure>

Next, following the test performed in [67], we let different solitary waves of nonlinearity $\epsilon = 0.3$, $0.45$, $0.6$ shoal on a slope $s = 1:15$, and plot the envelopes of the wave peaks along the shoal. The results are plotted on figure 3 (soliton nonlinearity increasing from left to right curves). The ted circles represent the experimental breaking point of each wave. Our results confirm only partly the ones by [67]. Indeed, for the equations of Nwogu, the breaking point is “missed” by the model, which predict a substantial over-shoaling of the waves. The more nonlinear the wave, the stronger the over prediction. However, the weakly nonlinear equations written in amplitude flux form provide results much closer to the fully nonlinear model, and a reasonable approximation of the wave peak at the breaking point. This is somewhat in contrast with the dogmatic view that all the weakly nonlinear models over predict wave heights, and that provide wrong values of the free surface elevation. To complete the picture, in figure 4 we report the solitary shape for the case $\epsilon = 0.3$ at breaking point. We can see that only the Nwogu model provides the nonphysical peak also obtained in [67], the other models giving a wave with similar topology.

228

229

230

231

232

233

234

235

236

237

238

239

240

<figure><img src="image_6.png" /><figcaption>Figure 4: Free surface for the different models for $\epsilon = 0.3$, $s = 1:15$ at breaking point $t = 11.32$.</figcaption></figure>

241

As a final remark, note that in [18] it has thoroughly shown that weakly nonlinear models written
starting from the Peregrine equations, thus with linear dispersive effects expressed as space-time derivatives of the velocity, have very close nonlinear behavior. Indeed, results similar to those of Nwogu's model can be obtained with other BT systems, e.g. Peregrine's [50], Beji-Nadaoka [4], etc. (see [18] and references therein).

To conclude, despite showing important differences in the nonlinear behavior of the models, this benchmark does not settle the question on whether full nonlinearity is necessary, the models in amplitude-flux form following relatively closely the fully nonlinear model up to the breaking point as shown clearly from Figure 3.

## 3.2 Undular bore propagation

The propagation of undular bores is also a classical problem, studied for example in [49, 67], and for which well known experiments also exist as e.g. [17, 65]. The initial setup is essentially the same as in [67]. We consider a computational domain of 300 m with wall conditions at both domain ends. A transition between uniform flow and still water is initially centered at $x = 0$ m, with velocity and surface elevation given by

$$u = \frac{1}{2} u_0 \left(1 - \tanh\left(\frac{x}{a}\right)\right), \quad \eta = u + \frac{1}{4} u^2$$

where $u_0$ is the velocity of uniform flow from the left boundary and $a$ set to 5 m as in [67]. We consider two initial states with jumps in the initial surface elevation $\eta_0 = 0.1$ and $0.3$, which also correspond to nonlinearity values $\epsilon = 0.1$ and $\epsilon = 0.3$, and to Froude numbers $Fr = 1.07$ and $1.22$ respectively. Both waves evolve in undulating bores propagating through the channel. The amplitude as well as the wavelength of the oscillations depend on the initial data.

<figure><img src="image_7.png" /><figcaption>Figure 5: Evolution of an undular bore. Up Left $\epsilon = 0.1$ Up right $\epsilon = 0.3$.</figcaption></figure>

As in [67] we plot the bore profiles at different (dimensionless) times $t' = 50, 60, \text{ and } 70$. The results are shown on figure 5. As one could expect, for the least nonlinear case the models provide nearly identical results, with only slight over amplification of the oscillations by the Nwogu model, and even less visible under amplifications for the MS and NA models. As the nonlinearity grows these differences become more visible. The wave heights predicted by Nwogu's model are substantially larger compared to the GN solution. Note that the latter has been shown to provide results almost identical to a fully nonlinear potential model [9]. On the contrary, both the MS and NA models underestimate the amplitude of the oscillation. Note that this difference in behavior is somewhat in line with the second order analysis done in [18] showing that amplitudes of Stokes-Waves are systematically under-predicted by BT models in amplitude-flux, while models in amplitude-velocity
272 underestimate very short waves, and over estimate longer ones (see also figure 1).

273

<figure><img src="image_8.png" /><figcaption>Figure 6: Undular bore propagation. +: Data of Treske and Favre, o: Numerical simulations</figcaption></figure>

274 To have a more quantitative view of these differences, we look at the dependence of the amplitudes of the secondary waves on the Froude number. Some experimental data are also available for this from Favre [17] and Treske [65]. So simulations are run by considering different values of $u_0$ corresponding to bore Froude numbers in between 1 and 1.3. Amplitudes are computed as the difference between the heights of the first peak and the first trough. On figure 6 we report the plots of both the amplitude, and the peak/trough elevations $A_{max}$ and $A_{min}$. The amplitude plot shows a growth with Froude which looks perfectly linear. For a Froude number below roughly 1.15 all the models provide similar predictions, all close to the experimental values. For more nonlinear situations, we see three family of models: the weakly nonlinear model of Nwogu over-predicting the amplitude; the models of MS and NA under-predicting the amplitudes; the fully nonlinear GN model in between the previous two, and somewhat closer to the experiments. It is worth mentioning that as the Froude increases the difference of the result for Nwogu's equations wrt the data (and to the GN results) is at least twice as much as that of the other weakly nonlinear models which remain in a range of amplitudes reasonably close to the data. From the right figure, we can see that the elevation of the first trough does not change much with the Froude number, so the behavior observed is really associated to the dynamics of the first peak which grows much faster with models in amplitude-velocity form, as Nwogu's, as remarked at the end of the previous section (see also [18]).

As the previous one, this benchmark does show the impact of the nonlinear behavior of the models. However, looking at the amplitudes on figure 6, the improvement of using the fully nonlinear GN model is not striking, compared to the MS and NA ones.

## 3.3 Monochromatic waves on a submerged bar

295 We consider not the laboratory experiments performed by Beji and Battjes [3] on monochromatic waves over a submerged bar (cf. sketch in figure 7). According to [3] the experimental set up was designed to investigate both frequency dispersion characteristics of monochromatic waves, and non-linear interactions. This is possible the most popular set benchmarks in the literature on numerical
models for dispersive free surface waves, widely used for the validation of new BT models as well as of their discretization. We refer to the original reference for the setup, or to one of the numerous papers using this test, see e.g. [26, 54, 27, 18, 2] by the present authors and references therein.

<figure><img src="image_9.png" /><figcaption>Figure 7: Monochromatic waves on a submerged bar, benchmark sketch</figcaption></figure>

We consider here two cases. The first is a mildly dispersive case with incoming wave is $A = 0.01m$, and the nonlinearity and dispersion parameters $\epsilon = 0.025$ and $\mu = 0.11$. The second, is a mildly nonlinear case with incoming waves of amplitude $A = 0.054m$ and nonlinearity and dispersion equal to $\epsilon = 0.125$ and $\mu = 0.083$ respectively. In the second case, the breaking point is reached and waves breaking is observed above the plateau of the submerged bar. Wave breaking is modeled for this case with the hybrid approach discussed in section §2.4 for all the BT models considered.

The first case involving incoming waves with very little nonlinearity, we expect all models to behave similarly. This is absolutely the case, as can be seen in the pictures of the left column of figure 8. The plots report a few periods of the water elevation time series in the gauge (a) right after the plateau (cf figure 7). The results of the GN, MS and Nwogu model are compared to the experimental data. This gauge allows to include effects of both the increased nonlinearity acquired at the end of the shoaling phase, and some dispersive effects during the propagation along the plateau. The numerical results obtained with the enhanced Green-Naghdi model (taken from [19]), and those obtained with the MS and Nwogu models (taken from [18]) are extremely close to one another.

For the mildly nonlinear case, with wave breaking, one might expect to see a substantial difference between different weakly nonlinear models, and wrt the fully nonlinear one. The right column on figure 8 reports the results for this case in terms of water elevation time series in gauge (b) of figure 8, in correspondence of which incipient breaking is observed. The results obtained with the enhanced Green-Naghdi model (taken from [19]), those obtained with the MS model (taken from [2]), and those provided by Nwogu's equations (taken from [27]) using the exact same wave breaking detection, and hybrid breaking closure, and the same mesh size of 4 cm are essentially the same.

Despite being one of the most commonly benchmarks, and while possibly relevant to assess other aspects, this test seems to be little sensitive to the amount of linearity accounted for in the BT model, unless of course the breaking closure is removed to reach nonphysical regimes as done in [18].

## 3.4 Propagation, breaking, and over-topping of a reef

We consider now a more complex case, involving highly nonlinear waves. It consists of the propagation, breaking and reef over-topping of a soliton, investigated numerically and experimentally in [56], and numerically in several other more recent works (see [64, 57, 30, 29] and references therein). The geometry of the problem is sketched in the left picture on figure 9: it consists of a reef with a constant fore slope, hiding on the lee side a shallow flat with a wall on its right end. Water height distributions at time series in 14 wave gauges have been measured in the flume experiments of [56].
<figure><img src="image_10.png" /><figcaption>Figure 8: Left column: Free surface elevation for GN, MS, and Nwogu's models for the case where initial wave's amplitude is A=0.01m (non breaking case). Right column: Free surface elevation for GN MS and Nwogu's models when A=0.054m (breaking case)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: Solitary wave over-topping a reef: geometry and gauges (left); initial propagation (right).</figcaption></figure>

A solitary wave of offshore nonlinearity $\epsilon = 0.3$ propagates on the left of the reef, as shown on the right on figure 9. Simulations with all the BT models of sections §2.1 and §2.2 are discussed. Breaking is modeled using the Turbulent Kinetic Energy (TKE) closure (see section §2.4 and [30]). We visualize the results via a few snapshots of the free surface distribution and water elevation series in 4 gauges denoted by red circles in the left picture on figure 9.

During propagation there is hardly any difference in the results as seen in the right picture on figure 9. Minor differences are observed at the end of the shoaling phase (top left on figure 10), and during over-topping (top right on figure 10) with the MS and NA models giving almost identical height distributions, and Nwogu's model giving the most strongly damped wave, possibly as a result of the slight over-shoaling observed in the top left picture triggering wave breaking a little earlier. After the over-topping, we see the generation of a very nonlinear bore propagating in the shallow area, mid-left picture on figure 10. Here again we can see some differences in the nature of the bore (some dispersion still appearing in the MS and GN results), or in its speed and magnitude (faster and stronger
<figure><img src="image_12.png" /><figcaption>Figure 10: Solitary wave over-topping a reef: free surface snapshots.</figcaption></figure>

for Nwogu's model). In the next snapshot (mid-right picture on figure 10) we see the propagation of the wall reflected bore for which all models are again in phase. Perhaps a main feature to remark is the fact that the first backwash on the front side is stronger with Nwogu's model. After the reflected bore overtops the reef (bottom-left picture on figure 10) we can see a weak dispersive jump on the front side of the reef, and a secondary reflected more in the shallow area. On the shallow side, the results of Nwogu's model now lag phase compared to the others. This is even more clear in the wall reflection of the bore in the bottom-right picture on figure 10. The same figure shows the second backwash on the front side of the reef, for which all models are very close to one another.

Going to the wave gauge signals, we can see that only minor differences between the GN, MS, and NA models. These are mostly small details of the undular bore generated by the second over-topping, and the first bore in the shallow area which shows a lot more dispersion with the MS model. Compared to the data (in red) one can however hardly judge one better than another. Nwogu's results are not far off but show high amplitude oscillations in the reflected wave from the first backwash, visible at dimensionless times ~ 45 – 50 in the first two wave gauges, and from the first over-topping, dimensionless times ~ 55 – 60 in the first two wave gauges. The phase lag of the secondary bores in the shallow area is also visible in the lagging arrival of these waves in the third and fourth gauge, at dimensionless times ~ 55 and ~ 68 approximately.

Despite its nonlinearity, and relevance in mimicking the inundation process, we cannot claim this test to allows to evaluate between the benefits of using fully or weakly nonlinear BT models. The differences observed, when not minor, are mostly originating in very shallow areas. This suggests that they may be more closely related to the parametrisation of wave breaking closure and of the wet/dry treatment, than to the choice of the propagation model. The interested reader car refer to e.g. [23] for a study on the sensitivity of breaking closures, or to [29] for a thorough study of this case using the hybrid closure. So far, we do seem to see some sort of trend, suggesting that weakly-nonlinear amplitude-flux models such as the MS and NA provide consistently results close to the fully nonlinear
<figure><img src="image_13.png" /><figcaption>Figure 11: Surface elevation on the wave gauges WG4-WG8-WG10-WG14</figcaption></figure>

375 GN model up to the breaking point. This will be contradicted in the following benchmarks.

376 **3.5 Extreme run-up on a vertical barrier**

377 Introduced in [66], this benchmark is relevant for wave structure applications, as e.g. the design of
378 coastal structures like sea walls and breakwaters. The core of the problem is to correctly predict the
379 wall run-up which can exceed six times their initial wave amplitude due to the interaction of steepening
380 and dispersion processes. The setup involves a wave tank with flat bottom of depth $d = 10$ m and
381 length $2L$. The actual physical domain is form 0 to $L$ while from $L$ to $2L$ a mirror image of the
382 initial wave is placed (cf figure 12). This symmetric arrangement allows to free the analysis of the
383 dependence of numerical wall boundary fluxes at the point $x = L$. In [66] the problem is studied using
384 the free surface Euler equations, and the classical (non enhanced) Green-Naghdi equations. Here we
385 perform the study with all the BT models discussed in section §2. Both the effects of changing the
386 nonlinearity and dispersion parameters are investigated.

<figure><img src="image_14.png" /><figcaption>Figure 12: Initial conditions and set up of the physical and computational domains</figcaption></figure>
The initial conditions are given by a linear wave multiplied by an envelope function i.e

$$ \eta(x) = W(x)\alpha_0 \sin k_0(x - x_0), \qquad (12) $$

$$ u(x) = W(x)\alpha_0 \frac{\omega_0}{dk_0} \sin k_0(x - x_0) \qquad (13) $$

where

$$ W(x) = \frac{1}{2} \tanh \left( \frac{x - x_0}{\delta} \right) - \frac{1}{2} \tanh \left( \frac{x - x_0 - N_w \lambda_0}{\delta} \right). \qquad (14) $$

We have denoted by $\lambda_0 = \frac{2\pi}{k}$ is the initial wavelength, and with $x_0$ and $\alpha_0$ the initial position and amplitude of the wave, respectively. The main parameters that characterize the initial conditions are the frequency $\omega_0$, which depends on the dispersion relation of the propagation model, the amplitude $\alpha_0$ and the initial wave's location $x_0$ or in other words, the wave's distance from the wall. The remaining parameters $N_w$ and $\delta$ denote respectively the number of waves in the physical domain, set to 3 like in [66], and a thickness parameter affecting the wave dynamics set to $0.2\lambda_0$ as in the original reference.

We first we examine the effect of the wavelength on the run-up. For a fixed initial wave amplitude $\alpha_0 = 0.05d$ we span the range $\lambda_0/d \in [30, 625]$, roughly corresponding to $\mu \in [0.0016, 0.033]$. In all simulations we set $L = 6\lambda_0$. The maximum run-up as a function of the initial wavelength is plotted on figure 13 for all the BT models considered. The reference solution of the Euler model is taken from [66]. For this case the situation is clearly reversed compared to the previous ones. The fully nonlinear models are quite close to the reference as expected, the enhanced model providing the best results. In this case, the amplitude-flux models, namely the MS and NA ones, provide a considerable underestimation of the run-up (up to 15%) for a wide range of wavelengths. The classical Nwogu equations instead perform quite well, as well as the enhanced Green-Naghdi model, for the longest waves, reversing the trend set up from the previous benchmarks. For short waves, the Nwogu model still shows a noticeable over-prediction of the run-up, compared to the Euler equations. Only the fully nonlinear model provides a correct prediction of these waves, the weakly nonlinear equations in amplitude-flux form still underestimating the run-up.

Next, as in [66], we examine the effect of the normalized initial amplitude, for a fixed wavelength set to $\lambda_0 = 125d$, roughly corresponding to a value $\mu \approx 8 \times 10^{-3}$. We consider nonlinearities in the range $\epsilon = \alpha_0/d \in [0.01, 0.15]$. The resulting maximum run-up is plotted for all the models in figure 14 (reference solution of the Euler model taken from [66]). For small nonlinearity $\epsilon \le 0.02$ all the models provide almost superimposed results, with a small underestimation of the run-up height, compared to the Euler equations. As in [66] we observe a sharp transition taking place in the range $0.02 < \epsilon < 0.05$. This is the region in which the run up regime changes from linear to fully nonlinear. Above $\epsilon = 0.05$ some sort of equilibrium is attained between all the processes, with the Euler equations providing run-up values all close to about $R_{max} = 6d$. Both the slopes in the second phase, as well as the run-up values and slopes in the third one show significant differences between the models. For both the MS and NA models the run-up increase in the second phase is slower, and reaches a value roughly 10% lower than all other models. In the last phase, we see a decrease in run-up, leading for the most nonlinear waves to an underestimation of the run-up of the order of 25% compared to the Euler equations. All the remaining models, including the weakly nonlinear Nwogu equations, provide a very close transition in the second phase, all reaching values of $R_{max}/d$ in between 5.5 and 5.75 for $\epsilon = 0.05$, and with values of for higher nonlinearities considerably closer to the Euler model than the MS and NA equations.
<figure><img src="image_15.png" /><figcaption>Figure 13: Normalized maximum run-up vs non-dimensional initial wavelength for different BT models and fixed $\alpha_0 = 0.05d$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: Normalized maximum run-up vs non-dimensional initial amplitude for different BT models and fixed $\lambda_0 = 125d$.</figcaption></figure>

Note that, for the range of $kd$ used in this test, the difference between the weakly nonlinear models is somewhat in line with the second order analysis recalled in section §2.3 (cf. rightmost picture in figure 1), and already invoked to interpret the results observed for the propagation of undular bores. Also it is important to stress that for all the waves considered the hypothesis $\epsilon = O(\mu^2)$ underlying the weakly nonlinear models is violated systematically, the nonlinearity exceeding this limit by at least one order of magnitude. Despite of this, we are fare from breaking conditions (we experimentally confirmed that the breaking sensor is far from breaking onset). So all the weakly nonlinear propagation models are used well outside their initial ansatz.

It is interesting to see in figure 15 the comparison of the time series of the run-up at the wall for the different BT models. We present the time series for two initial wave amplitudes, $\alpha_0 = 0.1d$ and
<figure><img src="image_17.png" /><figcaption>Figure 15: Run-up at the wall for two initial wave amplitudes $\alpha_0 = 0.1d$ (left) and $\alpha_0 = 0.15d$ (right) and fixed $\lambda_0 = 125d$.</figcaption></figure>

$\alpha_0 = 0.15d$. We focus on the first dispersive wave, were $16 < \omega_0 t < 20$ in order to best appreciate the differences. All the weakly nonlinear models written in an amplitude-flux form (MS and NA) produce results with a slight phase lag and, as seen before, an under-prediction of the amplitudes of the trailing waves with respect to the ones of the eGN equations. The later are really close to the Euler equations, see [66], while Nwogu's model shows a phase advance which grows as the nonlinearity of the initial wave grows.

## 3.6 Long distance propagation, shoaling and run-up of a rectangular wave

This test case was first proposed as a benchmark for tsunami-type surface wave simulation codes in [5] as a part of the French TANDEM project. The main idea was to compare the codes based on different mathematical models, from Euler to BT and shallow water equations, on the dynamics of wave train propagation and run up on the vertical walls of the domain. It involves the formation of an undular bore which propagates over a large distance, shoals on a near coastal bathymetry and finally runs up on a vertical wall. The computational domain consists of a numerical wave tank of 30km length and 50m depth. The sloping bed of 1:125 is placed at the right part of the domain where $x \in [25, 29km]$ followed by flat shelf where the water depth is constant and equals 18m. The initial conditions of the free surface elevation are:

$$ \eta(x) = \begin{cases} 5\text{km}, & \text{if } 0 \le x < 1900\text{m} \\ 2.5(1 - \tanh(229.756 \frac{x}{2000} - 1)), & \text{if } 1900\text{m} \le x \le 2100\text{m} \\ 0, & \text{if } 2100\text{m} < x \le 30000\text{m} \end{cases}. \quad (15) $$

Wall boundary conditions are imposed at both ends of the domain. Figures 16-17 show the first hour of life of the wave, including its propagation, shoaling, and reflection on the right and left wall. In the large scale the numerical solution provided by all models looks almost identical, at least up to the first reflection. Some differences in phase and secondary wave structure can be seen after the first reflection (bottom-right picture on figure 16. These differences are somewhat compensated during propagation toward the left wall, and the position and structure of the wave up to the 5th or 6th peak is quite similar for all models before and after the second reflection (middle- and bottom- right pictures on figure 17),
<figure><img src="image_18.png" /><figcaption>Figure 16: Propagation, shoaling, and wall reflection of the undular bore. A zoom on the wave trains is presented in the right column.</figcaption></figure>

aside perhaps from a larger phase advance and amplitude over-prediction for the Nwogu model. A closer look of the run-up of the waves on the vertical right wall is provided in figure 18. This plot allows to better visualize the differences between all models. We see again a clear under-estimation of the maximum run-up from the amplitude-flux weakly nonlinear models by 15%-20% compared to the GN equations. Note that the reference value of the maximum run-up [5] is of about 24.5 m, which means that the GN model already provides and underestimation of roughly 18.4%, and the MS and NA models have errors of 30%-35% wrt this value. These are much worse numbers than for the previous case. We suspect that the implementation of the wall boundary conditions may affect these values and we are investigating this point. The message is however that, as for the previous case, weakly nonlinear BT models which seem to perform well in other quite nonlinear cases, under estimate substantially the run-up on vertical walls. The tendency to over-predict amplitudes somehow leads to better results for the Nwogu model for this case as for the previous.

# 4 Two dimensional benchmarks

## 4.1 Bore propagation in a trapezoidal channel

We consider now a well known two dimensional variant of the bore propagation problem of section §3.2. The problem is similar, but involves topographic variations across the section which has a trapezoidal shape. As illustrated by the experimental work Treske [65], and by the numerical and theoretical study of [9], the bore dynamics are much more complex in this case. More precisely Treske already showed that the secondary wave field has strong variations in crest height across the
<figure><img src="image_19.png" /><figcaption>Figure 17: Reflected wave propagation, and reflection on the left hand wall. A zoom on the wave trains is presented in the right column.</figcaption></figure>

channel section, and may even have variations in frequency content, depending on the Froude number, a transition taking place around a critical value $Fr_t = 1.15$. According to [9] for Froude above $Fr_t$ the flow is truly multi-dimensional. The waves along the centerline of the channel are similar to the Favre waves and have higher amplitudes, while much longer and less steep waves are observed on the banks. Below $Fr_t$ the flow is quasi-one dimensional with longer waves of higher amplitudes on the banks of the channel. In this regime the same wavelength is observed all across the section. This benchmark, featured quite nonlinear waves, with Froude numbers approaching, and locally exceeding, the breaking limit.

We performed simulations for a range of Froude numbers from 1.05 to 1.25 using the GN, Nwogu and Nwogu-Abbot models. This range of numbers contains the transition, and remains below the transition to a fully breaking bore along the whole channel section. The computational domain of 27566 nodes, consist of a channel with $(x, y) \in [0, 25] \times [0, 1.5]$ in which the half of the channel is placed. More on the set up of the problem can be found in [9]. Free surface visualizations for Froude 1.05, 1.17 and 1.25 are reported in figures 19, 20 and 21 respectively. We can visually see the appearance of more oscillations along the same length of channel shown, and the presence of stronger cross-sectional waves for higher Froude. To confirm these observations, the distribution of the free surface elevation along the banks and centerline are reported in figure 22 for the three models.

Following [9] we post processed the results in order to provide a more quantitative assessment. In particular, we have estimated the height of the first peak and trough, and the wavelength. The results are compared in figure 23 with the experimental data of Treske [65]. The main difference observed among the models is related to the wave amplitudes after the transition, i.e for Froude numbers grater that 1.15. This is roughly the limit after which wave breaking starts taking place close to the banks.
<figure><img src="image_20.png" /><figcaption>Figure 18: Run up of the wave on the right end. A zoom of the peak is presented in the right column.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 19: 2 dimensional free surface elevation for GN (left) Nwogu (centered) Nwogu-Abott (right) at $t \approx 13.5$ sec. The Froude number is 1.05</figcaption></figure>

The computed wavelengths before and after the transition are almost identical for all the models. The only difference observed is for the model of Nwogu in which the transition seems to take place for slightly lower Froude. This may be related to the over shoal of the wave on the banks, leading to earlier breaking and possibly earlier transition.

Despite the strong nonlinearities involved, and Froude numbers close to (or higher than) the breaking point, the amount of nonlinearity embedded in the BT models doesn't seem to affect their capability to correctly predict both amplitudes, and wavelengths.

## 4.2 The seaside test case

The last test case presented is the laboratory benchmark test case for tsunami inundation through urban waterfront. The physical model was a 1:50 scale idealization of the town Seaside, Oregon, designed to observe the complex tsunami flow around the macro-roughness such as buildings idealized as impermeable, rectangular blocks. The laboratory benchmark took place in the Oregon State University. It has been served as a blind benchmark test case for the NTHMP Mapping and Modeling Benchmarking Workshop: Tsunami Currents. All the experimental data are available<fnref n="1" />. Numerical predictions using fully nonlinear using Boussinesq type equations can be found in [48, 28]. The numerical set up of

<footnote n="1">¹http://coastal.usc.edu/currents_workshop/index.html</footnote>
<figure><img src="image_22.png" /><figcaption>Figure 20: 2 dimensional free surface elevation for GN (left) Nwogu (centered) Nwogu-Abott (right) at $t \approx 13.5$ sec. The Froude number is 1.17</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 21: 2 dimensional free surface elevation for GN (left) Nwogu (centered) Nwogu-Abott (right) at $t \approx 13.5$ sec. The Froude number is 1.25</figcaption></figure>

the test case is as follows: the rectangular domain of 42m in $x$ and 14.5m in $y$ is discretized using an unstructured triangulation with 172854 nodes, refined in the region where the buildings are placed. The maxim side of the triangles is $h_{ref} \approx 0.7$m and the minimum one $h_{ref} \approx 0.03$m. The individual structures and buildings are approximated as steep-sided topography and a manning coefficient of 0.01 is used. The available experimental data involve the topography, the forcing signal and the time series of the waves measured at the wave gauges. The forcing signal is entering the computational domain through the left boundary and wall boundary conditions are placed elsewhere. The computational domain is shown in 24. The 31 wave gauges, measure the free surface and the velocity and are depicted in figure 24. The incoming solitary wave shoals and breaks close to the seawall and then floods an urban area which contains buildings of different height. A detailed description of behavior of the free surface elevation can be found in [48, 28]. Figure 25 shows snapshots of the free surface elevation, for each model, at $t \approx 30$sec during the flooding of the urban area.

Figures 26-29 show the free surface elevation for the 31 wave gauges. The experimental data are compared with the numerical results of the equations of Nwogu, Nwogu-Abbot and Green-Naghdi. We observe that in almost all the wave gauges the numerical results are very close to the experimental data and that, especially for the first wave are almost indistinguishable. A difference in the arrival time of the secondary wave at roughly 60sec, is observed in the wave gauges located in front of the buildings.

Overall, when comparing the free surface visualizations, and the time histories at gauges, for
<figure><img src="image_24.png" /><figcaption>Figure 22: Free surface profile for different Froude numbers. Continuous line: axis. Dashed line: banks.</figcaption></figure>

this challenging case, we cannot see extreme differences in the main features of the flow. One major
difference we see, which also affects some of the previous tests, is in very shallow areas, and especially
in initially dry flooding areas, where all weakly nonlinear models never include dispersion. This
explains in our opinion the details of the wave gauge histories for this tests which somewhat favor
the fully nonlinear model, especially maybe for Locations A and B. But the differences are in small
details, and also some of these results maybe could be improved by “tuning” the roughness.

# 5 Discussion and conclusions

In this paper we have reviewed a certain number of benchmarks for dispersive wave propagation
models. It is not simple to draw clear conclusions. One disappointing fact is perhaps that, when
considering both amplitude-flux and amplitude-velocity weakly-nonlinear models, many very well
known and common test cases are not adequate to reveal differences between fully nonlinear and
weakly nonlinear equations. This is the case even for quite nonlinear problems in which the breaking
point is reached and exceeded.

One exception is perhaps Favre waves. Indeed bore propagation on flat topography shows some
advantage of having full nonlinearity in a short window of high Froude numbers close to the breaking
limit. The relatively more recent benchmarks proposed in [66] and [5] allow to highlight more clearly
the need of full nonlinearity at least when considering wave run-up on vertical walls, already for mild
offshore nonlinearity, and as long as breaking conditions are not reached.

These are however relatively academic tests. In complex physical benchmarks with bathymetric
variations, wave breaking, and wet/dry transitions, both in one and two space dimensions it seems
much harder to draw clear conclusions. For models which are well known to over-estimate wave
amplitudes (as Nwogu's model) we do observe some early breaking affecting the overall wave dynamics.
This is however not the case with all weakly-nonlinear models. Some details of the secondary waves
appearing in complex problems as e.g. the Seaside case, also show some impact of not including
full-nonlinearity, which is certainly related to the inability of weakly nonlinear models to represent
some weak dispersive propagation as soon as it takes place above the initial shoreline. However the
differences are again not striking. Some genuinely 2D extensions of the 1D cases which reveal most
clearly the need of nonlinearity, while being of physical and engineering relevance, with possibly
openly available reference solutions or experimental data, seem still missing in literature.
<figure><img src="image_25.png" /></figure>

Figure 23: Wavelength and water level of undular bores in a channel with trapezoidal section. up-left: Water elevation (axis), up right: water elevation (banks), down wavelength. x: Data of Treske (1994) axis, +: data of Treske (1994) banks, circle: simulated values on the channel axis, squares: simulated values on the banks. Blue color: GN, Pink color: Nwogu, Green color: Nowgu-Abbot
<figure><img src="image_26.png" /><figcaption>Figure 24: Set up and location of the wave gauges for the seaside test case.</figcaption></figure>

We can thus draw two main conclusions. The first is that benchmarking habits need to be slightly revised. Authors, editors, and reviewers must be more careful in considering/suggesting appropriately chosen tests, so that when new models/schemes are proposed their capability for highly nonlinear propagation are clearly highlighted together with all other aspects (phase prediction, breaking location and strength, and so on). Also it should be always underlined when a benchmark, although involving highly nonlinear waves, is not sensitive to the fully/weakly nonlinear nature of the model. This work by no means aims at performing this task, which is possibly to be performed with a larger set of models/codes. The second conclusion we draw from our tests is that of course fully nonlinear models allow to safely approach problems of different nature and involving highly nonlinear waves. However, as a rule of thumb weakly nonlinear models can still be safely used in many problems of practical interest not involving planar non-breaking wave fronts reflecting on steep coastal structures in deep waters. This is true even for quite complex urban run-up of nonlinear waves, as seen in the Seaside example.

To answer the question in the title, the humble opinion of the authors is that full nonlinearity is neither luxury nor necessity. There is a large class of problems for which both appropriately chosen weakly and fully nonlinear models provide reasonable answers. So it is up to the user to be aware of the differences existing, and choose which of these models is most suited for his application.

## References

[1] Luca Arpaia and Mario Ricchiuto. r-adaptation for shallow water flows: conservation, well balancedness, efficiency. *Computers & Fluids*, 160:175–203, 2018.

[2] Paola Bacigaluppi, Mario Ricchiuto, and Philippe Bonneton. Implementation and evaluation of breaking detection criteria for a hybrid boussinesq model. *Water waves*, 2(2):207–241, 2020.

[3] S Beji and Jurjen A Battjes. Numerical simulation of nonlinear wave propagation over a bar. *Coastal Engineering*, 23(1-2):1–16, 1994.
<figure><img src="image_27.png" /><figcaption>Figure 25: 2 dimensional free surface elevation for GN (left) Nwogu (centered) Nwogu-Abott (right) at $t \approx 30$ sec.</figcaption></figure>

588

589

[4] S. Beji and K. Nadaoka. A formal derivation and numerical modeling of the improved boussinesq equations for varying depth. *Ocean Engineering*, 23, 1996.

590

591

592

[5] Michel Benoit, F Dias, J Herterich, and Yves-Marie Scolan. Un cas-test discriminant pour la simulation de la propagation et du run-up de trains de vagues de type tsunami. In *Actes des 16èmes Journées de l'Hydrodynamique*, 2018.

593

594

595

[6] Philippe Bonneton, Florent Chazel, David Lannes, Fabien Marche, and Marion Tissier. A splitting approach for the fully nonlinear and weakly dispersive green–naghdi model. *Journal of Computational Physics*, 230(4):1479–1498, 2011.

596

597

598

599

[7] Joseph Boussinesq. Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond. *Journal de mathématiques pures et appliquées*, pages 55–108, 1872.

600

601

602

[8] Maurizio Brocchini. A reasoned overview on boussinesq-type models: the interplay between physics, mathematics and numerics. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 469(2160):20130496, 2013.

603

604

605

[9] Rémi Chassagne, Andrea Gilberto Filippini, Mario Ricchiuto, and Philippe Bonneton. Dispersive and dispersive-like bores in channels with sloping banks. *Journal of Fluid Mechanics*, 870:595–616, 2019.

606

607

[10] F. Chazel, D. Lannes, and F. Marche. Numerical simulation of strongly nonlinear and dispersive waves using a green–naghdi model. *Journal of Scientific Computing*, 48(1):105–116, 2011.

608

609

610

[11] Didier Clamond, Denys Dutykh, and Dimitrios Mitsotakis. Conservative modified serre–green–naghdi equations with improved dispersion characteristics. *Communications in Nonlinear Science and Numerical Simulation*, 45:245–257, 2017.
611

612

613

[12] DHI. MIKE 21 Boussinesq wave module: scientific documentation.

https://manuals.mikepoweredbydhi.help/latest/Coast_and_Sea/MIKE21BW_Sci_Doc.php

614

615

616

[13] Arnaud Duran and Fabien Marche. A discontinuous galerkin method for a new class of green-naghdi equations on simplicial unstructured meshes. *Applied Mathematical Modelling*, 45:840-864, 2017.

617

618

[14] F. Shi et al. Funwave website: https://fengyanshi.github.io/build/html/index.html.

2023.

619

620

[15] P. Lynett et al. Coulwave website: http://isec.nacse.org/models/coulwave_description.php.

2004.

621

[16] S. Popinet et al. Basilisk website: http://basilisk.fr/Front%20Page. 2023.

622

623

[17] Henry Favre. Etude theorique et experimental des ondes de translation dans les canaux decouverts. *Dunod*, 150, 1935.

624

625

626

[18] Andrea Gilberto Filippini, Stevan Bellec, Mathieu Colin, and Mario Ricchiuto. On the nonlinear behaviour of boussinesq type models: Amplitude-velocity vs amplitude-flux forms. *Coastal Engineering*, 99:109-123, 2015.

627

628

629

[19] Andrea Gilberto Filippini, Maria Kazolea, and Mario Ricchiuto. A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up. *Journal of Computational Physics*, 310:381-417, 2016.

630

631

[20] Francesco Gallerano, Giovanni Cannata, and Mara Villani. An integral contravariant formulation of the fully non-linear boussinesq equations. *Coastal Engineering*, 83:119-136, 2014.

632

633

[21] Albert E Green and Paul M Naghdi. A derivation of equations for wave propagation in water of variable depth. *Journal of Fluid Mechanics*, 78(2):237-246, 1976.

634

635

[22] ST Grilli, R Subramanya, IA Svendsen, and J Veeramony. Shoaling of solitary waves on plane beaches. *Oceanographic Literature Review*, 8(42):608, 1995.

636

637

[23] S. M. Joshi, M. Kazolea, and M. Ricchiuto. Parameter sensitivity for wave-breaking closures in boussinesq-type models. *Water Waves*, 4(3):491-515, 2022.

638

639

[24] Maria Kazakova and Pascal Noble. Discrete transparent boundary conditions for the linearized green-naghdi system of equations. *SIAM Journal on Numerical Analysis*, 58(1):657-683, 2020.

640

641

642

[25] M Kazolea and AI Delis. A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1d boussinesq models. *Applied Numerical Mathematics*, 67:167-186, 2013.

643

644

[26] M Kazolea, AI Delis, IK Nikolos, and CE Synolakis. An unstructured finite volume numerical scheme for extended 2d boussinesq-type equations. *Coastal Engineering*, 69:42-66, 2012.

645

646

647

[27] M Kazolea, Argiris I Delis, and Costas E Synolakis. Numerical treatment of wave breaking on unstructured finite volume approximations for extended boussinesq-type equations. *Journal of computational Physics*, 271:281-305, 2014.
648

649

650

[28] M Kazolea, AG Filippini, and M Ricchiuto. Low dispersion finite volume/element discretization of the enhanced green–naghdi equations for wave propagation, breaking and runup on unstructured meshes. *Ocean Modelling*, 182:102157, 2023.

651

652

653

654

[29] Maria Kazolea, Andrea Filippini, Mario Ricchiuto, Stéphane Abadie, M Martin Medina, Denis Morichon, Camille Journeau, Richard Marcer, Kévin Pons, Sylvestre LeRoy, et al. Wave propagation, breaking, and overtopping on a 2d reef: A comparative evaluation of numerical codes for tsunami modelling. *European Journal of Mechanics-B/Fluids*, 73:122–131, 2019.

655

656

[30] Maria Kazolea and Mario Ricchiuto. On wave breaking for boussinesq-type models. *Ocean Modelling*, 123:16–39, 2018.

657

658

[31] James T Kirby. Boussinesq models and their application to coastal processes across a wide range of scales. *Journal of Waterway, Port, Coastal, and Ocean Engineering*, 142(6):03116005, 2016.

659

660

[32] D. Lannes and P. Bonneton. Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation. *Physics of Fluids*, 21, 2009. 016601 doi:10.1063/1.3053183.

661

662

[33] D Lannes and L Weynans. Generating boundary conditions for a boussinesq system*. *Nonlinearity*, 33(12):6868, 2020.

663

[34] David Lannes. Modeling shallow water waves. *Nonlinearity*, 33(5):R1, 2020.

664

665

666

[35] David Lannes and Fabien Marche. A new class of fully nonlinear and weakly dispersive green-naghdi models for efficient 2d simulations. *Journal of Computational Physics*, 282:238–268, 2015.

667

668

669

[36] Maojun Li, Liwei Xu, and Yongping Cheng. A cdg-fe method for the two-dimensional green-naghdi model with the enhanced dispersive property. *Journal of Computational Physics*, 399:108953, 2019.

670

671

672

[37] P. Lynett, P. Liu, K. Sitanggang, and D.-H. Kim. Modeling Wave Generation Evolution and Interaction with Depth-Integrated Dispersive Wave Equations COULWAVE Code Manual. Technical Report, Cornell University, 2004.

673

674

[38] Per A Madsen, Russel Murray, and Ole R Sørensen. A new form of the boussinesq equations with improved linear dispersion characteristics. *Coastal engineering*, 15(4):371–388, 1991.

675

676

677

[39] Per A Madsen and Ole R Sørensen. A new form of the boussinesq equations with improved linear dispersion characteristics. part 2. a slowly-varying bathymetry. *Coastal engineering*, 18(3-4):183–204, 1992.

678

679

680

[40] Per A Madsen, OR Sørensen, and HA Schäffer. Surf zone dynamics simulated by a boussinesq type model. part i. model description and cross-shore motion of regular waves. *Coastal Engineering*, 32(4):255–287, 1997.

681

682

683

[41] Christopher Michalak and Carl Ollivier-Gooch. Accuracy preserving limiter for the high-order accurate solution of the euler equations. *Journal of Computational Physics*, 228(23):8693–8711, 2009.
684
685

[42] Dimitrios Mitsotakis, Boaz Ilan, and Denys Dutykh. On the galerkin/finite-element method for the serre equations. *Journal of Scientific Computing*, 61(1):166–195, 2014.

686
687
688

[43] Sebastian Noelle, Martin Parisot, and Tabea Tscherpel. A class of boundary conditions for time-discrete green-naghdi equations with bathymetry. *SIAM Journal on Numerical Analysis*, 60(5):2681–2712, 2022.

689
690
691

[44] O.G. Nwogu and Z. Demirbilek. BOUSS-2D: A Boussinesq Wave Model for Coastal Regions and Harbors. US Army Corps of Engineers, Engineering Research and Development Center - Report TR-01-25, September 2001.

692
693

[45] Okey Nwogu. Alternative form of boussinesq equations for nearshore wave propagation. *Journal of waterway, port, coastal, and ocean engineering*, 119(6):618–638, 1993.

694
695
696
697

[46] Nishant Panda, Clint Dawson, Yao Zhang, Andrew B. Kennedy, Joannes J. Westerink, and Aaron S. Donahue. Discontinuous galerkin methods for solving boussinesq-green-naghdi equations in resolving non-linear and dispersive surface water waves. *Journal of Computational Physics*, 273:572–588, 2014.

698
699

[47] M. Parisot. Entropy-satisfying scheme for a hierarchy of dispersive reduced models of free surface flow. *International Journal for Numerical Methods in Fluids*, 91(10):509–531, 2019.

700
701
702

[48] Hyoungsu Park, Daniel T Cox, Patrick J Lynett, Dane M Wiebe, and Sungwon Shin. Tsunami inundation modeling in constructed environments: A physical and numerical comparison of free-surface elevation, velocity, and momentum flux. *Coastal Engineering*, 79:9–21, 2013.

703
704

[49] D Howell Peregrine. Calculations of the development of an undular bore. *Journal of Fluid Mechanics*, 25(2):321–330, 1966.

705

[50] D Howell Peregrine. Long waves on a beach. *Journal of fluid mechanics*, 27(4):815–827, 1967.

706
707

[51] Jordan P.A. Pitt, Christopher Zoppou, and Stephen G. Roberts. Numerical scheme for the generalised serre-green-naghdi model. *Wave Motion*, 115:103077, 2022.

708
709

[52] Stéphane Popinet. A quadtree-adaptive multigrid solver for the serre-green-naghdi equations. *Journal of Computational Physics*, 302:336–358, 2015.

710
711

[53] M. Ricchiuto. An explicit residual based approach for shallow water flows. *J.Comput.Phys.*, 80:306–344, 2015.

712
713
714

[54] M. Ricchiuto and A.G. Filippini. Upwind residual discretization of enhanced boussinesq equations for wave propagation over complex bathymetries. *J. Comput. Phys.*, 271:306–341, 2014.

715
716

[55] Mario Ricchiuto and Andreas Bollermann. Stabilized residual distribution for shallow water simulations. *Journal of Computational Physics*, 228(4):1071–1115, 2009.

717
718

[56] Volker Roeber. *Boussinesq-type model for nearshore wave processes in fringing reef environment*. PhD thesis, [Honolulu]:[University of Hawaii at Manoa], [December 2010], 2010.

719
720

[57] Volker Roeber and Kwok Fai Cheung. Boussinesq-type model for energetic breaking waves in fringing reef environments. *Coastal Engineering*, 70:1–20, 2012.
721
722

[58] François Serre. Contribution à l'étude des écoulements permanents et variables dans les canaux. *La Houille Blanche*, (6):830-872, 1953.

723
724
725

[59] Fengyan Shi, James T. Kirby, Jeffrey C. Harris, Joseph D. Geiman, and Stephan T. Grilli. A high-order adaptive time-stepping tvd solver for boussinesq modeling of breaking waves and coastal inundation. *Ocean Modelling*, 43-44:36-51, 2012.

726
727

[60] J. J. Stoker. *Water Waves: The Mathematical Theory with Applications*. Pure and Applied Mathematics. Interscience Publishers, New York, USA, 1957.

728
729
730

[61] Sasan Tavakkol and Patrick Lynett. Celeris: A gpu-accelerated open source software with a boussinesq-type wave solver for real-time interactive simulation and visualization. *Computer Physics Communications*, 217:117-127, 2017.

731
732

[62] Sasan Tavakkol and Patrick Lynett. Celeris base: An interactive and immersive boussinesq-type nearshore wave simulation software. *Computer Physics Communications*, 248:106966, 2020.

733
734
735

[63] Sergey Tkachenko, Sergey Gavrilyuk, and Jacques Massoni. Extended Lagrangian approach for the numerical study of multidimensional dispersive waves: Applications to the Serre-Green-Naghdi equations. *Journal of Computational Physics*, 477:111901, 2023.

736
737

[64] Mara Tonelli and Marco Petti. Simulation of wave breaking over complex bathymetries by a boussinesq model. *Journal of Hydraulic Research*, 49(4):473-486, 2011.

738
739

[65] Andreas Treske. Undular bores (favre-waves) in open channels-experimental studies. *Journal of Hydraulic Research*, 32(3):355-370, 1994.

740
741

[66] Claudio Viotti, Francesco Carbone, and Frédéric Dias. Conditions for extreme wave runup on a vertical barrier by nonlinear dispersion. *Journal of Fluid Mechanics*, 748:768-788, 2014.

742
743
744

[67] Ge Wei, James T Kirby, Stephan T Grilli, and Ravishankar Subramanya. A fully nonlinear boussinesq model for surface waves. part 1. highly nonlinear unsteady waves. *Journal of fluid mechanics*, 294:71-92, 1995.
<figure><img src="image_28.png" /><figcaption>Figure 26: Free surface elevation measured on wave gauges at location A</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Figure 27: Free surface elevation measured on wave gauges at location B</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 28: Free surface elevation measured on wave gauges at location C</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure 29: Free surface elevation measured on wave gauges at location D</figcaption></figure>