# On a high-order shallow-water wave model with canonical non-local Hamiltonian structure

Christos Papoutsellis, Michel Benoit

► To cite this version:

Christos Papoutsellis, Michel Benoit. On a high-order shallow-water wave model with canonical non-local Hamiltonian structure. Physica D: Nonlinear Phenomena, 2025, 479, pp.134691. 10.1016/j.physd.2025.134691. hal-05070684v2

HAL Id: hal-05070684

https://edf.hal.science/hal-05070684v2

Submitted on 22 May 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>

# On a high-order shallow-water wave model with canonical non-local Hamiltonian structure

Christos E. Papoutsellis<sup>a,∗</sup>, Michel Benoit<sup>b,c</sup>

## ARTICLE INFO

Communicated by Ye Zhou

**Keywords:**

Free-surface long gravity waves

Hamiltonian formulation

Solitary waves

## ABSTRACT

We derive and study a new family of non-local partial differential equations (PDEs) that model free-surface long gravity waves over a flat bottom. To derive the model equations we approximate the velocity potential as a series of vertical polynomials derived from the shallow-water expansion of the Dirichlet-to-Neumann problem in the Hamiltonian formulation of free-surface potential flow and invoke Luke's variational principle. The resulting evolution equations exhibit a non-local Hamiltonian structure being coupled with a system of linear elliptic spatial PDEs on the horizontal plane. A key advantage of this approach is that it directly yields canonical Hamiltonian equations, which are well-suited for numerical solutions using standard methods. This class of model equations offers high-order shallow-water approximations of the water-wave problem. It contains terms whose spatial derivatives are at most of order two, distinguishing it from asymptotic methods involving higher-order mixed spatio-temporal derivatives. We explore the first non-trivial member of this family, highlighting its connections to other mathematical models and emphasizing its practical utility. We then analyse and discuss its linear dispersive properties and demonstrate that it does not exhibit a specific type of instability known as wave-trough instability. Additionally, we demonstrate its effectiveness in simulating the long-distance steady propagation of strongly non-linear solitary waves and the head-on collision of two counter-propagating solitary waves. In the latter case, comparisons with experimental data confirm the model's ability to capture complex wave dynamics, including wave transformation in the presence of strong non-linearity and dispersion. The extension of this approach to accommodate variable bottom topography is briefly discussed.

## 1. Introduction

Equations modelling free-surface gravity waves are continuously evolving in connection to applications requiring extended spatial domains, longer simulation times and higher fidelity in the representation of physical processes. They aim to replace the original free-boundary hydrodynamic problem with one that is simpler to solve numerically while retaining the maximum amount of physics possible. Since the time of Boussinesq, several models have been proposed (see, e.g., Brocchini [1], Kirby [2], Lannes [3]), but the need for improvement still exists, particularly in capturing non-linear and dispersive effects. Improved wave models are crucial for tsunami prediction, extreme wave simulations, and coastal engineering applications, where balancing physical realism with numerical efficiency is essential. In this paper, we derive and study a new system of non-linear partial differential equations (PDEs) for modelling shallow water waves, that is convenient for numerical implementation while preserving key physical properties.

In finite depth, two important parameters to consider are the characteristic water depth $h_0$ and wavelength $L$. These parameters form the ratio $\mu \equiv (h_0/L)^2$, which is commonly referred to as the shallowness or dispersion parameter. Assuming that $\mu \ll 1$, we are placed into the shallow-water or long-wave regime, and we may invoke asymptotic expansions of the unknowns in terms of $\mu$. The original problem can then be simplified by neglecting terms above a certain order in $\mu$.

A well-known model of this class was originally derived by Serre [4] and later by many others in various equivalent forms [5–13], and it is predominantly called the Serre-Green–Naghdi (SGN) model nowadays. The SGN model is accurate to $O(\mu)$ or $O(\mu^2)$-accurate, i.e. it is derived by neglecting $O(\mu^2)$ terms in the invoked asymptotic expansion. It consists of two PDEs governing the evolution of the free-surface elevation and

<footnote><sup>a</sup> Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla, CA 92093, USA</footnote>

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

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

<footnote>* Corresponding author.
E-mail addresses: cpapoutsellis@ucsd.edu, cpapoutsellis@gmail.com (C.E. Papoutsellis).</footnote>
the depth-averaged horizontal velocity. This model has been extensively studied theoretically (see Lannes [14] and the references therein) and numerically [15–23]. The SGN equations involve space–time derivatives and their evolution in time requires the inversion of an elliptic operator on the horizontal plane. While SGN is a well-established model, it is important to note that it needs a linear dispersion improvement in strongly dispersive cases.

A natural improvement of $O(\mu)$-accurate models is to retain higher order terms in the associated asymptotic expansions, as in [24–28] for instance. This strategy requires special attention and the accuracy and applicability of the resulting equations depend on the choice of variables and/or the order of precision [29]. Recently, Matsuno [30], Matsuno [31] (see also Horikis et al. [32]) extended the SGN equations and derived their Hamiltonian structure. One of the findings in [31] is that the SGN equations accurate to $O(\mu^p)$ with $p$ even are linearly ill-posed — their dispersion relation is singular at a finite wave number. Thus, beyond $O(\mu)$ the first practically useful SGN model is the one accurate to $O(\mu^3)$ which involves seventh-order spatial derivatives and its numerical solution would be quite challenging. A similar situation occurs when one derives model equations in terms of the free-surface or bottom-surface potential by truncating the asymptotic expansion of the velocity potential. In this case, the equations accurate to $O(\mu^p)$ with $p$ odd are linearly ill-posed [29,33]. It is widely recognized that some models based on truncated asymptotic expansions may lead to non-physical high-frequency instabilities. We refer to Milder [34], Yoon and Liu [35], Radder [36], Madsen and Fuhrman [37], Kirby [38], Choi [33] for more details on this issue. A recent study of the properties and performance of several asymptotic models can be found in [39].

Recently, Choi [33] addressed the above challenges by exploring asymptotic expansions for both the velocity potential and its vertical derivative, developing models that are linearly well-posed at any desired order. These models are formulated in terms of the free-surface elevation and bottom surface potential and lead to new approximate high-order solitary-wave solutions. Furthermore, the study revealed that finite-difference schemes are unstable for this class of equations, prompting the development of a Fourier-spectral method for the stable numerical solution of solitary-wave problems. We note finally that asymptotic model derivation is not limited to one-layer fluids and several SGN-type models have also been developed for two-layer fluids (interfacial waves), see Choi and Camassa [40], Craig et al. [41], Guyenne [42], Choi [43] for instance.

A different strategy to derive model equations is to use convergent vertical series expansions for the velocity potential without an a priori smallness assumption. In this direction, we mention the coupled-mode systems obtained from Luke [44]'s variational principle [45,46] using an enhanced eigenfunction expansion [47], and the system obtained from a Chebyshev expansion and a Galerkin technique [48,49]. The evolution equations of these models retain the Hamiltonian structure of the original problem [50,51], known as the Zakharov/Craig-Sulem (ZCS) formulation. This approach replaces the Laplace equation for the velocity potential with a system of horizontal spatial PDEs. The number of PDEs depends on the order of truncation of the expansion and the accuracy increases with the order of truncation [46–48,52,53]. Provided a sufficient number of PDEs is kept (typically 4–7), these models can treat highly demanding cases with precision comparable to direct numerical solutions [48,49,52–58]. An advantage of this non-perturbative class of models compared to the asymptotic ones is that they contain up to second-order derivatives, which makes their numerical implementation straightforward. On the other hand, solving more than one or two horizontal spatial PDEs during the evolution may be considered expensive in some cases and acceleration techniques may be needed.

Another approach is to use a simple expression for the velocity potential or vertical series representations without an a priori assumption of convergence. This approach focuses on the development of simplified models, typically within the variational framework of the problem. For a detailed exposition, we refer to the recent monograph of Clamond et al. [59]. Isobe [60], Kakinuma [61] used the vertical polynomials of the classical Boussinesq-Rayleigh expansion [62–64]. The resulting system of equations has a non-trivial mathematical structure that attracted the attention of Iguchi and collaborators. In a series of works, they established its well-posedness and rigorously justified it as a high-order approximation of the full problem in terms of $\mu$ [65–67]. They also constructed solitary wave solutions [68] and established its relation to the ZCS formulation [69]. For some numerical results, we refer to Abohadima and Isobe [70]. Klopman et al. [71] proposed several other choices of vertical functions that lead directly to two evolution equations coupled with horizontal PDEs. A version of this class of models uses a hyperbolic vertical profile involving a free parameter. By choosing the parameter that optimizes the linear dispersion, this model shows very good properties [72]. Another possibility is to choose representations for the velocity variables and derive model equations from the free-surface Euler's equations as in the so-called Irrotational Green–Naghdi equations [73].

Despite significant progress in developing shallow-water wave models, existing approaches face notable challenges. Low-order asymptotic models often have a limited range of applicability, failing to capture strong dispersive effects. On the other hand, higher-order asymptotic models introduce high-order derivatives, which complicate numerical solutions and may be a source of instabilities [31,37,38]. The model recently derived by Choi [33] provides a promising alternative, however, it does not exactly conserve energy, which may pose issues for long-time simulations. Moreover, its instability when implemented with finite-difference schemes further limits its versatility. To address these challenges, this study aims to develop and analyse a new class of shallow-water wave models that provide high-order approximations of the original problem while avoiding high-order derivatives. To achieve this, we propose a new ansatz containing the vertical functions of an asymptotic expansion of the velocity potential in the context of the ZCS formulation. By invoking Luke's variational principle, we obtain a new family of equations in the form of two evolution PDEs in terms of the free-surface elevation and potential coupled with a system of elliptic spatial second-order PDEs on the unknown coefficients of the ansatz. This new family of equations preserves the Hamiltonian structure of the original problem, as it is formulated in canonical variables, while also maintaining the inherent non-locality of the ZCS formulation.

The remainder of this paper is organized as follows. In Section 2, we review the formulations of the water wave problem with an emphasis on the Hamiltonian one. In Section 3, we give a systematic approach to derive approximate Euler–Lagrange equations from Luke's functional. Then, we construct our ansatz and derive the corresponding family of equations. In Section 4, we show the Hamiltonian structure of this new family of equations. In Section 5, we show the equivalence of the new family with the Isobe–Kakinuma (IK) family. In Section 6, we consider the first non-trivial member of this family of equations and discuss its relation to several other existing asymptotic models. In Section 7, we discuss the linear dispersive properties of the model. We also demonstrate that, unlike some other wave models, it is not subject to a specific wave trough instability. In Section 8, we present two numerical examples to demonstrate the model's performance. The first examines the long-distance propagation of strongly non-linear solitary waves. The second focuses on the head-on collision of two counter-propagating solitary waves, comparing the results with experimental data.
## 2. The irrotational water-wave problem

We consider waves on the surface of a homogeneous, incompressible and inviscid fluid driven by gravity (surface tension is neglected). In a coordinate system $OXz$ with horizontal coordinates $X = (x_1, x_2) \in \mathbb{R}^2$ and the vertical axis, $z$, pointing upwards (with $z = 0$ at the still water level), the fluid domain is bounded by the free surface $z = \eta(X, t)$, and a fixed bottom surface $z = -h(X)$, with $h(X) > 0$. We assume the atmospheric pressure at the free surface is constant and homogeneous, and taken here to be 0 without loss of generality. For an irrotational flow, the motion is described in terms of the velocity potential $\Phi = \Phi(X, z, t)$ and the free surface elevation $\eta = \eta(X, t)$ by the following set of equations, see e.g. [14,64,74]

$$ \Delta_{X,z} \Phi = 0, \quad -h \le z \le \eta, \qquad (1a) $$

$$ N_h \cdot \nabla \Phi = 0, \quad z = -h, \qquad (1b) $$

$$ \partial_t \eta - N_\eta \cdot \nabla \Phi = 0, \quad z = \eta, \qquad (1c) $$

$$ \partial_t \Phi + \frac{1}{2}(\nabla \Phi)^2 + g\eta = 0, \quad z = \eta, \qquad (1d) $$

where $g$ is the acceleration of gravity and we denote the three-dimensional gradient and Laplacian operators by $\nabla = (\partial/\partial x_1, \partial/\partial x_2, \partial/\partial z) \equiv (\partial_1, \partial_2, \partial_z)$ and $\Delta_{X,z} = \partial_1^2 + \partial_2^2 + \partial_z^2$ respectively. We also introduce the horizontal gradient and Laplacian operators, denoted respectively by $\nabla$ and $\Delta$, so that $\nabla = (\nabla, \partial_z)$ and $\Delta_{X,z} = \Delta + \partial_z^2$, and the (non-unit) outward normal vectors at the free surface and seabed, respectively:

$$ N_\eta = (-\nabla\eta, 1), \quad N_h = (-\nabla h, -1). \qquad (2) $$

We also introduce the total fluid depth $H = H(X, t) = \eta(X, t) + h(X) > 0$ and assume vanishing at infinity conditions

$$ \lim_{|X| \to \infty} (|\eta| + |\nabla \Phi|) = 0, \qquad (3) $$

which we will regularly use to neglect terms at infinity without explicit mention.

### 2.1. Hamiltonian formulation

As was originally observed by Zakharov [50], the above equations allow for a neat reformulation:

$$ \partial_t \eta = \delta_\psi H = G[\eta]\psi, \qquad (4a) $$

$$ \partial_t \psi = -\delta_\eta H = -g\eta - \frac{1}{2} |\nabla \psi|^2 + \frac{(G[\eta]\psi + \nabla\eta \cdot \nabla\psi)^2}{2(1 + |\nabla\eta|^2)}, \qquad (4b) $$

where $\psi(X, t) \stackrel{\text{def}}{=} \Phi(X, \eta(X, t), t)$ is the trace of the velocity potential on the free surface and $G[\eta]\psi$ is the Dirichlet-to-Neumann (DtN) operator introduced by Craig and Sulem [51],

$$ G[\eta]\psi \stackrel{\text{def}}{=} N_\eta \cdot [\nabla(\Phi[\eta]\psi)]_\eta = -\nabla\eta \cdot \nabla\psi + ((\nabla\eta)^2 + 1) [\partial_z(\Phi[\eta]\psi)]_\eta, \qquad (5) $$

with $\Phi = \Phi[\eta]\psi$ the solution of the boundary value problem

$$ \Delta_{X,z} \Phi = 0, \quad N_h \cdot [\nabla \Phi]_{-h} = 0, \quad [\Phi]_\eta = \psi, \qquad (6a, b, c) $$

where we have used the notation $[u]_s = u(X, s(X, t), t)$ to denote the trace of a function $u(X, z, t)$ at $z = s(X, t)$. The system (4) is known as the ZAKHAROV/CRAIG-SULEM (ZCS) formulation. The Hamiltonian $H = H(\eta, \psi)$ represents the total mechanical energy of the fluid, composed of kinetic energy $K = K(\eta, \psi)$ and potential energy $V = V(\eta)$,

$$ H = K + V, \quad \text{with} \quad \begin{cases} K = \frac{1}{2} \int \int_{-h}^{\eta} (\nabla \Phi[\eta] \psi)^2 dz \, dX = \frac{1}{2} \int \int \psi G[\eta] \psi \, dX \, dt. \\ V = \frac{1}{2} g \int \eta^2 \, dX. \end{cases} \qquad (7) $$

The operator $G$ is non-local unlike a local differential operator. This means that determining $G[\eta]\psi$ for some $t$ at a given position $X$ requires knowledge of the entire fields $(\eta, \psi)$ and the solution of (6a,b,c), rather than their local values around $X$. The variational derivatives of the non-local Hamiltonian $H$, $\delta_\eta H \equiv \delta H/\delta\eta$ and $\delta_\psi H \equiv \delta H/\delta\psi$, appearing in (4), are calculated in Zakharov [50], Benjamin and Olver [75], see also Craig et al. [76].

### 2.2. Luke's variational principle

Luke [44] noted that the equations in (1) are the EULER-LAGRANGE (EL) equations of the pressure-type functional

$$ S(\eta, \Phi) = \int \int \int_{-h}^{\eta} \left[ \partial_t \Phi + \frac{1}{2} (\nabla \Phi)^2 + gz \right] dz \, dX \, dt. \qquad (8) $$

To derive (1) from (8), we require the first variation of $S$ to vanish for arbitrary $(\delta\eta, \delta\Phi)$ that vanish at the boundaries of the horizontal domain and temporal interval:

$$ D_\eta S(\eta, \Phi)(\delta\eta) + D_\Phi S(\eta, \Phi)(\delta\Phi) = 0. \qquad (9) $$

Here, $D_\eta S(\eta, \Phi)(\delta\eta)$ and $D_\Phi S(\eta, \Phi)(\delta\Phi)$ are given by

$$ D_\Phi S(\delta\Phi) = \int \int \left\{ - \int_{-h}^{\eta} \Delta_{X,z} \Phi \, \delta\Phi \, dz - (\partial_t \eta - N_\eta \cdot [\nabla\Phi]_\eta) [\delta\Phi]_\eta + N_h \cdot [\nabla\Phi \, \delta\Phi]_{-h} \right\} dX \, dt, \qquad (10) $$
and

$$D_{\eta} S(\delta\eta) = \iint \left[ \partial_t \Phi + \frac{1}{2} (\nabla \Phi)^2 + gz \right]_{\eta} \delta\eta \, dX \, dt. \qquad (11)$$

### 3. Variational derivation of model equations

Our goal in this section is to derive model equations for (1) by invoking (9) with a vertical series representation for $\Phi$. A general form of such a representation is

$$\Phi^a(X, z, t) = \sum_{n=0}^{K} \phi_n(X, t) Z_n(z; \eta(X, t), h(X)). \qquad (12)$$

In (12), we assume that $Z_n = Z_n(z; \eta(X, t), h(X))$ are known smooth functions of $z, h$ and $\eta$, and that $\eta$ and $\phi_n = \phi_n(X, t)$ are unknown smooth functions for which equations are to be derived. We also assume that $\Phi^a$ satisfies the vanishing at infinity conditions (3). By starting the derivation with (12), we can better comprehend the structure of the resulting equations in relation to the original problem (1). We will specify the functions $Z_n$ subsequently.

#### 3.1. Euler-Lagrange equations

Substituting $\Phi^a$ into (8), we obtain a new functional $S^a(\eta, \phi_0, \dots, \phi_K) = S(\eta, \Phi^a)$. The EL equations of $S^a$ can be derived using the chain rule [46,52,77],

$$D_{\phi_m} S^a(\eta, \phi_0, \dots, \phi_K)(\delta\phi_m) = D_{\phi} S(\eta, \Phi^a)(Z_m \delta\phi_m) = 0, \quad m = 0, \dots, K, \qquad (13a)$$

$$D_{\eta} S^a(\eta, \phi_0, \dots, \phi_K)(\delta\eta) = D_{\eta} S(\eta, \Phi^a)(\delta\eta) + D_{\phi} S(\eta, \Phi^a)(\partial_{\eta} \Phi^a \delta\eta) = 0, \qquad (13b)$$

where $\partial_{\eta} \Phi^a \equiv \partial \Phi^a / \partial \eta = \sum_n \phi_n \partial_{\eta} Z_n$. Eqs. (13a)-(13b) with (10)-(11), (12) and the relation $\partial_t \Phi^a = \sum_n \partial_t \phi_n Z_n + \phi_n \partial_{\eta} Z_n \partial_{\eta} \eta$ lead to the following equations on $(\eta, \phi_0, \dots, \phi_K)$

$$(m = 0, \dots, K) \quad \left( \partial_t \eta - N_{\eta} \cdot [\nabla(\phi_n Z_n)]_{\eta} \right) [Z_m]_{\eta} + L_{mn} \phi_n = 0, \qquad (14)$$

$$\partial_t \phi_n [Z_n]_{\eta} + g\eta + \frac{1}{2} [\nabla(\phi_n Z_n)]^2_{\eta} + N_{\eta} \cdot [\nabla(\phi_n Z_n)(\phi_n \partial_{\eta} Z_n)]_{\eta} - \Lambda_{mn} \phi_n \phi_m = 0, \qquad (15)$$

where we used the Einstein summation convention

$$a_n b_n = \sum_n a_n b_n. \qquad (16)$$

In (14) and (15) we have introduced the operators $L_{mn} \equiv L_{mn}(\eta, h)$ and $\Lambda_{mn} \equiv \Lambda_{mn}(\eta, h)$ given by

$$L_{mn} = A_{mn} \Delta + \underline{B}_{mn} \cdot \nabla + C_{mn}, \quad \Lambda_{mn} = a_{mn} \Delta + \underline{b}_{mn} \cdot \nabla + c_{mn}, \qquad (17a,b)$$

where $A_{mn}, C_{mn}, a_{mn}, c_{mn} \in \mathbb{R}$ and $\underline{B}_{mn}, \underline{b}_{mn} \in \mathbb{R}^2$ are given by

$$A_{mn} = (Z_n, Z_m) = A_{nm}, \qquad a_{mn} = (Z_n, \partial_{\eta} Z_m), \qquad (18a)$$

$$\underline{B}_{mn} = 2(\nabla Z_n, Z_m) + \nabla h [Z_n Z_m]_{-h}, \qquad \underline{b}_{mn} = 2(\nabla Z_n, \partial_{\eta} Z_m) + \nabla h [Z_n \partial_{\eta} Z_m]_{-h} \qquad (18b)$$

$$C_{mn} = (\Delta_{X,z} Z_n, Z_m) - N_h \cdot [\nabla Z_n Z_m]_{-h}, \qquad c_{mn} = (\Delta_{X,z} Z_n, \partial_{\eta} Z_m) - N_h \cdot [\nabla Z_n \partial_{\eta} Z_m]_{-h} \qquad (18c)$$

with $(Z, W) \stackrel{\text{def}}{=} \int_{-h}^{h} ZW \, dz$. These operators stem from the first and last terms of (10) and satisfy the identities

$$\begin{pmatrix} L_{mn} \\ \Lambda_{mn} \end{pmatrix} \phi_n = \int_{-h}^{h} \Delta_{X,z}(\phi_n Z_n) \begin{pmatrix} Z_m \\ \partial_{\eta} Z_m \end{pmatrix} dz - N_h \cdot [\nabla(\phi_n Z_n)]_{-h} \left[ \begin{pmatrix} Z_m \\ \partial_{\eta} Z_m \end{pmatrix} \right]_{-h}, \qquad (19)$$

which hold for $m = 0, \dots, K$ and will be used in the sequel.

#### 3.2. Shallow-water ansatz

We can now specify the vertical functions $Z_n$ in $\Phi^a$. Our choice is based on an asymptotic expansion of $\Phi$ in terms of the shallowness parameter $\mu = (h_0/L)^2$, where $h_0$ and $L$ are the typical depth and wavelength respectively, due to Lannes and Bonneton [9].

**Proposition 3.1.** Consider the dimensionless version of the DtN problem (6a,b,c) for a flat bottom, ($N_h = (0, 0, -1)$):

$$\mu \Delta \Phi + \partial_z^2 \Phi = 0, \quad -1 \le z \le \eta, \quad [\partial_z \Phi]_{-1} = 0, \quad [\Phi]_{\eta} = \psi, \qquad (20)$$

where $\mu$ is the shallowness parameter, see Appendix A. If $\mu \ll 1$, then the solution of (20) is approximated by

$$\Phi = \psi + \mu \left\{ \Delta\psi + \mu \left[ \frac{1}{2} \Delta(H^2) + \nabla(H^2) \cdot \nabla + \frac{H^2}{3} \Delta \right] \Delta\psi \right\} \tilde{Z} + \mu^2 \left\{ \frac{1}{6} \Delta^2 \psi \right\} \tilde{Z}^2 + O(\mu^3), \qquad (21)$$

where the vertical dependence is given by the function

$$\tilde{Z} \equiv \tilde{Z}(z; \eta) = -\frac{1}{2}(\eta - z)^2 + (\eta + 1)(\eta - z). \qquad (22)$$
The proof of Proposition 3.1 follows Lannes [14, Lemma 3.42] and is given in Appendix B. An equivalent expression appears in Iguchi [78, p. 958] which is (21) transformed to the flat strip, see also Khorbatly et al. [79]. The expression (21) is distinct in its ability to highlight the vertical structure of $\Phi$ in a concise way that suits our purposes. Boussinesq [62], Rayleigh [63] (see also Whitham [64]) derived another expansion by using the first two equations in (20). This expansion reads

$$ \Phi = \sum_{n=0}^{K} (-1)^n \mu^n \frac{1}{(2n)!} \Delta^n f (z+1)^{2n} + O(\mu^{K+1}), \quad (23) $$

where $f = f(X,t)$ denotes the bottom-surface potential. As mentioned in the introduction, constructing an ansatz based on (23) leads to the IK model. These two approaches are reconciled in Section 5.

Here, we construct $\Phi^a$ by picking the vertical functions in (21) in dimensional form. Slightly changing the notation to distinguish this $\Phi^a$ from the general expression (12), we write

$$ \Phi^a(X, z, t) = \psi(X, t)Z_0 + \sum_{n=1}^{K} \varphi_n(X, t) Z_n(z; \eta(X, t), h(X)), \quad (24) $$

$$ \text{with } Z_0 = 1, Z_n \equiv Z^n, n \ge 1 \text{ and } Z \equiv Z(z; \eta) = -\frac{1}{2}(\eta - z)^2 + (\eta + h)(\eta - z). $$

In other words, we extremize the functional (8) for functions of the form (24). The vertical polynomials in (24) satisfy

$$ [Z_n]_{\eta} = 0, \quad n = 1, \dots, K, \quad (25) $$

which is the condition considered by Klopman et al. [71]. This is a key property because it turns (14) into one evolution equation for $\eta$ (for $m=0$) and a system of spatial PDEs on $\varphi_n$ (for $m=1, \dots, K$), and (15) into an evolution equation on $\psi$. Moreover, from (25) we obtain the relations

$$ [\partial_z Z_n]_{\eta} + [\partial_{\eta} Z_n]_{\eta} = 0, \quad [\partial_h Z_n]_{\eta} = 0, \quad (26) $$

which combined with (24) and (25) give

$$ [\nabla \Phi^a]_{\eta} = \nabla \psi + \varphi_n [\partial_{\eta} Z_n]_{\eta} \quad \nabla \eta = \nabla \psi - \varphi_n [\partial_z Z_n]_{\eta} \nabla \eta. \quad (27) $$

Our vertical functions also satisfy

$$ [\partial_z Z_n]_{-h} = 0, \quad (28) $$

which implies that $\Phi^a$ satisfies exactly the flat bottom condition $[\partial_z \Phi^a]_{-h} = 0$ for any $h$. For the sake of consistency, we assume a flat bottom ($\partial_t h = 0, i = 1, 2$) from this point onward. It is important to note that this assumption is not mandatory for deriving the model equations from (14)–(15); minor modifications will be needed should this assumption be lifted. See the discussion at the end of Section 9 regarding how $\Phi^a$ can be consistently modified in this case.

### 3.3. Derivation of the model equations

Using the above properties, the EL Eqs. (14)–(15) become

$$ \partial_t \eta - N_{\eta} \cdot [\nabla(\psi + \varphi_n Z_n)]_{\eta} + H \Delta \psi + L_{0n} \varphi_n = 0, \quad (29) $$

$$ (m = 1, \dots, K) \quad (Z_m, 1) \Delta \psi + L_{mn} \varphi_n = 0, \quad (30) $$

$$ \partial_t \psi + g \eta + \frac{(\nabla \psi)^2}{2} - \frac{(\nabla \eta)^2 + 1}{2} (\varphi_n [\partial_z Z_n]_{\eta})^2 - \Delta \psi (\partial_{\eta} Z_m, 1) \varphi_m - \Lambda_{mn} \varphi_n \varphi_m = 0, \quad (31) $$

where the summation is performed for $m, n \ge 1$ and the identities in (19) take the form

$$ L_{mn} \varphi_n = \int_{-h}^{\eta} \Delta_{X,z}(\varphi_n Z_m) Z_m dz, \quad \Lambda_{mn} \varphi_n = \int_{-h}^{\eta} \Delta_{X,z}(\varphi_n Z_n) \partial_{\eta} Z_m dz, \quad (32a,b) $$

for $m = 0, \dots, K$. Eqs. (29)–(31) are valid for any $\Phi^a$ of the form (24) with (25) and (28) satisfied. For the present choice (24) we also have the following

**Lemma 3.1.** If $\varphi_n$, $n = 1, \dots, K$ satisfy (30) with $Z_n$ as in (24), then

$$ \Delta \psi (\partial_{\eta} Z_m, 1) \varphi_m + \Lambda_{mn} \varphi_n \varphi_m = (H \Delta \psi + L_{0n} \varphi_n) H \varphi_1. \quad (33) $$

The proof of Lemma 3.1 is detailed in Appendix C. In view of Lemma 3.1, it remains to compute the coefficients of $L_{mn}$ for $m \ge 0, n \ge 1$ from the integrals in the left column of (18). The derivation is further simplified by expressing the integrals in $B_{mn}$ and $C_{mn}$ in terms of $A_{mn}$ (see Appendix D). The final result can be written as

$$ \partial_t \eta + \nabla \eta \cdot \nabla \psi + (\nabla \eta)^2 H \varphi_1 + H \varphi_1 + H \Delta \psi + L_{0n} \varphi_n = 0, \quad (34) $$

$$ (m = 1, \dots, K) \quad L_{mn} \varphi_n = \zeta_m H^{2m+1} \Delta \psi, \quad (35) $$

$$ \partial_t \psi + g \eta + \frac{1}{2}(\nabla \psi)^2 - \frac{1}{2}((\nabla \eta)^2 + 1)(\varphi_1 H)^2 - (H \Delta \psi + L_{0n} \varphi_n) H \varphi_1 = 0, \quad (36) $$

where for $m = 0$ and $n = 1, \dots, K$,

$$ L_{0n} = 2^n B(n+1) H^{2n+1} \Delta + 2^{n+1} B(2n+1) H^{2n} \nabla \eta \cdot \nabla + 2^n (2n+1) B(n+1) (2n(\nabla \eta)^2 + H \Delta \eta) - \delta_{1n} n H - n H \delta_{1n} (\nabla \eta)^2, $$
where $\delta_{11} = 1$, $\delta_{1n} = 0$ for $n > 1$ is the Kronecker delta, and

$$ B(k) \equiv B(k, k) = \frac{(k-1)!(k-1)!}{(2k-1)!} \qquad (37) $$

is the Beta function for equal positive integers. For $m, n = 1, \dots, K$,

$$ L_{mn} = \alpha_{mn} H^{2(m+n)+1} \Delta + \beta_{mn} H^{2(m+n)} \nabla\eta \cdot \nabla + \kappa_{mn} H^{2(m+n)-1} ((\nabla\eta)^2 + \lambda_{mn} H \Delta\eta + \mu_{mn}), \qquad (38) $$

where the Greek-letter coefficients are quotient numbers depending only on $m$ and $n$:

$$ \alpha_{mn} = 2^{m+n} B(m+n+1), \quad \beta_{mn} = 2^{m+n} n B(m+n), $$

$$ \kappa_{mn} = 2^{m+n-1} n \left( \frac{m}{2} B(m+n-1) + 2(m+n)B(m+n) \right), $$

$$ \lambda_{mn} = 2^{m+n-1} n B(m+n), $$

$$ \mu_{mn} = \frac{2^{m+n-2} mn}{1-2(m+n)} B(m+n-1), $$

$$ \zeta_m = -2^m B(m+1). $$

For $K=0$, Eqs. (34)–(36) are just the non-dispersive non-linear shallow water (or Saint-Venant) system in terms of $(\eta, \psi)$:

$$ \partial_t \eta + \nabla \cdot (H \nabla \psi) = 0, \quad \partial_t \psi + g\eta + \frac{1}{2} (\nabla \psi)^2 = 0. $$

The first non-trivial member of (34)–(36) corresponds to $K=1$ in (24) and, after straightforward manipulations, is written as

$$ \partial_t \eta + \nabla \cdot \left[ H \nabla \psi + \frac{1}{3} \nabla (H^3 \varphi_1) \right] = 0, \qquad (39a) $$

$$ \partial_t \psi + g\eta + \frac{1}{2} (\nabla \psi)^2 + \frac{1}{2} (1 + (\nabla \eta)^2) (H \varphi_1)^2 - \left[ H \Delta \psi + \frac{1}{3} \Delta (H^3 \varphi_1) \right] H \varphi_1 = 0, \qquad (39b) $$

$$ \frac{2}{5} H^2 \Delta \varphi_1 + 2 H \nabla \eta \cdot \nabla \varphi_1 + (H \Delta H + (\nabla \eta)^2 - 1) \varphi_1 = -\Delta \psi. \qquad (39c) $$

The structure of the K-family (34)–(36) is reminiscent of the ZCS formulation; $(\eta, \psi)$ are the evolving variables and $\varphi_n, n = 1, \dots, K$ are determined from $(\eta, \psi)$ through the system of spatial PDEs (35). In fact, these equations possess a Hamiltonian structure which we will address in the next section. For $K=1$, we will refer to the (K = 1)-model, and more generally, the (K = K*)-model for arbitrary $K^* \ge 1$. The (K = 1)-model, Eqs. (39a)–(39c), will be used for comparisons with asymptotic models in Section 6 and for numerical simulations in Section 8.

## 4. Hamiltonian structure of the model equations

We show here formally that the evolution equations obtained from Luke's variational principle using an ansatz of the form (24) satisfying the conditions (25), i.e. Eqs. (29) and (31), are obtained from a Hamiltonian functional defined in terms of $(\eta, \psi)$ through the system of linear spatial PDEs (30). In other words, the Hamiltonian structure of the variationally derived model equations remains non-local. However, rather than solving a DtN problem for the Laplace equation (Eq. (6a,b,c)), the evolution requires solving a system of linear PDEs in the horizontal plane, Eqs. (30). Let L be the matrix differential operator $L = -(L_{mn})_{m,n=1,...,K}$, associated with (30), acting on column vector functions $u = (u_1, \dots, u_K)^t$. Then, in view of (32a,b), Green's identity and (25), we have

$$ \begin{aligned} \iint L u u^t \, dX &= - \iint_{-h}^{\eta} \Delta_{X,z}(u_n Z_n) (u_m Z_m) \, dz \, dX \\ &= \iint_{-h}^{\eta} [\nabla(u_n Z_n)]^2 \, dz \, dX \\ &= \iint_{-h}^{\eta} [ (\nabla(u_n Z_n))^2 + (u_n \partial_z Z_n)^2 ] \, dz \, dX \\ &= \iint_{-h}^{\eta} [ (\nabla u_n Z_n + u_n \partial_\eta Z_n \nabla \eta)^2 + (u_n \partial_z Z_n)^2 ] \, dz \, dX \end{aligned} \qquad (40) $$

Following Nemoto and Iguchi [67, Lemma 3.1] (see also Duchêne and Iguchi [80, Lemma 6.1]) it can be shown that L is coercive for sufficiently smooth bounded η provided that $\{Z_n, \partial_\eta Z_n\}_{n=1}^K$ and $\{\partial_z Z_n\}_{n=1}^K$ are linearly independent. This requirement is satisfied by the present choice of $Z_n$ in (24) as well as by the ones appearing in [71] for instance. Standard elliptic PDE theory ensures the classical solvability of the equation $Lu = g$ for sufficiently smooth right-hand sides (RHS) g. In the case of Eq. (30), the right-hand side takes the form $g^t = ((Z_m, 1)\Delta\psi)_{m=1}^K$. Consequently, given smooth $(\eta, \psi)$ the system (30) is solvable and in analogy to Section 2.1 we may introduce the mapping $\psi \to (f_n[\eta]\psi)_{n=1,...,K} = (\varphi_n)_{n=1,...,K}$, that maps $\psi$ to the solution of (30), and the corresponding approximate velocity potential

$$ (\Phi^a[\eta]\psi)(X, z, t) = \Phi^a(X, z, t) = \psi(X, t) + (f_n[\eta]\psi)(X, t) Z_n(z; \eta(X, t)). \qquad (41) $$

Using the identity (32a), the system (30) takes the form

$$ \int_{-h}^{\eta} \Delta_{X,z}(\Phi^a[\eta]\psi) Z_m \, dz = 0, \quad m = 1, \dots, K. \qquad (42) $$

It is also convenient to rewrite Eqs. (29) and (31) as

$$ \partial_t \eta = N_\eta \cdot [\nabla(\Phi^a[\eta]\psi)]_\eta - \int_{-h}^{\eta} \Delta_{X,z}(\Phi^a[\eta]\psi) \, dz \stackrel{\text{def}}{=} G^a[\eta]\psi, \qquad (43a) $$

$$ \partial_t \psi = -g\eta - \frac{(\nabla\psi)^2}{2} + \frac{(\nabla\eta)^2 + 1}{2} [\partial_z\Phi^a[\eta]\psi]_\eta^2 + \left( \int_{-h}^{\eta} \Delta_{X,z}(\Phi^a[\eta]\psi) \partial_\eta Z_n \, dz \right) f_n, \qquad (43b) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1. Schematic representation of the derivation process and summary of important equations.</figcaption></figure>

where we have introduced the approximate “DtN operator” $G^a[\eta]$. The Hamiltonian is obtained by replacing $\Phi$ with $\Phi^a$ in Zakharov’s Hamiltonian (7),

$$H^a(\eta, \psi) \stackrel{\text{def}}{=} \frac{1}{2} \int \left\{ \int_{-h}^{\eta} (\nabla \Phi^a[\eta] \psi)^2 dz + g\eta^2 \right\} dX = \frac{1}{2} \int (\psi G^a[\eta] \psi + g\eta^2) dX, \quad (44)$$

where the second equality is obtained by Green’s identity, Eq. (42), the fact that $[\Phi^a[\eta]\psi]_\eta = \psi$ and Eq. (43a). In a similar way, we see that

$$\int G^a[\eta] \psi_1 \psi_2 \, dX = \int \int_{-h}^{\eta} (\nabla \Phi^a[\eta] \psi_1) \cdot (\nabla \Phi^a[\eta] \psi_2) \, dz \, dX = \int \psi_1 G^a[\eta] \psi_2 \, dX, \quad (45)$$

for any $\psi_i$, $i = 1, 2$, where $\Phi^a[\eta]\psi_i$ satisfies (42) with $\psi = \psi_i$. Hence, $G^a[\eta]$ retains the symmetry and the positive definiteness properties of the original DtN operator, as expected. We can now compute $\delta_\psi H^a$ and $\delta_\eta H^a$ to show that the system (43) (or equivalently Eqs. (29) and (31)) is expressed in the equivalent Hamiltonian form

$$\partial_t \eta = \delta_\psi H^a, \quad (46a)$$

$$\partial_t \psi = -\delta_\eta H^a, \quad (46b)$$

where $H^a$ is given by (44) with $\Phi^a = \psi + \phi_n Z_n$, where $\phi_n$, $n = 1, \dots, K$ solve (30) and $Z_n$ satisfy (25). These calculations are provided in the Appendix E to maintain the flow of the paper.

Obviously, the solution of (46) conserves the Hamiltonian $H^a$, i.e. $dH^a/dt = 0$. Moreover, it conserves the (excess of) mass and the impulse in the direction $x_i$, $i = 1, 2$,

$$\frac{dM^a}{dt} = 0, \quad M^a(t) = \int \eta \, dX, \quad (47)$$

$$\frac{dI_i^a}{dt} = 0, \quad I_i^a(t) = \int \eta \, \partial_i \psi \, dX. \quad (48)$$

To verify (47), we note that $G^a[\eta]1 = 0$ which together with (46a) and the $L^2$-symmetry of $G^a[\eta]$ gives (47). To verify (48), we differentiate the Hamiltonian density (the bracketed term in (44)) and use Green’s identity and Eqs. (25), (26) and (42) to find

$$\begin{aligned} \int \frac{1}{2} \partial_i \left\{ \int_{-h}^{\eta} (\nabla \Phi^a[\eta] \psi)^2 dz + g\eta^2 \right\} dX &= \int \left\{ \frac{1}{2} \partial_i \eta [\nabla \Phi^a]^2_{\eta} + \int_{-h}^{\eta} \nabla \Phi^a \cdot \nabla (\partial_i \Phi^a) dz + g\eta \partial_i \eta \right\} dX \\ &= \int \left\{ \partial_\psi H^a \partial_i \psi + \partial_\eta H^a \partial_i \eta \right\} dX. \end{aligned}$$

Moreover, from (46) we have

$$\frac{dI_i^a}{dt} = \int \partial_i \eta \partial_i \psi \, dX + \int \eta \partial_i \partial_i \psi \, dX = \int (\partial_\psi H^a \partial_i \psi + \partial_\eta H^a \partial_i \eta) \, dX,$$

and $dI_i/dt = 0$ follows from the vanishing at infinity conditions.

Note that in the derivation above, the explicit calculation of the Hamiltonian (44) was not required. However, for future reference and for the evaluation of our numerical solutions, we give here the Hamiltonian of the system (39a)–(39b),

$$H^a(\eta, \psi) = \frac{1}{2} \int \left\{ \left[ H \nabla \psi + \frac{1}{3} \nabla (H^3 f_1[\eta] \psi) \right] \cdot \nabla \psi + g\eta^2 \right\} dX. \quad (49)$$

where $f_1[\eta]\psi$ solves (39c). To summarize and for easy reference, Fig. 1 presents a schematic outlining our derivation process and the most important equations.

## 5. Relation to the Isobe-Kakinuma model

As mentioned in the introduction, Isobe [60] used the vertical polynomials of the classical Boussinesq-Rayleigh expansion to derive model equations from Luke’s variational principle:

$$\Phi^a_*(X, z, t) = \sum_{n=0}^{K} \phi_n(X, t) Z_n^*(z; h), \quad \text{with} \quad Z_n^*(z; h) = (z+h)^{2n}. \quad (50)$$

In this case, the EL Eqs. (14)–(15) lead to $K + 1$ evolution equations on $\eta$ and one evolution equation on $(\phi_0, \phi_1, \dots, \phi^K)$. These equations are known as the Isobe-Kakinuma (IK) model. The representations (24) and (50) are inspired by two different asymptotic expressions, (21) and (23), respectively. Below, we establish the equivalence of these two approaches and explore their differences.
To simplify the presentation, we proceed with a direct comparison of the evolution equations in the case $K = 1$ and defer the calculations for the case $K > 1$ to Appendix F. With $Z_n = Z_n^*$ as in (50), the EL Eqs. (14)–(15) yield

$$ \partial_t \eta + \nabla \eta \cdot (\nabla \phi_0 + H^2 \nabla \phi_1) + H \Delta \phi_0 + \frac{H^3}{3} \Delta \phi_1 = 0, \qquad (51a) $$

$$ \partial_t \eta + \nabla \eta \cdot (\nabla \phi_0 + H^2 \nabla \phi_1) + \frac{H}{3} \Delta \phi_0 + \frac{H^3}{5} \Delta \phi_1 - \frac{4}{3} H \phi_1 = 0, \qquad (51b) $$

$$ \partial_t \phi_0 + \partial_t \phi_1 H^2 + g\eta + \frac{1}{2}(\nabla \phi_0)^2 + H^2 \nabla \phi_0 \cdot \nabla \phi_1 + \frac{1}{2} H^4 (\nabla \phi_1)^2 + 2H^2 \phi_1^2 = 0. \qquad (51c) $$

We will show that this system is equivalent to (39a)–(39c). Eqs. (51a)–(51c) can be written as a Hamiltonian system after eliminating $\partial_t \eta$ from (51a)–(51b) and introducing the trace of (50) on $z = \eta$, $\psi = \phi_0 + H^2 \phi_1$ [69]. This leads to two evolution equations on $(\eta, \psi)$ coupled with a system of equations on $(\phi_0, \phi_1)$,

$$ \partial_t \eta + \nabla \eta \cdot (\nabla \phi_0 + H^2 \nabla \phi_1) + H \Delta \phi_0 + \frac{H^3}{3} \Delta \phi_1 = 0, \qquad (52a) $$

$$ \partial_t \psi - 2H \partial_t \eta \phi_1 + g\eta + \frac{1}{2}(\nabla \phi_0)^2 + H^2 \nabla \phi_0 \cdot \nabla \phi_1 + \frac{1}{2} H^4 (\nabla \phi_1)^2 + 2H^2 \phi_1^2 = 0, \qquad (52b) $$

$$ \frac{2}{3} \Delta \phi_0 + \frac{2}{15} H^2 \Delta \phi_1 + \frac{4}{3} \phi_1 = 0, \qquad (52c) $$

$$ \phi_0 + H^2 \phi_1 = \psi. \qquad (52d) $$

Substituting $\phi_0$ from (52d) into (52c), we obtain (39c) with $\phi_1 = -2\phi_1$. Then, from (52d) we obtain that $\phi_0 = \psi + (H^2/2)\phi_1$. Thus, we have expressed $\phi_0$ and $\phi_1$ in terms of $\psi$ and $\phi_1$. Plugging these expressions into (52a), we easily obtain (39a). Similarly, from (52b) we obtain (39b) using additionally (52a). In a similar way, we can derive Eqs. (52a)–(52d) from Eqs. (39a)–(39c). Therefore, Eqs. (51a)–(39a) and (51c)–(39c) are equivalent.

From a numerical point of view, the formulation (51a)–(51c) is the less straightforward to implement. We refer to Abohadima and Isobe [70] for an implicit second-order scheme that has been proposed (see also Kakinuma [61]). In contrast, for the Hamiltonian formulations (52a)–(52d) or (39a)–(39c), a scheme for their solution can be readily devised. Indeed, starting from an initial condition $(\eta_0, \psi_0)$ the evolution Eqs. (52a)–(52b) (resp. (39a)–(39b)) can be marched in time provided that the spatial linear PDE (52d) (resp. (39c)) is solved at every time step to calculate $(\phi_0, \phi_1)$ (resp. $\phi_1$) in terms of $(\eta, \psi)$ in the evolution equations. In Section 8, we solve (39a)–(39c).

Note however that (39c) contains second-order derivatives on $\eta$ which are absent in (52d). Moreover, in contrast to (51a)–(51c), the evolution Eqs. (39a)–(39b) are fully non-linear and non-local in the sense that they contain non-linear terms with second-order derivatives of the solution of (39c) and this would obstruct their theoretical analysis.

From this equivalence, the local-in-time well-posedness proved for the initial value problem for the IK model [65,67] and its justification as a high-order shallow-water model [67,78] could be transferred to (34)–(36). The interested reader can find detailed and rigorous statements and conditions in these works. Here, we only highlight that the IK model with $K$ vertical functions is a $O(\mu^{2K+1})$ approximation of the water wave equations in the case of a flat bottom, and that its dispersion relation is a $[2K, 2K]$-Padé approximant of the dispersion relation of the linearized water wave equations (see Section 7.1). Our ($K = 1$)–model encompasses any of the reformulations (39a)–(39c) or (51a)–(51c) or (52a)–(52d); for the subsequent discussion and applications, we shall use (39a)–(39c).

# 6. Comparison with asymptotic models

In this section, we show that the ($K = 1$)–model, (39a)–(39c), represents a higher-order shallow-water approximation of the original problem and we compare it with models derived from asymptotic approximations.

## 6.1. Dimensionless expression of the ($K = 1$)–model

We begin by introducing the dimensionless version of (39). Since $\Phi^a$ and $\psi$ scale as $L\sqrt{gh_0}$ and $Z$ as $h_0^2$, we use the scalings $\Phi^a = \Phi^{a'} L\sqrt{gh_0}$, $\psi = \psi' L\sqrt{gh_0}$ and $\phi_1 = \phi_1' \sqrt{gh_0}/L$ together with (A.1) so that, after dropping primes,

$$ \Phi^a = \psi + \mu \phi_1 \tilde{Z}, \qquad (53) $$

where $\tilde{Z}$ is given in Proposition 3.1. Similarly, the equations in (39) become

$$ \partial_t \eta + \nabla \cdot \left[ H \nabla \psi + \mu \frac{1}{3} \nabla (H^3 \phi_1) \right] = 0, \qquad (54a) $$

$$ \partial_t \psi + \eta + \frac{(\nabla \psi)^2}{2} + \frac{\mu}{2} (\mu (\nabla \eta)^2 + 1) (H \phi_1)^2 - \mu \left[ H \Delta \psi + \frac{1}{3} \mu \Delta (H^3 \phi_1) \right] H \phi_1 = 0, \qquad (54b) $$

$$ \mu \mathcal{L} \phi_1 - \phi_1 + \Delta \psi = 0, \qquad (54c) $$

with $H = \eta + 1$ and

$$ \mathcal{L} = \mathcal{L}(\eta) = \frac{2}{5} H^2 \Delta + 2H \nabla \eta \cdot \nabla + \frac{1}{2} \Delta(H^2). \qquad (55) $$

The associated Hamiltonian (49) takes the form

$$ \mathcal{H}^a = \mathcal{H}^a(\eta, \psi) = \frac{1}{2} \int \left[ \int_{-1}^{\eta} \frac{1}{\mu} (\nabla_{\mu} \Phi^a)^2 dz + \eta^2 \right] dX, \quad \text{with} \quad \nabla_{\mu} \stackrel{\text{def}}{=} (\sqrt{\mu} \nabla, \partial_z), \qquad (56) $$

where $\Phi^a$ is given by (53) with $\phi_1$ solving the PDE (54c) with data $(\eta, \psi)$. Zakharov’s dimensionless Hamiltonian is given by

$$ \mathcal{H} = \mathcal{H}(\eta, \psi) = \frac{1}{2} \int \left[ \int_{-1}^{\eta} \frac{1}{\mu} (\nabla_{\mu} \Phi)^2 dz + \eta^2 \right] dX, \qquad (57) $$
where $\Phi$ solves (20) with data $(\eta, \psi)$. We show below that $|H - H^a| = O(\mu^3)$. While this result follows from Duchêne and Iguchi [69]'s work and the equivalence demonstrated in Section 5, we provide a brief formal proof within the present framework for completeness. Moreover, this approach yields a useful expression that facilitates comparisons with other asymptotic models. Using (54c) we can express $\varphi_1$ in terms of $(\eta, \psi)$ asymptotically as

$$ \varphi_1 = \Delta\psi + \mu\mathcal{L}\Delta\psi + O(\mu^2). \qquad (58) $$

From (53) and (58) we have that

$$ \Phi^a = \psi + (\mu\Delta\psi + \mu^2\mathcal{L}\Delta\psi)\tilde{Z} + O(\mu^3), \qquad (59) $$

consequently, it follows from (21) in Proposition 3.1 that

$$ \Phi - \Phi^a = O(\mu^2). \qquad (60) $$

Then,

$$ \begin{aligned} H - H^a &= \frac{1}{2\mu} \int \int_{-1}^{\eta} (\nabla_{\mu} \Phi)^2 \mathrm{d}z \mathrm{d}X - \frac{1}{2\mu} \int \int_{-1}^{\eta} (\nabla_{\mu} \Phi^a + \nabla_{\mu} \Phi - \nabla_{\mu} \Phi)^2 \mathrm{d}z \mathrm{d}X \\ &= -\frac{1}{2\mu} \int \int_{-1}^{\eta} (\nabla_{\mu} (\Phi^a - \Phi))^2 \mathrm{d}z \mathrm{d}X - \frac{1}{\mu} \int \int_{-1}^{\eta} \nabla_{\mu} \Phi \cdot \nabla_{\mu} (\Phi^a - \Phi) \mathrm{d}z \mathrm{d}X. \end{aligned} $$

The second term in the last RHS expression vanishes because $\Phi$ solves (20) and $\Phi^a - \Phi$ vanishes on the free surface. From the remaining first term and (60) we obtain that $|H - H^a| = O(\mu^3)$.

We will now use (58) to reformulate the dimensionless version of the (K = 1)-model, Eqs. (54a)-(54b), to explore its relation to asymptotically derived models.

## 6.2. Relation to a model from Whitham [81]

Plugging (58) into (54a)-(54b), the (K = 1)-model transforms to

$$ \partial_t \eta + \nabla \cdot \left[ H \nabla \psi + \mu \frac{1}{3} \nabla (H^3 \Delta \psi) + \mu^2 \frac{1}{3} \nabla (H^3 \mathcal{L} \Delta \psi) \right] = 0, \qquad (61a) $$

$$ \partial_t \psi + \eta + \frac{1}{2}(\nabla \psi)^2 - \mu \frac{1}{2}(H \Delta \psi)^2 + \mu^2 \left[ \frac{1}{2}(\nabla H)^2 H^2 \Delta \psi - \frac{1}{3}\Delta(H^3 \Delta \psi)H \Delta \psi \right] = 0, \qquad (61b) $$

after neglecting the $O(\mu^3)$ terms. This is a $O(\mu^2)$-accurate system that is also derived from the dimensionless version of the ZCS formulation (see Lannes and Bonneton [9] for instance) and the asymptotic expansion of the dimensionless DtN operator $G_{\mu}[\eta]\psi$ obtained from the relation $G_{\mu}[\eta]\psi = -\mu\nabla \cdot (\int_{-1}^{\eta} \nabla\Phi dz)$ and (21) after a simple computation:

$$ G_{\mu}[\eta]\psi = -\mu\nabla \cdot \left[ H\nabla\psi + \mu\frac{1}{3}\nabla(H^3\Delta\psi) + \mu^2\frac{1}{3}\nabla(H^3\mathcal{L}\Delta\psi) \right] + O(\mu^4). \qquad (62) $$

The $O(\mu^2)$-accurate component of the above expression is the flat-bottom counterpart to the variable-bottom expression derived in [14, §3.6.1] (see also Clamond [82] for an alternative derivation). Eq. (62) also follows by using (59) into the definition of the approximate DtN operator $G^a[\eta]\psi$, Eq. (43a). The system (61) has the favourable property of involving only spatial differential operators. In contrast, Eqs. (54a)-(54b) contain the non-local term $\varphi_1 = \varphi_1[\eta]\psi$ which is determined by the solution of the PDE (54c). However, this advantage does not come without a cost, as we shall see.

First, let us remark that neglecting the $O(\mu^2)$ terms in (61) leads to the $O(\mu)$-accurate model originally derived by Whitham [81] and later by Miles and Salmon [7] (see also Clamond and Dutykh [10, Section 3.4]). The Hamiltonian of this system is obtained from Zakharov's Hamiltonian (57) with $\Phi$ as in (21) after neglecting $O(\mu^2)$ terms,

$$ H^{(1)} = \frac{1}{2} \int (\psi T(\eta)\psi + \eta^2) \mathrm{d}X \quad \text{with} \quad T(\eta)\psi = -\nabla \cdot \left[ H \nabla \psi + \mu \frac{1}{3} \nabla (H^3 \Delta \psi) \right]. \qquad (63) $$

The truncated kinetic energy $\frac{1}{2} \int \psi T(\eta)\psi \mathrm{d}X$ in (63) is not the integral of a square and although $H^{(1)}$ is formally a $O(\mu^2)$ approximation of $H$ (i.e. $|H - H^{(1)}| = O(\mu^2)$), it is not positive definite in general, which is a nonphysical feature. In particular, it can be shown that the system (61) without the $O(\mu^2)$ terms is linearly ill-posed in the sense that the frequency of its plane waves of wave number $k$ over a flat bottom $h$ becomes imaginary. This occurs for $kh > \sqrt{3}$, resulting in an instability that is absent in the original problem. Such an effect has been documented in earlier works in relation to truncated asymptotic expansions of the ZCS Hamiltonian [34–36]. On the other hand, the full $O(\mu^2)$ system (61) leads to a real linear frequency that quickly diverges from the exact one away from $kh = 0$ and tends to infinity as $kh$ tends to infinity (not shown here). Thus, asymptotically derived models of this type may be of limited applicability (see Fig. 2 in Section 7.1). It is important to note that the Hamiltonian (44) in the present variational approach is positive definite for any $(\eta, \psi)$. Additionally, the associated linear dispersion relation is non-singular [67,71].

## 6.3. Relation to the SGN model

A useful and well-known $O(\mu)$-accurate model is obtained from (61) by introducing the asymptotic expansion of the depth-averaged horizontal velocity using (21),

$$ \bar{V} = \frac{1}{H} \int_{-1}^{\eta} \nabla (\psi + \mu \Delta \psi \tilde{Z}) \mathrm{d}X + O(\mu^2) = \nabla \psi + \mu \frac{1}{3H} \nabla (H^3 \nabla \cdot \nabla \psi) + O(\mu^2), \qquad (64) $$

and inverting (64) to find [9]

$$ \nabla \psi = \bar{V} - \mu \frac{1}{3H} \nabla (H^3 \nabla \cdot \bar{V}) + O(\mu^2). \qquad (65) $$
Then, the equations in (61) are transformed into (after taking the gradient of (61b))

$$ \partial_t \eta + \nabla \cdot (H \bar{V}) = 0, \qquad (66a) $$

$$ \partial_t \bar{V}^* + \nabla \left[ \eta + \frac{1}{2} \bar{V}^2 - \mu \left( \frac{1}{3H} \bar{V} \cdot \nabla (H^3 \nabla \cdot \bar{V}) + \frac{1}{2} H^2 (\nabla \cdot \bar{V})^2 \right) \right] = 0, \qquad (66b) $$

$$ \bar{V}^* = \bar{V} - \mu \frac{1}{3H} \nabla(H^3 \nabla \cdot \bar{V}), \qquad (66c) $$

upon dropping the $O(\mu^2)$ terms. The above system is a variant of the SGN model appearing in Choi [33], it is asymptotically equivalent to the other variants of the $O(\mu)$-accurate SGN model [9,33], and is linearly well-posed. As in the (K = 1)-model (54a)-(54b), the time evolution of (66a)-(66b) requires the inversion of the elliptic operator in (66c). Therefore, Eqs. (54a)-(54b), being $O(\mu^2)$-accurate, exhibit better precision at a similar computational complexity. The next in order SGN model, which exhibits the same precision as Eqs. (54a)-(54b), is linearly ill-posed and will not be further discussed here. We refer to Matsuno [30], Matsuno [31] for its detailed derivation and study of its Hamiltonian structure, and to Fig. 2 in the next section for a visual comparison of its dispersion relation with other models.

An additional noteworthy distinction between the SGN model and the (K = 1)-model (54a)-(54c) pertains to the irrotationality and incompressibility of the modelled velocity fields. Suppose that we have obtained $\eta$ and $\bar{V}$ by solving (66a)-(66c). Then, the reconstructed velocity field

$$ U = (\sqrt{\mu} V, w) \text{ with } V(X, z) = \bar{V} \text{ and } w(X, z) = -\mu(z+1)\nabla \cdot \bar{V}, $$

satisfies $\nabla_\mu \cdot U = 0$ and $|\nabla_\mu \times U| = O(\mu)$. Thus, the flow is exactly incompressible but weakly rotational. On the other hand, the velocity field reconstructed from the solution $(\eta, \psi)$ of the (K = 1)-model, Eqs. (54a)-(54c), is

$$ U = (\sqrt{\mu} V, w) \text{ with } V(X, z) = \nabla\psi + \mu\nabla(\phi_1 \bar{Z}) \text{ and } w(X, z) = -\mu(z+1)\phi_1. $$

From the above equation, we have $|\nabla_\mu \times U| = 0$ and, using (58), we see that $\nabla_\mu \cdot U = O(\mu^2)$. Thus the flow is weakly compressible but exactly irrotational.

## 6.4. Relation to Choi [33]'s model

In a recent paper, Choi [33] derived an interesting new class of asymptotic models that are linearly well-posed at any order of $\mu$. In this class, the evolving variables are $\eta$ and the bottom-surface potential, $f$. To see how this class relates to the (K = 1)-model (54a)-(54b) we can use the relation

$$ \psi = f - \mu \frac{H^2}{2} \Delta f + \mu^2 \frac{H^4}{24} \Delta^2 f + O(\mu^3), \qquad (67) $$

obtained from (23), to express $\phi_1$ in terms of $f$. Plugging (67) into (58) we obtain after a simple calculation that

$$ \phi_1 = \Delta f - \mu \frac{1}{10} H^2 \Delta f + O(\mu^2). \qquad (68) $$

Using (67) and (68), the evolution Eqs. (54a)-(54b) take the form

$$ \partial_t \eta + \nabla \cdot \left[ H \nabla f - \mu \frac{1}{6} H^3 \nabla \Delta f + \mu^2 \frac{1}{120} H^5 \nabla^2 f \right] = 0, \qquad (69a) $$

$$ \partial_t f^* + \eta + \frac{1}{2}(\nabla f)^2 - \mu \nabla \cdot \left( \frac{H^2}{2} \Delta f \nabla f \right) + \mu^2 \nabla \cdot \left[ \frac{H^4}{24} \Delta^2 f \nabla f + \frac{H^4}{16} \nabla (\Delta f)^2 \right] = 0, \qquad (69b) $$

$$ f^* = f - \mu \frac{H^2}{2} \Delta f + \mu^2 \frac{H^4}{24} \Delta^2 f, \qquad (69c) $$

after neglecting the $O(\mu^3)$ terms. This is the $O(\mu^2)$-accurate model derived and studied in [33] (see also Choi [83]). Both systems (54a)-(54c) and (69a)-(69c) require the solution of an elliptic spatial PDE throughout the time evolution. The maximum order of spatial derivatives is two for (54a)-(54c) and six for (69a)-(69b). Moreover, the elliptic PDE (54c) is of second order while (69c) is of fourth order. Increasing the precision of the (K = 1)-model (54a)-(54c) by adding more vertical functions will increase the unknowns in (54c) while the order of derivatives will remain two (see (34)-(36)). On the other hand, increasing the precision in (69a)-(69c) by keeping more terms in the relevant asymptotic expansions, (69c) will remain one PDE but the order of derivatives in (69a)-(69c) will increase. Finally, as pointed out in Choi [33] the truncated asymptotic expansion of the total energy associated to (69a)-(69c) is asymptotically conserved, whereas Eqs. (54a)-(54c) conserve the approximate Hamiltonian (49) exactly.

# 7. Discussion of linear properties of the model

## 7.1. Dispersion relations

In this section, we consider linear plane waves. Nemoto and Iguchi [67] proved that the phase speed of the IK model, $C_{IK}$, is a high-order approximation of the phase speed of the linearized water wave (WW) equations, $C_{WW}$,

$$ \frac{C_{WW}^2}{gh} = \frac{(\omega/k)^2}{gh} = \frac{\tanh(kh)}{kh}. \qquad (70) $$

In particular, when even functions are used for the approximation of $\Phi$, $C_{IK}^2$ is the [2K, 2K]-Padé approximant of $C_{WW}^2$. This fact remains true for Eqs. (34)-(36) due to the equivalence of the two formulations established in Section 5. To verify this for K = 1, we linearize the (K = 1)-model, Eqs. (39a)-(39c), around $(\eta, \psi) = (0, 0)$, to get

$$ \partial_t \eta + \nabla \cdot \left( h \nabla \psi + \frac{h^3}{3} \nabla \phi_1 \right) = 0, \qquad (71a) $$
$$ \partial_t \psi + g\eta = 0, \qquad (71b) $$

$$ \frac{2}{5}h^2\Delta\varphi_1 - \varphi_1 = -\Delta\psi. \qquad (71c) $$

Assuming without loss of generality, plane waves propagating along the x axis of the form $(\eta, \psi) = (\bar{\eta}, \bar{\psi})e^{i(kx-\omega t)}$, then we can easily solve Eq. (71c),

$$ \varphi_1 = - \frac{k^2}{1 + \frac{2}{5}(kh)^2} \bar{\psi} e^{i(kx-\omega t)}, \qquad (72) $$

and obtain from the evolution Eqs. (71a)–(71b) the dispersion relation

$$ \frac{C^2}{gh} = \frac{1 + \frac{1}{15}(kh)^2}{1 + \frac{2}{5}(kh)^2}. \qquad (73) $$

Eq. (73) is indeed the [2, 2]–Padé approximant of $C_{WW}^2$ and coincides with the dispersion relation of the improved SGN equations in [18, Section 2.6] with $a = 6/5$, and in [84] with $\beta = 2/15$. It also coincides with the dispersion relation of the parabolic model in [71] which has a very similar structure with (39a)–(39c) concerning their non-linear terms. The key difference is that the vertical functions in the present approach are chosen from an asymptotic expansion; therefore the inclusion of more functions and the refinement of the approximation is straightforward.

For the linearized ($K=2$)–model in (34)–(36), one can verify that its dispersion relation is the [4, 4]–Padé approximant of $C_{WW}^2$

$$ \frac{C^2}{gh} = \frac{1 + \frac{1}{9}(kh)^2 + \frac{1}{945}(kh)^4}{1 + \frac{4}{9}(kh)^2 + \frac{1}{63}(kh)^4}, \qquad (74) $$

which coincides with the dispersion relation of the model in [26] (abbreviated as GK99) obtained by using the $O(\mu^2)$-accurate Boussinesq expansion.<fnref n="1" /> GK99 has the same precision as the ($K=1$)–model but has better linear dispersion. This is due to the use of an evolving variable defined in terms of the velocity potential at two different elevations and a free parameter. GK99 has the same linear dispersion as the ($K=2$)–model which is, however, $O(\mu^4)$ accurate. The non-linear terms in GK99 involve up to fourth-order mixed spatial derivatives and mixed space-time derivatives that require special treatment, we refer to Gobbi and Kirby [26], Gobbi et al. [85], Wei et al. [86] for more details and various numerical results. Another model that embodies the [4, 4]–Padé dispersion relation (74) is the one developed in [27]. This model is based on the time evolution of the ZCS formulation, namely Eq. (4a) and the gradient of (4b). Throughout the time evolution, a system of linear fourth-order PDEs, resulting from a perturbative solution of the DtN problem, needs to be solved to evaluate the vertical and horizontal derivatives of the velocity potential on the free surface. We recall that variationally derived models (Eqs. (14)–(15)) do not contain mixed spatial derivatives and their maximum order is two.

In Fig. 2, we compare (73) and (74) with the phase speed of the $O(\mu)$-accurate SGN system,

$$ \frac{C^2}{gh} = \frac{1}{1 + \frac{1}{3}(kh)^2}, \qquad (75) $$

the phase speeds of the $O(\mu^2)$ and $O(\mu^4)$ accurate systems of Choi (see Choi [33, Eq. 2.10]),

$$ \frac{C^2}{gh} = \frac{1 + \frac{1}{3!}(kh)^2 + \frac{1}{5!}(kh)^4}{1 + \frac{1}{2!}(kh)^2 + \frac{1}{4!}(kh)^4}, \quad \frac{C^2}{gh} = \frac{1 + \frac{1}{3!}(kh)^2 + \frac{1}{5!}(kh)^4 + \frac{1}{7!}(kh)^6 + \frac{1}{9!}(kh)^8}{1 + \frac{1}{2!}(kh)^2 + \frac{1}{4!}(kh)^4 + \frac{1}{6!}(kh)^6 + \frac{1}{8!}(kh)^8}, \qquad (76) $$

and the exact phase speed (70). For completeness, we have also included the phase speed of the $O(\mu^2)$-accurate SGN model and the ones of Eqs. (61) at $O(\mu)$ and $O(\mu^2)$ accuracy derived and discussed in Section 6.

Fig. 2(a) clearly shows that Choi's $O(\mu^2)$-model and the ($K=1$)-model provide a better approximation of $C_{WW}$ compared to the SGN systems and Eqs. (61). To make more detailed comparisons we consider the relative difference $C/C_{WW}-1$ in Fig. 2(b) and set arbitrarily a 2% relative error on $C_{WW}$ as the limit of applicability of a given model. First, we remark that the improvement achieved by increasing the precision in Eqs. (61) or in the SGN model is not significant and is limited at low wave numbers. The ($K=1$)-model has slightly better dispersion in comparison to Choi's $O(\mu^2)$ model; the ($K=1$)-model would be applicable up to $kh \approx 2.84$ while Choi's $O(\mu^2)$ model would be applicable up to $kh \approx 2.46$. At $O(\mu^4)$, these limits are significantly improved but their respective distance grows; $kh \approx 7.39$ for the ($K=2$)-model and $kh \approx 5.49$ for Choi's $O(\mu^4)$ model. It is evident that both the present approach and Choi's method show promise in deriving high-order models, particularly in terms of their linear dispersive properties.

## 7.2. Trough instability

Recently, Madsen and Fuhrman [37] unveiled a new type of instability, called “trough instability”. They showed that several fully non-linear Boussinesq-type models, including the model of Agnon et al. [27] and GK99, are prone to this instability, and they proposed techniques for its suppression. To expose this instability they introduced into the model equations a small-amplitude periodic wave with a below-zero mean free-surface elevation and showed that the resulting dispersion relation may become singular. Following this approach for the ($K=1$)-model in (39a)–(39c), we assume waves of the form

$$ (\eta, \psi) = (\delta h, -\delta h g t) + \epsilon(\bar{\eta}, \bar{\psi})e^{i(kx-\omega t)}, \quad \text{with } \epsilon \ll 1 \text{ and } \delta \in (-1, 0). \qquad (77) $$

Here, $\delta h$ represents the trough of the wave. The term $-\delta h g t$ is introduced in the expression for $\psi$ to account for a non-zero Bernoulli constant and does not alter the result; we could equivalently work with the gradient of (39b). Substituting (77) into (39c) and assuming that $\varphi_1 = \epsilon\varphi_1^{(1)}$, we easily find that

$$ \frac{2}{5}(h + \delta h)^2 \partial_x^2 \varphi_1^{(1)} - \varphi_1^{(1)} = -\partial_x^2 \psi. \qquad (78) $$

<footnote n="1">¹ Note that GK99 uses the notation $\mu = 2\pi h_0/L$ instead of $\mu = (h_0/L)^2$ used in the present paper</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 2. (a) Dimensionless phase speed $C/\sqrt{gh}$ as a function of $kh$ for Eqs. (61) at $O(\mu)$ (↓), Eqs. (61) at $O(\mu^2)$ (↓), the $O(\mu)$-accurate SGN model (←), the $O(\mu^2)$ (←) and $O(\mu^4)$ (←) accurate Choi's models, the $(K=1)$-model (←), the $(K=2)$-model (↓) and the exact linearized water wave (WW) equations, given by (70) (---). (b) Relative difference of the above wave speeds with respect to $C_{WW}$, defined as $C/C_{WW} - 1$.</figcaption></figure>

Using Eq. (72), it follows that

$$ \varphi_1 = \varepsilon \varphi_1^{(1)} = -\varepsilon \frac{k^2}{1 + \frac{2}{5}(k(h + \delta h))^2} \tilde{\psi} e^{i(kx - \omega t)}. \quad (79) $$

Using (77) and (79) in (39a)–(39b), we obtain at leading order the relation

$$ \frac{C^2}{gh} = \frac{1 + \frac{1}{15}(k(\delta h + h))^2}{1 + \frac{2}{5}(k(\delta h + h))^2}, $$

which shows that $C^2$ is always strictly positive and finite, thus proving that Eqs. (39a)–(39c) are free of the trough instability. A similar calculation, which we carried out for Choi [33]'s model (69a)–(69c), shows that such an instability does not occur for this model either. The same holds for the SGN model, as pointed out by Kirby [38].

The theoretical analysis presented so far demonstrates that the $(K=1)$-model possesses desirable properties, including high-order dispersion, exact Hamiltonian structure and absence of the wave trough instability. These characteristics establish a strong foundation for the model's applicability. To further assess its performance, we examine how the model behaves in non-linear regimes where analytical solutions are not available. To this end, the following section presents numerical simulations in the presence of strong non-linearities and complex wave interactions.

## 8. Numerical examples

In this section, we present some illustrative numerical solutions of the $(K=1)$-model (54a)–(54c) with $\mu=1$ in one spatial dimension (i.e. $\partial_1 \equiv \partial_x$ and $\partial_2 \equiv 0$) and periodic lateral boundary conditions. Our objective is to showcase that the $(K=1)$-model accurately propagates strongly non-linear solitary waves. This is possible using either Fourier or finite-difference (FD) methods, unlike the high-order model of Choi [33] which lacks this flexibility. Furthermore, we show that our model provides a more accurate representation of wave dynamics during solitary wave collisions compared to the SGN model. While these results highlight the advantages of the proposed model, a comprehensive numerical analysis is beyond the scope of this study and is left for future work.

### 8.1. Numerical methods

We test and compare two numerical methods for spatial discretization (hereafter abbreviated as A and B). For the time integration we use the classical four-step fourth-order explicit Runge-Kutta scheme with constant time-step $\delta t \sqrt{g/h_0} = 0.1$. The initial solitary waves are interpolated onto our computational grid covering $x \in [-L, L]$ using $L/h_0 = 200$, with their crests located at $x=0$ at $t=0$.

In method A, we solve the system (54a)–(54c) using fourth-order FD approximations for spatial derivatives on 5-node centred stencils. For a solitary wave in $[-L, L]$, the function $\psi$ (and $\phi_1$) is anti-periodic, thus we consider the periodic extension of $(\eta, \psi)$ in a domain of length $4L$. The spatial discretization is uniform, with $\delta x/h_0 = 2(L/h_0)/2048 = 0.1953125 \approx 0.2$. The discretization of the one-dimensional version of Eq. (54c) results in a cyclic pentadiagonal linear system which is solved efficiently with the algorithm of Neossi Nguetchue and Abelman [87]. No filtering nor smoothing is applied for all cases presented here.

The method B is a Fourier pseudo-spectral method [88] in the domain $[-L, L]$ for the system obtained from Eqs. (54a)–(54c) after introducing $v \equiv \psi_x$ and taking the $x$-derivative of Eq. (54b):

$$ \partial_t \eta = -\partial_x \left[ H \nu + \frac{1}{3} \partial_x (H^3 \varphi_1) \right], \quad (80a) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Evolution of solitary waves of the free-surface Euler equations. The initial amplitude and phase speed increase from left to right, $(a/h_0, C/\sqrt{gh_0}) = (0.2, 1.0943)$, $(0.4, 1.1780)$ and $(0.6, 1.2500)$. (a)-(c): Space-time plot of the free-surface elevation $\eta/h_0$ for $(x/h_0, t\sqrt{gh_0}) \in [-200, 200] \times [0, 120]$ using method A. (d)-(f) Free-surface elevation computed using method A (—) and method B (-), and initial condition shifted by $120C/\sqrt{gh_0}$ (e) The insets in panel (f) correspond to the boxed regions. (g)-(i) Relative errors on the conservation of the excess of mass $M^a$ and energy $H^a$, calculated by Eqs. (82).</figcaption></figure>

$$ \partial_t v = -\partial_x \eta - v \partial_x v - \partial_x \left\{ \frac{1}{2} (1 + (\partial_x \eta)^2) (H \varphi_1)^2 - \left[ H \partial_x v + \frac{1}{3} \partial_x^2 (H^3 \varphi_1) \right] H \varphi_1 \right\}, \quad (80b) $$

$$ \left[ \frac{2}{5} H^2 \partial_x^2 + 2H \partial_x \eta \partial_x + H \partial_x^2 \eta + (\partial_x \eta)^2 - 1 \right] \varphi_1 = -\partial_x v. \quad (80c) $$

Note that $v$ is not the horizontal velocity at the free surface, which is given in our model by $u(\eta) = \psi_x + H\varphi_1\partial_x\eta$. We decompose $\eta$ and $v$ into $N = 2048$ Fourier modes, e.g., $\eta(x,t) = \sum_{n=-N/2+1}^{N/2} \eta_n(t)e^{i\pi nx/L}$ and consider $N$ uniformly distributed grid points $x_i \in [-L, L]$, $i = 1, \dots, N$. To solve Eq. (80c) at each time step, we collocate the solution at $x_i$. The derivatives in the variable coefficients, $\partial_x \eta$ and $\partial_x^2 \eta$, appearing on the right-hand side of Eq. (80c), as well as $\partial_x v$ on the right-hand side, are computed in Fourier space. The first and second derivatives of the operator in the brackets in Eq. (80c) are replaced by the corresponding spectral differentiation matrices. This leads to a linear system for the values $\varphi_1(x_i)$, $i = 1, \dots, N$, which we solve by a LU decomposition method. Having established the calculation of $\varphi_1$ at the grid points, we evaluate the RHS of Eqs. (80a)–(80b), by computing the derivatives in Fourier space and the products in physical space.

## 8.2. Solitary wave propagation

Following Mitsotakis et al. [89], we examine how well the ($K = 1$)–model propagates solitary waves of the free-surface Euler equations calculated using the code provided in [90]. We consider three cases of increasing dimensionless amplitude $a/h_0$ and phase speed $C/\sqrt{gh_0}$, namely $(a/h_0, C/\sqrt{gh_0}) = (0.2, 1.0943)$, $(0.4, 1.1780)$ and $(0.6, 1.2500)$, and run the model up to dimensionless time $T\sqrt{gh_0} = 120$. Depending upon $a/h_0$, the solitary wave propagates over a distance ranging from 131.3 to 150 times the still water depth $h_0$.

In the extreme case $a/h_0 = 0.6$, it was found necessary to apply a filtering scheme in order to reduce the aliasing errors and obtain a stable simulation. For the non-linear terms, we discard at every time step the Fourier modes corresponding to $|n| > 3N/10$ (we have also tested the 2/3 rule obtaining similar results).

We plot our results in Fig. 3. In panels (a)–(c), we show the space-time evolution by using method B (the plots corresponding to method A are indistinguishable and are thus omitted). We see that the initial solitary waves stably propagate towards the positive of the $x$-axis in all three cases. For $a/h_0 = 0.6$, the amplitude of the solitary wave slightly increases initially and quickly reaches a steady value. A very small-amplitude wave trailing the main solitary wave is also visible. This is due to the fact that the initial conditions are not exact solutions of the travelling-wave version of Eqs. (54a)–(54c).

In panels (d)–(f), we compare our solutions at the end of the simulation with the exact solutions for the same travelling duration. For $a/h_0 = 0.2$ and $0.4$ the agreement between methods A and B and the exact solution is excellent. For the case $a/h_0 = 0.6$, the agreement remains very good.
<caption>**Table 1**</caption>

Relative errors on the amplitude and the travelling wave celerity at $T\sqrt{g/h_0} = 120$, calculated by Eqs. (81).

<table><thead><tr><th rowspan="2">Method</th><th colspan="2">0.2</th><th colspan="2">0.4</th><th colspan="2">0.6</th></tr><tr><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>ε[a]</td><td>1.2 × 10<sup>-3</sup></td><td>1.2 × 10<sup>-3</sup></td><td>6.6 × 10<sup>-3</sup></td><td>6.9 × 10<sup>-3</sup></td><td>3.1 × 10<sup>-2</sup></td><td>3.6 × 10<sup>-2</sup></td></tr><tr><td>ε[C]</td><td>3.7 × 10<sup>-5</sup></td><td>4.0 × 10<sup>-5</sup></td><td>2.7 × 10<sup>-4</sup></td><td>2.8 × 10<sup>-4</sup></td><td>1.2 × 10<sup>-3</sup></td><td>1.5 × 10<sup>-3</sup></td></tr></tbody></table>

Zooming in around the crest reveals that both numerical methods produce a solitary wave with slightly higher amplitude (see the right inset of panel (f)). In the left inset of panel (f), we focus on the region behind the solitary wave, where a small trailing wave is observed. Its maximum amplitude remains very small, not exceeding $2 \times 10^{-3}$, for method A and $2.5 \times 10^{-3}$ for method B indicating only a minor discrepancy. Note that when the SGN system is used in this test, the amplitude becomes significantly smaller than the exact one. Moreover, the numerically computed solitary wave deviates significantly from the initial Euler solitary wave [89].

In Table 1, we show the signed relative errors on the amplitude, $ε[a]$, and the travelling wave celerity, $ε[C]$, at the end of the simulation defined by

$$ε[a] = \frac{\max\{\eta(x, T)\} - a}{a}, \quad ε[C] = \frac{x_c(T)/T - C}{C}. \qquad (81)$$

To calculate $\max\{\eta(x, T)\}$, we interpolate the numerical solutions using cubic interpolation on a very fine grid and then take the maximum. The abscissa of the maximum is denoted $x_c(T)$. The differences between the two methods are extremely small, with method A performing slightly better. In both methods, the errors increase with $a/h_0$ but maintain a consistently small magnitude.

Finally, in panels (g)–(h) we show the relative errors in the conservation of excess of mass and energy,

$$E[M] = \left| \frac{M^{\alpha}(t) - M^{\alpha}(0)}{M^{\alpha}(0)} \right|, \quad E[H] = \left| \frac{H^{\alpha}(t) - H^{\alpha}(0)}{H^{\alpha}(0)} \right|, \qquad (82)$$

in a semi-log plot. $M^{\alpha}(t)$ and $H^{\alpha}(t)$ are calculated from the integrals in Eqs. (47) and (49), respectively. The error on the impulse follows the same trend as the energy, being slightly smaller and is not shown for the clarity of the figure. Considering $E[M]$, we see that method B results in an error below $10^{-14}$ in all three cases. Method A gives higher $E[M]$, which nevertheless can be considered as very low (e.g. at the end of the simulation for $a/h_0 = 0.6$, $E[M] = 2.4 \times 10^{-4}$). The error $E[H]$ reaches an almost steady state after $t\sqrt{g/h_0} = 20$. It is almost identical in both methods and clearly increases with $a/h_0$. At the end of the simulation for $a/h_0 = 0.6$, $E[H] \approx 6.3 \times 10^{-4}$ for method A and $E[H] \approx 10^{-3}$ for method B.

These numerical results showcase that the present model consistently approximates permanent form solutions of the free-surface Euler equations in the strongly non-linear regime.

## 8.3. Head-on collision of solitary waves

In this section, we evaluate the performance of the ($K=1$)-model in capturing the significant transformation of solitary waves during head-on collision. To compare with high-precision experimental data, we consider the configuration of two counter-propagating solitary waves as in [91]. To ensure a close match between the simulation and the experimental results before the collision, the relative amplitudes of the solitary waves are set to $a/h_0 = 0.4$ (wave heading to the right) and $0.39$ (wave heading to the left), as in [33]. This slight asymmetry reflects the experimental conditions observed in the first available snapshot of the free-surface elevation. We solve the ($K=1$)-model in the domain $x/h_0 \in [-80, 80]$ using both Methods A and B with $N = 1024$ and $\delta t\sqrt{g/h_0} = 0.01$. However, as the results are indistinguishable in the figures, we only plot results from Method B for clarity. We also consider the Fourier pseudo-spectral solution of the SGN model using the same numerical parameters and RK4 time-stepping, implemented with the code provided in [16]. Note that the initial conditions for the SGN model are constructed using its solitary wave solution, while for the ($K=1$) model, we use the Euler solitary wave solution, as in the previous section. The numerical results for Choi's $O(\mu^2)$ model are digitized from Choi [33]. These solutions are compared in Fig. 4 with experimental data, digitized from Chen and Yeh [91], at eight time instants before and after the collision. An animation illustrating the evolution of the free surface elevation is provided as Supplementary material.

All three models agree well with the experiment during the initial, highly non-linear stages of the collision process, up to the point of maximum elevation shown in panel (d). The maximum normalized elevation is 0.9 in the experiment, while the SGN, Choi's model, and the present ($K=1$)-model yield values of 0.88, 0.9, and 0.91, respectively. Visual differences become apparent in the two subsequent post-collision snapshots (e) and (f). Specifically, the SGN model overestimates the surface elevation at $x/h_0 = 0$ and fails to capture the characteristic dip observed in panel (f). In contrast, Choi's model and the ($K=1$)-model accurately replicate the experimental results at these times. This indicates the presence of higher-order effects that are not accounted for in the first-order SGN model. In the final two snapshots, (g) and (h), the model results show overall good agreement with each other and with the experiment. To quantify the model performance, we also include the Root Mean Square Error (RMSE) in Table 2, comparing the model results with the experimental data for the available snapshots in Fig. 4. The RMSE is calculated as

$$\text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (\eta_{\text{exp}} - \eta_{\text{model}})^2},$$

where $n$ is the number of data points. The results indicate that the ($K=1$)-model provides the best overall agreement with the experiment, with Choi's model yielding very close results. Finally, we note that the maximum instantaneous relative error in the conservation of energy in our solution with Method B is $3.3 \times 10^{-10}$, demonstrating a high level of numerical accuracy and precise conservation of the Hamiltonian during the strong interaction of solitary waves.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Head-on collision of two solitary waves of amplitude $a/h_0 = 0.39$ and 0.4. Snapshots of the free-surface elevation computed using the SGN model (dot-dashed black line), the Choi $O(\mu^2)$ model (dashed green line) and the $(K = 1)$-model (red line). Experimental data are represented by black circles. $t_s$ denotes the dimensionless time relative to the instant of maximum free surface-elevation ($t_s = 0$) shown in panel (d).</figcaption></figure>

<table><caption>**Table 2**
Root Mean Square Error (RMSE) between model results and experimental data for the snapshots (a)-(h) shown in figure 4.</caption><thead><tr><th></th><th>(a)</th><th>(b)</th><th>(c)</th><th>(d)</th><th>(e)</th><th>(f)</th><th>(g)</th><th>(h)</th></tr></thead><tbody><tr><td>SGN</td><td>0.024</td><td>0.029</td><td>0.044</td><td>0.056</td><td>0.037</td><td>0.033</td><td>0.021</td><td>0.010</td></tr><tr><td>Choi</td><td>0.015</td><td>0.026</td><td>0.027</td><td>0.033</td><td>0.026</td><td>0.024</td><td>0.036</td><td>0.019</td></tr><tr><td>(K = 1)-model</td><td>0.013</td><td>0.024</td><td>0.024</td><td>0.030</td><td>0.022</td><td>0.014</td><td>0.008</td><td>0.007</td></tr></tbody></table>

## 9. Concluding remarks

We have derived a new family of model equations for free-surface long gravity waves over flat bottom based on a vertical series representation for the velocity potential and Luke's variational principle. This representation contains the first $K + 1$ vertical polynomials of the shallow-water expansion of the solution to the Laplace boundary value problem in the ZCS formulation (DtN problem). The model equations derived this way possess a Hamiltonian structure analogous to the ZCS formulation. They are derived by a non-local Hamiltonian defined in terms of the solution of a system of $K$ elliptic spatial PDEs with variable coefficients on the horizontal plane. Additionally, we have shown that the new family is equivalent
to the IK family also produced by Luke's variational principle but with a different series representation containing the vertical polynomials of the Boussinesq expansion of the velocity potential. The usefulness of the present approach is that the equations are derived directly in a canonical Hamiltonian form that is convenient for numerical calculations.

An important feature of this class of model equations is that they are $O(\mu^{2K})$-accurate and, at the same time, the order of spatial derivatives in the evolution equations never exceeds two for all K. This is particularly relevant in comparison to asymptotic models which inevitably introduce higher-order derivatives as their accuracy improves, often leading to numerical challenges. We have demonstrated that our (K = 1)-model (accurate to $O(\mu^2)$) reduces to the classical $O(\mu)$-accurate SGN model that contains third-order spatial derivatives and mixed space-time derivatives. The time evolution of both models requires solving a single elliptic PDE in the case of one horizontal dimension and a flat bottom. In the case of two horizontal dimensions, the (K = 1)-model formulated in the canonical variables $(\eta, \psi)$ still requires solving only one elliptic PDE along with the time integration of two evolution PDEs. In contrast, the SGN model, expressed in velocity variables, requires solving a system of elliptic PDEs and advancing three evolution PDEs in time. Consequently, the (K = 1)-model provides a more accurate approximation of the water wave equations while maintaining comparable computational complexity. It is important to recall that, while the higher-order $O(\mu^2)$ SGN model has the same accuracy as the (K = 1)-model, it exhibits a singularity in its linear dispersion relation [31]. Additionally, the (K = 1)-model transforms into the $O(\mu^2)$-accurate model of Choi [33] which is formulated in terms of the bottom-surface potential and involves fourth-order spatial derivatives and mixed space-time derivatives. Note that, unlike the (K = 1)-model, Choi's model conserves the energy only asymptotically. The key conclusion from this analytical comparison is that the precision of the (K = 1)-model matches or surpasses that of asymptotic models with comparable mathematical complexity. Moreover, the (K = 1)-model has better linear dispersive properties and is free of the wave trough instability which is present in several other fully non-linear shallow-water models considered in [37].

We have also presented numerical solutions of the (K = 1)-model using either a FD method or a Fourier spectral method to discretize the spatial derivatives. Our results demonstrate that the model accurately propagates highly non-linear solitary waves of the free-surface Euler equations over long distances. Furthermore, we simulated the head-on collision of two solitary waves and compared the results with experimental measurements, demonstrating the model's ability to capture strong wave transformations, notably outperforming the lower-order SGN model. The (K = 1)-model and Choi's $O(\mu^2)$ model yield comparable results in this experiment. However, Choi's model is unstable with FD approximations due to its high-order spatial derivatives, requiring a Fourier pseudo-spectral method combined with a filtering technique for stable solutions [33]. In contrast, the FD approximation of the (K = 1)-model yields stable long-distance simulations without the need for numerical filters. The same is true for the Fourier-spectral method, provided the non-linearity remains moderate; a filter is required for cases with higher non-linearity though. These results underscore the versatility of the present model, highlighting its potential as a valuable alternative for studying strongly non-linear and dispersive water wave dynamics.

To achieve further improvement in dispersive properties, one may consider $K > 1$, which entails solving a system of K linear PDEs throughout the evolution, see Eqs. (34)–(36). We emphasize that, irrespective of the choice of K, the highest order of derivatives in the system remains two. For large-scale and long-term simulations, one of the challenges lies in efficiently solving the linear system along the evolution, which can become computationally demanding despite the advantage of low-order derivatives. In such cases, iterative methods will be typically required. In the context of Fourier-spectral methods, recent examples of such techniques include [16, 22, 23, 33].

The results presented in this work lay the foundation for ongoing and future analysis. Having established that our new model accurately simulates strongly non-linear solitary waves and their collision, one might now explore its behaviour in less idealized cases where more complicated wave dynamics are expected (e.g multi-wave interactions). A key direction for future work is to examine the limits of this modelling approach as wave conditions extend into deeper water regimes. A limitation of the present approach is its restriction to a flat bottom, and our goal is to extend it to variable bottom cases to enhance its practical utility. As already mentioned, one way to do this is by using the same $\Phi^a$, Eq. (24), but allow the still water depth $h$ to vary in $(x_1, x_2)$. In this case, we still have $[\partial_x\Phi^a]_{-h} = 0$, which would remain a good approximation only in a mild slope setting [71]. We remark, however, that [72] employed such an approach and, despite this inconsistency, achieved good agreement for the free surface elevation in variable depth experiments. A more consistent approach would be to incorporate appropriate vertical functions into Eq. (24). Following the asymptotic solution of the DtN problem in this context, these vertical functions can be represented as odd polynomials, such as $Z_n^b = (z - n)^{2n+1}$, $n = 0, \dots, K - 1$. We have already verified that the EL equations obtained from Eqs. (14)–(15) in this case have a Hamiltonian structure with a positive definite Hamiltonian and a non-singular dispersion relation. The additional complexities in this case arise from the fact that the two evolution Hamiltonian PDEs become coupled with a system of two linear second-order PDEs. Another limitation concerns the vanishing-at-infinity assumptions. This can be addressed by deriving lateral boundary conditions for the functions $\varphi_n$ from Luke's variational principle. Implementing such boundary conditions will require modifications to the FD scheme, whereas the Fourier-spectral method, being inherently periodic, is not directly applicable. The details and a numerical study of these extensions will be reported in a forthcoming article. Another natural extension of the present approach is the consideration of internal (interfacial) waves in multi-layer fluids of different (constant) densities. Although significant work remains to be done, the results of this study establish the basis for a practical model that can be valuable in coastal engineering applications, particularly for long-duration simulations relevant to wave energy and coastal hazards investigations.

### CRediT authorship contribution statement

**Christos E. Papoutsellis:** Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization. **Michel Benoit:** Writing – review & editing, Validation, Software, Methodology, Investigation, Formal analysis.

### Declaration of competing interest

We have no competing interests to declare.

### Acknowledgements

The authors thank the anonymous reviewers of an earlier version of this article for their insightful comments and suggestions. The authors also acknowledge Prof. Harry Yeh (Oregon State University) for discussions on experimental data of solitary wave collision, and Guillaume Coulaud (EDF R&D LNHE and LHSV) for fruitful discussions. Ch. P also acknowledges Prof. Gerassimos Athanassoulis (National Technical University of Athens) for valuable discussions during the early stage of this work.
## Appendix A. Dimensionless water wave problem

In the case of a flat bottom, two important parameters are the typical horizontal wavelength L and the constant depth h₀. The non-dimensional variables (denoted with primes) are defined as

$$x_1 = x'_1 L, \quad x_2 = x'_2 L, \quad z = z'_0, \quad t = t'_0 L/\sqrt{gh_0}, \quad \eta = \eta'_0 h_0, \quad \Phi = (L\sqrt{gh_0})\Phi'. \quad (A.1)$$

Introducing the shallowness parameter $\mu = (h_0/L)^2$, the dimensionless version of (1) reads (after dropping the primes)

$$\mu \Delta \Phi + \partial_z^2 \Phi = 0, \quad -1 \le z \le \eta, \quad (A.2)$$

$$\partial_z \Phi = 0, \quad z = -1, \quad (A.3)$$

$$\partial_t \eta - \nabla \eta \cdot \nabla \Phi + \frac{1}{\mu} \partial_z \Phi = 0, \quad z = \eta, \quad (A.4)$$

$$\partial_t \Phi + \frac{1}{2}(\nabla \Phi)^2 + \frac{1}{2\mu}(\partial_z \Phi)^2 + \eta = 0, \quad z = \eta. \quad (A.5)$$

## Appendix B. Shallow-water solution of the DtN problem

For $\mu \ll 1$, we look for an approximate solution of Eq. (20) (Eqs. (A.2)–(A.3) and $[\Phi]_\eta = \psi$) in the form of a series expansion in powers of $\mu$, $\Phi^{\text{app}}(X, z) = \sum_{j=0}^{K} \mu^j \Phi^{(j)}(X, z)$. Plugging $\Phi^{\text{app}}$ into the first equation of (20) we get

$$\sum_{j=0}^{K} \mu^j (\Delta\Phi^{(j-1)} + \partial_z^2 \Phi^{(j)}) = O(\mu^{K+1}), \quad (B.1)$$

with $\Phi^{(-1)} \equiv 0$ by convention. By cancelling the $O(\mu^{K+1})$ residual, we first find for the term of order $j = 0$ the BVP,

$$\partial_z^2 \Phi^{(0)} = 0, \quad -1 \le z \le \eta, \quad \Phi^{(0)}(X, \eta) = \psi, \quad \partial_z \Phi^{(0)}(X, -1) = 0, \quad (B.2)$$

whose solution is $\Phi^{(0)}(X, z) = \psi(X)$, independent of $z$. For orders $1 \le j \le K$, we obtain the recurrence relation:

$$\partial_z^2 \Phi^{(j)} = -\Delta \Phi^{(j-1)}, \quad -1 \le z \le \eta, \quad \Phi^{(j)}(X, \eta) = 0, \quad \partial_z \Phi^{(j)}(X, -1) = 0. \quad (B.3)$$

Integrating the above differential equation twice and using the boundary conditions we find

$$\Phi^{(j)}(X, z) = \int_z^\eta \int_{-1}^{z'} \Delta\Phi^{(j-1)}(X, z'') dz'' dz', \quad 1 \le j \le K. \quad (B.4)$$

For $j = 1$, as $\Phi^{(0)}(X, z) = \psi(X)$, Eq. (B.4) yields

$$\Phi^{(1)}(X, z) = \Delta\psi(X)\tilde{Z}(z; \eta) \quad (B.5)$$

where the vertical dependence $\tilde{Z}(z; \eta) \equiv \tilde{Z}$ is obtained explicitly as a quadratic polynomial in $z$,

$$\tilde{Z} \equiv \int_z^\eta \int_{-1}^{z'} dz'' dz' = \frac{(1+\eta)^2 - (1+z)^2}{2} = -\frac{1}{2}(z-\eta)^2 - (\eta+1)(z-\eta), \quad (B.6)$$

We see from Eq. (B.6) that $\tilde{Z}(\eta; \eta) = 0$ and $\tilde{Z}(z; \eta) > 0$ for $-1 \le z < \eta$, with $\tilde{Z}(-1; \eta) = (1+\eta)^2/2 = H^2/2$. To compute $\Phi^{(2)}$, we first differentiate Eq. (B.5) to get

$$\Delta\Phi^{(1)} = \Delta^2\psi\tilde{Z} + 2\nabla(\Delta\psi) \cdot \nabla\tilde{Z} + \Delta\psi\Delta\tilde{Z}, \quad (B.7)$$

where $\Delta^2 \equiv \Delta\Delta$ is the horizontal bi-Laplacian operator and

$$\nabla\tilde{Z} = \partial_\eta\tilde{Z}\nabla\eta, \quad \Delta\tilde{Z} = \partial_{\eta\eta}\tilde{Z}(\nabla\eta)^2 + \partial_\eta\tilde{Z}\Delta\eta, \quad \partial_\eta\tilde{Z} = 1 + \eta = H. \quad (B.8)$$

Next, by direct integration we have

$$\int_z^\eta \int_{-1}^{z'} \tilde{Z}(z'') dz'' dz' = \frac{1}{6} \tilde{Z}^2 + \frac{H^2}{3} \tilde{Z}, \quad (B.9)$$

$$\int_z^\eta \int_{-1}^{z'} \nabla \tilde{Z}(z'') dz'' dz' = H \nabla \eta \tilde{Z} = \frac{1}{2} \nabla(H^2) \tilde{Z}, \quad (B.10)$$

$$\int_z^\eta \int_{-1}^{z'} \Delta \tilde{Z}(z'') dz'' dz' = [(\nabla \eta)^2 + H \Delta \eta] \tilde{Z} = \frac{1}{2} \Delta(H^2) \tilde{Z}. \quad (B.11)$$

Then, Eq. (B.4) for $j = 2$ gives

$$\Phi^{(2)}(X, z) = \left(\frac{1}{6}\tilde{Z}^2 + \frac{H^2}{3}\tilde{Z}\right) \Delta^2\psi + \tilde{Z} \left(\nabla(H^2) \cdot \nabla(\Delta\psi) + \frac{1}{2}\Delta(H^2)\Delta\psi\right), \quad (B.12)$$

and the expression in Eq. (B.21) follows from a simple rearrangement.
## Appendix C. Proof of Lemma 1

Using (32b), (24) and the notation (16) we obtain

$$ 
\begin{aligned}
\Delta\psi(\partial_{\eta}Z_m, 1)\varphi_m + \Lambda_{mn}\varphi_n\varphi_m &= \sum_{m=1}^{K} \varphi_m \int_{-h}^{\eta} \Delta_{X,z}\Phi^a\partial_{\eta}Z_m dz \\
&= \sum_{m=1}^{K} mH\varphi_m \int_{-h}^{\eta} \Delta_{X,z}\Phi^a Z_{m-1} dz \\
&= H\varphi_1 \int_{-h}^{\eta} \Delta_{X,z}\Phi^a dz \\
&= (H\Delta\psi + L_{0n}\varphi_n) H\varphi_1
\end{aligned}
$$

where we have used that

$$ \int_{-h}^{\eta} \Delta_{X,z} \Phi^a Z_{m-1} dz = 0 \quad \text{for } m = 2, \dots, K, $$

which holds in view of (30) and identity (32a).

## Appendix D. Vertical integrals appearing in Eq. (18)

Let $Z_m \equiv Z^m$ for $m \ge 0$ as in (24). Consider the integral $(\nabla Z_n, Z_m)$. By taking the gradient of $A_{mn}$ and using Eq. (25) and the relation $\nabla Z_n = \nabla\eta\partial_\eta Z_n = n\nabla\eta H Z_{n-1}$, we obtain

$$ (\nabla Z_n, Z_m) = \nabla A_{mn} - mH\nabla\eta A_{m-1,n}, \quad m \ge 0, n \ge 1. \qquad (D.1) $$

Similarly, taking the divergence of $(\nabla Z_n, Z_m)$ we find for $m \ge 0, n \ge 1$

$$ (\Delta Z_n, Z_m) = \nabla \cdot (\nabla Z_n, Z_m) - n\delta_{1n}\delta_{0m}(\nabla\eta)^2 H - mnH^2(\nabla\eta)^2 A_{m-1,n-1}. \qquad (D.2) $$

where $\delta$ is the Kronecker delta. For $(\partial_z^2 Z_n, Z_m)$, we use integration by parts together with $\partial_z Z_n = -n(z+h)Z^{n-1} = n(z+h)Z_{n-1}$ and

$$ -2Z + H^2 = (z + h)^2, \qquad (D.3) $$

to obtain

$$ (\partial_z^2 Z_n, Z_m) = -\delta_{1n}\delta_{0m} nH + mnA_{m,n-1} - mnH^2 A_{m-1,n-1}. \qquad (D.4) $$

## Appendix E. Calculation of $\delta_{\psi H}$ and $\delta_{\eta H}$

We start with $\delta_{\psi H}$. From (44) and (45), we easily obtain

$$ D_{\psi} H^{\alpha}(\eta, \psi)(\delta\psi) = \int G^{\alpha}[\eta]\psi \, \delta\psi \, dX, \qquad (E.1) $$

hence, $\delta_{\psi} H^{\alpha} = G^{\alpha}[\eta]\psi$. To compute $\delta_{\eta} H^{\alpha}$ we need the derivatives of $\Phi^{\alpha}[\eta]$ and $f_{n}[\eta]$ at $\eta$ denoted by $D_{\eta}\Phi^{\alpha}[\eta](\delta\eta)\psi$ and $D_{\eta}f_{n}[\eta](\delta\eta)\psi$ respectively. From (41), we have

$$ D_{\eta}\Phi^{\alpha}[\eta](\delta\eta)\psi = (D_{\eta}f_{n}[\eta](\delta\eta)\psi) Z_{n} + (f_{n}[\eta]\psi) \partial_{\eta}Z_{n}\delta\eta. \qquad (E.2) $$

Differentiating (42) with respect to $\eta$ using (25) we see that $D_{\eta}\Phi^{\alpha}[\eta](\delta\eta)\psi$ satisfies

$$ \int_{-h}^{\eta} \Delta_{X,z} (D_{\eta}\Phi^{\alpha}[\eta](\delta\eta)\psi) Z_m dz + \int_{-h}^{\eta} \Delta_{X,z} (\Phi^{\alpha}[\eta]\psi) \partial_{\eta} Z_m \delta\eta dz = 0, \quad m = 1, \dots, K. \qquad (E.3) $$

Using (E.2), (25) and ((17a,b)a) (with $\partial_l h = 0$ and (28)), Eq. (E.3) is equivalently written as $L\mathbf{u} = g\mathbf{a}$ with

$$ \mathbf{u} = (D_{\eta}f_{1}[\eta](\delta\eta)\psi, \dots, D_{\eta}f_{K}[\eta](\delta\eta)\psi)^{t}, \quad g\mathbf{a} = (\partial_{\eta}L_{1n}(f_{n}[\eta]\psi), \dots, \partial_{\eta}L_{Kn}(f_{n}[\eta]\psi))^{t} \delta\eta, \qquad (E.4) $$

so that $\mathbf{u}$ and consequently $D_{\eta}\Phi^{\alpha}[\eta](\delta\eta)\psi$ are well defined. Posing $\Phi^{\alpha} \equiv \Phi^{\alpha}[\eta]\psi$ and $D_{\eta}\Phi^{\alpha}(\delta\eta) \equiv D_{\eta}\Phi^{\alpha}[\eta](\delta\eta)\psi$ for brevity and differentiating $H^{\alpha}$ we get

$$ D_{\eta}H^{\alpha}(\delta\eta) = \frac{1}{2} \int [\nabla\Phi^{\alpha}]_{\eta}^{2} \delta\eta \, dX + \int \int_{-h}^{\eta} \nabla\Phi^{\alpha} \cdot \nabla (D_{\eta}\Phi^{\alpha}(\delta\eta)) \, dz \, dX + g\eta \, \delta\eta. \qquad (E.5) $$

Using (26), we see that

$$ [\nabla\Phi^{\alpha}]_{\eta}^{2} = (\nabla\psi)^{2} - 2\nabla\eta \cdot \nabla\psi [\partial_{z}\Phi^{\alpha}]_{\eta} + ((\nabla\eta)^{2} + 1) [\partial_{z}\Phi^{\alpha}]_{\eta}^{2}. \qquad (E.6) $$

Using Green's identity we obtain

$$ \int \int_{-h}^{\eta} \nabla\Phi^{\alpha} \cdot \nabla (D_{\eta}\Phi^{\alpha}(\delta\eta)) dz dX = \int \left\{ \mathbf{N}_{\eta} \cdot [\nabla\Phi^{\alpha} D_{\eta}\Phi^{\alpha}(\delta\eta)]_{\eta} - \int_{-h}^{\eta} \Delta_{X,z}\Phi^{\alpha} D_{\eta}\Phi^{\alpha}(\delta\eta) dz \right\} dX. \qquad (E.7) $$

For the first term in the RHS of (E.7) we use (E.2) together with (25) and (26):

$$ \mathbf{N}_{\eta} \cdot [\nabla\Phi^{\alpha} D_{\eta}\Phi^{\alpha}(\delta\eta)]_{\eta} = (\nabla\eta \cdot \nabla\psi [\partial_{z}\Phi^{\alpha}]_{\eta} - ((\nabla\eta)^{2} + 1) [\partial_{z}\Phi^{\alpha}]_{\eta}^{2}) \delta\eta. \qquad (E.8) $$
For the second term in the RHS of (E.7) we use (42) and (E.2) to get

$$ \int_{-h}^{n} \Delta_{X,z} \Phi^a D_{\eta} \Phi^a (\delta\eta) dz = \left( \int_{-h}^{\eta} \Delta_{X,z} \Phi^a \partial_{\eta} Z_n dz \right) f_n \delta\eta. \quad (E.9) $$

Substituting (E.6) and (E.7) (with (E.8)–(E.9)) into (E.5) we find $\delta_{\eta}H^a$.

## Appendix F. Equivalence for arbitrary K

We establish here that the EL Eqs. (14)–(15) obtained by using the representations (24) and (50) are equivalent. Using Eq. (D.3), we can express each $Z_n$ in Eq. (24) in terms of $\{Z_n^*\}_{n=0}^K = \{(z+h)^{2n}\}_{n=0}^K$ in Eq. (50) as

$$ Z_n = \frac{1}{2^n} (H^2 - (z+h)^2)^n = \frac{1}{2^n} \sum_{k=0}^{n} \binom{n}{k} H^{2(n-k)} (-1)^k (z+h)^{2k}. \quad (F.1) $$

Taking into account Eq. (F.1), the ansatz (24) can be equivalently reformulated as:

$$ \begin{aligned} \Phi^a &= \varphi_0 + \frac{\varphi_1}{2} (H^2 - (z+h)^2) + \frac{\varphi_2}{4} (H^2 - (z+h)^2)^2 + \dots + \frac{\varphi_K}{2^K} (H^2 - (z+h)^2)^K \\ &= \varphi_0 + \frac{\varphi_1}{2} (H^2 - (z+h)^2) + \frac{\varphi_2}{4} (H^4 - 2(z+h)^2 H^2 + (z+h)^4) \\ &\quad + \dots + \frac{\varphi_K}{2^K} \sum_{k=0}^{K} \binom{K}{k} (-1)^k H^{2(K-k)} (z+h)^{2k} \\ &= \varphi_0 + \frac{\varphi_1}{2} H^2 + \frac{\varphi_2}{4} H^4 + \frac{\varphi_3}{4} H^6 + \dots + \frac{\varphi_K}{2^K} \binom{K}{0} H^{2K} \\ &\quad - \left[ \frac{\varphi_1}{2} + \frac{\varphi_2}{4}(2H^2) + \frac{\varphi_3}{8}(3H^4) + \dots + \frac{\varphi_K}{2^K} \binom{K}{1} H^{2(K-1)} \right] (z+h)^2 \\ &\quad + \left[ \frac{\varphi_2}{4} + \frac{\varphi_3}{8}(3H^2) + \dots + \frac{\varphi_K}{2^K} \binom{K}{2} H^{2(K-2)} \right] (z+h)^4 \\ &\quad + \dots + (-1)^K \left[ \frac{\varphi_K}{2^K} \binom{K}{K} H^{2(K-K)} \right] (z+h)^{2K} \\ &= \sum_{n=0}^{K} \left( (-1)^n \sum_{k=n}^{K} \frac{1}{2^k} H^{2(k-n)} \binom{k}{n} \varphi_k \right) (z+h)^{2n}. \end{aligned} \quad (F.2) $$

where we have used the notation $\varphi_0 = \psi$. By comparing (50) and (F.2) we see that

$$ \varphi_n = (-1)^n \sum_{k=n}^{K} \frac{1}{2^k} H^{2(k-n)} \binom{k}{n} \varphi_k, \quad n = 0, \dots, K. \quad (F.3) $$

The relations (F.3) can be inverted through binomial inversion [92]:

$$ \varphi_n = (-2)^n \sum_{k=n}^{K} H^{2(k-n)} \binom{k}{n} \varphi_k, \quad n = 0, \dots, K. \quad (F.4) $$

Let us denote $\varphi = (\varphi_0, \dots, \varphi_K)^T$ and $\varphi = (\varphi_0, \dots, \varphi_K)^T$, and introduce the functionals $S^a(\eta, \varphi) = S(\eta, \Phi_*^a)$ and $S^a(\eta, \varphi) = S(\eta, \Phi^a)$, where we recall that $S$ is Luke's functional (8). $S_*^a(\eta, \varphi)$ corresponds to the IK model and our functional $S^a(\eta, \varphi)$ is related to $S_*^a(\eta, \varphi)$ by

$$ S_*^a(\eta, \varphi) = S^a(\eta, \varphi M(H)), \quad (F.5) $$

where $M(H)$ is a $(K+1) \times (K+1)$ invertible upper triangular matrix of polynomials in $H$ whose entries are the coefficients of $\varphi_k$ in the $K+1$ relations (F.4). From the above relation, we can deduce that $(\eta, \varphi)$ is a critical point of $S_*^a$ if and only if $(\eta, \varphi)$ is a critical point of $S^a$. In other words, one can derive EL equations using either $S_*^a$ or $S^a$ and switch between formulations using (F.3) or (F.4).

## Appendix G. Supplementary data

Supplementary material related to this article can be found online at https://www.sciencedirect.com/science/article/pii/S016727892500168X via%3Dihub#appSH.

## Data availability

No data was used for the research described in the article.

## References

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

[2] J.T. Kirby, Boussinesq models and their application to coastal processes across a wide range of scales, J. Waterw. Port Coast. Ocean. Eng. 142 (6) (2016) 03116005, http://dx.doi.org/10.1061/(ASCE)WW.1943-5460.0000350.

[3] D. Lannes, Modeling shallow water waves, Nonlinearity 33 (5) (2020) R1–R57, http://dx.doi.org/10.1088/1361-6544/ab6c7c.

[4] F. Serre, Contribution à l'étude des écoulements permanents et variables dans les canaux, Houille Blanche 39 (6) (1953) 830–872, http://dx.doi.org/10.1051/lhb/1953058.
[5] C.H. Su, C.S. Gardner, Korteweg de Vries equation and generalizations. III. Derivation of the Korteweg de Vries equation and Burgers equation, J. Math. Phys. 10 (3) (1969) 536-539, http://dx.doi.org/10.1063/1.1664873.

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

[7] J. Miles, R. Salmon, Weakly dispersive nonlinear gravity waves, J. Fluid Mech. 157 (1985) 519-531, http://dx.doi.org/10.1017/S0022112085002488.

[8] S.V. Bazdenkov, N.N. Morozov, O.P. Pogutse, Dispersive effects in two-dimensional hydrodynamics, Dokl. Akad. Nauk SSSR 293 (1987) 818-822, www.mathnet.ru/eng/dan47874.

[9] D. Lannes, P. Bonneton, Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation, Phys. Fluids 21 (2009) 016601, http://dx.doi.org/10.1063/1.3053183.

[10] D. Clamond, D. Dutykh, Practical use of variational principles for modeling water waves, Phys. D: Nonlinear Phenom. 241 (1) (2012) 25-36, http://dx.doi.org/10.1016/j.physd.2011.09.015.

[11] D. Dutykh, D. Clamond, M. Chhay, Serre-type equations in deep water, Math. Model. Nat. Phenom. 12 (1) (2017) 23-40, http://dx.doi.org/10.1051/mmnpp/201712103.

[12] V. Duchêne, S. Israwi, Well-posedness of the Green-Naghdi and Boussinesq-Peregrine systems, Ann. Math. Blaise Pascal 25 (1) (2018) 21-74, http://dx.doi.org/10.5802/ambp.372.

[13] G. Khakimzyanov, D. Dutykh, Z. Fedotova, O. Gusev, Dispersive shallow water waves: Theory, modeling, and numerical methods, in: Lecture Notes in Geosystems Mathematics and Computing, Springer, 2020, http://dx.doi.org/10.1007/978-3-030-46267-3, URL: https://link.springer.com/book/10.1007/978-3-030-46267-3.

[14] D. Lannes, The Water Waves Problem Mathematical Analysis and Asymptotics, in: Mathematical Surveys and Monographs, vol. 188, American Mathematical Society, 2013.

[15] O. Le Métayer, S. Gavrilyuk, S. Hank, A numerical scheme for the Green-Naghdi model, J. Comput. Phys. 229 (6) (2010) 2034-2045, http://dx.doi.org/10.1016/j.jcp.2009.11.021.

[16] D. Dutykh, D. Clamond, P. Milewski, D. Mitsotakis, Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations, Eur. J. Appl. Math. 24 (5) (2013) 761-787, http://dx.doi.org/10.1017/S0956792513000168.

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

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

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

[20] D. Mitsotakis, C. Synolakis, M. McGuinness, A modified Galerkin/finite element method for the numerical solution of the Serre-Green-Naghdi system, Int. J. Numer. Methods Fluids 83 (10) (2017) 755-778, http://dx.doi.org/10.1002/fld.4293.

[21] N. Favrie, S. Gavrilyuk, A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves, Nonlinearity 30 (7) (2017) 2718, http://dx.doi.org/10.1088/1361-6544/aa712d.

[22] V. Duchêne, C. Klein, Numerical study of the serre-green-naghdi equations and a fully dispersive counterpart, Discret. Contin. Dyn. Syst. B 27 (10) (2022) 5905-5933, http://dx.doi.org/10.3934/dcdsb.2021300, URL: https://www.aimsciences.org/article/article/id/61ca74fe2d80b7246ec61e55.

[23] S. Gavrilyuk, C. Klein, Numerical study of the Serre-Green-Naghdi equations in 2D\*, Nonlinearity 37 (4) (2024) 045014.

[24] T.Y. Wu, Nonlinear waves and solitons in water, Phys. D: Nonlinear Phenom. 123 (1) (1998) 48-63, http://dx.doi.org/10.1016/S0167-2789(98)00111-0.

[25] P.A. Madsen, H.A. Schäffer, Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis, Philos. Trans. R. Soc. A 356 (1749) (1998) 3123-3181, http://dx.doi.org/10.1098/rsta.1998.0309.

[26] M.F. Gobbi, J.T. Kirby, Wave evolution over submerged sills: tests of a high-order Boussinesq model, Coast. Eng. 37 (1) (1999) 57-96, http://dx.doi.org/10.1016/S0378-3839(99)00015-0.

[27] Y. Agnon, P.A. Madsen, H.A. Schäffer, A new approach to high-order Boussinesq models, J. Fluid Mech. 399 (1999) 319-333, http://dx.doi.org/10.1017/S0022112099006394.

[28] P.A. Madsen, H.B. Bingham, H. Liu, A new Boussinesq method for fully nonlinear waves from shallow to deep water, J. Fluid Mech. 462 (2002) 1-30, http://dx.doi.org/10.1017/S0022112002008467.

[29] P.A. Madsen, Y. Agnon, Accuracy and convergence of velocity formulations for water waves in the framework of Boussinesq theory, J. Fluid Mech. 477 (2003) 285-319, http://dx.doi.org/10.1017/S0022112002003257.

[30] Y. Matsuno, Hamiltonian formulation of the extended Green-Naghdi equations, Phys. D: Nonlinear Phenom. 301-302 (2015) 1-7, http://dx.doi.org/10.1016/j.physd.2015.03.001.

[31] Y. Matsuno, Hamiltonian structure for two-dimensional extended Green-Naghdi equations, Proc. R. Soc. A: Math. Phys. Eng. Sci. 472 (2016) 2190, http://dx.doi.org/10.1098/rspa.2016.0127.

[32] T.P. Horikis, D.J. Frantzeskakis, N.F. Smyth, Extended shallow water wave equations, Wave Motion 112 (2022) 102934, http://dx.doi.org/10.1016/j.wavemoti.2022.102934.

[33] W. Choi, High-order strongly nonlinear long wave approximation and solitary wave solution, J. Fluid Mech. 945 (2022) A15, http://dx.doi.org/10.1017/jfm.2022.544.

[34] D.M. Milder, The effects of truncation on surface-wave Hamiltonians, J. Fluid Mech. 217 (1990) 249-262, http://dx.doi.org/10.1017/S0022112090000714.

[35] S.B. Yoon, P.L.F. Liu, A note on Hamiltonian for long water waves in varying depth, Wave Motion 20 (4) (1994) 359-370, http://dx.doi.org/10.1016/0165-2125(94)90019-1.

[36] A.C. Radder, Hamiltonian dynamics of water waves, in: Advances in Coastal and Ocean Engineering: Volume 4, World Scientific, 1999, pp. 21-59.

[37] P.A. Madsen, D.R. Fuhrman, Trough instabilities in Boussinesq formulations for water waves, J. Fluid Mech. 889 (2020) A38, http://dx.doi.org/10.1017/jfm.2020.76.

[38] J.T. Kirby, A new instability for Boussinesq-type equations, J. Fluid Mech. 894 (2020) F1, http://dx.doi.org/10.1017/jfm.2020.257.

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

[40] W. Choi, R. Camassa, Fully nonlinear internal waves in a two-fluid system, J. Fluid Mech. 396 (1999) 1-36, http://dx.doi.org/10.1017/S0022112099005820.

[41] W. Craig, P. Guyenne, H. Kalisch, A new model for large amplitude long internal waves, C. R. Méc. 332 (7) (2004) 525-530.

[42] P. Guyenne, Large-amplitude internal solitary waves in a two-fluid model, C. R. Méc. 334 (6) (2006) 341-346.

[43] W. Choi, High-order strongly nonlinear long wave approximation and solitary wave solution. Part 2. Internal waves, J. Fluid Mech. 952 (2022) A41, http://dx.doi.org/10.1017/jfm.2022.950.

[44] J. Luke, A variational principle for a fluid with a free surface, J. Fluid Mech. 27 (1967) 395-397, http://dx.doi.org/10.1017/S0022112067000412.

[45] K.A. Belibassakis, G.A. Athanassoulis, A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions, Coast. Eng. 58 (2011) 337-350, http://dx.doi.org/10.1016/j.coastaleng.2010.11.007.

[46] C.E. Papoutsellis, A.G. Charalamopoulos, G.A. Athanassoulis, Implementation of a fully nonlinear Hamiltonian Coupled-Mode Theory, and application to solitary wave problems over bathymetry, Eur. J. Mech. B Fluids 72 (2018) 199-224, http://dx.doi.org/10.1016/j.euromechflu.2018.04.015.

[47] G.A. Athanassoulis, C.E. Papoutsellis, Exact semi-separation of variables in waveguides with non-planar boundaries, Proc. R. Soc. A: Math. Phys. Eng. Sci. 473 (2201) (2017) 20170017, http://dx.doi.org/10.1098/rspa.2017.0017.

[48] M.L. Yates, M. Benoit, Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves, Internat. J. Numer. Methods Fluids 77 (2015) 616-640, http://dx.doi.org/10.1002/fld.3992.

[49] J. Zhang, M. Benoit, Wave-bottom interaction and extreme wave statistics due to shoaling and de-shoaling of irregular long-crested wave trains over steep seabed changes, J. Fluid Mech. 912 (2021) A28, http://dx.doi.org/10.1017/jfm.2020.1125.

[50] V.E. Zakharov, Stability of periodic waves of finite amplitude on the surface of a deep fluid, J. Appl. Mech. Tech. Phys. 9 (1968) 86-94, http://dx.doi.org/10.1007/BF00913182.

[51] W. Craig, C. Sulem, Numerical simulation of gravity waves, J. Comput. Phys. 108 (1993) 73-83, http://dx.doi.org/10.1006/jcph.1993.1164.

[52] C.E. Papoutsellis, G.A. Athanassoulis, A new efficient Hamiltonian approach to the nonlinear water-wave problem over arbitrary bathymetry, 2017, https://arxiv.org/abs/1704.03276.

[53] C. Raoult, M. Benoit, M.L. Yates, Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments, Coast. Eng. 114 (2016) 194-207, http://dx.doi.org/10.1016/j.coastaleng.2016.04.003.
[54] T.K. Papathanasiou, C.E. Papoutsellis, G.A. Athanassoulis, Semi-explicit solutions to the water-wave dispersion relation and their role in the non-linear Hamiltonian coupled-mode theory, J. Engrg. Math. 114 (1) (2018) 87–114, http://dx.doi.org/10.1007/s10665-018-09983-1.

[55] G.A. Athanassoulis, C.P. Mavroeidis, P.E. Koutsogiannakis, C.E. Papoutsellis, A numerical study of the run-up and the force exerted on a vertical wall by a solitary wave propagating over two tandem trenches, J. Ocean. Eng. Mar. Energy 5 (2019) 311–331, http://dx.doi.org/10.1007/s40722-019-00148-5.

[56] C. Raoult, M. Benoit, M.L. Yates, Development and validation of a 3D RBF-spectral model for coastal wave simulation, J. Comput. Phys. 378 (2019) 278–302, http://dx.doi.org/10.1016/j.jcp.2018.11.002.

[57] C.E. Papoutsellis, B. Simon, M.L. Yates, M. Benoit, Modelling of depth-induced wave breaking in a fully nonlinear free-surface potential flow model, Coast. Eng. 157 (2019) 103579, http://dx.doi.org/10.1016/j.coastaleng.2019.103579.

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

[59] D. Clamond, D. Dutykh, D. Mitsotakis, Variational Approach to Water Wave Modelling, International Association for Hydro-Environment Engineering and Research (IAHR), 2024, http://dx.doi.org/10.3850/978-90-833476-6-0_IAHR_Watermonograph_003.

[60] M. Isobe, Time-dependent mild-slope equations for random waves, in: Proc. 24th Int. Conf. Coastal Eng., ICCE'1994, 1994, pp. 285–299, http://dx.doi.org/10.1061/9780784400890.023.

[61] T. Kakinuma, A set of fully nonlinear equations for surface and internal gravity waves, WIT Trans. Built Environ. 58 (2001) http://dx.doi.org/10.2495/CE010201.

[62] J. 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, J. Math. Pures Appl. 17 (1872) 55–108, http://eudml.org/doc/234248.

[63] L. Rayleigh, XXXII. On waves, Lond. Edinb. Dublin Philos. Mag. J. Sci. 1 (4) (1876) 257–279, http://dx.doi.org/10.1080/14786447608639037.

[64] G.B. Whitham, Linear and Non Linear Waves, Wiley-Interscience, New-York, 1974.

[65] Y. Murakami, T. Iguchi, Solvability of the initial value problem to a model system for water waves, Kodai Math. J. 38 (2015) 470–491, http://dx.doi.org/10.2996/kmj/1436403901.

[66] T. Iguchi, A mathematical justification of the Isobe–Kakinuma model for water waves with and without bottom topography, J. Math. Fluid Mech. 20 (2018) 1985–2018, http://dx.doi.org/10.1007/s00021-018-0398-x.

[67] R. Nemoto, T. Iguchi, Solvability of the Initial Value Problem to the Isobe–Kakinuma model for water waves, J. Math. Fluid Mech. 20 (2018) 632–653, http://dx.doi.org/10.1007/s00021-017-0338-1.

[68] M. Colin, T. Iguchi, Solitary wave solutions to the Isobe-Kakinuma model for water waves, Stud. Appl. Math. 145 (1) (2020) 52–80, http://dx.doi.org/10.1111/sapm.12310.

[69] V. Duchêne, T. Iguchi, A Hamiltonian structure of the Isobe–Kakinuma model for water waves, Water Waves 3 (1) (2021) 193–211, http://dx.doi.org/10.1007/s42286-020-00025-x.

[70] S. Abohadima, M. Isobe, Linear and nonlinear wave diffraction using the nonlinear time dependent mild slope equations, Coast. Eng. 37 (2) (1999) 175–192, http://dx.doi.org/10.1016/S0378-3839(99)00020-4.

[71] G. Klopman, B. van Groesen, M.W. Dingemans, A variational approach to Boussinesq modelling of fully nonlinear water waves, J. Fluid Mech. 657 (2010) 36–63, http://dx.doi.org/10.1017/S0022112010001345.

[72] C. Lawrence, D. Adytia, E. van Groesen, Variational Boussinesq model for strongly nonlinear dispersive waves, Wave Motion 76 (2018) 78–102, http://dx.doi.org/10.1016/j.wavemoti.2017.10.009.

[73] J.W. Kim, K.J. Bai, R.C. Ertekin, W.C. Webster, A derivation of the Green-Naghdi equations for irrotational flows, J. Engrg. Math. 40 (1) (2001/05/01) 17–42, http://dx.doi.org/10.1023/A:1017541206391.

[74] J. Stoker, Water Waves, Interscience, New-York, 1957.

[75] T.B. Benjamin, P.J. Olver, Hamiltonian structure, symmetries and conservation laws for water waves, J. Fluid Mech. 125 (1982) 137–185, http://dx.doi.org/10.1017/S0022112082003292.

[76] W. Craig, P. Guyenne, C. Sulem, The water wave problem and Hamiltonian transformation theory, in: T. Bodnár, G.P. Galdi, v. Nečasová (Eds.), Springer International Publishing, Cham, 2021, pp. 113–196.

[77] C. Papoutsellis, Nonlinear Water Waves Over Varying Bathymetry: Theoretical and Numerical Study Using Variational Methods (Ph.D. thesis), National Technical University of Athens, 2016, http://dx.doi.org/10.26240/heal.ntua.2614.

[78] T. Iguchi, Isobe–Kakinuma model for water waves as a higher order shallow water approximation, J. Differential Equations 265 (3) (2018) 935–962, http://dx.doi.org/10.1016/j.jde.2018.03.019.

[79] B. Khorbatly, I. Zaiter, S. Isrwal, Derivation and well-posedness of the extended Green-Naghdi equations for flat bottoms with surface tension, J. Math. Phys. 59 (7) (2018) 071501, http://dx.doi.org/10.1063/1.5020601.

[80] V. Duchêne, T. Iguchi, A mathematical analysis of the kakinuma model for interfacial gravity waves. Part I: Structures and well-posedness, Ann. Inst. H. Poincaré, Anal. Non Linéaire 41 (2) (2023) 257–315, http://dx.doi.org/10.4171/AIHPC/82.

[81] G.B. Whitham, Variational methods and applications to water waves, Proc. R. Soc. A: Math. Phys. Eng. Sci. 299 (1456) (1967) 6–25, http://dx.doi.org/10.1098/rspa.1967.0119.

[82] D. Clamond, Explicit Dirichlet–Neumann operator for water waves, J. Fluid Mech. 950 (2022) A33, http://dx.doi.org/10.1017/jfm.2022.830.

[83] W. Choi, High-order strongly nonlinear long wave approximation and solitary wave solution 2013 CORRIGENDUM, J. Fluid Mech. 952 (2022) E2, http://dx.doi.org/10.1017/jfm.2022.928.

[84] D. Clamond, D. Dutykh, D. Mitsotakis, Conservative modified Serre–Green–Naghdi equations with improved dispersion characteristics, Commun. Nonlinear Sci. Numer. Simul. 45 (2017) 245–257, http://dx.doi.org/10.1016/j.cnsns.2016.10.009.

[85] M.F. Gobbi, J.T. Kirby, G. Wei, A fully nonlinear Boussinesq model for surface waves. Part 2. Extension to $O(kh)^4$, J. Fluid Mech. 405 (2000) 181–210, http://dx.doi.org/10.1017/S0022112099007247.

[86] G. Wei, J.T. Kirby, S.T. Grilli, R. Subramanya, A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves, J. Fluid Mech. 294 (1995) 71–92, http://dx.doi.org/10.1017/S0022112095002813.

[87] S. Neossi Nguetchue, S. Abelman, A computational algorithm for solving nearly penta-diagonal linear systems, Appl. Math. Comp. 203 (2) (2008) 629–634, http://dx.doi.org/10.1016/j.amc.2008.05.012.

[88] L.N. Trefethen, Spectral Methods in MATLAB, Society for Industrial and Applied Mathematics, 2000.

[89] D. Mitsotakis, D. Dutykh, J. Carter, On the nonlinear dynamics of the traveling-wave solutions of the Serre equations, Wave Motion 70 (2017) 166–182, http://dx.doi.org/10.1016/j.wavemoti.2016.09.008.

[90] D. Dutykh, D. Clamond, Efficient computation of steady solitary gravity waves, Wave Motion 51 (1) (2014) 86–99, http://dx.doi.org/10.1016/j.wavemoti.2013.06.007.

[91] Y. Chen, H. Yeh, Laboratory experiments on counter-propagating collisions of solitary waves. Part 1. Wave interactions, J. Fluid Mech. 749 (2014) 577–596, http://dx.doi.org/10.1017/jfm.2014.231.

[92] J. Riordan, Combinatorial Identities, John Wiley & Sons, 1968, reprint 1979.