# A new open boundary condition for Boussinesq-type models, applied to irregular wave fields

Mathieu Rigal, Philippe Bonneton, David Lannes

► To cite this version:

Mathieu Rigal, Philippe Bonneton, David Lannes. A new open boundary condition for Boussinesq-type models, applied to irregular wave fields. 2025. hal-05068186v2

HAL Id: hal-05068186

https://hal.science/hal-05068186v2

Preprint submitted on 1 Sep 2025

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

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

<figure><img src="image_1.png" /></figure>
A new open boundary condition for Boussinesq-type
models, applied to irregular wave fields*

Mathieu Rigal<sup>a,*</sup>, Philippe Bonneton<sup>b</sup>, David Lannes<sup>c</sup>

<sup>a</sup>Sorbonne Université, LJLL, UMR 7598, 4 place Jussieu, 75005, Paris, France

<sup>b</sup>Université de Bordeaux, CNRS, UMR 5805 EPOC, Allée Geoffroy Saint-Hilaire - CS
50023, 33615, Pessac, France

<sup>c</sup>Université de Bordeaux, CNRS, UMR 5251, 351 Cours de la
libération, 33405, Talence, France

Abstract

We present a novel approach to handle open boundary conditions for a Boussinesq-
type wave model coupled with the nonlinear shallow water equations. Traditional
methods for managing open boundaries — such as sponge layers and source
functions — are computationally intensive and require ad hoc calibration. To
address this, we reformulate the Boussinesq equations as a system of conservation
laws with nonlocal flux and a rapidly decaying source term. This reformulation is
adapted to generate waves at the boundary of the numerical domain, from surface
elevation data in situations where both incoming and outgoing waves are present.
The proposed numerical scheme employs a MacCormack prediction-correction
strategy combined with finite volume and finite difference methods, preserving
key physical properties and ensuring stability. Comparison with laboratory
experiments demonstrates that our approach avoids boundary reflection issues.
In particular, it is able to accurately reproduce infragravity waves associated
with a random wave field propagating over a sloping beach. This work opens
important perspectives for improving phase-resolving coastal wave models, with
the aim of forecasting complex random wave conditions in natural environments.

Keywords: Boussinesq-type equations, open boundary condition, nonlocal reformulation, irregular wave fields, hybrid coupling, MacCormack scheme

*This material is based upon work supported by ANR BOURGEONS (ANR-23-CE40-0014) and ERC BLOC.

*Corresponding author

Email addresses: mathieu.rigal@math.u-bordeaux.fr (Mathieu Rigal),
philippe.bonneton@u-bordeaux.fr (Philippe Bonneton),
david.lannes@math.u-bordeaux.fr (David Lannes)
# 1. Introduction

When considering Boussinesq-type models, the issue of boundary conditions is a significant challenge. To this day the behavior of solutions near the extremities of the domain is not well understood, hindering the design of a robust handling of boundary conditions that avoids stability issues and boundary reflections. Yet this is a prerequisite to perform realistic simulations of transient and irregular wave fields in natural coastal environments, by accurately generating the incoming waves and absorbing the outgoing ones. More precisely, the presence of high-order derivatives in dispersive models prevents the use of standard hyperbolic boundary treatment through Riemann invariants, and other strategies have to be used to deal with open boundaries.

Several techniques have been proposed to address this issue. Among the most widely used are sponge layers (Israeli and Orszag [8], 1981), which attenuate outgoing waves via artificial damping regions, and the source function method (Wei, Kirby, and Sinha [25], 1999), which generates waves internally using an added source term. When used in combination, these methods offer a practical framework for managing open boundary conditions. However, each approach has inherent limitations. Extending the spatial domain to include a sponge layer can lead to a significant increase in computational cost. On the other hand, calibrating the source function requires prior knowledge of the solution, which restricts its applicability to predefined scenarios.

In this paper, we propose a new approach for handling open boundary conditions at the entrance of the numerical domain in Boussinesq-type models that avoids these limitations. It is built on a recent strategy [12] developed by Lannes and Weynans (2020), where the Boussinesq equations restricted to a flat bathymetry are recast as a system of conservation laws with nonlocal flux and a rapidly decaying source term. Our main contribution is to extend this method to a Boussinesq-Abbott model with an optimized dispersion relation, coupled with the nonlinear shallow water equations (NSW) to handle wave breaking, and applicable over general bathymetries — therefore enabling the simulation of natural coastal environments beyond simple academic test cases. The open ocean boundary condition considered here is a prescription of the surface elevation at the entrance of the numerical domain. Such data are typically obtained from field measurements, which inherently contain both incoming and outgoing wave signals. To assess the performance of the model, we conduct wave simulations and compare the results with experimental observations.

The remainder of the paper is structured as follows. Section 2 introduces the Boussinesq model under consideration and its coupling with the NSW equations. Section 3 details the reformulation of the model into a nonlocal form suitable for enforcing open boundary conditions. Section 4 presents a numerical discretization of the proposed model and discusses implementation aspects. Finally, in Section 5, we validate the approach through wave simulations, comparing the results with experimental data. Additional technical details and mathematical proofs are provided in the appendices for completeness.
## 2. Physical model

In the nearshore area, water waves experience complex transformations, including shoaling, breaking, and run-up — processes that can be particularly challenging to model accurately. Capturing these phenomena requires accounting for both dispersive and nonlinear effects, as they are key drivers of the transformations mentioned. Boussinesq-type models offer a compelling balance between physical accuracy and computational efficiency in this context. However, their use becomes problematic in the surf zone, where steep gradients can trigger nonphysical oscillations due to the dispersive terms. In contrast, the nonlinear shallow water (NSW) model naturally incorporates dissipation, making it better suited for simulating wave breaking. To leverage the strengths of both approaches, many studies have adopted a hybrid modeling strategy, coupling a Boussinesq-type model with the NSW equations to capture the full spectrum of nearshore dynamics, see for instance Bonneton et al. [4] (2011). In this section, we begin by introducing the Boussinesq-type model under consideration, and we then detail the coupling strategy.

<figure><img src="image_2.png" /><figcaption>Figure 1: Configuration of the fluid domain with variable bathymetry and open ocean boundary condition on the left.</figcaption></figure>

### 2.1. Improved dispersion over a flat bottom

As a starting point we consider the following Boussinesq-type model with improved dispersion over a flat bathymetry:

$$ (1) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ (1 - \frac{\beta}{3} h_0^2 \partial_{xx}^2) \partial_t q + \partial_x (f_{\text{NSW}} - \frac{\beta-1}{3} g h_0^3 \partial_{xx}^2 \zeta) = 0, \end{cases} $$

with $\zeta$ the elevation, $h_0$ the depth at rest, $h = h_0 + \zeta$ the water height, $q$ the depth-averaged horizontal discharge, and $f_{\text{NSW}}(\zeta, q)$ the flux defined by

$$ (2) \qquad f_{\text{NSW}}(\zeta, q) = \frac{1}{2} g (h^2 - h_0^2) + \frac{1}{h} q^2. $$
The optimization parameter $\beta \ge 1$ is used to tune the associated dispersion relation which is given by

$$ (3) \qquad \omega^2 = gh_0k^2 \frac{1 + \frac{\beta-1}{3}h_0^2k^2}{1 + \frac{\beta}{3}h_0^2k^2}. $$

In [15], Madsen et al. derived a new form of the Boussinesq equations in order to improve the linear dispersion characteristics, and in [16] a dispersion correction is introduced. In our case we can minimize the error in terms of the phase and group velocity between (1) and the linear Stokes equations — whose dispersion relation is given by $\omega^2 = gk \tanh(kh_0)$ — by setting $\beta$ to 1.159.

In the next section we discuss the coupling between system (1) and the NSW equations.

## 2.2. Coupling with NSW

In this study, we adopt a simplified approach to surf zone wave modeling based on the nonlinear shallow water (NSW) equations. As demonstrated by Bonneton [3], the NSW system is capable of capturing the nonlinear transformation and energy dissipation of broken waves in the inner surf zone by representing wave fronts as shock waves. This motivates a hybrid modeling strategy in which dispersive effects are locally deactivated, and the NSW equations are employed in regions where broken waves are detected. This hybrid coupling technique was initially introduced by Tonelli and Petti [22] in the context of Boussinesq-type models. Their approach relied on a simple switching criterion based on the wave amplitude-to-depth ratio $A/h_0$, where values exceeding 0.8 triggered the suppression of dispersive terms. While effective for idealized test cases, this criterion lacks robustness for irregular wave fields. To address this limitation, more sophisticated detection and termination criteria were later proposed by Tissier et al. [21], enabling the individual tracking of broken waves with improved accuracy. These advancements have extended the applicability of hybrid approaches to realistic coastal environments. Further developments in the numerical treatment of wave breaking within weakly nonlinear Boussinesq models can be found in [22, 23, 9, 18], while the extension to the fully nonlinear Boussinesq framework (i.e. the Serre-Green-Naghdi equations) was explored in [4, 21, 19].

For simplicity, we consider the surf zone to be fixed in time, and thus the transition area between Boussinesq and NSW models is also fixed (see Figure 1). However, the strategy outlined below is readily extendable to a time-dependent transition zone, for example to account for tidal variability over long durations. In the aforementioned studies, the hybrid strategy was implemented numerically using a time-splitting approach, where the contributions from the NSW and dispersive components are treated separately. In regions where dispersive effects are to be suppressed, the dispersive step is simply omitted. In contrast, our approach does not rely on time splitting. Instead, we incorporate a smooth, decreasing step function $\sigma(x)$ directly into the model. This function satisfies $0 \le \sigma \le 1$, and is such that $\sigma = 0$ fully deactivates dispersive effects, while $\sigma = 1$ retains them entirely. A smooth transition between dispersive and non-dispersive
dynamics is important since it mitigates numerical artifacts such as spurious reflections and oscillations near the coupling interface — a point further discussed in Section 4.3. The use of a cut-off function to combine various models (including dispersive ones) has been considered by Parisot [17] through the framework of thick interface coupling.

Over a flat bathymetry, our modified model reads

$$ (4) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ (1 - \frac{\beta}{3} h_0^2 \partial_x (\sigma \partial_x \cdot)) \partial_t q + \partial_x (f_{\text{NSW}} - \frac{\beta-1}{3} \sigma g h_0^3 \partial_{xx}^2 \zeta) = 0. \end{cases} $$

Note that the dispersive operator $1 - \frac{\beta}{3}h_0^2\partial_x(\sigma\partial_x\cdot)$ from above is invertible. Setting $\sigma \equiv 1$, one obviously recovers (1); whereas setting $\sigma \equiv 0$, one recovers the NSW system as expected. Next we extend this coupled model to the framework of a varying bathymetry with vanishing water depth at the shoreline.

### 2.3. Case of a varying bottom with possibly vanishing depth

To account for a varying bathymetry parameterized by $b(x)$, we introduce the depth at rest $h_b = h_0 - b$. We consider the following set of PDEs:

$$ (5) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ [1 + \beta h_b \mathcal{T}_\sigma] \partial_t q + \partial_x \tilde{f}_\beta = -g \zeta \partial_x \tilde{b}, \end{cases} $$

where the second order operator $\mathcal{T}_\sigma(\cdot)$ is defined as

$$ (6) \qquad \mathcal{T}_\sigma(\cdot) = S^*(\sigma h_b S(\cdot)) + \frac{\sigma (\partial_x b)^2}{4} h_b(\cdot), $$

with $h_b S$ and the adjoint operator $S^*$ given by

$$ (7) \qquad h_b S(\cdot) = -\frac{h_b^2}{\sqrt{3}} \partial_x \left(\frac{\cdot}{h_b}\right) + \sqrt{3} \frac{\cdot}{2} \partial_x b, \quad S^*(\cdot) = \frac{1}{\sqrt{3}h_b} \partial_x (h_b \cdot) + \sqrt{3} \frac{\cdot}{2} \frac{\partial_x b}{h_b}. $$

Since the definition of $\mathcal{T}_\sigma(\cdot)$ involves a division by $h_b$ and that $h_b$ may vanish, it is necessary to assume that $\sigma = 0$ when $h_b$ vanishes, that is, that (5) then reduces to the NSW equations. We therefore always assume that

$$ (8) \qquad h_b > 0 \quad \text{on} \quad \text{supp } \sigma. $$

The optimization parameter $\beta$ and the smooth step function $\sigma$ introduced earlier are featured in the modified flux and bathymetry terms in the following way:

$$ (9) \qquad \begin{cases} \tilde{f}_\beta = f_{\text{NSW}}(\zeta, q) - f_{\text{NSW}}(0, 0) + (1-\beta) \frac{gh_b^2}{3} \left( \sigma h_b \partial_{xx}^2 \zeta - \frac{3}{2} \zeta \partial_x (\sigma \partial_x b) \right), \\ \tilde{b} = b + \frac{1-\beta}{2} h_b^2 \partial_x (\sigma \partial_x b). \end{cases} $$

Note that in (5) the topography source term is expressed in terms of $\zeta$ and not $h$, which is possible since we subtract $f_{\text{NSW}}(0, 0)$ to $f_{\text{NSW}}(\zeta, q)$ in the definition (9). The reason for using this pre-balanced writing is to later obtain a numerical scheme that will naturally preserve the hydrostatic equilibrium (lake-at-rest steady state $\zeta = q = 0$). We refer to Appendix A for a formal justification of the model (5). In the next section we will investigate how the equivalent nonlocal rewriting from [12] can be adapted to (5).
### 3. Boundary conditions and coupling with the NSW model

In practice, we will consider the coupled system (5) over some bounded domain in the configuration shown in Figure 1. The domain is delimited on the right by the shoreline which is numerically handled by the NSW equations ($\sigma = 0$ in (5)). On the left extremity of the domain, located for instance at $x = 0$, we impose the boundary condition

$$ (10) \qquad \zeta(t, x = 0) = g(t), $$

where $g$ is known. This is an open boundary condition because this time series results both from incoming and outgoing waves.

A key difficulty arises when dispersive effects are active at the left boundary (i.e., $\sigma(0) = 1$), as the inversion of the dispersive operator $1 + h_b T$ requires knowledge of the time derivative of the discharge, $\partial_t q(t, 0)$. However, the boundary condition (10) provides information only about the free surface elevation, $\zeta(t, 0)$. As a result, the discharge $q(t, 0)$ constitutes missing data, for which no straightforward recovery method exists. Importantly, a boundary condition is not directly imposed on the discharge because, in many practical applications, available measurements — typically obtained from buoys or wave gauges — are limited to surface elevation, with no direct access to the discharge.

It was shown in [11] how to solve the Boussinesq-Abbott system on a finite interval, with boundary conditions at both extremities. The essential step is based on the method of [12] which provides a way to reconstruct the “missing data” $q$ at the boundary. The situation here is slightly different, for at least two reasons. Firstly, we work with a Boussinesq model with improved frequency dispersion; when compared with $f_{NSW}$, the modified flux term $\tilde{f}_{\beta}$ provided by (9) contains an additional second derivative of $\zeta$. Because of this we cannot directly implement the same steps as in [11] to reconstruct the missing data; however, performing approximations consistent with the weakly nonlinear and weakly dispersive regime, we are able to provide a generalization of this method in the present case. The second difference with respect to [11] is that there is no boundary condition at the right extremity of the numerical domain. Instead, the Boussinesq model is coupled to the NSW equations that are solved until the physical boundary (the shoreline). This coupling is performed using a smooth cut-off function $\sigma$ that cancels the dispersive terms in the surf zone; the presence of this cut-off function must be taken into account in the treatment of the nonlocal terms involved in the analysis of the boundary condition at the entrance of the numerical domain.

#### 3.1. The Boussinesq-Abbott model over a flat bathymetry

We recall that the original Boussinesq-Abbott system for a flat bathymetry is expressed as

$$ (11) \qquad \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0 \\ (1 - \frac{1}{3} h_0^2 \partial_{xx}^2) \partial_t q + \partial_x f_{NSW}(\zeta, q) = 0 \quad \text{in } \mathbb{R}_+ \end{array} \right. $$
We are interested in a configuration where the surface elevation is imposed at the entrance of the numerical domain, located at $x = 0$:

$$ (12) \qquad \zeta(t, 0) = g(t), $$

and we also need to provide an initial condition:

$$ (13) \qquad (\zeta, q)(0, x) = (\zeta^{\text{in}}, q^{\text{in}})(x), $$

where $g, \zeta^{\text{in}}, q^{\text{in}}$ are some given datum. We briefly explain how to deal with this boundary condition following the procedure of [12] on the half-line $x \in \mathbb{R}_+$. To this end we introduce $R_0$ and $R_1$ defined as the inverses of $(1 - \frac{1}{3}h_0^2\partial_{xx}^2)$ with homogeneous Dirichlet and homogeneous Neumann boundary conditions respectively:

$$ (14) \qquad R_0 : f \mapsto u \quad \text{such that} \quad \left\{ \begin{array}{l} (1 - \frac{1}{3}h_0^2\partial_{xx}^2)u = f \quad \text{in } \mathbb{R}_+, \\ u(0) = 0, \end{array} \right. $$

$$ (15) \qquad R_1 : f \mapsto v \quad \text{such that} \quad \left\{ \begin{array}{l} (1 - \frac{1}{3}h_0^2\partial_{xx}^2)v = f \quad \text{in } \mathbb{R}_+, \\ \partial_x v(0) = 0. \end{array} \right. $$

It is then possible to recast the momentum equation from (11) as

$$ \partial_t q + R_0 \partial_x f_{\text{NSW}} = \dot{q}|_{x=0} \mathfrak{s}, $$

rapidly decaying function $\mathfrak{s}(x)$, that can be understood as a dispersive boundary layer, defined through

$$ (16) \qquad \left\{ \begin{array}{l} (1 - \frac{1}{3}h_0^2\partial_{xx}^2)\mathfrak{s} = 0 \\ \mathfrak{s}(0) = 1 \end{array} \right. \quad \Longrightarrow \quad \mathfrak{s}(x) = \exp\left(-\frac{\sqrt{3}}{h_0}x\right); $$

using the relation $R_0\partial_x = \partial_x R_1$ (see Lemma 1 from [12] for a justification) we recover a conservative writing of the interior momentum equation with nonlocal flux $R_1 f_{\text{NSW}}$:

$$ (17) \qquad \partial_t q + \partial_x R_1 f_{\text{NSW}} = \dot{q}|_{x=0} \mathfrak{s}. $$

Applying $\partial_x$ to (17), using $\partial_{xt}^2 q = -\partial_{tt}^2 \zeta$ and $\partial_{xx}^2 R_1 = \frac{3}{h_0^2}(R_1 - \text{id})$, taking the trace at $x = 0$ and dividing by $\mathfrak{s}'(0)$, one also finds the following equation for the missing data $q|_{x=0}$:

$$ (18) \qquad \begin{aligned} \dot{q}|_{x=0} &= \frac{1}{\mathfrak{s}'(0)} \left( -\ddot{g} + \frac{3}{h_0^2} \left( (R_1 f_{\text{NSW}})|_{x=0}(\zeta, q) - f_{\text{NSW}}(g, q|_{x=0}) \right) \right) \\ &=: \mathcal{Q}(q|_{x=0}, g, \ddot{g}, (R_1 f_{\text{NSW}})|_{x=0}). \end{aligned} $$

The important thing here is that the right-hand side depends on $\zeta$ and $q$ only through $(R_1 f_{\text{NSW}})|_{x=0}$; this latter quantity can be numerically evaluated with
interior values of $\zeta$ and $q$. To conclude, it can be shown that if the compatibility conditions below are satisfied

$$ egin{cases} \zeta_{|t=0}(0) = g(0), \\ -\partial_x q_{|t=0}(0) = \dot{g}(0), \end{cases} $$

then the initial boundary value problem (11), (12) and (13) is equivalent to the initial value problem consisting in the same initial condition and the following reformulated interior equations:

$$ (19) \quad \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0 \\ \partial_t q + \partial_x R_1 f_{\text{NSW}} = Q(q_{|_{x=0}}, g, \ddot{g}, (R_1 f_{\text{NSW}})_{|_{x=0}}) \mathfrak{s} \end{array} \right. \quad \text{in } \mathbb{R}_+ $$

This nonlocal formulation has several advantages. First, it provides knowledge on the time evolution of the missing data $q|_{x=0}$ through equation (18); it constitutes an ODE on $(\zeta, q)$ in some functional space, which is locally well-posed as a consequence of the Cauchy-Lipschitz theorem [12]; at last it allows for an efficient numerical implementation of the boundary conditions without resorting to sponge layers or the source function method [11]. This feature is particularly advantageous to deal with natural coastal environments featuring both incoming and outgoing waves, as it eliminates the need for calibration and doesn't require to extend the computational domain — thereby significantly reducing computational costs. In our context, this efficiency is important for analyzing the energy spectra of irregular wave fields, which require long-duration simulations to yield meaningful results.

## 3.2. Adaptation for the improved model over a flat bottom

We proceed to derive a nonlocal formulation of the improved model (4) by adapting the strategy presented in Section 3.1. To keep things simple, we assume that $\sigma$ is a decreasing function equal to one in the left part of the numerical domain (i.e. the shoaling zone, where the Boussinesq model is used) and to zero in the right part (i.e. the surf zone up to the shoreline, where the NSW equations are used); the area where $0 < \sigma < 1$ corresponds to the transition zone between both models, see Figure 1. We assume that the entrance of the numerical domain is located at $x = 0$.

The first step is to generalize the definition of the operator $R_0$ given in (14) to take into account the transition function $\sigma$ present in the dispersive operator and the new terms related to the optimization coefficient $\beta$. We now have:

$$ (20) \quad R_0 : f \mapsto u \quad \text{such that} \quad \left\{ \begin{array}{l} (1 - \frac{\beta}{3} h_0^2 \partial_x (\sigma \partial_x \cdot)) u = f \quad \text{in } \mathbb{R}_+, \\ u(0) = 0, \end{array} \right. $$

and one can check that the original definition (14) of $R_0$ is recovered by setting $\sigma \equiv 1$ and $\beta = 1$. Similarly, we shall extend the definition (16) of the dispersive boundary layer function $\mathfrak{s}(x)$ so as to now verify

$$ (21) \quad \left\{ \begin{array}{l} (1 - \frac{\beta}{3} h_0^2 \partial_x (\sigma \partial_x \cdot)) \mathfrak{s} = 0, \\ \mathfrak{s}(0) = 1, \end{array} \right. $$
and which does not admit an analytical expression in general. We have the following formula which we prove in Appendix B.1:

$$ (22) \qquad \partial_t q = -R_0 \partial_x \tilde{f}_\beta(\zeta, q) + \dot{q}|_{x=0} \mathfrak{s}. $$

Next we show that the operator $R_0 \partial_x$ can be recast under a conservative form $\partial_x R_1$ where we define $R_1$ as

$$ (23) \qquad R_1 : f \mapsto v \quad \text{such that} \quad \left\{ \begin{array}{l} (1 - \frac{\beta}{3} h_0^2 \sigma \partial_{xx}^2) v = f \quad \text{in } \mathbb{R}_+, \\ \partial_x v(0) = 0. \end{array} \right. $$

Note that the expression $(1 - \frac{\beta}{3} h_0^2 \sigma \partial_{xx}^2)$ appearing above is not self-adjoint; nevertheless this problem can be made well posed under a reasonable assumption on $\sigma$, namely, $|\sigma'| \le \sqrt{3/(\beta h_0^2)} \sqrt{\sigma}$, see Appendix B.2. We then show the identity

$$ (24) \qquad \partial_x R_1 = R_0 \partial_x $$

remains true for the generalized definitions (20) and (23) of $R_0$ and $R_1$. Let $v = R_1 f$ for some function $f$, then by construction of $R_1$ there holds

$$ \left\{ \begin{array}{l} (1 - \frac{\beta}{3} h_0^2 \sigma \partial_{xx}^2) v = f \\ \partial_x v(0) = 0 \end{array} \right. \implies \left\{ \begin{array}{l} (1 - \frac{\beta}{3} h_0^2 \partial_x (\sigma \partial_{x\cdot})) \partial_x v = \partial_x f \\ \partial_x v(0) = 0 \end{array} \right. $$

where we applied $\partial_x$ to the first equation and used $\partial_x(\sigma \partial_{xx}^2 \cdot) = (\partial_x(\sigma \partial_{x\cdot})) \partial_x$. Therefore $\partial_x v = R_0 \partial_x f$ and thus $\partial_x R_1 f = R_0 \partial_x f$. Using identity (24) in the momentum equation (22) we get:

$$ (25) \qquad \partial_t q + \partial_x R_1 \tilde{f}_\beta = \dot{q}|_{x=0} \mathfrak{s}. $$

The final step is to obtain an equation on the missing data $q|_{x=0}$. As for the standard Boussinesq-Abbott equations above, we apply $\partial_x$ to the reformulated momentum equation (25) and substitute $\partial_{xx}^2 q$ by $-\partial_{tt}^2 \zeta$:

$$ -\partial_{tt}^2 \zeta + \partial_{xx}^2 R_1 \tilde{f}_\beta = \dot{q}|_{x=0} \mathfrak{s}'. $$

Multiplying by $\sigma$ and remarking that (23) implies the identity $\frac{\beta}{3} h_0^2 \sigma \partial_{xx}^2 R_1 = R_1 - \text{id}$ we get:

$$ -\sigma \partial_{tt}^2 \zeta + \frac{3}{\beta h_0^2} (R_1 - \text{id}) \tilde{f}_\beta = \sigma \dot{q}|_{x=0} \mathfrak{s}'. $$

Taking the trace at $x=0$ and using the assumption $\sigma(0) = 1$ we have:

$$ -\ddot{g} + \frac{3}{\beta h_0^2} (R_1 \tilde{f}_\beta)|_{x=0}(\zeta, q) - \frac{3}{\beta h_0^2} \tilde{f}_\beta(g, q|_{x=0}) = \dot{q}|_{x=0} \mathfrak{s}'(0). $$

Recall that the definition (9) of $\tilde{f}_\beta$ involves the quantity $gh_0 \partial_{xx}^2 \zeta$. Such a term was not present in the previous work [12, 11], and to be able to conduct a similar
derivation of an ODE on the missing data, we suggest to approximate it by $-\partial_{tt}^2\zeta$, since at first order in $A/h_0$ and $h_0^2/\lambda^2$, where $A$ and $\lambda$ are the characteristic amplitude and wavelength of the waves, $\zeta$ solves a wave equation of speed $\sqrt{gh_0}$. With this in mind, the approximated trace equation at $x=0$ reads:

$$ -\ddot{g} + \frac{3}{\beta h_0^2} (R_1 \tilde{f}_\beta)|_{x=0}(\zeta, q) - \frac{3}{\beta h_0^2} \left( f_{\text{NSW}}(g, q|_{x=0}) + \frac{\beta - 1}{3} h_0^2 \ddot{g} \right) = \dot{q}|_{x=0} \mathfrak{s}'(0), $$

hence

$$ (26) \quad \begin{aligned} \dot{q}|_{x=0} &= \frac{1}{\mathfrak{s}'(0)} \left( \left( \frac{\beta-1}{\beta} - 1 \right) \ddot{g} + \frac{3}{\beta h_0^2} (R_1 \tilde{f}_\beta)|_{x=0} - \frac{3}{\beta h_0^2} f_{\text{NSW}}(g, q|_{x=0}) \right) \\ &=: Q^\dagger(q|_{x=0}, g, \ddot{g}, (R_1 \tilde{f}_\beta)|_{x=0}). \end{aligned} $$

As before, the right-hand side depends on $\zeta$ and $q$ only through $(R_1 \tilde{f}_\beta)|_{x=0}$; which is evaluated numerically given the interior values of $\zeta$ and $q$. Gathering the previous steps, we conclude that an approximate nonlocal reformulation of the improved model (1) is

$$ (27) \quad \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \partial_x R_1 \tilde{f}_\beta = Q^\dagger(q|_{x=0}, g, \ddot{g}, (R_1 \tilde{f}_\beta)|_{x=0}) \mathfrak{s}. \end{array} \right. $$

**Remark 3.1.** Taking the trace of the momentum equation from (27) at $x=0$, one recovers the ODE (26) on the missing data $q|_{x=0}$ since $\partial_x R_1 \tilde{f}_\beta$ vanishes at that point by construction of $R_1$ and $\mathfrak{s}(0) = 1$. A similar property will be satisfied in the variable bathymetry case treated below.

### 3.3. Extension to a varying bathymetry

We now aim to provide an approximate nonlocal reformulation of the improved model (5) with varying bathymetry. As above the first step is to introduce $R_0$ the inverse of the dispersive operator $1 + \beta h_b \mathcal{T}_\sigma$ with homogeneous Dirichlet condition at $x=0$:

$$ (28) \quad R_0 : f \mapsto u \quad \text{such that} \quad \left\{ \begin{array}{l} [1 + \beta h_b \mathcal{T}_\sigma] u = f \quad \text{in } \mathbb{R}_+, \\ u(0) = 0, \end{array} \right. $$

where $\mathcal{T}_\sigma$ has been defined in (6). We also extend the definition of $R_1$ which maps any function $f$ to the solution $u$ of the problem below:

$$ (29) \quad \left\{ \begin{array}{l} \left[ 1 + \frac{\beta}{h_b} \sigma(h_b S) \left( \frac{h_b}{\alpha_b} (h_b S)^* (\cdot) \right) \right] u = f, \quad \text{in } \mathbb{R}_+ \\ (h_b S)^* u(0) = 0, \end{array} \right. $$

where $S(\cdot)$ and $(h_b S)^*(\cdot) = S^*(h_b \cdot)$ have been defined through (7), and where

$$ (30) \quad \alpha_b = 1 + \frac{\sigma\beta}{4} (\partial_x b)^2; $$
As above in the case of a flat bathymetry, the elliptic operator (29) used to define $R_1$ is not self-adjoint; we show however in Proposition C.1 that it can be inverted under the condition $|\sigma'| \le \sqrt{3\alpha_b/(\beta h_b^2)}\sqrt{\sigma}$. The latter definition also allows one to recover a commuting identity similar to (24) even when the bathymetry is not flat (see Proposition C.5).

Finally we introduce the following quantity

$$ (31) \qquad B(\zeta, q) = -g\zeta\partial_x\tilde{b} + \frac{3}{2}\frac{\partial_x b}{h_b}\tilde{f}_\beta(\zeta, q), $$

and the nonlocal terms

$$ (32) \qquad f_b = \frac{h_b^2}{\alpha_b} R_1 \left( \frac{\tilde{f}_\beta}{h_b^2} \right), \qquad j_b = \left( \frac{\partial_x \alpha_b}{\alpha_b} + \frac{3}{2} \frac{\partial_x b}{h_b} \right) f_b, $$

with $\tilde{f}_\beta$ and $\tilde{b}$ given in (9). It is then shown in Appendix C.2 that the momentum equation from (5) can be rewritten under the following form:

$$ (33) \qquad \partial_t q + \partial_x f_b = R_0 B(\zeta, q) - j_b + \dot{q}|_{\zeta=0} \mathfrak{s}, $$

where $\mathfrak{s}(x)$ is the solution of the problem

$$ (34) \qquad \begin{cases} [1 + \beta h_b \mathcal{T}_\sigma]\mathfrak{s} = 0, & \text{in } \mathbb{R}_+, \\ \mathfrak{s}(0) = 1. & \end{cases} $$

To conclude we propose to recover the missing boundary data $q|_{x=0}(t)$ by the following equation, which we explain how to obtain in Appendix C.3:

$$ (35) \qquad \begin{aligned} \dot{q}|_{0} &= \frac{1}{\mathfrak{s}'(0)} \left[ -\frac{1}{\beta}\ddot{g}_0 + \frac{3}{\beta}\left(R_1\left(\frac{\tilde{f}_\beta}{h_b^2}\right)\right) - (\partial_x R_0 B) \right. \\ & \qquad \left. - \frac{3}{\beta h_b(0)^2}\left(f_{\text{NSW}}(\zeta, q) - f_{\text{NSW}}(0,0) - (1-\beta)\frac{g h_b^2}{2}\zeta \partial_x(\sigma \partial_x b)\right) \right] \Big|_{x=0} \\ &=: Q^*(q|_{x=0}, g, \ddot{g}, (\partial_x R_0 B)|_{x=0}, (R_1(\frac{\tilde{f}_\beta}{h_b^2}))|_{x=0}). \end{aligned} $$

The final problem is then

$$ (36) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \partial_x f_b = R_0 B(\zeta, q) - j_b + Q^*(q|_{x=0}, g, \ddot{g}, (\partial_x R_0 B)|_{x=0}, (R_1(\frac{\tilde{f}_\beta}{h_b^2})))|_{x=0} \mathfrak{s}. \end{cases} $$

# 4. Numerical scheme

## 4.1. Discrete operators

We assume that there exists $\ell > 0$ such that the solution is a dry state to the right of this point: $x \ge \ell$ implies that for all times $t$ there holds $h(t, x) =$
$q(t, x) = 0$. Hence it is sufficient to discretize the domain $(0, \ell)$; we denote $N$ the number of cells all of size $\Delta x = \ell/(N - 1)$ and whose center is $x_i = (i - 1)\Delta x$, $1 \le i \le N$. The approximation of the quantities of interest in cell $i$ and at time $t^n = n\Delta t$ shall be written $U_i^n = (\zeta_i^n, q_i^n)$.

We introduce $\delta_x$ the second order accurate discretization of $\partial_x$ such that for any $v \in \mathbb{R}^N$,

$$ (37) \quad (\delta_x v)_i = \begin{cases} \frac{1}{2\Delta x}(v_{i+1} - v_{i-1}) & \text{if } 2 \le i \le N-1, \\ \frac{1}{2\Delta x}(-3v_1 + 4v_2 - v_3) & \text{if } i = 1, \\ \frac{1}{2\Delta x}(v_{N-2} - 4v_{N-1} + 3v_N) & \text{if } i = N. \end{cases} $$

We set $b_i = b(x_i)$, $h_{b,i} = h_0 - b_i$ and we approximate $S(\cdot)$ defined in (7) by

$$ \underline{S}(v)_i = -\frac{h_{b,i}}{\sqrt{3}} \delta_x \left(\frac{v}{h_b}\right)_i + \frac{\sqrt{3}}{2} \delta_x(b)_i v_i, $$

which remains valid in the cells where $h_{b,i} > 0$. To discretize $(h_b S)^* = S^*(h_b)$ we remark that its definition (7) is equivalent to

$$ (h_b S)^*(\cdot) = \frac{\alpha_b}{\sqrt{3}h_b} \partial_x \left(\frac{h_b^2 v}{\alpha_b}\right) + \left(\frac{h_b}{\sqrt{3}} \frac{\partial_x \alpha_b}{\alpha_b} + \sqrt{3} \frac{\partial_x b}{2}\right)(\cdot), $$

the advantage of this form being that it makes (C.3) conservative due to the coefficient $\frac{\sqrt{3}h_b}{\alpha_b}$ in front of $(h_b S)^*$. Therefore we propose the following discretization:

$$ (\underline{h_b S})^*(v)_i = \frac{\alpha_{b,i}}{\sqrt{3}h_{b,i}} \delta_x \left(\frac{h_b^2 v}{\alpha_b}\right)_i + \left(\frac{h_{b,i}}{\sqrt{3}} \frac{\delta_x \alpha_b}{\alpha_{b,i}} + \frac{\sqrt{3}}{2}(\delta_x b)_i\right) v_i, $$

where we use

$$ \alpha_{b,i} = 1 + \sigma(x_i) \frac{\beta}{4} ((\delta_x b)_i)^2. $$

Finally we can define $\underline{R}_0, \underline{R}_1$ the discretization of the nonlocal operators $R_0, R_1$:

$$ (38) \quad \underline{R}_0 : f \in \mathbb{R}^N \mapsto v = \underline{R}_0 f \in \mathbb{R}^N \\ \text{solving} \quad \begin{cases} \alpha_{b,i} v_i + \beta h_{b,i} (\underline{h_b S})^* (\sigma \underline{S}(v))_i = f_i, & \forall i \in [2, N] \quad \text{s.t. } \sigma_i > 0 \\ v_i = f_i, & \forall i \in [2, N] \quad \text{s.t. } \sigma_i = 0, \\ v_1 = 0, & \end{cases} $$

$$ (39) \quad \underline{R}_1 : f \in \mathbb{R}^N \mapsto v = \underline{R}_1 f \in \mathbb{R}^N \\ \text{solving} \quad \begin{cases} v_i + \beta \sigma_i \underline{S} \left( \frac{h_b}{\alpha_b} (\underline{h_b S})^*(v) \right)_i = f_i, & \forall i \in [2, N] \quad \text{s.t. } \sigma_i > 0, \\ v_i = f_i, & \forall i \in [2, N] \quad \text{s.t. } \sigma_i = 0, \\ (\underline{h_b S})^*(v)_1 = 0; & \end{cases} $$

we recall that under the assumption (8), $h_{b,i} > 0$ when $\sigma_i > 0$, so that the division by $h_b$ in the above expressions is always well defined.
## 4.2. MacCormack second order scheme

In this section we propose to discretize our model (36) by a stable second order accurate MacCormack scheme adapted from the work of Beck, Lannes and Weynans [1] (2023). This scheme consists in a prediction-correction strategy that combines a Heun time integrator with an upwinded flux whose direction of upwinding alternates at every sub-step. More precisely, if we denote by $\mathcal{L}^-$, $\mathcal{L}^+$ the spatial discretizations respectively with a left and right upwinding of the nonlocal flux $f_b$ and centered source term, the overall procedure reads:

$$ \frac{U_i^{n+1,p} - U_i^n}{\Delta t} + \mathcal{L}^{-}(U)_i^n = 0, \qquad (\text{prediction step}) $$

$$ \frac{U_i^{n+1,c} - U_i^n}{\Delta t} + \mathcal{L}^{+}(U)_i^{n+1,p} = 0, \qquad (\text{correction step}) $$

$$ U_i^{n+1} = \frac{1}{2}(U_i^{n+1,p} + U_i^{n+1,c}). \qquad (\text{final averaging}) $$

In the following we review each one of these steps in detail. We shall make use of various strategies found in the literature to deal with critical issues such as: (i) avoiding a blow up of the velocity as the water height $h$ vanishes, (ii) keep the water height positive, (iii) make the MacCormack procedure stable in the NSW area through a TVD correction in the final averaging step, and (iv) preserve the lake-at-rest equilibrium with the pre-balanced formulation.

**Prediction.** For every interior cell $2 \le i \le N - 1$ the predicted elevation is of the following form

$$ \frac{\zeta_i^{n+1,p} - \zeta_i^n}{\Delta t} + \frac{1}{\Delta x} \left( \tau_{i+1/2} q_{i+1/2}^n - \tau_{i-1/2} q_{i-1/2}^n \right) = 0, $$

where $q_{i \pm 1/2}^n$ is the numerical flux at interfaces $x_i \pm \Delta x/2$, and where $\tau_{i \pm 1/2}$ is a flux limiter used to ensure the positivity of the water height. More specifically we use a left upwinded flux

$$ (40) \qquad q_{i+1/2}^n = h_i^n \mathcal{V}(h_i^n, q_i^n), \qquad \mathcal{V}(h, q) = \frac{\sqrt{2hq}}{\sqrt{h^4 + \max(h^4, \Delta x^4)}}, $$

with $\mathcal{V}$ the desingularized velocity associated to the state $(h, q)$ introduced by Kurganov and Petrova [10]. The advantage of computing the velocity in this way is that it remains robust even when the water height converges to zero, in the sense that there is no blow up of the velocity; in fact the above formula sets the velocity to zero in dry cells, is equivalent to $q/h$ in the case where $h \ge \Delta x$, and connects these two states in a smooth way as $h$ ranges from 0 to $\Delta x$.

The quantity $\tau_{i \pm 1/2}$ acts as a flux limiter by constraining the outflow of each cell so that the predicted water heights $(h_0 + \zeta_i^{n+1,p} - b_i)_i$ are positive. This limiting procedure is based on a draining time step first introduced by
Bollermann et al. [2] and later used by Liu [14]. For each interface $i + 1/2$ we define the limiter $\tau_{i+1/2}$ as:

$$ \tau_{i+1/2} = \begin{cases} \min\left(1, \frac{\delta t_i}{\Delta t}\right) & \text{if } q_{i+1/2}^n > 0, \\ \min\left(1, \frac{\delta t_{i+1}}{\Delta t}\right) & \text{otherwise,} \end{cases} \qquad \delta t_i = \frac{h_i^n \Delta x}{[q_{i-1/2}^n]_- + [q_{i+1/2}^n]_+}, $$

where $[\cdot]_-$ and $[\cdot]_+$ are respectively the negative and positive parts, and where $\delta t_i$ can be interpreted as the time required to reach a dry state in cell $i$ given the outflow $-([q_{i-1/2}^n]_- + [q_{i+1/2}^n]_+)$ and neglecting the inflow defined by $[q_{i+1/2}^n]_- + [q_{i-1/2}^n]_+$.

In order to compute the predicted discharge we first need to evaluate the approximation of $\dot{q}_{|x=0}(t^n)$. To this end we discretize Equation (35) as follows

$$ (\delta_x \mathfrak{s})_1 \delta_t q_0^n = -\frac{1}{\beta} \frac{g_0^{n-1} - 2g_0^n + g_0^{n+1}}{\Delta t^2} + \frac{3\alpha_{b,1}}{\beta h_{b,1}^2} f_{b,1}^n - (\delta_x \underline{R}_0 B^n)_1 \\ - \frac{3}{\beta h_{b,1}^2} \left( f_{\text{NSW}}(\zeta, q) - f_{\text{NSW}}(0,0) - (1-\beta) \frac{gh_b}{2} \zeta \delta_x(\sigma \delta_x b) \right)_1^n. $$

The predicted discharge is then obtained through a left upwinding:

$$ \frac{q_i^{n+1,p} - q_i^n}{\Delta t} + \frac{1}{\Delta x} (f_{b,i}^n - f_{b,i-1}^n) = (\underline{R}_0 B^n)_{i-1} - j_{b,i-}^n + \delta_t q_0^n \mathfrak{s}_{i-1}, $$

with

$$ \begin{aligned} f_{b,i}^n &= \begin{cases} \frac{(h_{b,i})^2}{\alpha_{b,i}} \left( \underline{R}_1 \left( \frac{\tilde{f}_\beta}{h_b^2} \right)_i \right)^n & \sigma_i > 0, \\ \tilde{f}_{\text{NSW},i}^n & \sigma_i = 0, \end{cases} \\ j_{b,i-}^n &= \begin{cases} \left( \frac{(\delta_x \alpha_b)_i}{\alpha_{b,i}} + \frac{3}{2} \frac{(\delta_x b)_i}{h_{b,i}} \right) f_{b,i-1}^n & \sigma_i > 0, \\ 0 & \sigma_i = 0, \end{cases} \\ B_i^n &= -g \zeta_i^n (\delta_x \tilde{b})_i + \begin{cases} \frac{3}{2} \frac{(\delta_x b)_i}{h_{b,i}} \tilde{f}_\beta(\zeta_i^n, q_i^n) & \sigma_i > 0, \\ 0 & \sigma_i = 0. \end{cases} \end{aligned} $$

The boundary cells are set according to:

$$ \begin{cases} \zeta_1^{n+1,p} = g_0(t^{n+1}), \\ \zeta_N^{n+1,p} = [h_{b,N}]_- \end{cases}, \qquad \begin{cases} q_1^{n+1,p} = q_1^n + \Delta t \delta_t q_0^n, \\ q_N^{n+1,p} = 0. \end{cases} $$

*Correction.* This step is similar to the first one but involves the predicted state $(h, q)^{n+1,p}$ previously obtained in the numerical flux

$$ \frac{\zeta_i^{n+1,c} - \zeta_i^n}{\Delta t} + \frac{1}{\Delta x} \left( \tau_{i+1/2} q_{i+1/2}^{n+1,p} - \tau_{i-1/2} q_{i-1/2}^{n+1,p} \right) = 0, $$
where

$$q_{i+1/2}^{n+1,p} = h_{i+1}^{n+1,p} V(h_{i+1}^{n+1,p}, q_{i+1}^{n+1,p}),$$

and where $\tau_{i\pm 1/2}$ is computed as before, the only difference being that $\delta t_i$ is now computed using the predicted outflow:

$$\delta t_i = \frac{h_i^n \Delta x}{[q_{i-1/2}^{n+1,p}]_- + [q_{i+1/2}^{n+1,p}]_+}.$$

We then solve

$$\begin{aligned} (\delta_x \mathfrak{s})_1 \delta_t q_0^{n+1,p} = & -\frac{1}{\beta} \frac{g_0^{n-1} - 2g_0^n + g_0^{n+1}}{\Delta t^2} + \frac{3\alpha_{b,1}}{\beta h_{b,1}^2} f_{b,1}^{n+1,p} - (\delta_x \underline{R}_0 B^{n+1,p})_1 \\ & - \frac{3}{\beta h_{b,1}^2} \left( f_{\text{NSW}}(\zeta, q) - f_{\text{NSW}}(0,0) - (1-\beta) \frac{gh_b}{2} \zeta \delta_x(\sigma \delta_x b) \right)_1^{n+1,p}. \end{aligned}$$

and the formula for the corrected discharge is also right upwinded:

$$\frac{q_i^{n+1,c} - q_i^n}{\Delta t} + \frac{1}{\Delta x} \left( f_{b,i}^{n+1,p} - f_{b,i-1}^{n+1,p} \right) = (\underline{R}_0 B)_{i+1}^{n+1,p} - j_{b,i+}^{n+1,p} + \delta_t q_0^{n+1,p} \mathfrak{s}_{i+1},$$

with $f_{b,i}^{n+1,p}$, $j_{b,i+}^{n+1,p}$ and $B_i^{n+1,p}$ defined as before but using the predicted values computed in the previous step. Finally the corrected boundary values are given by

$$\begin{cases} \zeta_1^{n+1,c} = g_0(t^{n+1}), \\ \zeta_N^{n+1,c} = [h_{b,N}]_- \end{cases}, \quad \begin{cases} q_1^{n+1,c} = q_1^n + \Delta t \delta_t q_0^{n+1,p}, \\ q_N^{n+1,c} = 0. \end{cases}$$

**Averaging with TVD correction.** In its standard version, the final step of the MacCormack scheme would consist to define the final update $U^{n+1}$ by averaging the predicted and corrected states $U^{n+1,p}$ and $U^{n+1,c}$ in every cell $1 \le i \le N$. Although this produces stable results in the dispersive area ($\sigma > 0$), it fails to do so in the shallow water part of the domain ($\sigma = 0$). The explanation is that brutal variations and even discontinuous shocks appear in the solution, which are in turn transformed into undesirable spurious oscillations by the scheme. To circumvent this, a stabilization technique has been introduced by Davis [5], and has been adapted for the shallow water system by Liang, Falconer and Lin in [13]. The method consists to add a correction term making the scheme TVD in the hyperbolic area of the domain. The averaging step reads:

$$\begin{aligned} U_i^{n+1} = \frac{1}{2} (U_i^{n+1,p} + U_i^{n+1,c}) &+ (G(r_i^+) + G(r_{i+1}^-)) \delta U_{i+1/2}^n \\ &- (G(r_{i-1}^+) + G(r_i^-)) \delta U_{i-1/2}^n, \end{aligned}$$
where $\delta U_{i\pm 1/2}^n \in \mathbb{R}^2$ denotes the discrete variation of $U$ at the interfaces

$$ \delta U_{i+1/2}^n = U_{i+1}^n - U_i^n, \quad \delta U_{i-1/2}^n = U_i^n - U_{i-1}^n, $$

and where $r_i^{\pm}$ are based on the scalar products between these vectors

$$ r_i^+ = \frac{(\delta U_{i-1/2}^n, \delta U_{i+1/2}^n)_{\mathbb{R}^2}}{(\delta U_{i+1/2}^n, \delta U_{i+1/2}^n)_{\mathbb{R}^2}}, \quad r_i^- = \frac{(\delta U_{i-1/2}^n, \delta U_{i+1/2}^n)_{\mathbb{R}^2}}{(\delta U_{i-1/2}^n, \delta U_{i-1/2}^n)_{\mathbb{R}^2}}. $$

The function $G$ is given by

$$ G(x) = \frac{C}{2}(1 - \varphi(x)), \quad \varphi(x) = \max(0, \min(2x, 1)), \\ C = \begin{cases} Cr(1 - Cr), & Cr \le 1/2 \\ 1/4, & Cr > 1/2 \end{cases} \quad Cr = \max_i (|u_i| + \sqrt{gh_i}) \frac{\Delta t}{\Delta x}. $$

### 4.3. Practical implementation

**Stability.** To ensure the stability of the above scheme a Courant-Friedrichs-Lewy condition should be satisfied, which in the hyperbolic case consists to restrict the time step according to $\Delta t \le \text{CFL} \times \frac{\Delta x}{v}$ where CFL is the Courant constant and $v$ the maximum speed at which the information is propagated in the solution. However the Boussinesq-Abbott is non hyperbolic, in fact it features nonlocal effects arising from the inversion of the elliptic operator, and as a consequence the maximum wave speed cannot in theory be deduced from the eigenvalues of the flux Jacobian. Nevertheless, it is a common practice to simply neglect the dispersive effects and to estimate $v$ from the hyperbolic NSW model, that is to say $v = \max_{i,n}|u_i^n| + \sqrt{gh_i^n}$. In actual use this strategy seem to yield stable results assuming the CFL constant is taken small enough; in our case we shall set it to 0.45.

**Coupling with NSW.** For performing numerical simulations presented in the next section, we shall use the following smooth transition function $\sigma$:

$$ (41) \qquad \begin{cases} \sigma(x) = \phi(\pm(x - x_c)/w), \\ \phi(x) = \mathbb{1}_{x \le -1/2} + \frac{1}{2}(1 + \tanh(-4x/(1 - 4x^2))) \cdot \mathbb{1}_{|x| < 1/2}, \end{cases} $$

with $x_c$ the center point and $w$ the width of the transition that will be set depending upon the test case, and with the sign $\pm$ depending on whether the open boundary is on the left (+) or the right (-) of the domain (note that the proposed scheme can be adapted to the second setting by symmetry). In particular, we should avoid positioning the transition too far away from the beach, in which case the shoaling might be underestimated; neither should we put it too close, since then oscillations might appear due to the growing steepness of the waves. As for the choice of the transition width, we know from Proposition C.1 that operator $R_1$ given in (29) is well-defined under the
condition $|\sigma'| \le \sqrt{3\alpha_b/(\beta h_b^2)}\sqrt{\sigma}$. This condition can be made true by setting $w$ large enough. More precisely, since the functions $\sigma, \phi$ defined in (41) satisfy

$$ \sup_{x \in \mathbb{R}} |\phi'(x)| = 2 \quad \text{and} \quad \sup_{x \in \mathbb{R}} |\sigma'(x)| = \frac{2}{w}, $$

a sufficient condition is to verify

$$ (42) \qquad w \ge 2 \cdot \sup_{|x-x_c|<w} \sqrt{\beta h_b^2/(3\alpha_b)}. $$

While inequality (42) ensures that $R_1$ is well-defined, it does not prevent the generation of reflections whenever a wave passes through from the Boussinesq domain to the NSW one. This behavior has been observed numerically in the past for various Boussinesq-type models and coupling strategies, see for instance the work of J. D. Galaz Mora [7]. A way to attenuate the reflections is to further increase the width of the transition; in this regard we suggest to set it at least to the peak wavelength characterizing the solution (which can be approximated from the peak frequency of the input signal $g$ and using the dispersion relation provided in (3)).

**Well-balancedness.** The above scheme is well-balanced in the sense that it preserves the rest state $(\zeta, q) = (0, 0)$. However, due to the presence of a sloped beach with a dry area in the NSW subdomain, the elevation at rest corresponds in reality to $\zeta = \max(0, b - h_0) = [h_b]_-$ which is not constant anymore. To recover a well-balanced scheme near the water line, a solution is to substitute the modified flux $\tilde{f}_\beta$ and source term $-g\zeta\partial_x\tilde{b}$ respectively by

$$ f_{\text{NSW}}(\zeta, q) - f_{\text{NSW}}([h_b]_-, 0) \quad \text{and} \quad -g(\zeta - [h_b]_-)\partial_x b $$

in every cell within the NSW area (i.e. outside of the support of $\sigma$). This change does not affect the scheme where dispersive effects are active ($\sigma_i > 0$).

# 5. Wave simulations

In this part, we aim to validate the coupled Boussinesq-NSW model (36) and its discretization provided in Section 4 against real-world experiments performed in a wave flume featuring a sloping beach. More specifically, the test cases considered are that of Synolakis [20] (solitary wave) in Section 5.1 and van Noorloos [24] (irregular wave field on a plane slope) in Section 5.2. The primary objective is to evaluate the model's ability to correctly handle boundary conditions; in this respect, the two test cases are particularly relevant as they involve both incoming and outgoing waves.

## 5.1. Solitary wave breaking over a sloped beach

This experiment was conducted by Synolakis in 1987, and consists in a solitary wave propagating and breaking on a beach with a 1:19.85 slope. The
non dimensional nonlinearity parameter is set to $A/h_0 = 0.28$ and $h_0 = 0.25$ m corresponds to the still water depth in the offshore flat bathymetry area of the flume. We consider a bounded computational domain $x_{dry} \le x \le x_{open}$ where $x$ increases seaward. The origin $x = 0$ corresponds to the water line at rest, the boundary condition is enforced offshore at $x_{open} = 4.775$ m, while the onshore boundary of the domain is located at $x_{dry} = -2.5$ m. The cell size is $\Delta x = 0.025$ m and the time step is chosen according to the CFL condition presented in Section 4.3. We also stress that in the literature, numerical simulations of this test case are performed by setting the solitary wave through the initial condition and not the open boundary. In our case, the initial condition is a lake-at-rest and the solitary wave is generated through the open boundary:

$$ (43) \qquad \zeta(t, x_{open}) = g(t), $$

where $g(t)$ is interpolated from the time series of the rightmost gauge.

In the definition (41) of the smooth transition function we set $\sigma(x) = \phi(-(x - x_c)/w)$ with $x_c = 7.5h_0$ and $w \approx 4.26h_0$. The comparison between the experimental and simulated elevation can be seen in Figure 2 at various dimensionless times $t^* = t\sqrt{g/h_0}$. In our data sample the first gauge used to enforce the boundary condition (43) was seemingly out of synchronization with the others, and we had to shift the associated timestamps by +0.1 s. After this change, the simulation remains in phase with the experimental data of the other gauges at later times. Until $t^* = 15$ included, the solitary wave propagates towards the shore and develops asymmetry due to the non flat bathymetry; the breaking hasn't occurred yet and we see a good agreement with the data. Around time $t^* = 20$, the wave begins to break and we see that the wavefront location is correctly predicted and remains so at times $t^* = 21$ and $t^* = 22$, however the shoaling is underestimated. This is a known drawback of Boussinesq-type models written in $(\zeta, q)$-coordinates (see [6]). At time $t^* = 23$, the simulated wave reaches the shore and starts collapsing with a good matching compared to the experiment. We then witness a run-up reaching its maximum at time $t^* = 45$ and receding afterwards; we see that our scheme produces good results overall, however with a run-down slower to that of the experiment ($t^* = 55$).

To conclude, these results attest the good prediction of the main transformations undergone by the wave (shoaling, breaking, and run-up) thanks to our strategy, whose fidelity is comparable to that of other Boussinesq-based methods in the literature. Though it doesn't make use of sponge layers nor the source function method, our model provides a robust treatment of the open boundary condition able to generate the incoming solitary wave properly. Evidence that outgoing waves are also properly handled is given in the next section.

## 5.2. Irregular waves on a sloping bottom

We now focus on the experiment performed by van Noorloos in 2003 involving irregular waves propagating over a sloped plane. The waves are generated by an offshore paddle ($x = 0$ m), such that their energy distribution in frequency reproduces a Jonswap spectrum. More precisely, we focus on experiment reference
<figure><img src="image_3.png" /><figcaption>Figure 2: Synolakis test case for $A/h_0 = 0.28$ and $h_0 = 0.25$ m at various dimensionless times $t^* = t\sqrt{g/h_0}$. Light blue stars represent the experiment gauges, while the dark blue line corresponds to the numerical approximation.</figcaption></figure>
C3 described in [24] which corresponds to a peak frequency of $f_p = 0.5$ Hz and a significant wave height of 0.1 m. Waves propagate onshore over a flat bottom in the first 8.5 meters of the flume, where the depth at rest is $h_0 = 0.7$ m; then they start evolving over a plane with a 1:35 slope. The toe of the bed and the shoreline at rest are respectively located at $x = 8.5$ m and $x = 33$ m. There is a total of 80 gauges recording the elevation throughout the channel with a sampling frequency of $f_s = 25$ Hz. The most offshore gauge is found at $x = 6$ m, however we noticed some inconsistencies in the associated data, so we preferred to use the second one at $x_{open} = 6.5$ m to generate the boundary data. The onshore boundary of the computational domain is located at $x_{dry} = 36.2$ m which always remains in the dry area. Hence the numerical approximation is computed on $[x_{open}, x_{dry}]$.

In the experiment, the waves break around $x_c = 25.1$ m, although individual events may occur slightly earlier or later. This variability helps explain why the simulation occasionally underestimates the elevation when the breaking occurs to the right of $x_c$, while it may overestimate the elevation when the breaking takes place to the left of $x_c$. We hypothesize that this issue could be mitigated by allowing the smooth transition function $\sigma$ to evolve dynamically in time and track several wave fronts simultaneously, guided by a wave-breaking criterion as in [21]. Such an approach should integrate seamlessly with our current coupling strategy between the Boussinesq-Abbott and NSW models, and is left for future work.

We define the smooth transition function by $\sigma(x) = \phi((x - x_c)/w)$ with $\phi$ given in (41) and $w = (x_{dry} - x_{open})/5$. The spatial discretization uses a refinement of $\Delta x = 0.025$ m and $\Delta t$ is set according to the same CFL condition as before. The duration of the simulation is $T_f = 1400$ s.

In what follows, we analyze the significant wave height, distinguishing between short waves ($H_{s,short}$) and infragravity waves ($H_{s,infra}$), defined respectively as

$$H_{s,short} = 4\sqrt{\int_{f_{cut}}^{f_{max}} E(f) \, df}, \quad H_{s,infra} = 4\sqrt{\int_{f_{min}}^{f_{cut}} E(f) \, df},$$

where $E(f)$ is the power spectral density (m$^2$/Hz), where $f_{min} = 10^{-2}$ Hz, $f_{cut} = 0.25$ Hz and $f_{max} = 1$ Hz. Note that to filter the short waves, we do not use Nyquist's frequency as the upper bound $f_{max}$, since Boussinesq-type models are not intended to capture well the very high frequencies.

Good agreement is observed in Figure 3a between the computed and experimental short-wave significant wave heights. However, $H_{s,short}$ is slightly underestimated in the first part of the surf zone. This may be partly explained by the simplified modeling of the surf zone that we have used. Figure 3b shows very good agreement between computed and experimental infragravity significant wave heights. This is an important result, as infragravity waves are weakly dissipated in the surf zone and are therefore reflected at the shoreline. Figure 3b thus highlights the capability of our open boundary condition to manage both incoming and reflected waves. Finally, the mean surface wave elevation is accurately reproduced by our Boussinesq model, as shown in Figure 3c.
<figure><img src="image_4.png" /><figcaption>Figure 3: Spatial evolution of computed (black line) and experimental data (light blue stars) for van Noorloos' experiment C3. (a) short wave significant wave height; (b) infragravity significant wave height; (c) mean surface elevation.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Energy spectrum $E(f)$ [m$^2$/Hz] at several locations. (a) $x = 7$ m; (b) $x = 20.5$ m; (c) $x = 25.7$ m; (d) $x = 30.2$ m. The black line represents the numerical approximation, while the light blue line corresponds to the experimental data. Note that the y-axis has been magnified in plot (d).</figcaption></figure>
Characteristic energy spectra at various locations in the wave flume are shown in Figure 4: (a) initial flat-bottom region; (b) shoaling zone; (c) just after the onset of wave breaking; and (d) inner surf zone. Overall, the spectra are very well reproduced by the model throughout the entire flume. A slight underestimation of energy is observed at frequencies higher than about twice the peak frequency. The shoreward growth of infragravity wave energy is accurately captured by our Boussinesq model, even within the surf zone. The model also successfully reproduces the spectral structure of the infragravity wave energy and the location of its extrema. This highlights the ability of our open boundary condition to manage both incident and reflected waves.

## 6. Conclusion

This study addresses the challenge of implementing open boundary conditions in a Boussinesq-type wave model coupled with the nonlinear shallow water (NSW) equations. This coupling framework is particularly well-suited to capture the full range of nearshore wave transformations — including shoaling, breaking, and run-up. Accurate treatment of boundary conditions is critical for realistic simulations of irregular wave fields. Traditional approaches such as sponge layers and source functions, though commonly used, are computationally expensive and difficult to calibrate. To overcome these limitations, we extended a recent nonlocal reformulation strategy to our model. This approach reformulates the Boussinesq-type model as a system of conservation laws with nonlocal flux and a rapidly decaying source term, enabling wave generation and absorption at open boundaries without auxiliary damping techniques. A key novelty lies in adapting this reformulation to the Boussinesq-Abbott model with an optimized dispersion parameter and its coupling to the NSW system over variable bathymetry.

We developed a numerical discretization based on a MacCormack-type prediction-correction scheme, involving finite-difference approximations of the nonlocal terms and a finite volume treatment with nonlocal fluxes for the time-stepping. The scheme preserves critical properties such as positivity of water depth and the lake-at-rest equilibrium, while ensuring stability through a TVD correction. The proposed method was validated through comparisons of numerical simulations with experimental data, including a solitary wave breaking on a sloped beach and the propagation of irregular waves. In these test cases, the model successfully generated incoming waves at the offshore boundary, captured key nearshore dynamics, and was able to properly evacuate reflected waves. Moreover, it accurately reproduced the spectral energy distribution across different frequencies and locations, demonstrating robustness despite the complexity of the wave field.

Future work could focus on enhancing the hybrid Boussinesq/NSW coupling by incorporating a criterion for identifying and tracking individual broken waves. This addition would likely improve modeling accuracy within the surf zone. Furthermore, unlike traditional source function methods, our approach does not rely on calibration based on a priori knowledge of the solution. This fundamental
advantage opens many perspectives, particularly in scenarios involving short and infragravity outgoing waves arising in natural coastal environments.

A. Derivation of the momentum equation with improved dispersion

We explain below how to obtain the modified flux and bathymetry terms $\tilde{f}_{\beta}$, $\tilde{b}$ appearing in (5) and defined in (9). The starting point is to consider the momentum equation with varying bathymetry, to which we subtract the identity $\partial_x f_{NSW}(0, 0) = -gh_b \partial_x b$ to get the so called pre-balanced formulation:

$$ (A.1) \qquad [1 + h_b \mathcal{T}_\sigma] \partial_t q + \partial_x (f_{NSW}(\zeta, q) - f_{NSW}(0, 0)) = -g\zeta \partial_x b. $$

Note that $f_{NSW}(0, 0) = \frac{1}{2}gh_b^2$ depends on $x$ through $h_b$. Next we consider the following approximation:

$$ (A.2) \qquad [1 + h_b \mathcal{T}_\sigma] \partial_t q \approx [1 + \beta h_b \mathcal{T}_\sigma] \partial_t q + (1 - \beta) h_b \mathcal{T}_\sigma(-gh_b \partial_x \zeta), $$

which comes from $\partial_t q = -gh_b \partial_x \zeta + O(A/h_0, h_0^2/\lambda^2)$ with $A$ and $\lambda$ the characteristic wave amplitude and wavelength. Using the definitions (6) and (7) we can develop $h_b \mathcal{T}_\sigma$ and check that

$$ h_b \mathcal{T}_\sigma(\cdot) = -\frac{1}{3} \partial_x \left( \sigma h_b^2 \partial_x(\cdot) + h_b(\partial_x b)(\cdot) \right) + \frac{h_b}{2} \partial_x (\sigma \partial_x b)(\cdot). $$

Furthermore, since there holds

$$ \left\{ \begin{aligned} & -\frac{1}{3} \partial_x \left( \sigma h_b^2 \partial_x(\cdot) + h_b(\partial_x b)(\cdot) \right) (gh_b \partial_x \zeta) = -\frac{1}{3} \partial_x (\sigma g h_b^3 \partial_{xx}^2 \zeta), \\ & \frac{h_b}{2} \partial_x (\sigma \partial_x b) (gh_b \partial_x \zeta) = \partial_x \left( \frac{g}{2} h_b^2 \zeta \partial_x (\sigma \partial_x b) \right) - g \zeta \partial_x \left( \frac{1}{2} h_b^2 \partial_x (\sigma \partial_x b) \right), \end{aligned} \right. $$

one can rewrite the approximation (A.2) as

$$ [1 + h_b \mathcal{T}_\sigma] \partial_t q \approx [1 + \beta h_b \mathcal{T}_\sigma] \partial_t q + \frac{1 - \beta}{3} \partial_x \left( \sigma g h_b^3 \partial_{xx}^2 \zeta - \frac{3}{2} g h_b^2 \zeta \partial_x (\sigma \partial_x b) \right) \\ + (1 - \beta) g \zeta \partial_x \left( \frac{1}{2} h_b^2 \partial_x (\sigma \partial_x b) \right). $$

Injecting this in (A.1) we finally get

$$ \begin{aligned} & [1 + h_b \mathcal{T}_\sigma] \partial_t q \\ & \qquad + \partial_x \left( \underbrace{f_{NSW}(\zeta, q) - f_{NSW}(0, 0) + (1 - \beta) \frac{gh_b^2}{3} \left( \sigma h_b \partial_{xx}^2 \zeta - \frac{3}{2} \zeta \partial_x (\sigma \partial_x b) \right)}_{\tilde{f}_\beta(\zeta, q)} \right) \\ & \approx -g\zeta \partial_x \left( \underbrace{b + \frac{1 - \beta}{2} h_b^2 \partial_x (\sigma \partial_x b)}_{\tilde{b}} \right), \end{aligned} $$

which is exactly the momentum equation from (5) after dropping the error terms arising from the above approximation.
## B. Nonlocal reformulation over a flat bottom

We show here how to reformulate the equation for the momentum derived above as a conservation law with nonlocal flux and a rapidly decaying source term. We consider here the case of a flat bottom.

### B.1. Reformulation with $R_0$ and dispersive boundary layer

The first step is to invert the dispersive operator acting on $\partial_t q$.

**Lemma B.1.** The pair $(\zeta, q)$ verifies the momentum equation of (1) if and only if it satisfies the equation (22) recalled below:

$$ \partial_t q = -R_0 \partial_x \tilde{f}_\beta(\zeta, q) + \dot{q}|_{x=0} \mathfrak{s}, $$

where $\tilde{f}_\beta$, $R_0$ and $\mathfrak{s}(x)$ and are defined in (9), (20) and (21) respectively.

*Proof.* Let $\phi \in L^2(\mathbb{R}_+)$. As a consequence of Lax-Milgram's theorem, the ODE

$$ (1 - \frac{\beta}{3} h_0^2 \partial_x (\sigma \partial_x \cdot)) Y = \phi, \quad Y(0) = \dot{q}|_{x=0} $$

admits a unique solution in the weighted Hilbert space

$$ (B.1) \qquad \begin{cases} V_{\sigma} = \{v \in L^2(\mathbb{R}_+), \|v\|_{V_{\sigma}} := \sqrt{(v, v)_{V_{\sigma}}} < \infty\}, \\ (u, v)_{V_{\sigma}} = (u, v)_{L^2(\mathbb{R}_+)} + \frac{\beta}{3} h_0^2 (\sigma \partial_x u, \partial_x v)_{L^2(\mathbb{R}_+)} \end{cases} $$

By linearity of the dispersive operator, this solution can be decomposed as $Y = Y_h + Y_b$ where

$$ \begin{cases} [1 + h_b T_\sigma] Y_h = 0, \\ Y_h(0) = \dot{q}|_{x=0}, \end{cases} \quad \text{and} \quad \begin{cases} [1 + h_b T_\sigma] Y_b = \phi, \\ Y_b(0) = 0, \end{cases} $$

and we remark that $Y_b = R_0\phi$ and $Y_h = \dot{q}|_{x=0}\mathfrak{s}$ by construction of $R_0$ and $\mathfrak{s}$. Setting $\phi = -\partial_x \tilde{f}_\beta(t, \cdot)$ implies $Y = \partial_t q(t, \cdot)$ for any arbitrary time $t$, and in that case the decomposition $Y = Y_h + Y_b$ is exactly (22). $\square$

### B.2. Operator $R_1$ is well defined over flat bottoms

We prove here that the operator $R_1$ is well defined on $L^2(\mathbb{R}_+)$, that is, that for every $f \in L^2(\mathbb{R}_+)$ there exists $v \in V_\sigma$ such that $(1 - \frac{\beta}{3}h_0^2\partial_x^2)v = f$ and $\partial_x v(0) = 0$.

**Proposition B.2.** Assume that there holds $|\sigma'| \le \sqrt{3/(\beta h_0^2)}\sqrt{\sigma}$. Then for any $f \in L^2(\mathbb{R}_+)$ the quantity $R_1 f$ is well defined in $V_\sigma$.

*Proof.* The variational formulation associated to the definition (23) of $R_1$ consists in finding $v \in V_\sigma$ defined in (B.1) such that for all $w \in V_\sigma$ there holds $a(v, w) = l(w)$, where

$$ a(v, w) = (v, w)_{L^2(\mathbb{R}_+)} + \frac{\beta}{3} h_0^2 \left[ (\sigma \partial_x v, \partial_x w)_{L^2(\mathbb{R}_+)} + (\sigma' \partial_x v, w)_{L^2(\mathbb{R}_+)} \right], $$
$$l(w) = (f, w)_{L^2(\mathbb{R}_+)}.$$

The linear form $l$ is continuous on $L^2(\mathbb{R}_+)$ and therefore on $V_\sigma$. The first two terms of the bilinear form $a$ are continuous on $V_\sigma$ by Cauchy-Schwarz inequality, while for the third component, we observe that if $|\sigma'| \le M\sqrt{\sigma}$ for some $M > 0$, then

$$ \begin{aligned} |(\sigma' \partial_x v, w)_{L^2(\mathbb{R}_+)}| & \le M|\sqrt{\sigma} \partial_x v|_{L^2(\mathbb{R}_+)}|w|_{L^2(\mathbb{R}_+)} \\ & \le M(h_0 \sqrt{\frac{\beta}{3}})^{-1}|v|_{V_\sigma}|w|_{V_\sigma}, \end{aligned} $$

so that we can conclude that $a$ is continuous on $V_\sigma$. Setting $w = v$ we also have:

$$(B.2) \quad a(v, v) = \|v\|_{L^2(\mathbb{R}_+)}^2 + \frac{\beta}{3} h_0^2 \|\sqrt{\sigma} \partial_x v\|_{L^2(\mathbb{R}_+)}^2 + \left( \frac{\beta}{3} h_0^2 \sigma' \partial_x v, v \right)_{L^2(\mathbb{R}_+)}$$

and the third term can be bounded from below using Young's inequality:

$$ \begin{aligned} \left(\frac{\beta}{3} h_0^2 \sigma' \partial_x v, v\right)_{L^2(\mathbb{R}_+)} & \ge -\frac{1}{2\epsilon} \|\frac{\beta}{3} h_0^2 \sigma' \partial_x v\|_{L^2(\mathbb{R}_+)}^2 - \frac{\epsilon}{2} \|v\|_{L^2(\mathbb{R}_+)}^2 \\ & \ge -\frac{1}{2\epsilon} \frac{\beta}{3} h_0^2 \left(M^2 \frac{\beta}{3} h_0^2\right) \|\sqrt{\sigma} \partial_x v\|_{L^2(\mathbb{R}_+)}^2 - \frac{\epsilon}{2} \|v\|_{L^2(\mathbb{R}_+)}^2 \end{aligned} $$

for any $\epsilon > 0$. Under the assumption $M^2 \frac{\beta}{3} h_0^2 \le 1$, we then get

$$ \left(\frac{\beta}{3} h_0^2 \sigma' \partial_x v, v\right)_{L^2(\mathbb{R}_+)} \ge -\frac{1}{2\epsilon} \frac{\beta}{3} h_0^2 \|\sqrt{\sigma} \partial_x v\|_{L^2(\mathbb{R}_+)}^2 - \frac{\epsilon}{2} \|v\|_{L^2(\mathbb{R}_+)}^2; $$

injecting this inequality in (B.2) yields:

$$ a(v, v) \ge \left(1 - \max\left(\frac{\epsilon}{2}, \frac{1}{2\epsilon}\right)\right) \|v\|_{V_\sigma}^2. $$

We obtain the coercivity of $a$ by taking $\frac{1}{2} < \epsilon < 2$, and operator $R_1$ is well defined by application of Lax-Milgram's theorem. $\square$

### C. Nonlocal reformulation with varying bottom

We extend here the results of the previous section to the case of varying bathymetry.

#### C.1. Operator $R_1$ is well defined over varying bottoms

The following propositions establishes a condition on $\sigma$ that ensures that the definition of $R_1$ makes sense.

**Proposition C.1.** Assume that there holds $|\sigma'| \le \sqrt{3\alpha_b/(\beta h_b^2)}\sqrt{\sigma}$. Then for any $f \in L^2(\mathbb{R}_+)$ the problem (29) admits a unique solution so that $R_1 f$ is well defined.
Proof. Using (7) we can show the identity

$$ \sigma h_b S(\cdot) = h_b S(\sigma \cdot) + \frac{\sigma'}{\sqrt{3}} h_b(\cdot), $$

therefore definition (29) of the operator $R_1$ is equivalent to solve the following problem for $u$:

$$ (C.1) \quad \left\{ \begin{array}{l} \left[ 1 + \frac{\beta}{h_b}(h_b S)\left(\frac{\sigma}{h_b}\alpha_b(h_b S)^*(\cdot)\right) + \beta \frac{h_b}{\alpha_b \sqrt{3}}(h_b S)^*(\cdot) \right] v = f, \\ (h_b S)^* v(0) = 0. \end{array} \right. $$

Let us introduce the Hilbert space $(V_\sigma, (\cdot, \cdot)|_{V_\sigma})$ defined as

$$ (C.2) \quad \left\{ \begin{array}{l} V_\sigma = \{v \in L^2(\mathbb{R}_+) : |v|_{V_\sigma} := \sqrt{(v, v)|_{V_\sigma}|} < \infty\}, \\ (u, v)|_{V_\sigma} = (h_b u, v)|_{L^2(\mathbb{R}_+)} + \beta \left( \frac{\sigma h_b}{\alpha_b} (h_b S)^* u, (h_b S)^* v \right)|_{L^2(\mathbb{R}_+)} \end{array} \right. $$

The variational formulation associated to (C.1) consists in finding $u \in V_\sigma$ such that for all $v \in V_\sigma$ there holds $a(u, v) = l(v)$, where the bilinear and linear forms $a, l$ are given by

$$ \begin{gathered} a(u, v) = (h_b u, v)|_{L^2} + \beta \left( \frac{\sigma h_b}{\alpha_b} (h_b S)^* u, (h_b S)^* v \right)|_{L^2} + \beta \left( \frac{h_b^2}{\alpha_b \sqrt{3}} \frac{\sigma'}{\sqrt{3}} (h_b S)^* u, v \right)|_{L^2}, \\ l(v) = (h_b f, v)|_{L^2}. \end{gathered} $$

As in the proof of Proposition B.2, we just need to show that $a$ is coercive and continuous in $V_\sigma \times V_\sigma$, and that $l$ is continuous in $V_\sigma$; the result then follows from the Lax-Milgram theorem. The proof being a straightforward adaptation of the proof of Proposition B.2, it is omitted. $\square$

### C.2. Reformulation of the interior equations

We recall that our coupled Boussinesq-NSW model (5) reads

$$ \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0, \\ [1 + \beta h_b \mathcal{T}_\sigma] \partial_t q + \partial_x \tilde{f}_\beta = -g\zeta \partial_x \tilde{b}, \end{array} \right. $$

with $\mathcal{T}_\sigma$ given in (6) and $\tilde{f}_\beta$ and $\tilde{b}$ defined in (9). We want to prove that the momentum equation from this system is equivalent to the form (33), which is given by:

$$ \partial_t q + \partial_x \tilde{f}_b = R_0 B(\zeta, q) - j_b + \dot{q}|_{\partial \Omega} \mathfrak{s}. $$

We start with the following proposition.

**Proposition C.2.** The pair $(\zeta, q)$ verifies the momentum equation of (5) if and only if it satisfies the equation below:

$$ (C.3) \quad \partial_t q + \sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* \left( R_1 \left( \frac{\tilde{f}_\beta}{h_b^2} \right) \right) = R_0 B(\zeta, q) + \dot{q}|_{\partial \Omega} \mathfrak{s}, $$

where $(h_b S)^*(\cdot) = S^*(h_b \cdot)$, $\alpha_b$ and $B(\zeta, q)$ are respectively defined in (7), (30) and (31), and where $\mathfrak{s}(x)$ is the unique solution of the problem (34).
The desired form (33) then results from the following Corollary.

**Corollary C.3.** Owing to the identity

$$ (h_b S)^*(\cdot) = \frac{\alpha_b}{\sqrt{3}h_b} \partial_x \left( \frac{h_b^2}{\alpha_b} \cdot \right) + \left( \frac{h_b}{\sqrt{3}} \frac{\partial_x \alpha_b}{\alpha_b} + \sqrt{3} \frac{\partial_x b}{2} \right) (\cdot), $$

we have the equivalent rewriting of (C.3) under conservative form:

$$ \partial_t q + \partial_x \mathfrak{f}_b = R_0 B(\zeta, q) - \mathfrak{j}_b + \dot{q}|_{\zeta=0} \mathfrak{s}, $$

with the nonlocal flux $\mathfrak{f}_b$ and nonlocal quantity $\mathfrak{j}_b$ defined in (32) as

$$ \mathfrak{f}_b = \frac{h_b^2}{\alpha_b} R_1 \left( \frac{\tilde{f}_\beta}{h_b^2} \right), \quad \mathfrak{j}_b = \left( \frac{\partial_x \alpha_b}{\alpha_b} + \frac{3}{2} \frac{\partial_x b}{h_b} \right) \mathfrak{f}_b. $$

It is this version of the discharge equation that we discretize in Section 4.

To prove the result from Proposition C.2 we follow the ideas of the proof of Proposition 4 from [11], and we need the following Lemma.

**Lemma C.4.** Let $\beta \ge 1$. The pair $(\zeta, q)$ verifies the momentum equation of (5) if and only if it satisfies the equation below:

$$ (C.4) \qquad \partial_t q = -R_0 \left( \partial_x \tilde{f}_\beta + g \zeta \partial_x \tilde{b} \right) + \dot{q}|_{\zeta=0} \mathfrak{s}, $$

where $\mathfrak{s}(x)$ is solution of the elliptic problem (34).

*Proof.* Since $\mathcal{T}_\sigma$ is self-adjoint, we get by Lax-Milgram's theorem that $1 + \beta h_b \mathcal{T}_\sigma$ is invertible, and the proof is similar to that of Lemma B.1, where we set $\phi = -\partial_x \tilde{f}_\beta - g \zeta \partial_x \tilde{b}$. $\square$

Next we show that the operator $R_0 \partial_x$ can be recasted under a conservative form up to a nonlocal operator with lower order by using $R_1$ defined in (29). More specifically, the following result holds.

**Proposition C.5.** Let $\beta \ge 1$; we have the identity:

$$ R_0 \partial_x f = \sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* \left( R_1 \left( \frac{f}{h_b^2} \right) \right) - \frac{3}{2} R_0 \left( \frac{\partial_x b}{h_b} f \right). $$

*Proof.* We begin by remarking that definition (7) of $S^*$ yields

$$ (C.5) \qquad \partial_x(\cdot) = h_b S^* \left( \frac{\sqrt{3}}{h_b} \cdot \right) - \frac{3}{2} \frac{\partial_x b}{h_b} (\cdot). $$

Hence we have

$$ R_0 \partial_x f = R_0 \left( h_b S^* \left( \frac{\sqrt{3}}{h_b} f \right) \right) - \frac{3}{2} R_0 \left( \frac{\partial_x b}{h_b} f \right), $$
so that it remains to show that the following identity holds:

$$ (C.6) \qquad R_0(h_b S^* \tilde{f}) = \frac{h_b}{\alpha_b} (h_b S)^* \left( R_1 \left( \frac{\tilde{f}}{h_b} \right) \right); $$

and the proof will be done by setting $\tilde{f} = \frac{\sqrt{3}}{h_b} f$. For that matter, let us write $v = R_1(\frac{\tilde{f}}{h_b})$. By definition (29) of $R_1$, one has

$$ \left\{ \begin{array}{l} \left[ 1 + \frac{\beta}{h_b} (\sigma h_b S) \left( \frac{h_b}{\alpha_b} (h_b S)^*(\cdot) \right) \right] v = \frac{1}{h_b} \tilde{f}, \\ (h_b S)^* v(0) = 0. \end{array} \right. $$

Applying $(h_b S)^*(\cdot) = S^*(h_b \cdot)$ to the equation, one finds that

$$ [(h_b S)^* + \beta S^*(\sigma h_b S) \left( \frac{h_b}{\alpha_b} (h_b S)^*(\cdot) \right)] v = S^*(\tilde{f}), $$

or equivalently multiplying by $h_b$

$$ [\alpha_b + \beta h_b S^*(\sigma h_b S)] \left( \frac{h_b}{\alpha_b} (h_b S)^*(v) \right) = h_b S^*(\tilde{f}); $$

using that $1 + \beta h_b T_\sigma = \alpha_b + \beta h_b S^*(\sigma h_b S(\cdot))$, we deduce that

$$ [1 + \beta h_b T_\sigma] \left( \frac{h_b}{\alpha_b} (h_b S)^*(v) \right) = h_b S^*(\tilde{f}). $$

Since moreover $(\frac{h_b}{\alpha_b}(\sigma h_b S)^*(v))$ vanishes at $x=0$, we use the definition (28) of $R_0$ to conclude that

$$ \frac{h_b}{\alpha_b} (\sigma h_b S)^*(v) = R_0(h_b S^*(\sigma \tilde{f})), $$

and recalling that $v = R_1(\frac{\tilde{f}}{h_b})$, we just proved (C.6) which concludes the proof. $\square$

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

We can now prove Proposition C.5

*Proof.* Using the result of Proposition C.5 with $f = \tilde{f}_\beta$ we get that (C.4) is equivalent to

$$ \begin{align*} \partial_t q &= -\sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* \left( R_1 \left( \frac{\tilde{f}_\beta}{h_b^2} \right) \right) + R_0 \left( \frac{3}{2} \frac{\partial_x b}{h_b} \tilde{f}_\beta - g \zeta \partial_x \tilde{b} \right) + \dot{q}_{|_0} \mathfrak{s} \\ &= -\sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* \left( R_1 \left( \frac{\tilde{f}_\beta}{h_b^2} \right) \right) + R_0 B(\zeta, q) + \dot{q}_{|_0} \mathfrak{s}. \end{align*} $$

<figure><img src="image_7.png" /></figure>
### C.3. Boundary equation on the missing data

In this section we give a justification of equation (35) on the missing boundary data $q|_{x=0}$. Deriving (C.3) with respect to $x$ and using that the mass conservation implies $\partial_{tx}^2 q = -\partial_{tt}^2 \zeta$ we have

$$ (C.7) \qquad -\partial_{tt}^2 \zeta + \sqrt{3}\partial_x \left[ \frac{h_b}{\alpha_b} (h_b S)^* \left( R_1 \left( \frac{\tilde{f}_\beta}{h_b^2} \right) \right) \right] = \partial_x R_0 B + \dot{q}|_{0} \mathfrak{s}'. $$

We will make use of the following Lemma.

**Lemma C.6.** The following identity holds:

$$ \sqrt{3} \left( \sigma \partial_x \left[ \frac{h_b}{\alpha_b} (h_b S)^* (R_1(f)) \right] \right)_{|_{x=0}} = - \frac{3}{\beta} [f - R_1 f]_{|_{x=0}}. $$

*Proof.* We remark that $\sqrt{3}\partial_x = -3S + \psi$ with the operator $S$ defined through (7) and where $\psi$ is a function whose definition does not matter. We can then develop:

$$ \sigma \sqrt{3} \partial_x \left[ \frac{h_b}{\alpha_b} (h_b S)^* (R_1 f) \right] = - \frac{3\sigma}{h_b} (h_b S) \left[ \frac{h_b}{\alpha_b} (h_b S)^* (R_1 f) \right] + \sigma \psi \frac{h_b}{\alpha_b} (h_b S)^* (R_1 f), $$

and by definition (29) of $R_1$ one gets the identity $\sigma(h_b S)[\frac{h_b}{\alpha_b}(h_b S)^*(R_1 f)] = h_b f - h_b R_1 f$ so that

$$ \sqrt{3} \left( \sigma \partial_x \left[ \frac{h_b}{\alpha_b} (h_b S)^* (R_1(f)) \right] \right) = - \frac{3}{\beta} (f - R_1 f) + \sigma \psi \frac{h_b}{\alpha_b} (h_b S)^* (R_1 f). $$

But by construction of $R_1$ the boundary condition $(h_b S)^*(R_1 f)|_{x=0} = 0$ is satisfied, and the result follows after taking the trace of the above equality at the boundary $x = 0$. $\square$

Taking the trace of equation (C.7) at $x = 0$ and applying Lemma C.6 with $\tilde{f} = \frac{1}{h_b^2} \tilde{f}_\beta$ we get

$$ (C.8) \qquad \dot{q}|_{0} \mathfrak{s}'(0) + (\partial_x R_0 B)|_{x=0} = -\ddot{g}_0 - \frac{3}{\beta h_b(0)^2} \tilde{f}_\beta(g_0, q_0) + \frac{3}{\beta} \left( R_1 \left( \frac{\tilde{f}_\beta}{h_b^2} \right) \right)_{|_{x=0}}. $$

It is important to remember that the modified flux $\tilde{f}_\beta$ contains the second order derivative $\partial_{xx}^2 \zeta$ which can be approximated by $(gh_b)^{-1} \partial_{tt}^2 \zeta$ similarly to what was done in Section 3.2. Hence if we introduce

$$ f_{\beta}^{\dagger} := f_{\text{NSW}}(\zeta, q) - f_{\text{NSW}}(0, 0) - (1 - \beta) \frac{gh_b^2}{2} \zeta \partial_x (\sigma \partial_x b), $$

one has

$$ \tilde{f}_{\beta}(\zeta, q) = f_{\beta}^{\dagger} - \frac{\beta - 1}{3} \sigma g h_b^3 \partial_{xx}^2 \zeta = f_{\beta}^{\dagger} - \frac{\beta - 1}{3} \sigma h_b^2 \partial_{tt}^2 \zeta + O((1 - \beta) \frac{A}{\lambda}). $$
We propose to replace (C.8) with

$$ \dot{q}|_{0} \mathbf{s}'(0) + (\partial_{x} R_{0} B)|_{x=0} = \left(\frac{\beta - 1}{\beta} \sigma(0) - 1\right) \ddot{g}_{0} - \frac{3}{\beta h_{b}(0)^{2}} f_{\beta}^{\dagger}(g_{0}, q_{0}) \\ + \frac{3}{\beta} \left(R_{1}\left(\frac{\tilde{f}_{\beta}}{h_{b}^{2}}\right)\right)|_{x=0}, $$

which is exactly (35) after using the assumption $\sigma(0) = 1$ and dividing by $\mathbf{s}'(0)$.

## References

[1] G. Beck, D. Lannes, and L. Weynans. “A numerical method for wavestructure interactions in the Boussinesq regime”. In: *arXiv preprint 2307.01749* (2023).

[2] A. Bollermann, S. Noelle, M. Lukáčová-Medvid'ová. “Finite volume evolution Galerkin methods for the shallow water equations with dry beds”. In: *Communications in Computational Physics* 10 (2011), pp. 371–404.

[3] P. Bonneton. “Modelling of periodic wave transformation in the inner surf zone”. In: *Ocean Engineering* 34.10 (2007), pp. 1459–1471. ISSN: 0029-8018.

[4] P. Bonneton et al. “A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model”. In: *Journal of Computational Physics* 230.4 (2011), pp. 1479–1498. ISSN: 0021-9991.

[5] S. F. Davis. “TVD finite difference schemes and artificial viscosity”. In: *ICASE Report* No. 84-20 (1984).

[6] A. G. Filippini et al. “On the nonlinear behaviour of Boussinesq type models: Amplitude-velocity vs amplitude-flux forms”. In: *Coastal Engineering* 99 (2015), pp. 109–123.

[7] J. D. Galaz Mora. “Coupling methods of phase-resolving coastal wave models”. Theses. Université de Montpellier, June 2024.

[8] M. Israeli and S. A. Orszag. “Approximation of radiation boundary conditions”. In: *Journal of Computational Physics* 41.1 (1981), pp. 115–135. ISSN: 0021-9991.

[9] M. Kazolea and A. I. Delis. “A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models”. In: *Applied Numerical Mathematics* 67 (2013). NUMAN 2010, pp. 167–186. ISSN: 0168-9274.

[10] A. Kurganov and G. Petrova. “A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system”. In: *Commun. Math. Sci.* Vol. 5, No. 1 (2007), pp. 133–160.
[11] D. Lannes and M. Rigal. "General boundary conditions for a Boussinesq model with varying bathymetry". In: *Studies in Applied Mathematics* 153.4 (2024), e12751.

[12] D. Lannes and L. Weynans. "Generating boundary conditions for a Boussinesq system". In: *Nonlinearity* 33.12 (Oct. 2020), pp. 6868–6889.

[13] D. Liang, R. A. Falconer, and B. Lin. "Comparison between TVD-MacCormack and ADI-type solvers of the shallow water equations". In: *Advances in Water Resources* 29.12 (2006), pp. 1833–1845. ISSN: 0309-1708.

[14] X. Liu. "A well-balanced and positivity-preserving numerical model for shallow water flows in channels with wet-dry fronts". *J. Sci. Comput.* 85, 60 (2020).

[15] P. A. Madsen, R. Murray, and O. R. Sørensen. "A new form of the Boussinesq equations with improved linear dispersion characteristics". In: *Coastal Engineering* 15.4 (1991), pp. 371–388. ISSN: 0378-3839.

[16] P. A. Madsen and H. A. Schäffer. "Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis". In: *Phil. Trans. R. Soc. A.* 356 (1998), pp. 3123–3181.

[17] M. Parisot. "Thick interface coupling technique for weakly dispersive models of waves". In: *ESAIM: M2AN* 58, 4 (2024), pp. 1497–1522.

[18] V. Roeber and K. F. Cheung. "Boussinesq-type model for energetic breaking waves in fringing reef environments". In: *Coastal Engineering* 70 (2012), pp. 1–20. ISSN: 0378-3839.

[19] F. Shi et al. "A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation". In: *Ocean Modelling* 43-44 (2012), pp. 36–51. ISSN: 1463-5003.

[20] C. E. Synolakis. "The runup of solitary waves". In: *Journal of Fluid Mechanics* 185 (1987), pp. 523–545.

[21] M. Tissier et al. "A new approach to handle wave breaking in fully non-linear Boussinesq models". In: *Coastal Engineering* 67 (2012), pp. 54–66. ISSN: 0378-3839.

[22] M. Tonelli and M. Petti. "Hybrid finite volume-finite difference scheme for 2DH improved Boussinesq equations". In: *Coastal Engineering* 56.5 (2009), pp. 609–620. ISSN: 0378-3839.

[23] M. Tonelli and M. Petti. "Shock-capturing Boussinesq model for irregular wave propagation". In: *Coastal Engineering* 61 (2012), pp. 8–19. ISSN: 0378-3839.
[24] J. C. Van Noorloos. "Energy transfer between short wave groups and bound long waves on a plane slope". MA thesis. Delft University of Technology, 2003.

[25] G. Wei, J. T. Kirby, and A. Sinha. "Generation of waves in Boussinesq models using a source function method". In: *Coastal Engineering* 36.4 (1999), pp. 271-299. ISSN: 0378-3839.