# Semi-implicit strategies for the Serre-Green-Naghdi equations in hyperbolic form. Is hyperbolic relaxation really a good idea?

Emanuele Macca, Walter Boscheri, Mario Ricchiuto

► To cite this version:

Emanuele Macca, Walter Boscheri, Mario Ricchiuto. Semi-implicit strategies for the Serre-Green-Naghdi equations in hyperbolic form. Is hyperbolic relaxation really a good idea?. 2025. hal-05425479

HAL Id: hal-05425479

https://inria.hal.science/hal-05425479v1

Preprint submitted on 19 Dec 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
# Semi-implicit strategies for the Serre-Green-Naghdi equations in hyperbolic form. Is hyperbolic relaxation really a good idea?

Emanuele Macca<sup>××</sup>, Walter Boscheri<sup>†</sup> and Mario Ricchiuto<sup>+</sup>

<sup>×</sup> Università degli Studi di Catania, Catania, Italy <sup>*</sup> Corresponding author

<sup>†</sup> Laboratoire de Mathématiques UMR 5127 CNRS, Université Savoie Mont Blanc, France

<sup>+</sup> Centre Inria de l'Université de Bordeaux, CNRS UMR 5251, Talence, France

## Abstract

The Serre-Green-Naghdi (SGN) equations provide a valuable framework for modelling fully nonlinear and weakly dispersive shallow-water flows. However, their elliptic formulation can considerably increase the computational cost compared to the Saint-Venant equations. To overcome this difficulty, hyperbolic models (hSGN) have been proposed that replace the elliptic operators with first-order hyperbolic formulations augmented by relaxation terms, which recover the original elliptic formulation in the stiff limit. Yet, as the relaxation parameter $\lambda$ increases, explicit schemes face restrictive stability constraints that may offset these advantages.

To mitigate this limitation, we introduce a semi-implicit (SI) integration strategy for the hSGN system, where the stiff acoustic terms are treated implicitly within an IMEX Runge-Kutta framework, while the advective components remain explicit. The proposed approach mitigates the CFL stability restriction and maintains dispersive accuracy at a moderate computational cost. Numerical results confirm that the combination of hyperbolization and semi-implicit time integration provides an efficient and accurate alternative to both classical SGN and fully explicit hSGN solvers.

**Keywords:** Serre-Green-Naghdi equations; dispersive shallow-water models; hyperbolic relaxation; semi-implicit methods; IMEX Runge-Kutta schemes.

## 1 Introduction

Dispersive shallow-water models are fundamental tools in the mathematical and numerical study of free-surface flows, with applications that range from coastal engineering and tsunami propagation to laboratory experiments and reduced-order modelling [23, 28]. Among these, the Serre-Green-Naghdi (SGN) equations stand out as a reliable framework for describing fully nonlinear and weakly dispersive wave motion in shallow to intermediate water depths. Unlike the classical shallow-water (Saint-Venant) equations, the SGN system incorporates dispersive corrections, enabling a more accurate representation of wave dynamics in regimes where dispersion plays a significant role [22, 13, 1].

From a numerical perspective, two alternative approaches can be found in litteraure. One is based on the classical formulation of the equations which embeds mixed space and time derivatives [23]. More precisely such formulation embeds an elliptic nonlienar operator whose inversion is necessary to recover the depth averaged velocity [36, 2, 20, 17]. This leads to a
dramatic increase in computational cost, especially compared to the underlying shallow water equations for which efficient fully explicit conservative methods can be used.

One of the most common strategies to circumvent the elliptic step in the SGN equations is to resort to relaxation systems, often referred to as “hyperbolization”. In this case, rather than solving directly the dispersive equations one solves an augmented system of first order PDEs, which includes several auxiliary variables and appropriately designed relaxation sources which allow to approximate the dispersive operators when a relaxation parameter $\epsilon = 1/\lambda$ tends to zero. Several such approaches have been proposed over the years for dispersive equations [27, 15, 14, 6, 31, 29], as well as for other models involving high order derivatives [11, 9, 7, 10]. The first order part of the system is usually designed to be hiperbolic, so we can speak of a hyperbolic Serre-Green-Naghdi model (by some abuse of language). This hSGN system, is thus a hyperbolic set of balance laws, with stiff relaxation terms. In the asymptotic relaxation limit, the original parabolic equations are retrieved. In this way, the difficulties associated with nonlocal elliptic operators are removed, and the resulting system can be tackled with standard high-order methods for hyperbolic balanced laws. This advantage, however, comes with a significant trade-off: the augmented wave speeds typically grow with $\lambda$, so that for large relaxation parameters they can far exceed the natural shallow-water celerity $\sqrt{gh}$, leading to a much more restrictive CFL condition in explicit schemes.

The core question that this work addresses is therefore practical and quantitative: when accuracy demands push the model towards large $\lambda$, does hyperbolization still pay off? In other words, although hyperbolization removes elliptic inversions, it also increases system size and possibly the stiffness of acoustic modes, so the net computational gain is not automatically guaranteed. To answer this question we pursue a complementary strategy: instead of treating the entire hyperbolic system explicitly, we design semi-implicit (SI) time-stepping schemes that treat the stiff ($\lambda$-dependent) acoustic coupling implicitly while handling advective and non-stiff terms explicitly. This splitting aims to combine the structural advantages of the hSGN model with a time discretization that relaxes the severe hyperbolic CFL constraint and leads, in practice, to an implicit solver that is affordable and subject to a material-type CFL rather than to the overly restrictive acoustic CFL.

In this paper we present a new splitting of the hSGN equations tailored to SI integration, and we couple this splitting with a second-order IMEX Runge-Kutta time stepping technique. The combination yields a method that damps acoustic modes through the implicit sub-step and produces a linear or mildly nonlinear elliptic-type problem for the depth update that can be solved efficiently. We discuss stability and consistency properties of the splitting and we provide a careful assessment of the cost/accuracy trade-off between (i) the classical elliptic-hyperbolic SGN solver, (ii) fully explicit discretizations of the hyperbolic relaxation, and (iii) the proposed SI-hSGN approach.

Numerical experiments illustrate the practical performance of the method and delimit the regimes where SI-hSGN is advantageous. Benchmark tests include solitary-wave propagation, dispersive wave breaking scenarios and other standard tests that highlight dispersive accuracy, nonlinear interactions and long-time behaviour. Based on these results we provide pragmatic recommendations for choosing the stiffness parameter $\lambda$ and tuning the IMEX time integrator in applications of interest.

The remainder of the paper is organized as follows. We begin by recalling the standard SGN formulation and the hyperbolic relaxation variants used in our study. Then we introduce the semi-implicit splitting and detail the IMEX time discretization, together with stability considerations and an asymptotic analysis. Spatial discretization options for both explicit and SI-hSGN
solvers are discussed next. Afterwards we present numerical results that assess accuracy, stability and computational cost, and we conclude with final remarks and practical guidelines for users interested in adopting the proposed approach.

# 2 Serre-Green-Naghdi equations in hyperbolic and standard form

In this section we discuss the different forms of the Serre-Green-Naghdi (SGN) system used in this paper. The interested reader is referred to [15, 23, 14, 6, 21, 29] and references therein for more details on the SGN equations, and more generally on dispersive wave models.

## 2.1 Standard SGN equations

Let us consider a one-dimensional domain $\Omega$ defined by the position $x \in \mathbb{R}$ and let $t \in \mathbb{R}_+$ be the time coordinate. The propagation of non-linear dispersive waves on flat bathymetry can be modeled by the Serre-Green-Naghdi equations which write

$$ 
\begin{aligned}
h_t + (hu)_x &= 0, \\
(hu)_t + (hu^2 + gh^2/2 + p_{nh})_x &= 0,
\end{aligned} 
\quad (2.1) $$

where $h$ denotes the water depth, $u$ the depth averaged velocity, $g$ the gravity acceleration, and with $p_{nh}$ the non-hydrostatic pressure

$$ p_{nh} = -\frac{h^3}{3} \left(\frac{Du}{Dt}\right)_x + 2\frac{h^3}{3}(u_x)^2, \quad (2.2) $$

with $D/Dt$ denoting the material or Lagrangian derivative:

$$ \frac{D}{Dt}(\cdot) = (\cdot)_t + u(\cdot)_x. \quad (2.3) $$

The above system admits several equivalent formulations, more or less suitable for numerical approximation. To update the velocity, all of these formulations require inverting some equivalent form of the operator

$$ (hu)_t - \left(\frac{h^3}{3}u_{tx}\right)_x = \text{RHS} \quad (2.4) $$

which is embedded in (2.1)-(2.2). In this work, we investigate hyperbolic reformulations which allow in principle to avoid this inversion. These are discussed in the next sections.

## 2.2 Hyperbolic SGN (hSGN) equations without bathymetry

The hyperbolized SGN model (hSGN) obtained using an extended Lagrangian formulation can be written in the flat bathymetry case as (see [15, 6, 32])

$$ 
\left\{ 
\begin{aligned}
h_t + (hu)_x &= 0, \\
(hu)_t + (hu^2 + p)_x &= 0, \\
(h\eta)_t + (hu\eta)_x &= hw, \\
(hw)_t + (huw)_x &= -\lambda \left(\frac{\eta}{h} - 1\right),
\end{aligned} 
\right. 
\quad (2.5) $$
with

$$p := g \frac{h^2}{2} + h \frac{\lambda \eta}{3 h} \left(1 - \frac{\eta}{h}\right). \qquad (2.6)$$

In addition to the physical variables already defined, $w$ and $\eta$ are here relaxation variables and $\lambda \gg 1$ a relaxation parameter. The analysis of [12] shows that for large values of $\lambda$ then the above system provides an order $O(1/\lambda)$ approximation of the SGN model, with $\eta = h + O(1/\lambda)$.

The above system can be shown to be fully hyperbolic, with eigenvalues

$$\lambda_1 = u - c, \quad \lambda_2 = \lambda_3 = u, \quad \lambda_4 = u + c, \qquad (2.7)$$

and linearly independent eigenvectors. This allows in principle to use standard high order explicit numerical methods avoiding the inversion of second order operators.

Note that the extended celerity $c$ in the eigenvalues is defined by

$$c^2 = c_h^2 + \frac{\lambda \eta^2}{3 h^2}, \quad c_h^2 = gh, \qquad (2.8)$$

where $c_h$ is the usual shallow water celerity. Clearly, since $\eta/h = 1 + O(1/\lambda)$, one has $c \gg c_h$ in general, and especially for large $\lambda$.

For later use, it is quite interesting to note that the pressure derivative can be written as

$$\begin{aligned} p_x &= (c^2 - \lambda\eta\alpha)h_x + \lambda\alpha(h\eta)_x \\ &= a^2 h_x + \lambda\alpha(h\eta)_x \end{aligned} \qquad (2.9)$$

where

$$a^2 := c^2 - \lambda\eta\alpha, \quad \alpha := -\frac{1}{3h} \left(2\frac{\eta}{h} - 1\right). \qquad (2.10)$$

Because for $\lambda \gg 1$ we have $\eta/h = 1 + O(1/\lambda)$, it is reasonable to assume $\alpha < 0$.

## 2.3 Including bathymetric effects

The effects of bathymetry in the hyperbolic system can be accounted for via two simple modifications. A hyperbolized SGN model with bathymetry can be written according to [6, 29] as

$$\left\{ \begin{array}{l} h_t + (hu)_x = 0, \\ (hu)_t + \left(hu^2 + \frac{1}{2}gh^2 + \eta\tilde{p}\right)_x + (gh + \frac{3}{2}\tilde{p})b_x = 0, \\ (h\tilde{\eta})_t + (hu\tilde{\eta})_x = hw, \\ (hw)_t + (huw)_x = 3\tilde{p}, \\ \tilde{p} = \frac{\lambda}{3}\left(1 - \frac{\eta}{h}\right), \\ \tilde{\eta} = \eta + \frac{3}{2}b, \end{array} \right. \qquad (2.11)$$

where $b = b(x)$ is the topography. Note that, as remarked in [29], the asymptotic limit of the above system does not coincide with the usual SGN model but with a slightly different
approximation which can be obtained under classical mild-slope hypotheses. In particular, following the reasoning of [29], one can show that the limit of system (2.11) writes

$$ \left\{ \begin{array}{l} h_t + (hu)_x = 0, \\ (hu)_t + \left(hu^2 + \frac{1}{2}gh^2 + h\tilde{p}_{nh}\right)_x + (gh + \frac{3}{2}\tilde{p}_{nh})b_x = 0, \end{array} \right. \qquad (2.12) $$

where now

$$ \tilde{p}_{nh} = -\frac{h^2}{3} \left(\frac{Du}{Dt}\right)_x + \frac{hb_x}{2} \frac{Du}{Dt} + \frac{2}{3}h^2 u_x^2 + h^2 u^2 b_{xx}. \qquad (2.13) $$

Comparison with the expressions of the original SGN model, reported e.g. in [16, 23, 17], exhibits some differences. Numerical tests on classical benchmarks show that the original SGN and the above mild-slope approximation provide identical results even for rapidly varying bathymetries [29]. We will use thus (2.11) and (2.12) with no further discussion on the topic, and we refer to this system with the abbreviation hSGN. The interested reader can refer to [6, 29] for more details.

## 3 Semi-implicit time integration for the hSGN model

As previously stated in the Introduction, we aim at investigating whether hyperbolic relaxation is actually an effective modeling approach. Indeed, the main rationale behind the development of hyperbolic approximations is to avoid the inversion of the elliptic operators appearing in (2.1)-(2.2) and (2.12)-(2.13). Indeed, this is true. However, there are two points which may attenuate critically the advantage of using a hyperbolic formulation. One is the number of equations which is substantially larger (double in one dimension, and larger of a factor 5/3 in two dimensions). The second is the time step limitation. Several previous works show that dispersive terms do not necessarily affect the stability of the schemes. This is for example the case when using centered approximations, and certain residual based approaches (see e.g. [30, 16, 20] for spectral stability results). In this case, even in the presence of dispersive terms, the relevant time step limitation is the usual shallow water one, associated to the speed $|u| + c_h$. Conversely, the hyperbolic reformulation will require time-steps inversely proportional to $|u| + c$. As we saw, we expect $c$ to be much larger than $c_h$, and this is indeed the case, especially in intermediate to shallow depths where weakly dispersive corrections become important, as shown in Figure 1. We can see that, for large values of $\lambda$, this ratio may exceed one order of magnitude, and even beyond that. This means that, if the overhead in the inversion of the elliptic operator is less then a factor 10 or 20, the hyperbolic reformulation may actually be slower, as reported in some works e.g. [29, 18]. A natural question is whether, by some efficient semi-implicit (SI) time discretization, we may offset this issue in favor of the hyperbolic formulation. Indeed, designing and implementing SI and IMplicit-EXplicit (IMEX) schemes can be done in a very efficient way for first-order hyperbolic systems (see e.g. [4, 5] and references therein). This is what we propose to study here. For the ODE

$$ \frac{dQ}{dt} = R(Q, x), \qquad (3.1) $$

the simplest form of SI schemes involves some splitting of the right hand side evolution operator $R$ in a stiff part $R_I$, treated implicitly, and a non-stiff part $R_E$, treated explicitly:

$$ \frac{Q^{n+1} - Q^n}{\Delta t} = R_I(Q^{n+1}, x) + R_E(Q^n, x), \quad R_I(Q, x) + R_E(Q, x) = R(Q, x). \qquad (3.2) $$

Our first task is to derive a suitable splitting for the hSGN systems (2.5) and (2.11).
<figure><img src="image_2.png" /><figcaption>Figure 1: Celerity ratio $c/c_h$ as a function of the water depth for different values of $\lambda$.</figcaption></figure>

## 3.1 Splitting for the hSGN equations

Let us start from the model (2.5). The acoustic part of the eigenstructure is associated to the coupling between the pressure, mass conservation, and the wave equation defined by the last two equations. Following the general guideline (3.2), we thus split the system as follows:

$$ \frac{\partial}{\partial t} \begin{pmatrix} h \\ hu \\ h\eta \\ hw \end{pmatrix} = R_E^{\text{hSGN}} + R_I^{\text{hSGN}}, \qquad (3.3) $$

where

$$ R_E^{\text{hSGN}} := - \begin{pmatrix} 0 \\ hu^2 \\ hu\eta \\ huw \end{pmatrix}_x, \quad R_I^{\text{hSGN}} := - \begin{pmatrix} hu \\ p \\ 0 \\ 0 \end{pmatrix}_x - \begin{pmatrix} 0 \\ 0 \\ -hw \\ \lambda \left(\frac{\eta}{h} - 1\right) \end{pmatrix}. \qquad (3.4) $$

The above splitting can be further developed to show how its solution implies the inversion of a second order operator, allowing to dump all acoustic modes and increasing the stability range. To begin with, we introduce the intermediate update

$$ \begin{aligned} h^* &= h^n, \\ (hu)^* &= (hu)^n - \Delta t (hu^2)_x^n, \\ (h\eta)^* &= (h\eta)^n - \Delta t (hu\eta)_x^n, \\ (hw)^* &= (hw)^n - \Delta t (huw)_x^n, \end{aligned} \qquad (3.5) $$

that essentially solves the convective contribution. This allows us to write the implicit sub-system as

$$ \begin{aligned} h^{n+1} &= h^* - \Delta t (hu)_x^{n+1}, \\ (hu)^{n+1} &= (hu)^* - \Delta t p_x^{n+1}, \\ (h\eta)^{n+1} &= (h\eta)^* + \Delta t (hw)^{n+1}, \\ (hw)^{n+1} &= (hw)^* - \lambda \Delta t \left(\frac{\eta}{h} - 1\right)^{n+1}. \end{aligned} \qquad (3.6) $$
Combining the last two equations above, we can easily find that

$$ (h\eta)^{n+1} = (\widetilde{h\eta}) b_1(h^{n+1}), \quad b_1(h) := \frac{h^2}{h^2 + \lambda\Delta t^2}, \qquad (3.7) $$

having set

$$ (\widetilde{h\eta}) := (h\eta)^* + \Delta t (hw)^* + \lambda \Delta t^2. \qquad (3.8) $$

We now move on to the first equation in (3.6), which we can combine with the second one to obtain

$$ h^{n+1} = h^* - \Delta t q_x^* + \Delta t^2 p_{xx}^{n+1}. \qquad (3.9) $$

We now use (2.9) to evaluate the last term, hence obtaining

$$ h^{n+1} - \Delta t^2 (a^2 h_x^{n+1})_x = h^* - \Delta t q_x^* + \lambda \Delta t^2 (\alpha(h\eta)_x^{n+1})_x. \qquad (3.10) $$

To close the problem, we have to use (3.7). In particular, we can write

$$ (h\eta)_x^{n+1} = b_1(\widetilde{h\eta})_x + (\widetilde{h\eta}) \partial_h b_1 h_x^{n+1}. \qquad (3.11) $$

Simple computations show that

$$ \partial_h b_1 = 2\lambda\Delta t^2 \frac{h}{(h^2 + \lambda\Delta t^2)^2} = 2\frac{\lambda\Delta t^2}{h^3} b_1^2 \geq 0, \qquad (3.12) $$

so we can write

$$ (\alpha(h\eta)_x^{n+1})_x = [\alpha b_1(\widetilde{h\eta})_x + 2\alpha(\widetilde{h\eta}) \frac{\lambda\Delta t^2}{h^3} b_1^2 h_x^{n+1}]_x. \qquad (3.13) $$

We now propose to use the following linearized second order problem to evaluate the water depth at the next time step:

$$ \begin{gathered} h^{n+1} - \Delta t^2 (\kappa(h^n, u^n, \eta^n, w^n) h_x^{n+1})_x = h^* - \Delta t q_x^* + \lambda \Delta t^2 (\alpha b_1(\widetilde{h\eta})_x)_x, \quad (3.14) \\ \kappa(h^n, u^n, \eta^n, w^n) := a^2(h^n, \eta^n) + 2\alpha(h^n, \eta^n)(\widetilde{h\eta}) \frac{\lambda^2 \Delta t^2}{(h^n)^3} b_1^2(h^n) \end{gathered} $$

The nature of the last equation is not clear due to the undecided sign of the second term in $\kappa$. A more detailed study allows us to prove what follows.

**Proposition 1** (Asymptotic coercivity). *Provided $\alpha^n \le 0$, $h^n > 0$, and $\eta^n > 0$, and provided that the timestep $\Delta t$ is independent of $\lambda$, the second order equation (3.14) is asymptotically coercive in the sense that*

$$ \Delta t^2 \kappa = \Delta t^2 g h^n + \Delta t^2 \lambda (\Lambda_0 + \mathcal{O}(1/\lambda)), \quad \Lambda_0 = \frac{1}{3} \frac{(\eta^n)^2}{(h^n)^2} - \alpha^n h^n > 0. \qquad (3.15) $$

*Proof.* Using the definition of $a^2$ (eq. (2.10)) we can easily write (we omit the superscript $n$ to lighten the notation)

$$ \Delta t^2 \kappa = \Delta t^2 g h + \Delta t^2 \frac{\lambda}{3} \frac{\eta^2}{h^2} - \alpha \lambda \eta \Gamma, \qquad (3.16) $$

where

$$ \Gamma := \Delta t^2 - \frac{(\widetilde{h\eta})}{h\eta} \frac{2h^2}{\lambda} \frac{\lambda^2 \Delta t^4}{(h^2 + \lambda \Delta t^2)^2}. \qquad (3.17) $$
A Taylor development in terms of $1/\lambda$ shows that

$$ \Gamma = \Delta t^2 - \frac{1}{\lambda} 2h^2 \frac{\langle h\eta \rangle}{h\eta} \left( 1 - \frac{2h^2}{\Delta t^2} \frac{1}{\lambda} + \frac{6h^4}{\Delta t^4} \frac{1}{\lambda^2} + \mathcal{O}\left(\frac{1}{\lambda^3}\right) \right). \quad (3.18) $$

Since $\Delta t$ is independent on $\lambda$, we have $\Gamma/\Delta t^2 = 1 + \mathcal{O}(1/\lambda)$ from which we can deduce the estimate on $\Delta t^2\kappa$.

**Remark 2.** Note that for the semi-implicit scheme we will use $\Delta t \approx \Delta x/\max(|u| + \sqrt{gh})$ which is indeed independent of $\lambda$. In practice, we have found $\kappa > 0$ on all meshes used, for all benchmarks, and for all the values of $\lambda$ tested.

## 3.2 High order generalization in time

As done in the first-order case and following the approach of [26], we rewrite system (2.5) in the ODE form (3.1)–(3.2). The semi-implicit semi-discrete scheme can be written as

$$ \frac{dU}{dt} = \tilde{H}(U_E, U_I), \quad (3.19) $$

with

$$ \tilde{H}(U_E, U_I) = \begin{bmatrix} -\partial_x((hu)_I) \\ -\partial_x((hu^2)_E) & -\partial_x(p_I) \\ -\partial_x((hu\eta)_E) & -(hw)_I \\ -\partial_x((huw)_E) & \lambda\left(\frac{h\eta}{h^2}\right)_I \end{bmatrix}. \quad (3.20) $$

Based on this formulation, we apply an IMEX scheme to system (3.20). The general procedure to advance the numerical solution from time $t^n$ to $t^{n+1}$ using an s-stage Runge-Kutta IMEX method is as follows:

* **Stage values:** For $i = 1, \dots, s$, compute

$$ U_E^{(i)} = U^n + \Delta t \sum_{j=1}^{i-1} a_{i,j}^E H \left( U_E^{(j)}, U_I^{(j)} \right), \quad (3.21) $$

$$ U_I^{(i)} = U^n + \Delta t \left( \sum_{j=1}^{i-1} a_{i,j}^I H \left( U_E^{(j)}, U_I^{(j)} \right) + a_{i,i}^I H \left( U_E^{(i)}, U_I^{(i)} \right) \right). \quad (3.22) $$

* **Numerical solution:**

$$ U^{n+1} = U_I^{(s)}. $$

**Remark 3.** Note that if a system of the form (3.19) is autonomous (i.e., the right-hand side does not explicitly depend on time), and the s-stage double Butcher tableau has identical b coefficients, then the evolution requires only s evaluations of the function H. Furthermore, if the last row of the matrix A coincides with the weights b (i.e., the implicit tableau defines a stiffly accurate scheme), then the numerical solution corresponds to the last stage of the implicit part [3].
Here, we consider the IMEX scheme defined by the following double Butcher tableau [24]:

$$ \begin{array}{c|cc} & 0 & \gamma \\ c & c & 0 \\ \hline & 1 - \gamma & \gamma \end{array} \qquad \begin{array}{c|cc} & \gamma & \gamma \\ 1 & 1 - \gamma & \gamma \\ \hline & 1 - \gamma & \gamma \end{array} \qquad (3.23) $$

where $\gamma = 1 - \frac{1}{\sqrt{2}}$ and $c = \frac{1}{2\gamma}$.

In our case, applying the scheme defined in (3.23), we obtain:

1. $U_E^{(1)} = U^n$;

2. $U_I^{(1)} = U^n + \Delta t \gamma H(U_E^{(1)}, U_I^{(1)})$;

3. $U_E^{(2)} = U^n + \Delta t c H(U_E^{(1)}, U_I^{(1)})$;

4. $U_I^{(2)} = U^n + \Delta t (1 - \gamma) H(U_E^{(1)}, U_I^{(1)}) + \Delta t \gamma H(U_E^{(2)}, U_I^{(2)})$;

5. $U^{n+1} = U_I^{(2)}$.

**Remark 4.** We observe that $U_E^{(2)}$, $U_I^{(2)}$, and $U_I^{(1)}$ share a common term. Hence, steps 3 and 4 can be equivalently rewritten in efficient form as:

$$ U_E^{(2)} = \left(1 - \frac{c}{\gamma}\right) U^n + \frac{c}{\gamma} U_I^{(1)}, \qquad (3.24) $$

$$ U_I^{(2)} = \left(1 - \frac{1-\gamma}{\gamma}\right) U^n + \frac{1-\gamma}{\gamma} U_I^{(1)} + \Delta t \gamma H(U_E^{(2)}, U_I^{(2)}). \qquad (3.25) $$

### 3.3 Stability conditions

For an explicit scheme, the time step $\Delta t$ is determined at stage $t^n$ by the Courant-Friedrichs-Lewy (CFL) condition as $\Delta t = \text{CFL } \Delta x / \max_{\Omega} \mu$, where the CFL restriction is governed by the maximum spectral radius $\mu$ of the matrix $A(U)$ defining the hyperbolic system

$$ \text{CFL} = \frac{\mu_{\max}^n \Delta t^n}{\Delta x} \le C_{\text{ex}}, \qquad (3.26) $$

being $\mu_{\max} = \max_{\Omega} (u + \sqrt{gh + \lambda/3\eta^2/h^2})$, and $C_{\text{ex}}$ is a constant close to one.

In the case of our semi-implicit scheme, we empirically derive the following stability condition based on the velocity and not the celerity:

$$ \text{MCFL} = \frac{u_{\max}^n \Delta t^n}{\Delta x} \le C_{\text{im}}, \qquad (3.27) $$

where $u_{\max} = \max_{\Omega} |u|$ and $C_{\text{im}} \approx 0.75$.

This condition is considerably less restrictive than (3.26), as the condition for the classical CFL can be expressed as

$$ \text{CFL} = \frac{\mu_{\max}^n}{u_{\max}^n} \text{MCFL} \le \frac{\mu_{\max}^n}{u_{\max}^n} C_{\text{im}}, \qquad (3.28) $$

and $\mu_{\max}^n / u_{\max}^n \gg 1$ for small Froude numbers.
In the classical Shallow Water equations, the CFL condition is typically governed by the propagation speed of gravity waves, leading to a constraint of the form

$$ \Delta t \leq \text{CFL} \frac{\Delta x}{u + \sqrt{gh_{\text{max}}}}. \qquad (3.29) $$

However, for the hSGN model, the presence of higher-order dispersive terms modifies the characteristic wave speeds, potentially leading to a more restrictive CFL condition. The additional dispersive effects introduce new characteristic velocities that must be taken into account when selecting the time step. Consequently, the stability constraint is generally stricter than in the classical Shallow Water framework. This difference is particularly evident in numerical experiments, where maintaining stability in hSGN formulations often requires a lower effective CFL number than in Shallow Water models. As a result, the semi-implicit approach provides a significant advantage by allowing larger time steps while preserving stability and accuracy.

# 4 Spatial discretization

## 4.1 Explicit hSGN

The first-order explicit numerical scheme for system (2.5) reads as follows:

$$ \left\{ \begin{aligned} h_i^{n+1} &= h_i^n - \Delta t D_x(q^n)_i, \\ q_i^{n+1} &= q_i^n - \Delta t \hat{D}_x(q^n u^n)_i - \Delta t a^2(h^n, \eta^n) D_x(h^n)_i - \Delta t \lambda \alpha(h^n, \eta^n) D_x(h^n \eta^n)_i, \\ (h\eta)_i^{n+1} &= (h\eta)_i^n - \Delta t \hat{D}_x(h^n u^n \eta^n)_i + \Delta t h_i^n w_i^n, \\ (hw)_i^{n+1} &= (hw)_i^n - \Delta t \hat{D}_x(h^n u^n w^n)_i - \Delta t \lambda \left( \frac{(h\eta)_i^n}{(h_i^n)^2} - 1 \right). \end{aligned} \right. \qquad (4.1) $$

The differential operators $D_x$ and $\hat{D}_x$ applied to a given flux function $F(U)$ are respectively defined as:

* $D_x(F)_i = \frac{F_{i+\frac{1}{2}} - F_{i-\frac{1}{2}}}{\Delta x}$, in which $F_{i\pm\frac{1}{2}} = \frac{1}{2}(F(U_{i\pm1}) + F(U_i))$;

* $\hat{D}_x(F)_i = \frac{\hat{F}_{i+\frac{1}{2}} - \hat{F}_{i-\frac{1}{2}}}{\Delta x}$, where $\hat{F}_{i+\frac{1}{2}} = \frac{1}{2}(F(U_{i+\frac{1}{2}}^-) + F(U_{i+\frac{1}{2}}^+) - \alpha_{i+\frac{1}{2}}(U_{i+\frac{1}{2}}^+ - U_{i+\frac{1}{2}}^-))$ is the Rusanov flux and $\alpha_{i+\frac{1}{2}}$ is related to the eigenvalues of the explicit sub system. In our case, $\alpha \approx |u|$.

The values $U_{i+\frac{1}{2}}^{\pm}$ are computed via a conservative component-wise linear reconstruction [34]:

$$ v_{i\pm\frac{1}{2}}^{\mp} = v_i \pm v'_i \frac{\Delta x}{2}, \qquad (4.2) $$

where the reconstructed slope $v'_i$ is given by

$$ v'_{i} = \frac{v_{i+1} - v_{i-1}}{2\Delta x}, \qquad (4.3) $$

and $v_i$ represents the cell average $\int_{I_i} v(x, t) dx$.

The second order in time has been obtained by Heun time procedure [34].
## 4.2 Spatial discretization of the SI-hSGN

In the semi-implicit case, we use the semi-discrete first-order in time system<fnref n="1" /> as:

$$ egin{cases} 
\bar{h} = h - \Delta t D_x(\bar{q}) \\
\bar{q} = q - \Delta t \hat{D}_x(qu) - \Delta t a^2 D_x(\bar{h}) - \Delta t \lambda \alpha D_x(\bar{h} \eta) \\
\bar{h} \eta = h \eta - \Delta t \hat{D}_x(h u \eta) + \Delta t \bar{h} \bar{w} \\
\bar{h} \bar{w} = h w - \Delta t \hat{D}_x(h u w) - \Delta t \lambda \left( \frac{\bar{h} \eta}{\bar{h}^2} - 1 \right),
\end{cases} 
\quad (4.4) $$

where the differential operators $D_x$ and $\hat{D}_x$ have been defined above.

The fully-discrete second order scheme can be derived combining the previous Sec.3.1-3.2. In particular, the fully-discrete first order scheme is:

$$ egin{cases} 
\displaystyle hw_i^\dagger = hw_i - \frac{\Delta t}{\Delta x} (F_{i+\frac{1}{2}}^{huw} - F_{i-\frac{1}{2}}^{huw}) + \lambda \Delta t \\ 
\displaystyle h\eta_i^\dagger = h\eta_i - \frac{\Delta t}{\Delta x} (F_{i+\frac{1}{2}}^{h\eta} - F_{i-\frac{1}{2}}^{h\eta}) + \Delta t h w_i^\dagger \\ 
\displaystyle q_i^\dagger = q_i - \frac{\Delta t}{\Delta x} (F_{i+\frac{1}{2}}^{q\mu} - F_{i-\frac{1}{2}}^{q\mu}) \\ 
\beta_1 := 1 + \frac{\lambda \Delta t^2}{h^2} \\ 
\beta_2 := \frac{2\lambda \Delta t^2}{\beta_1^2 h^3} \\ 
\beta := a^2 + \lambda \alpha \beta_2 h \eta^\dagger \\ 
\delta := \frac{\alpha}{\beta_1} \\ 
\Delta[h]_i := h_i \left(1 + \frac{\Delta t^2}{\Delta x^2}(\beta_{i+\frac{1}{2}} + \beta_{i-\frac{1}{2}})\right) - \frac{\Delta t^2}{\Delta x^2}\beta_{i+\frac{1}{2}}h_{i+1} - \frac{\Delta t^2}{\Delta x^2}\beta_{i-\frac{1}{2}}h_{i-1} \\ 
h_i^\dagger = h_i - \frac{\Delta t}{2\Delta x}(q_{i+1}^\dagger - q_{i-1}^\dagger) + \frac{\lambda \Delta t^2}{2\Delta x^2}(\delta_{i+\frac{1}{2}}(h\eta_{i+1}^\dagger - h\eta_i^\dagger) - \delta_{i-\frac{1}{2}}(h\eta_i^\dagger - h\eta_{i-1}^\dagger)) \\ 
\Delta[\bar{h}]_i = h_i^\dagger \\ 
\bar{q}_i = q_i^\dagger - \frac{\Delta t}{2\Delta x}\beta_i(\bar{h}_{i+1} - \bar{h}_{i-1}) - \frac{\lambda \Delta t}{2\Delta x}\delta_i(h\eta_{i+1}^\dagger - h\eta_{i-1}^\dagger) \\ 
\bar{h}\eta_i = \frac{h\eta_i^\dagger}{1 + \frac{\lambda \Delta t^2}{\bar{h}_i^2}} \\ 
\bar{h}w_i = hw_i^\dagger - \lambda \Delta t \frac{\bar{h}\eta_i}{\bar{h}_i^2}.
\end{cases} 
\quad (4.5) $$

## 4.3 Standard SGN equations

To approximate (2.1) numerically, we will exploit the elliptic-hyperbolic form initially due to [16], and used for the Green-Naghdi and other dispersive systems in [33, 8, 19, 20, 18]. For solution purposes, we recast the system as follows

$$ 
\begin{cases} 
h_t + (hu)_x = 0 \\ 
(hu)_t + (hu^2 + gh^2/2)_x = h\phi, 
\end{cases} 
\quad (4.6) $$

<footnote n="1">¹Since we are focused on the spatial discretization and in order to simplify the readability $\bar{U} = U(x, t^{n+1})$, while $U = U(x, t^n)$.</footnote>
and now we look at the constraint that $\phi$ must verify to satisfy the Serre-Green-Naghdi equation. We note first that, using the mass equation we have

$$ \phi = \frac{Du}{Dt} + gh_x, \qquad (4.7) $$

so using the last formula of the non-hydrostatic pressure we see that

$$ p_{nh} = -\frac{h^3}{3} (\phi - gh_x)_x + 2 \frac{h^3}{3} (u_x)^2. \qquad (4.8) $$

For the last system to be equivalent to the Serre-Green-Naghdi equations we must have

$$ h\phi = -(p_{nh})_x = \left(\frac{h^3}{3}(\phi - gh_x)_x - 2\frac{h^3}{3}(u_x)^2\right)_x, \qquad (4.9) $$

which is equivalent to the constraint

$$ h\phi - \left(\frac{h^3}{3}\phi_x\right)_x = -\left(\frac{h^3}{3}gh_{xx} + 2\frac{h^3}{3}(u_x)^2\right)_x. \qquad (4.10) $$

So we can now solve the problem as follows

$$ \begin{aligned} h\phi - \left(\frac{h^3}{3}\phi_x\right)_x &= -\left(\frac{h^3}{3}gh_{xx} + 2\frac{h^3}{3}(u_x)^2\right)_x, \\ h_t + (hu)_x &= 0, \\ (hu)_t + (hu^2 + gh^2/2)_x &= h\phi, \end{aligned} \qquad (4.11) $$

marching in time by solving first for $\phi$ and then solving the shallow water equations with the additional source.

The above procedure can be easily generalized to (2.12)-(2.13), and leads to the coupled system

$$ \begin{aligned} h\phi - \left(\frac{h^3}{3}\phi_x\right)_x + \kappa(h, b)\phi &= -\left(\frac{h^3}{3}g\zeta_{xx} + 2\frac{h^3}{3}(u_x)^2 + \rho(h, u, b)\right)_x - hQ(h, u, b)b_x, \\ h_t + (hu)_x &= 0, \\ (hu)_t + (hu^2 + gh^2/2)_x &= h\phi, \end{aligned} \qquad (4.12) $$

having set

$$ \begin{aligned} \kappa(h, b) &:= \left(\frac{h^2 b_x}{2}\right)_x + 3\frac{hb_x^2}{4}, \\ \rho(h, u, b) &:= -\frac{h^2 b_x}{2}g\zeta_x + h^3 u^2 b_{xx}, \\ Q(h, u, b) &:= \frac{h}{2}g\zeta_{xx} - 3\frac{gb_x}{4}\zeta_x + hu_x^2 + \frac{3}{2}hu^2 b_{xx}. \end{aligned} \qquad (4.13) $$

The above can be marched in time by first inverting the second order PDE for $\phi$, and then performing the shallow water updated enhanced with the $h\phi$ source. The interested reader is referred to [33, 8, 19, 20, 18] for more details.
## 4.4 Explicit Scheme for the standard SGN equations

The fully-discrete explicit scheme for the standard Serre-Green-Naghdi is given by

$$ egin{cases} h_i^{n+1} = h_i^n - \frac{\Delta t}{\Delta x} (\hat{F}_{i+\frac{1}{2}}^{hu} - \hat{F}_{i-\frac{1}{2}}^{hu}), \\ hu_i^{n+1} = hu_i^n + \frac{\Delta t}{\Delta x} (\hat{F}_{i+\frac{1}{2}}^{hu^2+gh^2/2} - \hat{F}_{i-\frac{1}{2}}^{hu^2+gh^2/2}) + \Delta t \phi_i^n, \\ h\phi_i^n - \frac{1}{3\Delta x^2} (h_{i+\frac{1}{2}}^3 (\phi_{i+1}^n - \phi_i^n) - h_{i-\frac{1}{2}}^3 (\phi_i^n - \phi_{i-1}^n)) = h_i^*, \end{cases} \quad (4.14) $$

where

$$ h_i^* = -\frac{1}{6\Delta x^3} \left( gh_{i+1}^n (h_{i+2}^n - 2h_{i+1}^n + h_i^n) + \frac{(h^3)_{i+1}^n}{2} (u_{i+2}^n - u_i^n)^2 + \right. \quad (4.15) $$

$$ \left. - gh_{i-1}^n (h_i^n - 2h_{i-1}^n + h_{i-2}^n) - \frac{(h^3)_{i-1}^n}{2} (u_i^n - u_{i-2}^n)^2 \right). \quad (4.16) $$

In order to achieve second order in time, the Heun method has been applied.

# 5 Numerical tests

This section presents numerical tests to assess the reliability of the proposed model and its semi-implicit numerical treatment, highlighting its accuracy, ability to capture dispersive effects, and overall performance. First, we verify accuracy using a solitary wave (soliton) with a known exact solution, allowing us to quantify numerical errors. Then, we examine the model's ability to capture dispersive effects, a crucial aspect in wave propagation problems. Finally, we consider the evolution of Favre waves to further evaluate the method's effectiveness in complex scenarios. These experiments demonstrate the advantages of the proposed approach in terms of accuracy, stability, and computational efficiency.

The time step $\Delta t^n$ is determined by

$$ \Delta t^n = \text{CFL} \frac{\Delta x}{\mu_{\max}^n}, \quad (5.1) $$

where $\mu_{\max}^n$ denotes the maximum spectral radius of the Jacobian matrices $\partial F/\partial U$ across cells. However, $\Delta t^n$ must also satisfy the material CFL condition (MCFL, see (3.27)), which we have empirically determined to be not greater than 0.5 [25]. This introduces variability in the CFL condition, allowing it to differ from one time step to another.

The time step is thus computed in accordance with both CFL conditions, one associated with the free surface waves and the other with the material wave. Specifically, for each $n$,

$$ \bar{\Delta} t^n = \text{CFL} \frac{\Delta x}{\lambda_{\max}^n}, \quad (5.2) $$

after which the $\bar{\Delta} t$-material CFL condition is evaluated as $MCFL^* = \frac{u_{\max} \bar{\Delta} t^n}{\Delta x}$. If $MCFL^* \le MCFL$, then $\Delta t = \bar{\Delta} t$; otherwise, $\Delta t = MCFL \frac{\Delta x}{u_{\max}}$. This approach means that the CFL condition may vary within each time step. Consequently, when the CFL condition is not constant across iterations, we include a figure to display its evolution when significant.
Furthermore, we introduce distinct CFL numbers tailored to the characteristic speeds of each numerical scheme. For both the semi-implicit and explicit treatment of the hyperbolic Serre-Green-Naghdi model, the CFL condition is expressed as the ratio $\Delta t/\Delta x$ multiplied by the maximal eigenvalue $\nu_{max}$, which itself depends on the parameter $\lambda$. Concretely, we set

$$CFL_{IMEX} = \frac{\Delta t}{\Delta x} \nu_{max}, \quad CFL_{EX} = \frac{\Delta t}{\Delta x} \nu_{max}, \qquad (5.3)$$

with

$$\nu_{max} = u + \sqrt{gh + \frac{\lambda}{3} \left( \frac{\eta}{h} - 1 \right)}. \qquad (5.4)$$

Although the definitions of CFL<sub>IMEX</sub> and CFL<sub>EX</sub> are formally identical, they are subject to different stability constraints. The explicit CFL, $CFL_{EX}$, must satisfy the classical condition $CFL_{EX} < 1$ to ensure numerical stability. In contrast, the IMEX CFL, $CFL_{IMEX}$, is governed by the material stability condition and can typically exceed one. This distinction reflects the enhanced stability properties of IMEX schemes in the presence of stiff terms.

This formulation captures how the maximum wave speed in the hyperbolic system varies with $\lambda$ and the flow state, ensuring that both semi-implicit and explicit schemes respect stability constraints dictated by the most restrictive local propagation speed.

By contrast, for the standard Serre-Green-Naghdi equations we employ a more classical CFL number defined through the simpler characteristic speed $\sigma_{max}$. Specifically,

$$CFL_{SGN} = \frac{\Delta t}{\Delta x} \sigma_{max}, \quad \sigma_{max} = u + \sqrt{gh}. \qquad (5.5)$$

Here, $\sigma_{max}$ corresponds to the maximal signal velocity in the absence of the extra dispersive correction associated with $\lambda$. By distinguishing these CFL definitions, we align the time-step restrictions with the physics and numerical treatment of each variant of the Serre-Green-Naghdi model, thereby promoting both stability and efficiency across the different schemes.

In all the numerical tests, $g$ is the gravitational constant $g = 9.81 \text{ m}^2/\text{s}$ and the numerical errors have been computed in $L^1$ norm as:

$$\text{Error} = \frac{\|u_{\text{num}} - u_{\text{exact}}\|_1}{\|u_{\text{num}}\|_1}. \qquad (5.6)$$

# 6 Breakdown of a Gaussian bell

## 6.1 Flat bathymetry

Dispersive waves play a crucial role in many physical systems, particularly in shallow water models where higher-order terms introduce significant dispersive effects. These small-amplitude oscillations, commonly referred to as dispersive ripples, typically emerge behind a primary wave due to the interaction between nonlinearity and dispersion. The numerical resolution of such phenomena presents several challenges: high-frequency components may lead to spurious oscillations if not properly captured, whereas excessive numerical dissipation can suppress these effects entirely, distorting the physical behavior of the solution.

In this test, we evaluate the ability of the first- and second-order semi-implicit schemes to accurately resolve dispersive effects. The initial condition consists of a localized perturbation over a constant background state, given by:

$$h(x, 0) = h_{\infty} + Ae^{-x^2/20}, \quad h_{\infty} = 1, \quad A = 1. \qquad (6.1)$$
<figure><img src="image_3.png" /><figcaption>Figure 2: Numerical solutions for dispersive effects using second-order semi-implicit and explicit schemes for different values of $\lambda$, respectively 100 and 500. The final time is $t = 35$. The reference solution has been obtained with the fourth-order energy preserving scheme for the standard SGN model(4.11) developed in [29].</figcaption></figure>

The computational domain is set to [-200, 200] with a uniform grid of $N = 5000$ points, and the final time of the simulation is $t = 35$. The numerical stability is ensured with a CFL<sub>IMEX</sub> number of 2 for $\lambda \in \{100, 500, 1000, 5000\}$. Simultaneously, the CFL<sub>SW</sub> is, respectively, 1.34, 0.66, 0.48 and 0.22.

Figures 2a-3b show the numerical solutions obtained with the second-order semi-implicit and explicit schemes for different values of $\lambda$ (500, 1000, 5000, 10000). The results highlight the improved accuracy of the second-order method in capturing the dispersive ripples while
<figure><img src="image_4.png" /><figcaption>Figure 3: Numerical solutions for dispersive effects using second-order semi-implicit and explicit schemes for different values of $\lambda$, respectively 1000 and 5000. The final time is $t = 35$. The reference solution has been obtained with the fourth-order energy preserving scheme for the standard SGN model(4.11) developed in [29].</figcaption></figure>

maintaining numerical stability. The reference solution has been obtained with a fourth order energy preserving scheme by [29] applied to the classical SGN equations on 15000 points.

As shown in Table 1 indicates that the semi-implicit IMEX scheme offers clear computational advantages in the presence of strong dispersive effects (large $\lambda$), while the explicit scheme becomes rapidly less efficient as the required time-step diminishes. The comparison with the standard SGN model at the baseline case highlights the trade-off between model complexity (inclusion of dispersive corrections) and computational cost: even though SGN avoids the parameter $\lambda$, its
<table><caption>Table 1: Computational cost (in seconds) of semi-implicit, explicit and Serre-Green-Naghdi scheme for $\lambda \in \{100, 500, 1000, 5000\}$ adopting $N = 5000$ uniform mesh. The CFL<sub>IMEX</sub> = 2 while the CFL<sub>Exp</sub> = 0.4 and the CFL<sub>SGN</sub> = 0.9.</caption><thead><tr><th colspan="5">CPU time for Dispersive effects</th></tr><tr><th>Scheme</th><th>λ = 100</th><th>λ = 500</th><th>λ = 1000</th><th>λ = 5000</th></tr></thead><tbody><tr><td>IMEX</td><td>7.68</td><td>16.18</td><td>25.59</td><td>63.34</td></tr><tr><td>Explicit</td><td>19.51</td><td>39.938</td><td>75.86</td><td>169.73</td></tr><tr><td>SGN</td><td>12.43</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

stability constraint leads to higher CPU time than IMEX at $\lambda = 100$, and it cannot adapt to scenarios with varying $\lambda$ in the same framework.

## 6.2 Smooth bathymetric hump

In this test, we assess the performance of the semi-implicit (IMEX) and explicit schemes for the hyperbolic Serre-Green-Naghdi model with smooth (non-flat) bathymetry. The initial water depth and discharge are prescribed by

$$h_0 = 1, \quad A = 1, \\ b(x) = \frac{1}{20}(h_0 + Ae^{-2(x+50)^2}), \quad h(x, 0) = h_0 + Ae^{-x^2} - b(x), \quad u(x, 0) = 10^{-2}.$$

Here, $b(x)$ defines a smooth bottom variation and $h(x, 0)$ the initial perturbed free-surface above this bathymetry, while $u$ is a small non-zero initial velocity. The computational domain and final time are chosen consistently with previous experiments (e.g., domain $[-150, 150]$, final time $t = 35$), ensuring that waves interact with the varying bed.

We consider three values of the dispersive parameter $\lambda \in \{100, 500, 1000\}$. For each $\lambda$, the CFL conditions are set analogously to the flat-bottom case:

$$CFL_{\text{IMEX}} = 2, \quad CFL_{\text{Exp}} = 0.4.$$

Figures 4a-4c display the depth profiles at $t = 35$ for semi-implicit and explicit methods. A reference solution obtained with a fourth order energy preserving scheme by [29] applied to the classical SGN equations on 15000 points is also reported. The results highlight how dispersive ripples evolve over uneven bottom and how each scheme captures small-scale oscillations without introducing spurious artifacts or excessive dissipation.

<table><caption>Table 2: Computational cost (in seconds) of semi-implicit (IMEX) and explicit schemes for the hyperbolic dispersive model over smooth bathymetry, with $\lambda \in \{100, 500, 1000\}$, adopting $N = 5000$ uniform mesh. The CFL numbers are CFL<sub>IMEX</sub> = 2 and CFL<sub>Exp</sub> = 0.4.</caption><thead><tr><th colspan="4">CPU time for Dispersive Bathymetry</th></tr><tr><th>Scheme</th><th>λ = 100</th><th>λ = 500</th><th>λ = 1000</th></tr></thead><tbody><tr><td>IMEX</td><td>12.73</td><td>24.78</td><td>34.63</td></tr><tr><td>Explicit</td><td>36.02</td><td>72.40</td><td>109.99</td></tr></tbody></table>

Table 2 reports the CPU times for semi-implicit and explicit schemes under increasing $\lambda$. Even with smooth bathymetry, the qualitative behavior parallels the flat-bottom case: for $\lambda = 100$,
<figure><img src="image_5.png" /><figcaption>Figure 4: Numerical solutions for dispersive effects using second-order semi-implicit and explicit schemes for different values of $\lambda$, respectively 100, 500 and 1000. The final time is $t = 35$. The reference solution has been obtained with the fourth-order energy preserving scheme for the standard SGN model(4.11) developed in [29].</figcaption></figure>
the semi-implicit scheme typically outperforms the explicit one due to the larger allowable time step under CFL<sub>IMEX</sub> = 2. As λ increases to 500 and 1000, the maximum eigenvalue ν<sub>max</sub> grows, tightening the explicit time-step restriction more severely; hence the explicit scheme's CPU time rises faster than semi-implicit. The bathymetry variation also affects local wave speeds, but the semi-implicit treatment remains advantageous by handling stiff dispersive terms while accommodating depth-induced speed variability.

## 6.3 Solitary waves

The GN equations admits exact solitary wave solutions of the form

$$ 
\begin{aligned}
h_{\text{ex}} &= h_0(1 + \varepsilon \operatorname{sech}^2(\kappa(x - Ct))) \\
u_{\text{ex}} &= C \left(1 - \frac{h_0}{h_{\text{ex}}}\right)
\end{aligned} 
\qquad (6.2) 
$$

with $h_0$ the depth at rest, $\varepsilon = A/h_0$ the nonlinearity coefficient (amplitude of depth at rest), $C = \sqrt{gh_0(1 + \varepsilon)}$ the soliton speed, and $\kappa$ a shape parameter given by

$$ \kappa = \sqrt{\frac{3\varepsilon}{4h_0^2(1 + \varepsilon)}} $$

### 6.3.1 Accuracy and CPU time

To assess the accuracy of the proposed semi-implicit scheme, we consider a solitary wave (soliton) for which an exact solution is available. This allows us to quantify the numerical error and verify the expected convergence rates for both the first- and second-order schemes. The numerical tests are performed for different values of $\lambda$ (500, 1000, 5000, and 10000), with results presented in Figures (a) and (b). The CFL<sub>IMEX</sub> number is chosen as approximately 2.7 for most cases, except for $\lambda = 10000$, where it is set to 2.0. The numerical scheme follows the second version of the implicit discretization, in which $h^2$ is treated implicitly in the fourth equation. The initial conditions for the solitary wave are given by

$$ h_{\infty} = 1, \quad A = 1, \quad \varepsilon = \frac{A}{h_{\infty}}, \qquad (6.3) $$

$$ k = \sqrt{\frac{3\varepsilon}{4h_{\infty}^2(1 + \varepsilon)}}, \quad c = \sqrt{gh_{\infty}(1 + \varepsilon)}, \qquad (6.4) $$

where $g = 9.81$ is the gravitational acceleration. The initial profiles for water height and velocity are defined as

$$ h(x, 0) = h_{\infty} \left(1 + \varepsilon \operatorname{sech}^2(kx)\right), \qquad (6.5) $$

$$ u(x, 0) = c \left(1 - \frac{h_{\infty}}{h(x, 0)}\right), \qquad (6.6) $$

meanwhile $\eta(x, 0) = h(x, 0)$ and $w = -h(x, 0)u_x(x, 0)$.

The numerical results, shown in Figures 5a-5d, confirm the expected accuracy of the scheme. The first-order method exhibits an error scaling consistent with $O(\Delta x)$, while the second-order scheme achieves an error proportional to $O(\Delta x^2)$. These behaviors remain consistent across
<figure><img src="image_6.png" /><figcaption>(a) First order grid convergence: for the depth $h$ (b) First order grid convergence: for the velocity $u$ ($N_{\text{Points}}$-Error, log scale)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) Second order xgrid convergence: for the depth (d) Second order grid convergence: for the velocity $h$ ($N_{\text{Points}}$-Error, log scale)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: Test 6.3.1 $L^1$ errors vs number of points, $L^1$ errors vs $\lambda$ and $L^1$ errors vs CPU time expressed in second. Numerical accuracy analysis: 5a-5c error comparison for first- and second-order schemes for $h$; 5b-5d error comparison for first- and second-order schemes for $u$; while 5e error reduction for increasing $\lambda$ and 5f error vs CPU time expressed in second. The final time is 1.</figcaption></figure>

different values of $\lambda$, demonstrating the robustness of the numerical approach. Additionally, Figure 5e illustrates the error reduction as $\lambda$ increases, and Figure 5f the computational cost comparison further validating the semi-implicit formulation. These results highlight the accuracy and efficiency of the proposed scheme for capturing solitary wave dynamics with high fidelity.
The exact solution is given by:

$$h(x, t) = h_{\infty} \left(1 + \varepsilon \operatorname{sech}^2(kx - ct)\right), \qquad (6.7)$$

$$u(x, t) = c \left(1 - \frac{h_{\infty}}{h(x, t)}\right). \qquad (6.8)$$

<table><caption>Table 3: Computational cost (in seconds) of first- and second-order semi-implicit and explicit schemes for $\lambda = 1000$ and with $N = 2000$. The CFL<sub>IMEX</sub> = 2.5 while the CFL<sub>Exp</sub> = 0.4.</caption><thead><tr><th colspan="3">hSGN - CPU time</th></tr><tr><th>Scheme</th><th>Order 1</th><th>Order 2</th></tr></thead><tbody><tr><td>IMEX</td><td>1.2502</td><td>4.0350</td></tr><tr><td>Explicit</td><td>3.2105</td><td>16.9138</td></tr></tbody></table>

The computational efficiency of the semi-implicit (IMEX) schemes is evident when compared to their explicit counterparts. While explicit methods require a significantly smaller CFL number (0.4) to ensure stability, the IMEX approach remains stable even with a much larger CFL value (2.5). This results in a dramatic reduction in computational cost, as seen in Table 3. In particular, the first-order semi-implicit method is approximately 2.5 times faster than its explicit counterpart, while the second-order one exhibits a 4-fold reduction in computational sense. These results highlight not only the superior stability properties of the semi-implicit approach but also its remarkable efficiency, making it a compelling choice for large-scale simulations.

## 6.4 Favre waves

We consider a smooth hydraulic jump problem on the interval $[x_a, x_b] = [-50, 50]$ under constant gravity $g = 9.81$. The still-water depth is set to $h_0 = 0.2$, and we perform three experiments with jump amplitudes $\varepsilon \in \{1.1, 1.2, 1.3\}$. The upstream depth is

$$h_1 = \varepsilon h_0,$$

and the corresponding jump velocity is chosen to satisfy the Rankine-Hugoniot condition,

$$u_1 = \sqrt{\frac{g(h_1 + h_0)}{2h_0 h_1}} (h_1 - h_0).$$

We introduce a smooth transition of width $\alpha = 5h_0$ by means of the auxiliary functions

$$g(x) = 1 - \tanh(\frac{x}{\alpha}),$$

so that the initial water depth and velocity profiles read

$$h(x) = h_0 + \frac{1}{2}(h_1 - h_0)g(x), \quad u(x) = \frac{1}{2}u_1 g(x), \quad u_x(x) = \frac{1}{2}u_1 g'(x).$$

The discontinuities in $h$ and $u$ across the jump are smoothed by the hyperbolic tangent profile. Figures 6-8 present the numerical solutions for the three values of $\varepsilon$. In each subplot, the second-order explicit and semi-implicit solutions are shown together with the numerical solution of the Serre-Green-Naghdi model (4.11) at different final time $t$ where $t(T) = T\sqrt{gh_0}$ with $T = [50, 60, 70]$ for different values of $\lambda \in \{100, 500\}$ adopting $N = 5000$ uniform mesh and
<figure><img src="image_9.png" /><figcaption>(a) $T = 50$, $\epsilon = 1.1$ and $\lambda = 100$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c) $T = 60$, $\epsilon = 1.1$ and $\lambda = 100$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(e) $T = 70$, $\epsilon = 1.1$ and $\lambda = 100$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) $T = 50$, $\epsilon = 1.1$ and $\lambda = 500$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(d) $T = 60$, $\epsilon = 1.1$ and $\lambda = 500$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(f) $T = 70$, $\epsilon = 1.1$ and $\lambda = 500$.</figcaption></figure>

Figure 6: Numerical solutions of Favre waves using second-order semi-implicit and explicit schemes, and second-order scheme for the Serre-Green-Naghdi model (4.11), at different final times with nonlinearity coefficient $\epsilon = 1.1$ and $\lambda \in \{100, 500\}$. The solutions are plotted with the dimensionless variables $x/h_0$ on the horizontal axis and $(h-h_0)/h_0$ on the vertical axis. Reference: nonlinear potential flow solution by [36]. The CFL numbers are CFL<sub>IMEX</sub> = 4, CFL<sub>EX</sub> = 0.4 and CFL<sub>SGN</sub> = 0.9.

the CFL numbers are CFL<sub>IMEX</sub> = 4, CFL<sub>EX</sub> = 0.4 and CFL<sub>SGN</sub> = 0.9. Both schemes, in all configuration, consistently exhibit comparable solution with the standard Serre-Green-Naghdi problem, allowing it to accurately maintain the plateau at $\epsilon - 1$, which is a distinctive feature of Favre waves. The results compare well with the reference solution by [36] obtained with a fully nonlinear potential flow solver.

The effect of the nonlinearity coefficient $\epsilon$ is evident in the results. As $\epsilon$ increases, the waves exhibit more pronounced steepening, leading to sharper gradients and a more defined plateau. The second-order schemes consistently captures these characteristics with high accuracy, highlighting its robustness and ability to resolve and preserve intricate wave structures.

### 6.4.1 CPU time benchmarking

In these configurations, we conduct a detailed assessment of the computational cost associated with each scheme and benchmark these costs against those of the canonical SGN model (4.11). Specifically, Table 4 and Figures 6-8 illustrate that, for sufficiently small values of the para-
<figure><img src="image_15.png" /><figcaption>(a) $T = 50$, $\epsilon = 1.2$ and $\lambda = 100$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) $T = 60$, $\epsilon = 1.2$ and $\lambda = 100$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(e) $T = 70$, $\epsilon = 1.2$ and $\lambda = 100$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b) $T = 50$, $\epsilon = 1.2$ and $\lambda = 500$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(d) $T = 60$, $\epsilon = 1.2$ and $\lambda = 500$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(f) $T = 70$, $\epsilon = 1.2$ and $\lambda = 500$.</figcaption></figure>

Figure 7: Numerical solutions of Favre waves using second-order semi-implicit and explicit schemes, and second-order scheme for the Serre-Green-Nagdhi model (4.11), at different final times with nonlinearity coefficient $\epsilon = 1.1$ and $\lambda \in \{100, 500\}$. The solutions are plotted with the dimensionless variables $x/h_0$ on the horizontal axis and $(h-h_0)/h_0$ on the vertical axis. Reference: nonlinear potential flow solution by [36]. The CFL numbers are $CFL_{IMEX} = 4$, $CFL_{EX} = 0.4$ and $CFL_{SGN} = 0.9$.

meter $\lambda$, the semi-implicit scheme formulated for the hyperbolic SGN model achieves solution accuracy on par with both the original SGN discretization and a high-resolution reference solution. Remarkably, despite this comparable level of fidelity, the semi-implicit approach delivers a substantial speed-up in terms of CPU time. These results underscore not only the robustness of the semi-implicit scheme in capturing the essential dynamics of the standard SGN equations, but also its enhanced computational efficiency, rendering it an attractive and practical alternative for large-scale or real-time applications.
<figure><img src="image_21.png" /><figcaption>(a) $T = 50$, $\epsilon = 1.3$ and $\lambda = 100$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(c) $T = 60$, $\epsilon = 1.3$ and $\lambda = 100$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(e) $T = 70$, $\epsilon = 1.3$ and $\lambda = 100$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b) $T = 50$, $\epsilon = 1.3$ and $\lambda = 500$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(d) $T = 60$, $\epsilon = 1.3$ and $\lambda = 500$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(f) $T = 70$, $\epsilon = 1.3$ and $\lambda = 500$.</figcaption></figure>

Figure 8: Numerical solutions of Favre waves using second-order semi-implicit and explicit schemes, and second-order scheme for the Serre-Green-Nagdhi model (4.11), at different final times with nonlinearity coefficient $\epsilon = 1.1$ and $\lambda \in \{100, 500\}$. The solutions are plotted with the dimensionless variables $x/h_0$ on the horizontal axis and $(h - h_0)/h_0$ on the vertical axis. Reference: nonlinear potential flow solution by [36]. The CFL numbers are $CFL_{IMEX} = 4$, $CFL_{EX} = 0.4$ and $CFL_{SGN} = 0.9$.

### 6.4.2 Comparison with experimental data

To validate the performance of the semi-implicit hyperbolic relaxation scheme against the fully explicit one we now consider extended propagation distances until the dimensionless final time

$$ \tilde{t} = t \sqrt{\frac{g}{h_0}} = 256, $$

representative of the experiments by [35].

The computations are now run on the spatial domain $x \in [0, 73.58]$, initial jump at $x_0 = 68.58$, and hyperbolic relaxation parameter $\lambda = 500$.

Figure 9 presents the free-surface profile at $\tilde{t} = 256$, overlaying numerical solutions with the
<table><caption>Table 4: CPU time (in second) for different values of $\lambda$ and $\varepsilon$ for the semi-implicit and explicit second-order schemes applied to hyperbolic SGN and explicit second-order scheme applied to original SGN model. The CFL numbers are CFL<sub>IMEX</sub> = 4, CFL<sub>EX</sub> = 0.4 and CFL<sub>SGN</sub> = 0.9.</caption><thead><tr><th colspan="6">CPU time for Favre waves</th></tr><tr><th colspan="3">λ = 500</th><th colspan="3">λ = 100</th></tr><tr><th>ε</th><th>SI</th><th>EX</th><th>SI</th><th>EX</th><th>SGN</th></tr></thead><tbody><tr><td rowspan="3">0.1</td><td>7,09</td><td>173,28</td><td>0,63</td><td>10,65</td><td>2,16</td></tr><tr><td>8,50</td><td>207,51</td><td>0,76</td><td>13,00</td><td>2,45</td></tr><tr><td>9,81</td><td>241,80</td><td>0,85</td><td>15,51</td><td>3,13</td></tr><tr><td rowspan="3">0.2</td><td>7,70</td><td>180,10</td><td>1,03</td><td>10,23</td><td>2,26</td></tr><tr><td>7,83</td><td>198,22</td><td>1,34</td><td>14,79</td><td>3,27</td></tr><tr><td>8,88</td><td>232,00</td><td>1,53</td><td>17,53</td><td>4,07</td></tr><tr><td rowspan="3">0.3</td><td>7,61</td><td>162,67</td><td>2,50</td><td>10,19</td><td>2,60</td></tr><tr><td>8,54</td><td>210,23</td><td>3,03</td><td>12,68</td><td>3,33</td></tr><tr><td>10,26</td><td>267,89</td><td>3,79</td><td>15,68</td><td>3,98</td></tr></tbody></table>

<figure><img src="image_27.png" /><figcaption>Figure 9: Favre waves: heights of the first peak (circles) and trough (squares) at $\tilde{t} = 256$ as a function of the Froude number. Solid line with circles: data by [35].</figcaption></figure>

laboratory measurements by [35]. The numerical solutions of both schemes agree closely with the experimental data, with maximum elevation errors below 2% and phase deviations within one grid cell. However, the semi-implicit relaxation method allows a tenfold increase in the stable time step compared to the explicit solver, leading to a reduction in total CPU time of approximately 70% for this test. This substantial gain in computational efficiency, combined with the high-fidelity agreement to laboratory measurements, underscores the advantage of the semi-implicit scheme for long-time simulations of weakly dispersive free-surface flows.
### 6.4.3 Wave propagation over a shelf

In this section we examine the evolution of a solitary wave of amplitude $A = 0.2\text{m}$ as it traverses a submerged shelf. The shelf causes the still-water depth to vary from the upstream depth $h_0 = 1\text{m}$ to a shallower downstream depth of $0.5\text{m}$ via a ramp of slope 1 : 20. The computational domain extends from $x = 0$ to $x = 280\,	ext{m}$, and initially at time $t = 0$ the wave crest is positioned at $x_0 = 80\text{m}$.

<figure><img src="image_28.png" /><figcaption>Figure 10: Evolution of the solitary wave as it propagates over the submerged shelf, showing reflection and splitting into three pulses for the explicit scheme.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 11: Evolution of the solitary wave as it propagates over the submerged shelf, showing reflection and splitting into three pulses for the semi-implicit scheme.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 12: Comparison of the solitary wave at the final time as it propagates over the submerged shelf, showing reflection and splitting into three pulses for the semi-implicit and explicit schemes.</figcaption></figure>

The bathymetric profile is described by a piecewise linear function:

$$b(x) = \begin{cases} 0, & x \le 130, \\ \frac{x-130}{20}, & 130 < x < 140, \\ 0.5, & x \ge 140, \end{cases} \qquad (6.9)$$

Accordingly, the still-water depth is given by $h(x, 0) = h_0 - b(x)$.

Upstream of the shelf ($x \le 130, m$), we impose an exact solitary wave solution of the enhanced Boussinesq system:

$$h(x, 0) = h_0 [1 + \varepsilon \operatorname{sech}^2(k(x - x_0))], \qquad (6.10)$$

$$u(x, 0) = c \left(1 - \frac{h_0}{h(x, 0)}\right). \qquad (6.11)$$

We discretize the domain using a uniform grid with $N = 5000$ and advance in time with CFL= 3 and $\lambda = 500$.

As the incident soliton encounters the shoal, a small reflected wave is generated upstream, and the primary pulse splits into three transmitted pulses downstream, each of increasing amplitude.

Figures 10,11 and 12 illustrate the free-surface elevation after interaction, clearly showing the separated pulses for the explicit and semi-implicit scheme showing similar qualitative results.

# 7 Dingemann's experiment

The Dingemans flume experiment was originally conceived to investigate nonlinear wave shoaling and transformation over a gently sloping submerged shoal under controlled laboratory conditions. In the physical setup, a monochromatic wavetrain of small amplitude $A = 0.02$ m propagates in a flume of constant still-water depth $h_0 = 0.8$ m. The shoal consists of a smooth upslope
rising to 0.6 m above the undisturbed bottom, followed by a short plateau and a symmetrical downslope returning to the original depth. Four capacitance gauges located at

$$x = 3.04, 9.44, 20.04, \text{ and } 26.04 \text{ m}$$

record the time series of free-surface elevation over a duration of 70 s, providing high-resolution

<figure><img src="image_31.png" /><figcaption>Figure 13: Time-series comparison at gauge locations. Semi-implicit SGN (solid), explicit SGN (dashed), experimental data (symbols).</figcaption></figure>

experimental data for model validation.

To numerically reproduce these measurements, we solve the one-dimensional hyperbolic SGN equations on the domain $x \in [-140, 100] \text{ m}$. The initial condition imposes a sinusoidal surface displacement

$$z(x, 0) = \begin{cases} h_0 + A \cos(kx), & -\frac{34.5\pi}{k} \le x \le -\frac{8.5\pi}{k}, \\ h_0, & \text{otherwise,} \end{cases}$$

with $k = 0.8406220896381442 \text{ rad/m}$ determined by solving $gk \tanh(kh_0) = \omega^2$ for the given wave frequency and $h_0 = 0.8 \text{ m}$. The corresponding velocity is initialized as

$$u(x, 0) = \sqrt{\frac{g}{k} \tanh(kh_0)} \frac{z(x, 0) - h_0}{h_0}.$$

The shoal bathymetry $b(x)$ is defined piecewise to match the laboratory profile, and the initial water depth is $h(x, 0) = z(x, 0) - b(x)$.

Spatial discretization uses a uniform mesh with $N = 15000$ points, and time integration proceeds until $t_{final} = 70 \text{ s}$ under a CFL = 3, with $\lambda = 100$.

Figure 13 presents the free-surface elevation time series at all six gauges. Both numerical solutions closely follow the experimental measurements in phase and amplitude. No systematic
bias is observed between the two schemes, confirming that the hyperbolic SGN model faithfully
captures the wave dynamics even with the moderate relaxation parameter.

Remarkably, despite advancing with time steps three times larger, the semi-implicit scheme retains the same high fidelity as the explicit method when reproducing the experimental wave trains.

Figure 13 presents the time series of the free-surface elevation recorded data. Both numerical solutions exhibit an excellent agreement with the experimental measurements, accurately reproducing both the amplitude and phase of the wave trains. No systematic bias can be identified between the two schemes, confirming that the hyperbolic SGN model is able to faithfully capture the main features of the wave dynamics, even when using a moderate value of the relaxation parameter.

Remarkably, even when advancing with time steps three times larger, the semi-implicit scheme maintains the same level of accuracy as the explicit method in reproducing the experimental signals. This result highlights the enhanced stability of the semi-implicit formulation, which allows for a substantial relaxation of the time-step constraint without any loss of physical fidelity. In addition, the computational cost remains comparable, with total CPU times of $CPU_{IMEX} = 5.1806 \times 10^3$ and $CPU_{EX} = 5.7232 \times 10^3$. From a performance perspective, the semi-implicit solver completes the simulation in slightly less time than the explicit scheme. This gain arises from the enlarged time-step $\Delta t$, which reduces the total number of updates by a factor of three while maintaining identical accuracy. Such a computational advantage is critical for large-scale or real-time wave-prediction applications, where both precision and efficiency are essential.

# 8 Conclusions

In this paper we have designed a semi-implicit finite volume scheme for the solution the Serre-Green-Naghdi (SGN) equations written in first-order hyperbolic form. The stiff part of the system, which includes relaxation terms and pressure fluxes, has undergone an implicit discretization, while the convection sub-system and the remaining terms have been discretized explicitly. The novel scheme has been compared in terms of both accuracy and efficiency against the classical mixed elliptic-hyperbolic solver for the SGN model as well as against a fully explicit method for the hyperbolic SGN system. Several benchmarks have been performed to verify accuracy and robustness, and performing simulations of solitary wave propagation and dispersive wave breaking phenomena. This ultimately allows us to provide some pragmatic considerations for the effective use of semi-implicit schemes in the context of dispersive wave propagation models. Further developments will consider the extension of the proposed methodology to multidimensional unstructured meshes in order to capture realistic landscape geometries and the design of high order schemes in space and time to increase even further the accuracy and to reduce the error related to the propagation of dispersive waves.

# Acknowledgements

This research has received funding from the European Union's NextGenerationUE - Project: Centro Nazionale HPC, Big Data e Quantum Computing, "Spoke 1" (No. CUP E63C22001000006). E. Macca was partially supported by: GNCS No. CUP E53C24001950001 Research Project "Soluzioni Innovative per Sistemi Complessi: Metodi Numerici e Approcci Multiscala"; PRIN 2022 PNRR "FIN4GEO: Forward and Inverse Numerical Modeling of hydrothermal systems in
volcanic regions with application to geothermal energy exploitation", No. P2022BNB97; PRIN
2022 "Efficient numerical schemes and optimal control methods for time-dependent partial dif-
ferential equations", No. 2022N9BM3N - Finanziato dall'Unione europea - Next Generation EU
– CUP: E53D23005830006. E. Macca would like to thank the Italian Ministry of Instruction,
University and Research (MIUR) to support this research with funds coming from PRIN Pro-
ject 2022 (2022KA3JBA, entitled "Advanced numerical methods for time dependent parametric
partial differential equations and applications"). E. Macca and W. Boscheri are members of the
INdAM Research group GNCS. MR is a member of the CARDAMOM research team, Inria at
University of Bordeaux research center.

**Conflict of interest**

The authors declare that they have no conflict of interest.

**References**

[1] P. Bonneton, E. Barthelemy, F. Chazel, R. Cienfuegos, D. Lannes, F. Marche, and M. Tissier. Recent advances in serre-green naghdi modelling for wave transformation, breaking and runup processes. *European Journal of Mechanics - B/Fluids*, 30(6):589-597, 2011. (Cited page 1.)

[2] P. Bonneton, F. Chazel, D. Lannes, F. Marche, and M. Tissier. A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model. *Journal of Computational Physics*, 230(4):1479-1498, 2011. (Cited page 1.)

[3] S. Boscarino, F. Filbet, and G. Russo. High order semi-implicit schemes for time dependent partial differential equations. *Journal of Scientific Computing*, 68(8):975-1001, 2016. (Cited page 8.)

[4] S. Boscarino and G. Russo. Asymptotic preserving methods for quasilinear hyperbolic systems with stiff relaxation: a review. *SeMA*, 81:3-49, 2024. (Cited page 5.)

[5] Sebastiano Boscarino, Lorenzo Pareschi, and Giovanni Russo. *Implicit-Explicit Methods for Evolutionary Partial Differential Equations*. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2024. (Cited page 5.)

[6] S. Busto, M. Dumbser, C. Escalante, N. Favrie, and S. Gavrilyuk. On High Order ADER Discontinuous Galerkin Schemes for First Order Hyperbolic Reformulations of Nonlinear Dispersive Systems. *Journal of Scientific Computing*, 87:1-47, 2021. (Cited page 2, 3, 4 et 5.)

[7] S. Busto, M. Dumbser, S. Gavrilyuk, and K. Ivanova. On Thermodynamically Compatible Finite Volume Methods and Path-Conservative ADER Discontinuous Galerkin Schemes for Turbulent Shallow Water Flows. *Journal of Scientific Computing*, 88, 2021. (Cited page 2.)

[8] A. Cauquis, M. Ricchiuto, and P. Heinrich. Lax-wendroff schemes with polynomial extrapolation and simplified lax-wendroff schemes for dispersive waves: A comparative study. *Water Waves*, pages 1-33, 2022. (Cited page 11 et 12.)
[9] F. Dhaouadi, M. Dumbser, and S. Gavrilyuk. A first-order hyperbolic reformulation of the cahn–hilliard equation. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 481(2312):20240606, 2025. (Cited page 2.)

[10] F. Dhaouadi, N. Favrie, and S. Gavrilyuk. Extended Lagrangian approach for the defocusing nonlinear Schrödinger equation. *Studies in Applied Mathematics*, 142(3):336–358, 2019. (Cited page 2.)

[11] F. Dhaouadi and S. Gavrilyuk. An eulerian hyperbolic model for heat transfer derived via hamilton’s principle: analytical and numerical study. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 480(2283):20230440, 2024. (Cited page 2.)

[12] Vincent Duchêne. Rigorous justification of the Favrie-Gavrilyuk approximation to the Serre-Green-Naghdi model. *Nonlinearity*, 32(10):3772, 2019. (Cited page 4.)

[13] Vincent Duchêne and Christian Klein. Numerical study of the serre-green-naghdi equations and a fully dispersive counterpart. *Discrete and Continuous Dynamical Systems - B*, 27(10):5905–5933, 2022. (Cited page 1.)

[14] C. Escalante and T.M. de Luna. A general non-hydrostatic hyperbolic formulation for Boussinesq dispersive shallow flows and its numerical approximation. *J.Sci.Comput.*, 83(62), 2020. (Cited page 2 et 3.)

[15] N Favrie and S Gavrilyuk. A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves. *Nonlinearity*, 30(7):2718, 2017. (Cited page 2 et 3.)

[16] A.G. Filippini, M. Kazolea, and M. Ricchiuto. A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up. *Journal of Computational Physics*, 310:381–417, 2016. (Cited page 5 et 11.)

[17] S. Gavrilyuk and K.-M. Shyue. 2D Serre-Green-Naghdi equations over topography: El- liptic operator inversion method. *Journal of Hydraulic Engineering*, 150(1):04023054, 2024. (Cited page 1 et 5.)

[18] Sergey Gavrilyuk and Mario Ricchiuto. A geometrical green-naghdi-type system for dispersive-like waves in prismatic channels. *Journal of Fluid Mechanics*, 1014:A5, 2025. (Cited page 5, 11 et 12.)

[19] B. Jouy, D. Violeau, M. Ricchiuto, and M. Le. One dimensional modelling of Favre waves in channels. *Applied Mathematical Modelling*, 133:170–194, 2024. (Cited page 11 et 12.)

[20] M. Kazolea, A.G. Filippini, and M. Ricchiuto. Low dispersion finite volume/element discretization of the enhanced green-naghdi equations for wave propagation, breaking and runup on unstructured meshes. *Ocean Modelling*, 182:102157, 2023. (Cited page 1, 5, 11 et 12.)

[21] Maria Kazolea and Mario Ricchiuto. Full nonlinearity in weakly dispersive Boussinesq models: Luxury or necessity. *Journal of Hydraulic Engineering*, 150(1):04023061, 2024. (Cited page 3.)
[22] C. Klein and J.C. Saut. *Nonlinear Dispersive Equations: Inverse Scattering and PDE Methods*. Springer Nature Link, 2021. (Cited page 1.)

[23] D Lannes. Modeling shallow water waves. *Nonlinearity*, 33(5):R1, 2020. (Cited page 1, 3 et 5.)

[24] E. Macca, S. Avgerinos, M.-J. Castro-Diaz, and G. Russo. A semi-implicit finite volume method for the Exner model of sediment transport. *Journal of Computational Physics*, 499, 2024. (Cited page 9.)

[25] E. Macca and S. Boscarino. Semi-implicit-Type Order-Adaptive CAT2 Schemes for Systems of Balance Laws with Relaxed Source Term. *Communications on Applied Mathematics and Computation*, 2024. (Cited page 13.)

[26] E. Macca and G. Russo. Boundary effects on wave trains in the exner model of sedimental transport. *Bollettino dell'Unione Matematica Italiana*, 17(2):417-433, 2024. (Cited page 8.)

[27] A. Mazaheri, M. Ricchiuto, and H. Nishikawa. A first-order hyperbolic system approach for dispersion. *Journal of Computational Physics*, 321:593-605, 2016. (Cited page 2.)

[28] Detlev Müller. Dispersive shallow waters. *Wave Motion*, 17(1):89-95, 1993. (Cited page 1.)

[29] H. Ranocha and M. Ricchiuto. Structure-preserving approximations of the Serre-Green-Naghdi equations in standard and hyperbolic form. *Numerical methods for partial differential equations*, 41(4), 2025. (Cited page 2, 3, 4, 5, 15, 16, 17 et 18.)

[30] Mario Ricchiuto and Andrea Gilberto Filippini. Upwind residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries. *Journal of Computational Physics*, 271:306-341, 2014. (Cited page 5.)

[31] S. Tkachenko, S. Gavrilyuk, and J. Massoni. Extended Lagrangian approach for the numerical study of multidimensional dispersive waves: Applications to the Serre-Green-Naghdi equations. *Journal of Computational Physics*, 477:111901, 2023. (Cited page 2.)

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

[33] D. Torlo and M. Ricchiuto. Model order reduction strategies for weakly dispersive waves. *Mathematics and Computers in Simulation*, 205:997-1028, 2023. (Cited page 11 et 12.)

[34] E.F. Toro. *Riemann Solvers and Numerical Methods for Fluid Dynamics*. Springer, third edition, 2009. (Cited page 10.)

[35] A. Treske. Undular bores (Favre-waves) in open channels - Experimental studies. *Journal of Hydraulic Research*, 32(3):355-370, May 1994. (Cited page 24 et 25.)

[36] G. Wei, J.T. Kirby, S.T. Grilli, and R. Subramanya. A fully nonlinear Boussinesq model for surface waves. part 1. highly nonlinear unsteady waves. *Journal of Fluid Mechanics*, 294:71-92, 1995. (Cited page 1, 22, 23 et 24.)