# General boundary conditions for a Boussinesq model with varying bathymetry

David Lannes, Mathieu Rigal

► To cite this version:

David Lannes, Mathieu Rigal. General boundary conditions for a Boussinesq model with varying bathymetry. Studies in Applied Mathematics, 2024. hal-04441332

HAL Id: hal-04441332

https://hal.science/hal-04441332v1

Submitted on 6 Feb 2024

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 - ShareAlike 4.0 International License
# General boundary conditions for a Boussinesq model with varying bathymetry

David Lannes¹ and Mathieu Rigal¹

¹Institut de Mathématiques de Bordeaux, Université de Bordeaux et CNRS UMR 5251, 351 Cours de la libération 33405 Talence Cedex, France

February 6, 2024

## Abstract

This paper is devoted to the theoretical and numerical investigation of the initial boundary value problem for a system of equations used for the description of waves in coastal areas, namely, the Boussinesq-Abbott system in the presence of topography. We propose a procedure that allows one to handle very general linear or nonlinear boundary conditions. It consists in reducing the problem to a system of conservation laws with nonlocal fluxes and coupled to an ODE. This reformulation is used to propose two hybrid finite volumes/finite differences schemes of first and second order respectively. The possibility to use many kinds of boundary conditions is used to investigate numerically the asymptotic stability of the boundary conditions, which is an issue of practical relevance in coastal oceanography since asymptotically stable boundary conditions would allow one to reconstruct a wave field from the knowledge of the boundary data only, even if the initial data is not known.

in memory of Vassilios A. Dougalis

This work was partially supported by the Institut des Mathématiques pour la Planète Terre and by the BOURGEONS project, grant ANR-23-CE40-0014-01 of the French National Research Agency (ANR).

*Mathematics Subject Classification.* 35B30, 35G61, 35Q35, 65M08, 76B15.

*Keywords.* Boussinesq-Abbott model, Initial boundary value problem, Inhomogeneous boundary conditions, Dispersive boundary layer, Nonlocal flux, Finite volumes.
# 1 Introduction

## 1.1 General motivation and objectives

The dynamics of water waves over a non flat bottom in the littoral area is a subject of great interest with far-reaching societal and economic stakes. Among them, one can mention safety aspects linked to the risk of flooding and sailing accidents, the impact of waves on the morphodynamics of the beach through sedimentation and erosion, or the issue of marine renewable energies, for which a good knowledge of the waves dynamics would be beneficial to increase the yield of offshore mooring systems or to choose pertinent wave farm locations. When considering coastal flows, it is not always possible to neglect dispersive effects which play an important role in the shoaling zone and in the generation of extreme waves. In order to accurately describe these phenomena, one therefore needs a model of shallow free surface flows that takes into account both nonlinear and dispersive effects. In this regard, the family of weakly nonlinear vertically-averaged Boussinesq-type systems is pertinent as it presents a good compromise between accuracy and reduced complexity compared for instance to the shallow water system, not accurate enough because it neglects dispersive effects, or to the fully nonlinear Serre-Green-Naghdi equations, more accurate but also more complex [25].

For practical applications, one usually considers a bounded spatial domain, and one has to solve an initial boundary value problem: the solution is computed at all times in terms of the initial data in the interior of the domain, as well as of the boundary data imposed at all times on its boundary. In the presence of dispersive terms, the question of boundary conditions is difficult and widely open, yet it is crucial for applications. The purpose of this work is to investigate this issue both theoretically and numerically in the case of the one dimensional Boussinesq-Abbott model with a non flat bottom, namely,

$$ \begin{cases} \partial_t \zeta + \partial_x q = 0 \\ (1 + h_b T_b) \partial_t q + \partial_x (\frac{1}{h} q^2) + gh \partial_x \zeta = 0 \end{cases} \quad t > 0, \quad x \in (0, \ell), $$

where $\zeta(t, x)$ denotes the elevation of the fluid at time $t$ and horizontal coordinate $x$ with respect to the rest state $z = 0$; $q(t, x)$ is the horizontal discharge of the fluid; $h(t, x) = h_0 + \zeta(t, x) - b(x)$ is the water height, with $h_0$ the characteristic depth and $z = -h_0 + b(x)$ a parametrization of the bottom topography; $g$ is the acceleration of gravity; and finally $T_b$ is a positive second order differential operator, whose exact expression is not important at this point. This model was proposed by Abbott [1] as a variant of the original Boussinesq-Peregrine model, which is written in terms of $\zeta$ and the average velocity $\bar{v} = q/h$; see [19] for a comparison of these two approaches.

More specifically, we propose a new method that enables to prescribe a nonlinear function of the unknowns $(\zeta, q)$ at the boundaries of the domain by a given function of time. The possibility to enforce such general non homogeneous boundary conditions is a prerequisite to be able to go beyond simple academic
test-cases such as solitary waves, and instead force realistic and complex wave fields in the domain, which typically feature swell and infra-gravity waves characterized by very disparate spatial scales and multi-directional propagation. This question is also related to that of asymptotic stability in the following sense: if we know the boundary data for all times (via buoys measurements for instance) but not the initial condition, can we expect the solution computed by solving the initial boundary value problem with an arbitrary initial data to converge in time to the real solution, and how does the choice of boundary conditions impact this convergence? This is a completely open theoretical problem, for which we propose here some numerical insight based on a new second order numerical scheme.

## 1.2 State of the art and contribution

Before describing the state of the art for the mathematical and numerical analysis of the initial boundary value problem associated with nonlinear dispersive models such as the Boussinesq-Abbott equations, let us recall some elements of the hyperbolic theory. Indeed, the Boussinesq-Abbott equations are a dispersive perturbation of the nonlinear shallow water equations that can be deduced from the Boussinesq-Abbott equations by neglecting the operator $\mathcal{T}_b$, namely,

$$ \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \partial_x (\frac{1}{h} q^2) + gh\partial_x \zeta = 0 \end{array} \right. \quad t > 0, \quad x \in (0, \ell). $$

This model belongs to the class of hyperbolic systems for which initial boundary value problems are well understood (see for instance [9] as well as [20] for a detailed analysis of the one dimensional case). As we shall see, despite the fact that the Boussinesq-Abbott system is a perturbation of the nonlinear shallow water equations, the behavior of the associated initial boundary value problems differ drastically. When applied to the nonlinear shallow water equations, the general theory of hyperbolic systems implies that under certain well identified conditions, it is possible to solve the equations on a finite domain if we know $\zeta$ and $q$ at $t = 0$, and if we impose at all times one scalar boundary data at $x = 0$ and $x = \ell$: typically, one can impose $\zeta$, $q$ or the incoming Riemann invariant (see in particular V. Dougalis' contribution for the numerical analysis of characteristic boundary conditions [2]). The missing boundary data is then recovered locally through the method of characteristics applied to the outgoing Riemann invariants. This is a fundamental difference compared to dispersive models, which do not admit Riemann invariants due to their non hyperbolic nature, and a consequence is that it is no longer possible to recover the missing boundary data using local arguments as pointed in [28].

V. Dougalis was a pioneer in the analysis of initial boundary value problems for Boussinesq-type equations with non trivial (i.e. non periodic) boundary conditions. In [4] he considered systems of the Bona-Smith family (for which a dispersive term is also present in the equation for the surface elevation) with non homogeneous boundary conditions on the surface elevation and on the velocity.
In the case of homogeneous boundary conditions, he then proposed a numerical Galerkin-finite element scheme [5]; he even addressed the two-dimensional case in [18, 17] where he proposed a mathematical and numerical analysis for various types of homogeneous boundary conditions, still for systems of the Bona-Smith family. In [3] he considered the "classical" Boussinesq system (corresponding to the Boussinesq-Peregrine system with a flat topography), with a homogeneous boundary condition on the velocity. He also considered the case of a variable topography [24] with homogeneous condition on the velocity or approximate transparent boundary conditions. The contribution of V. Dougalis and his coauthors is without doubt the most important to the theoretical and numerical analysis of the initial boundary value problem of Boussinesq-type equations. In the wake of these works, we propose here to investigate the case of non homogeneous boundary conditions.

Several authors also considered the related problem of transparent boundary conditions: for the linear KdV, BBM or Boussinesq equations, this was addressed in [11, 10, 30, 22] while for the nonlinear case approaches based on perfectly matched layer (PML) methods have been proposed for the KdV equation [12] and for a hyperbolic relaxation of the Serre-Green-Naghdi equations [21]. This approach can be adapted to handle the initial boundary value problem; instead of imposing a boundary condition, one imposes an initial data in the PML. This is also the idea behind the source function method [33]. These methods are robust and popular but have the drawback of being only approximations of the initial boundary value problems and to require to work on a considerably larger computational domain, at the cost of computational time. Besides let us mention [29] where a class of boundary conditions for time-discrete Serre-Green-Naghdi equations are considered.

It was shown recently that wave-structure interaction problems for the Boussinesq equations could be reduced to initial boundary value problems with non homogeneous boundary conditions. A new method was proposed to handle this problem theoretically [15, 7] and numerically [8]; it was also generalized in [28] to handle generating boundary conditions (the boundary condition is on the surface elevation $ζ$). All these references considered the case of the Boussinesq-Abbott model with a flat topography, and with boundary conditions imposed on $ζ$ or $q$. For applications to coastal oceanography it is however important to consider a non-flat topography, as well as other kinds of boundary conditions. The goal of this article is to generalize the above method to propose a solution to these two problems, and to design a general second order well-balanced scheme to approximate the solutions.

## 1.3 Organization of the paper

We first propose in Section 2 a theoretical analysis of the initial boundary value problem associated with the Boussinesq-Abbott when the boundary conditions are imposed on the surface elevation at $x = 0$ and $x = l$. In the case of a flat bottom, considered in §2.1 this is a generalization of [28] to the case of a finite interval; in the case of a non-flat topography studied in §2.2, the analysis is
significantly more involved because some crucial commutation properties used
in [28] are no longer true. Well-posedness is proved using a reformulation of
the original initial boundary value problem as a simple initial value problem
(no boundary condition); a well-balanced version of this formulation is proposed
in §2.3, on which our numerical scheme is based on.

We then consider the case of general boundary conditions in Section 3. We
work in a general framework where a (possibly nonlinear) function of ζ and
q is imposed at each endpoint of the domain; these functions are called *input*
functions. There exist also in this framework so called *output* functions and a
reconstruction mapping that allows to recover ζ and q from the knowledge of the
input and output functions, see §3.1. We show in §3.2 that the output functions
can be recovered through the resolution of a first or second order ODE; the
number of initial conditions that must be provided to solve this ODE depends
of course on its order; in addition compatibility conditions must also be imposed
on the data. The analysis of these two types of conditions is performed in §3.3.
The well-posedness of the original initial boundary value problem with general
boundary conditions can then be established in §3.4. We end this section by
commenting in §3.5 on an open theoretical problem of high practical interest: the
asymptotical stability of the boundary conditions. In other terms, the question
is to know whether the solution of the initial boundary value problem can be
asymptotically recovered after a transitory regime if we start from a different
initial data but impose the same boundary data.

We then present in Section 4 the numerical schemes we use in our numerical
simulations. We use a hybrid finite volumes/finite differences approach. We
propose both a first order Lax-Friedrichs scheme in §4.2 and a second order
MacCormack scheme in §4.3.

Numerical simulations are finally presented in Section 5. We first show in §5.1
the ability of our numerical schemes to solve non trivial initial boundary value
problems in the presence of topography. Several kinds of boundary conditions
are investigated. We finally investigate in §5.2 the issue of asymptotic stability.
It is in particular shown that numerical asymptotic stability is achieved by
imposing nonlinear boundary conditions (the Riemann invariants associated
with the nonlinear shallow water equations), but not if we impose the surface
elevation or the discharge at the boundaries. This shows the relevance of being
able to enforce generic nonlinear boundary conditions when dealing with complex
applications of initial boundary value problems.

**Acknowledgement.** The authors warmly thank Philippe Bonneton for many fruitful discussions and for his insights on the physical motivations of this work.

## 2 The Boussinesq-Abbott system with boundary conditions on the surface elevation

We propose here a theoretical analysis of the initial boundary value problem for
the Boussinesq-Abbott system with boundary conditions on the surface elevation.
We first treat in §2.1 the case of flat bathymetries; this problem was considered in [28] on the half-line; we extend it here to treat the case of a finite interval, with boundary conditions on the surface elevation at both ends. We then consider in §2.2 the Boussinesq-Abbott equations with topography. As in the case of flat topography, it is possible to reduce the problem to an initial boundary value problem, which we use to prove well-posedness, but the analysis is considerably more difficult because the operator $T_b$ does not commute with space derivatives. In order to prepare the ground for our numerical scheme, a well balanced version of this reformulation is proposed in §2.3.

## 2.1 The case of flat bathymetries

We consider here the case of a flat bathymetry; the Boussinesq-Abbott model then reduces to

$$ (1) \qquad \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0 \\ (1 - \frac{h_0^2}{3} \partial_x^2) \partial_t q + \partial_x f_{\text{NSW}} = 0 \quad \text{in } (0, \ell), \end{array} \right. $$

where we recall that

$$ f_{\text{NSW}} = f_{\text{NSW}}(\zeta, q) = \frac{1}{2}g(h^2 - h_0^2) + \frac{1}{h}q^2, $$

where $h = h_0 + \zeta$. The initial boundary value problem for this system was studied in [28] on the half-line $(0, \infty)$, with a prescribed boundary value on the surface elevation $\zeta$ at $x = 0$ (see also [8] for the case of a prescribed boundary value on $q$). We extend this result here to the case of a finite interval $(0, \ell)$ with prescribed boundary values on $\zeta$ at $x = 0$ and $x = \ell$.

### 2.1.1 Notations

We first need to introduce $R_0$, the inverse of $(1 - \frac{h_0^2}{3} \partial_x^2)$ with homogeneous Dirichlet boundary data at $x = 0$ and $x = \ell$; more precisely, we set

$$ (2) \qquad R^0 : f \in L^2(0, \ell) \mapsto u \in H^2(0, \ell) \quad \text{solving} \quad \left\{ \begin{array}{l} [1 - \frac{h_0^2}{3} \partial_x^2] u = f \\ u(0) = u(\ell) = 0 \end{array} \right. $$

(note that $R^0$ is also well defined on $H^{-1}(0, \ell)$ with values in $H^1(0, \ell)$). Since, contrary to [28], we work on a finite interval, we also need to introduce $\mathfrak{s}^{(0)}$ and $\mathfrak{s}^{(\ell)}$ the solutions of the equations of the homogeneous equation but with non homogeneous boundary conditions at the left and right boundaries respectively,

$$ (3) \qquad \left\{ \begin{array}{l} [1 - \frac{h_0^2}{3} \partial_x^2] \mathfrak{s}^{(0)} = 0 \\ \mathfrak{s}^{(0)}(0) = 1, \quad \mathfrak{s}^{(0)}(\ell) = 0 \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{l} [1 - \frac{h_0^2}{3} \partial_x^2] \mathfrak{s}^{(\ell)} = 0 \\ \mathfrak{s}^{(\ell)}(0) = 0, \quad \mathfrak{s}^{(\ell)}(\ell) = 1 \end{array} \right. . $$

We also denote by $\mathcal{S}'$ the matrix

$$ (4) \qquad \mathcal{S}' = \begin{pmatrix} (\mathfrak{s}^{(0)})'(0) & (\mathfrak{s}^{(\ell)})'(0) \\ (\mathfrak{s}^{(0)})'(\ell) & (\mathfrak{s}^{(\ell)})'(\ell) \end{pmatrix}. $$
**Remark 2.1.** One can compute $\mathfrak{s}^{(0)}$ and $\mathfrak{s}^{(\ell)}$ explicitly, namely,

$$ (5) \qquad \mathfrak{s}^{(0)}(x) = \frac{\sinh(\sqrt{3}\frac{(\ell-x)}{h_0})}{\sinh(\sqrt{3}\frac{\ell}{h_0})}, \qquad \mathfrak{s}^{(\ell)}(x) = \frac{\sinh(\sqrt{3}\frac{x}{h_0})}{\sinh(\sqrt{3}\frac{\ell}{h_0})}, $$

from which an explicit expression for $\mathfrak{s}'$ easily follows.

We finally introduce the inverse of $(1 - \frac{h_0^2}{3}\partial_x^2)$ with homogeneous Neumann boundary data at $x = 0$ and $x = \ell$;

$$ (6) \quad R^1 : f \in L^2(0, \ell) \mapsto u \in H^2(0, \ell) \quad \text{solving} \quad \left\{ \begin{array}{l} [1 - \frac{h_0^2}{3}\partial_x^2]u = f \\ \partial_x u(0) = \partial_x u(\ell) = 0. \end{array} \right. $$

An important observation is the following commutation property,

$$ (7) \qquad \forall f \in L^2(0, \ell), \qquad R_0 \partial_x f = \partial_x R_1 f. $$

### 2.1.2 Well-posedness of the initial boundary value problem

We prove here that the initial boundary value problem formed by (1) with boundary conditions

$$ (8) \qquad \zeta(\cdot, 0) = g_0 \quad \text{and} \quad \zeta(\cdot, \ell) = g_\ell, $$

and initial condition

$$ (9) \qquad (\zeta, q)(0, \cdot) = (\zeta^{\text{in}}, q^{\text{in}}) $$

is well posed. Remark that from the first equation of (15) one has $\frac{d}{dt}(\zeta(t, 0)) = -\partial_x q(t, 0)$, and proceeding similarly at $x = \ell$, a necessary compatibility condition on the initial and boundary data to allow the possibility of a solution which is of class $C^1$ at the origin is that

$$ (10) \qquad \left\{ \begin{array}{l} \zeta^{\text{in}}(0) = g_0(0), \\ -\partial_x q^{\text{in}}(0) = \dot{g}_0(0), \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{l} \zeta^{\text{in}}(\ell) = g_\ell(0), \\ -\partial_x q^{\text{in}}(\ell) = \dot{g}_\ell(0). \end{array} \right. $$

The key step, in the spirit of [28, 8] is to reformulate the problem as an initial value problem (no boundary condition).

**Proposition 2.2.** Assume that the initial and boundary data $(\zeta^{\text{in}}, q^{\text{in}})$ and $(g_0, g_\ell)$ satisfy the compatibility condition (10). Then the two following assertions are equivalent:

i. The couple $(\zeta, q)$ is a regular solution to (1) such that the depth $h$ never vanishes and with boundary conditions (8) and initial condition (9).

ii. The couple $(\zeta, q)$ is a regular solution such that the depth $h$ never vanishes to

$$ (11) \qquad \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \partial_x (R_1 f_{\text{NSW}}) = \dot{q}_0 \mathfrak{s}^{(0)} + \dot{q}_\ell \mathfrak{s}^{(\ell)}, \end{array} \right. \qquad \text{in } (0, \ell), $$
with initial condition (9), and where $q_0$ and $q_\ell$ solve the ODE

$$ (12) \qquad \mathcal{S}'\left(\frac{\dot{q}_0}{\dot{q}_\ell}\right) + \frac{3}{h_0^2}\left(\frac{f_{\text{NSW}}(g_0, q_0)}{f_{\text{NSW}}(g_\ell, q_\ell)}\right) = -\left(\ddot{g}_0\right) + \frac{3}{h_0^2}\left(\frac{(R_1 f_{\text{NSW}})_{|x=0}}{(R_1 f_{\text{NSW}})_{|x=\ell}}\right), $$

with $\mathcal{S}'$ defined in (4), and with initial condition

$$ (13) \qquad q_0(0) = q^{\text{in}}(0) \quad \text{and} \quad q_\ell(0) = q^{\text{in}}(\ell). $$

*Proof.* Let $(\zeta, q)$ be a regular solution to (1) with boundary conditions $\zeta(\cdot, 0) = g_0$ and $\zeta(\cdot, \ell) = g_\ell$ and initial condition $(\zeta, q)(0, \cdot) = (\zeta^{\text{in}}, q^{\text{in}})$. By definition of $R_0$, the second equation of (1) can be rewritten as

$$ \partial_t q + R_0 \partial_x f_{\text{NSW}} = \dot{q}_0 \mathbf{s}^{(0)} + \dot{q}_\ell \mathbf{s}^{(\ell)}, $$

where $q_0$ and $q_\ell$ denote the trace of $q$ at $x=0$ and $x=\ell$ respectively. Differentiating with respect to $x$ we get

$$ \partial_t \partial_x q + \partial_x R_0 \partial_x f_{\text{NSW}} = \dot{q}_0 \mathbf{s}^{(0)} + \dot{q}_\ell \mathbf{s}^{(\ell)}. $$

Using the first equation one can replace $\partial_t \partial_x q = -\partial_t^2 \zeta$ while the identities $R_0 \partial_x = \partial_x R_1$ and $(1 - \frac{h_0^2}{3} \partial_x^2) R_1 = \text{Id}$ allow us to deduce that

$$ -\partial_t^2 \zeta - \frac{3}{h_0^2}(1 - R_1)f_{\text{NSW}} = \dot{q}_0 \mathbf{s}'(0) + \dot{q}_\ell \mathbf{s}'(\ell). $$

Taking the trace of this equation at $x=0$ and $x=\ell$ respectively, we therefore obtain the following differential system on $q_0$ and $q_\ell$,

$$ \begin{cases} \mathbf{s}'(0)(0)\dot{q}_0 + \mathbf{s}'(\ell)(0)\dot{q}_\ell = -\frac{3}{h_0^2}f_{\text{NSW}}(g_0, q_0) + \frac{3}{h_0^2}(R_1 f_{\text{NSW}})_{|x=0} - \ddot{g}_0 \\ \mathbf{s}'(0)(\ell)\dot{q}_0 + \mathbf{s}'(\ell)(\ell)\dot{q}_\ell = -\frac{3}{h_0^2}f_{\text{NSW}}(g_\ell, q_\ell) + \frac{3}{h_0^2}(R_1 f_{\text{NSW}})_{|x=\ell} - \ddot{g}_\ell; \end{cases} $$

this shows that the first point of the proposition implies the second one. For the reverse implication, we can observe that if $(\zeta, q)$ and $(q_0, q_\ell)$ solve the above system then $(\zeta, q)$ solves the Boussinesq system (1). By taking the trace of the second equation of (11) at $x=0$ and $x=\ell$ respectively, we also obtain that $\frac{d}{dt}q(t, 0) = \dot{q}_0(t)$ and $\frac{d}{dt}q(t, \ell) = \dot{q}_\ell(t)$. Since moreover we take by assumption $q(0, 0) = q_0(0)$ and $q(0, \ell) = q_\ell(0)$, we deduce that $q(t, 0) = q_0(t)$ and $q(t, \ell) = q_\ell(t)$ for all times. In contrast, the boundary conditions (8) are not necessarily satisfied. However, as for the derivation of (12), we obtain that $(\zeta_0(t), \zeta_\ell(t)) := (\zeta(t, 0), \zeta(t, \ell))$ satisfies

$$ \mathcal{S}'\left(\frac{\dot{q}_0}{\dot{q}_\ell}\right) + \frac{3}{h_0^2}\left(\frac{f_{\text{NSW}}(\zeta_0, q_0)}{f_{\text{NSW}}(\zeta_\ell, q_\ell)}\right) = -\left(\frac{\ddot{\zeta}_0}{\ddot{\zeta}_\ell}\right) + \frac{3}{h_0^2}\left(\frac{(R_1 f_{\text{NSW}})_{|x=0}}{(R_1 f_{\text{NSW}})_{|x=\ell}}\right), $$

and it follows that $(\zeta_0, \zeta_\ell)$ and $(g_0, g_\ell)$ satisfy the same second order ODE. We can therefore conclude that $(\zeta_0, \zeta_\ell) = (g_0, g_\ell)$ if and only if

$$ (\zeta_0(0), \zeta_\ell(0)) = (g_0(0), g_\ell(0)) \quad \text{and} \quad (\dot{\zeta}_0(0), \dot{\zeta}_\ell(0)) = (\dot{g}_0(0), \dot{g}_\ell(0)). $$

Replacing $\partial_t \zeta$ by $-\partial_x q$, these conditions are equivalent to (10). It follows that the boundary conditions (8) are also satisfied. This concludes the proof of the proposition. $\square$
The initial value problem formed by (11) and (12) with initial conditions (9) and (13) is actually an ordinary differential equation on $H^n(0, \ell) \times H^{n+1}(0, \ell) \times \mathbb{R}^2$ for $(\zeta, q, q_0, q_\ell)$ if $n$ is a nonnegative integer. The existence of a local in time solution therefore follows immediately from Cauchy-Lipschitz theorem (we do not give details for the proof, which can easily be adapted from [28]). Note that the condition $\inf(h_0 + \zeta^{\text{in}}) > 0$ ensures that the shallow water flux $f_{\text{NSW}}$ is well defined. By Proposition 2.2, this solution furnishes a solution to the original initial boundary value problem (1), (8) and (9).

**Proposition 2.3.** Let $g_0, g_\ell \in C^\infty(\mathbb{R}^+)$. Let also $n \in \mathbb{N} \setminus \{0\}$ and $(\zeta^{\text{in}}, q^{\text{in}}) \in H^n(0, \ell) \times H^{n+1}(0, \ell)$ be such that $\inf(h_0 + \zeta^{\text{in}}) > 0$ and satisfying the compatibility conditions (10).

Then there exists a maximal existence time $T^* > 0$ and a unique solution $(\zeta, q, q_0, q_\ell) \in C^\infty([0, T^*); H^n(0, \ell) \times H^{n+1}(0, \ell) \times \mathbb{R}^2)$ to (11) and (12) with initial conditions (9) and (13), and for all $t \in [0, T^*)$, one has $q_0(t) = q(t, 0)$ and $q_\ell(t) = q(t, \ell)$.

## 2.2 The Boussinesq-Abbott model with topography

We now consider the full Boussinesq-Abbott model which contains additional topography terms compared to (1). Throughout this section, we assume that the bottom is parameterized by a function $-h_0 + b$, with $b$ a smooth function on $[0, \ell]$. We assume also that the depth of the fluid at rest never vanishes, namely,

$$ (14) \qquad \inf_{(0,\ell)} (h_0 - b) > 0. $$

Denoting $h_b = h_0 - b$, the Boussinesq-Abbott model is given by

$$ (15) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0 \\ (1 + h_b T_b) \partial_t q + \partial_x f_{\text{NSW}} = -gh \partial_x b \end{cases} \quad \text{in } (0, \ell), $$

where $h$ now contains an additional topography term, $h = h_0 + \zeta - b$, and where the shallow water flux is still given by

$$ f_{\text{NSW}} = f_{\text{NSW}}(\zeta, q) = \frac{1}{2}g(h^2 - h_0^2) + \frac{1}{h}q^2, $$

which is the same expression as in the case of a flat topography, except that $h$ now depends also on $b$. The second order operator $T_b$ is given by

$$ (16) \qquad T_b(\cdot) = -\frac{1}{3h_b} \partial_x \left( h_b^3 \partial_x \left( \frac{\cdot}{h_b} \right) \right) + \frac{1}{2} \partial_x^2 b $$

and can be alternatively written under the form (see Section 5.6 of [26])

$$ T_b(\cdot) = S^*(h_b S(\cdot)) + \frac{1}{4} \frac{(\partial_x b)^2}{h_b} $$

with

$$ (17) \qquad S(\cdot) = -\frac{1}{\sqrt{3}} h_b \partial_x \left( \frac{1}{h_b} \cdot \right) + \frac{\sqrt{3}}{2} \frac{\partial_x b}{h_b}. $$
### 2.2.1 Notations and preliminary results

We generalize here, in the presence of a non-flat topography, the concepts introduced in Section 2.1.1. We introduce the operator $R_b^0$ as the inverse of $(1 + h_b T_b)$ with homogeneous Dirichlet boundary conditions, that is to say

$$ (18) \qquad R_b^0 : f \in L^2(0, \ell) \mapsto u \in H^2(0, \ell) \quad \text{solving} \quad \begin{cases} [1 + h_b T_b] u = f \\ u(0) = u(\ell) = 0 \end{cases}; $$

here again, the operator $R_b^0$ can be also extended as an operator mapping $H^{-1}(0, \ell)$ to $H^1(0, \ell)$. We also introduce the solutions $\mathfrak{s}^{(b,0)}$ and $\mathfrak{s}^{(b,\ell)}$ of the homogeneous equation with non homogeneous boundary conditions, (19)

$$ \left\{ \begin{array}{l} [1 + h_b T_b] \mathfrak{s}^{(b,0)} = 0 \\ \mathfrak{s}^{(b,0)}(0) = 1, \quad \mathfrak{s}^{(b,0)}(\ell) = 0 \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{l} [1 + h_b T_b] \mathfrak{s}^{(b,\ell)} = 0 \\ \mathfrak{s}^{(b,\ell)}(0) = 0, \quad \mathfrak{s}^{(b,\ell)}(\ell) = 1 \end{array} \right. $$

and we also denote by $\mathcal{G}'_b$ the matrix

$$ (20) \qquad \mathcal{G}'_b = \begin{pmatrix} (\mathfrak{s}^{(b,0)})'(0) & (\mathfrak{s}^{(b,\ell)})'(0) \\ (\mathfrak{s}^{(b,0)})'(\ell) & (\mathfrak{s}^{(b,\ell)})'(\ell) \end{pmatrix}. $$

Note that contrary to the case of a flat topography, there is in general no explicit expression for $\mathfrak{s}^{(b,0)}$ and $\mathfrak{s}^{(b,\ell)}$. We can however prove that $\mathcal{G}'_b$ is invertible under an assumption which is satisfied by all the topography profiles we considered in this paper, as well as by typical beach profiles.

**Proposition 2.4.** Assume that (14) is satisfied and that in addition

$$ (21) \qquad \min_{(0,\ell)} \left(1 + \frac{1}{4}(\partial_x b)^2 + \frac{1}{6}h_b \partial_x^2 b\right) > 0. $$

Then the matrix $\mathcal{G}'_b$ is invertible and $(\mathfrak{s}^{(b,0)})'(0) \neq 0$ and $(\mathfrak{s}^{(b,\ell)})'(\ell) \neq 0$.

*Proof.* We can remark that

$$ 1 + h_b T_b = -\frac{1}{3} \partial_x (h_b^2 \partial_x \cdot) - \frac{1}{3} h_b (\partial_x b) \partial_x + \left(1 + \frac{1}{3} (\partial_x b)^2 + \frac{1}{6} h_b \partial_x^2 b\right). $$

For all $u \in H^1(0, \ell)$ such that $\partial_x u(0) = \partial_x u(\ell) = 0$, we obtain therefore that

$$ \int_0^\ell u(1+h_b T_b)u = \frac{1}{3} \int_0^\ell (h_b^2(\partial_x u)^2 - h_b(\partial_x b)u\partial_x u) + \int_0^\ell \left(1+\frac{1}{3}(\partial_x b)^2 + \frac{1}{6}h_b\partial_x^2 b\right)u^2. $$

Since for all $\epsilon > 0$, one has

$$ -h_b(\partial_x b)u\partial_x u \geq -\frac{\epsilon}{2}h_b^2(\partial_x u)^2 - \frac{1}{2\epsilon}(\partial_x b)^2 u^2, $$

we deduce that

$$ \int_0^\ell u(1+h_b T_b)u \geq \frac{1}{3}\left(1-\frac{\epsilon}{2}\right) \int_0^\ell h_b^2(\partial_x u)^2 + \int_0^\ell \left(1+\frac{1}{3}\left(1-\frac{1}{2\epsilon}\right)(\partial_x b)^2 + \frac{1}{6}h_b\partial_x^2 b\right)u^2. $$
Under the assumptions of the proposition, it is possible to find $\frac{1}{2} < \epsilon < 2$ and a constant $C_{\epsilon} > 0$ such that

$$ \int_{0}^{\ell} u(1 + h_b \mathcal{T}_b)u \geq C_{\epsilon}|u|_{H^1(0, \ell)}^2; $$

in particular, if $u$ solves $(1 + h_b \mathcal{T}_b)u = 0$ and satisfies $\partial_x u(0) = \partial_x u(\ell) = 0$, then one has $u \equiv 0$.

Now, if the matrix $\mathcal{S}'_b$ were not invertible, there would exist $(\lambda, \mu) \neq (0, 0)$ such that the function $\mathfrak{s} := \lambda(\mathfrak{s}^{(b,0)}) + \mu(\mathfrak{s}^{(b,\ell)})$ would satisfy $\mathfrak{s}'(0) = \mathfrak{s}'(\ell) = 0$. From the above considerations, one would have $\mathfrak{s} \equiv 0$, which is not possible since $\mathfrak{s}(0) = \lambda$ and $\mathfrak{s}(\ell) = \mu$ and $(\lambda, \mu) \neq (0, 0)$.

In order to prove that $(\mathfrak{s}^{(b,0)})'(0) \neq 0$, we can modify the arguments above to get that if $u \in H^1(0, \ell)$ is such that $u'(0) = 0$ and $u(\ell) = 0$ then $u \equiv 0$. If we had $(\mathfrak{s}^{(b,0)})'(0) = 0$ then $(\mathfrak{s}^{(b,0)})$ would be such a function and would identically vanish, which is absurd. We get similarly that $(\mathfrak{s}^{(b,\ell)})'(\ell) \neq 0$, which concludes the proof. $\square$

Another important difference with the case of flat topography is that the commutation property (7) is no longer true, that is, in general we have $R_b^0 \partial_x \neq R_b^1 \partial_x$ if we define $R_b^1$ as the inverse of $(1 + h_b \mathcal{T}_b)$ with homogeneous Neumann boundary conditions. We can however define $R_b^1$ as follows

$$ R_b^1 : f \in L^2(0, \ell) \mapsto u \in H^2(0, \ell) $$

where $u$ solves

$$ (22) \qquad \left\{ \begin{array}{l} [1 + \frac{1}{h_b}(h_b S)(\frac{h_b}{\alpha_b}(h_b S)^*)]u = f, \\ (h_b S)^*u(0) = (h_b S)^*u(\ell) = 0, \end{array} \right. $$

and with $\alpha_b = 1 + \frac{1}{4}(\partial_x b)^2$. The following proposition shows that there is a commutation property that generalizes the identity $R^0 \partial_x = \partial_x R^1$ for the Boussinesq-Abbott model with topography.

**Proposition 2.5.** For all $f \in H^1(0, \ell)$, the following identity holds

$$ R_b^0 \partial_x f = \sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1 (\frac{1}{h_b^2} f)) - \frac{3}{2} R_b^0 (\frac{\partial_x b}{h_b} f), $$

where $\alpha_b = 1 + \frac{1}{4}(\partial_x b)^2$, and where we recall that $S(\cdot)$ is defined in (17).

**Remark 2.6.** When $b \equiv 0$ then $h_b = h_0$, $R_b^0 = R^0$, $R_b^1 = R^1$ and $S^* = \frac{1}{\sqrt{3}}\partial_x$, so that the identity of the lemma coincides with the identity $R^0 \partial_x = \partial_x R^1$ obtained previously.

*Proof.* Let us first remark that one can write

$$ (23) \qquad 1 + h_b \mathcal{T}_b = \alpha_b + h_b S^* (h_b S(\cdot)); $$
it follows that the equation $(1+h_b\mathcal{T}_b)u = \partial_x f$ can be equivalently written under the form

$$ (24) \qquad \begin{aligned} & \left[ \frac{\alpha_b}{h_b} + S^*(h_b S(\cdot)) \right] u = \frac{1}{h_b} \partial_x f \\ & = S^*\left(\frac{\sqrt{3}}{h_b} f\right) - \frac{3}{2} \frac{1}{h_b^2} (\partial_x b) f. \end{aligned} $$

We now need the following lemma.

**Lemma 2.7.** For all $\tilde{f} \in L^2(0, \ell)$, the following identity holds,

$$ R_b^0(h_b S^* \tilde{f}) = \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1(\frac{\tilde{f}}{h_b})). $$

*Proof of the lemma.* Let us write $v = R_b^1(\frac{\tilde{f}}{h_b})$; by definition, one has

$$ \left\{ \begin{array}{l} \left[ 1 + \frac{1}{h_b} (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) = (h_b S)^* v(\ell) = 0; \end{array} \right. $$

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

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

or equivalently

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

using (23), we deduce that

$$ [1 + h_b \mathcal{T}_b] \left( \frac{h_b}{\alpha_b} (h_b S)^*(v) \right) = h_b S^* \tilde{f} $$

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

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

which proves the lemma. $\square$

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

Owing to (23), one can write (24) under the form

$$ R_b^0 \partial_x f = R_b^0 (h_b S^* (\frac{\sqrt{3}}{h_b} f)) - \frac{3}{2} R_b^0 (\frac{\partial_x b}{h_b} f), $$

so that the result follows by applying the lemma with $\tilde{f} = \frac{\sqrt{3}}{h_b} f$. $\square$

<figure><img src="image_3.png" /></figure>
### 2.2.2 Well-posedness of the initial boundary value problem

We prove here that the initial boundary value problem formed by (15) with boundary conditions

$$ (25) \qquad \zeta(\cdot, 0) = g_0 \quad \text{and} \quad \zeta(\cdot, \ell) = g_\ell, $$

and initial condition

$$ (26) \qquad (\zeta, q)(0, \cdot) = (\zeta^{\text{in}}, q^{\text{in}}) $$

is well posed. As seen previously, a necessary compatibility condition on the initial and boundary data to allow the possibility of a solution which is of class $C^1$ at the origin is that

$$ (27) \qquad \begin{cases} \zeta^{\text{in}}(0) = g_0(0), \\ -\partial_x q^{\text{in}}(0) = \dot{g}_0(0), \end{cases} \quad \text{and} \quad \begin{cases} \zeta^{\text{in}}(\ell) = g_\ell(0), \\ -\partial_x q^{\text{in}}(\ell) = \dot{g}_\ell(0). \end{cases} $$

We first state the following generalization of Proposition 2.2 to the Boussinesq-Abbott model, in which we use the notation

$$ (28) \qquad B(\zeta, q) = -gh\partial_x b + \frac{3}{2} \frac{\partial_x b}{h_b} f_{\text{NSW}}. $$

We also introduce the two dimensional vectors

$$ (29) \qquad V_{\text{bdry}}(\zeta_0, q_0, \zeta_\ell, q_\ell) = \left( \frac{3}{h_b(0)^2} f_{\text{NSW}}(\zeta_0, q_0) \quad \frac{3}{h_b(\ell)^2} f_{\text{NSW}}(\zeta_\ell, q_\ell) \right) $$

which depends only on the boundary values of $(\zeta, q)$, and

$$ (30) \qquad V_{\text{int}}[\zeta, q] = \left( \begin{array}{l} 3(R_b^1(\frac{1}{h_b^2} f_{\text{NSW}}))|_{x=0} - (\partial_x(R_b^0 B))|_{x=0} \\ 3(R_b^1(\frac{1}{h_b^2} f_{\text{NSW}}))|_{x=\ell} - (\partial_x(R_b^0 B))|_{x=\ell} \end{array} \right), $$

which is a nonlocal function of the interior values of $\zeta$ and $q$ on the whole interval $(0, \ell)$.

**Proposition 2.8.** Let the bottom parametrization $b$ satisfy (14) and (21). Assume that the initial and boundary data $(\zeta, q^{\text{in}})$ and $(g_0, g_\ell)$ satisfy the compatibility condition (27). Then the two following assertions are equivalent:

i. The couple $(\zeta, q)$ is a regular solution to (15) such that the depth $h$ never vanishes and with boundary conditions (25) and initial condition (26).

ii. The couple $(\zeta, q)$ is a regular solution such that the depth $h$ never vanishes to (31)

$$ \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1(\frac{1}{h_b^2} f_{\text{NSW}})) = R_b^0 B + \dot{q}_0 \mathfrak{s}^{(b,0)} + \dot{q}_\ell \mathfrak{s}^{(b,\ell)}, \quad \text{in } (0, \ell), \end{array} \right. $$
where we recall that $(h_b S)^* = \frac{1}{\sqrt{3}} \frac{1}{h_b} \partial_x (h_b^2) + \frac{\sqrt{3}}{2} \partial_x b$ and that $B(\zeta, q)$ is given by (28), and with initial condition (26), while $q_0$ and $q_\ell$ solve the ODE

$$ (32) \qquad \mathcal{S}'_b \begin{pmatrix} \dot{q}_0 \\ \dot{q}_\ell \end{pmatrix} + V_{\text{bdry}}(g_0, q_0, g_\ell, q_\ell) = V_{\text{int}}[\zeta, q] - \begin{pmatrix} \ddot{g}_0 \\ \ddot{g}_\ell \end{pmatrix}, $$

with $\mathcal{S}'_b$ defined in (20), and with initial condition

$$ (33) \qquad q_0(0) = q^{\text{in}}(0) \quad \text{and} \quad q_\ell(0) = q^{\text{in}}(\ell). $$

**Remark 2.9.** When $b = 0$ (flat topography), one can check that $f_b = R^1 f_{\text{NSW}}$ and $g_b = 0$, so that (31) coincides as expected with (11).

*Proof.* We proceed as in the proof of Proposition 2.2. Applying $R_b^0$ to the equation in $q$, and using Proposition 2.5, we first obtain

$$ \begin{aligned} & \partial_t q + \sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1 (\frac{1}{h_b^2} f_{\text{NSW}})) - \frac{3}{2} R_b^0 (\frac{\partial_x b}{h_b} f_{\text{NSW}}) \\ & = -g R_b^0 (h \partial_x b) + \dot{q}_0 \mathfrak{s}^{(b,0)} + \dot{q}_\ell \mathfrak{s}^{(b,\ell)}. \end{aligned} $$

Regrouping the terms involving $R_b^0$, one deduces the second equation of (31). Applying $\partial_x$ to this equation and using the equation on $\zeta$, we then get

$$ (34) \qquad \begin{aligned} & -\partial_t^2 \zeta + \sqrt{3} \partial_x \left[ \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1 (\frac{1}{h_b^2} f_{\text{NSW}})) \right] - \frac{3}{2} \partial_x \left[ R_b^0 (\frac{\partial_x b}{h_b} f_{\text{NSW}}) \right] \\ & = -g \partial_x R_b^0 (h \partial_x b) + \dot{q}_0 (\mathfrak{s}^{(b,0)})' + \dot{q}_\ell (\mathfrak{s}^{(b,\ell)})'. \end{aligned} $$

We now need the following lemma.

**Lemma 2.10.** For all $\tilde{f} \in L^2(0, \ell)$, one has

$$ \sqrt{3} \partial_x \left[ \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1 \tilde{f}) \right] \Big|_{x=0,\ell} = -3 \left[ (\tilde{f} - R_b^1 \tilde{f}) \right] \Big|_{x=0,\ell}. $$

*Proof.* One can notice that $\sqrt{3}\partial_x = -3S + \beta$, where $\beta$ is a function whose exact expression is of no importance here. We then have

$$ \sqrt{3} \partial_x \left[ \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1 \tilde{f}) \right] = - \frac{3}{h_b} (h_b S) \left[ \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1 \tilde{f}) \right] + \beta \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1 \tilde{f}); $$

by definition of $R_b^1$, one has $(h_b S) \left[ \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1 \tilde{f}) \right] = h_b \tilde{f} - h_b R_b^1 \tilde{f}$ so that

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

Since by construction $(h_b S)^* (R_b^1 \tilde{f})$ vanishes at $x = 0$ and $x = \ell$, the result follows upon taking the trace of the above identity at the two boundary points. $\square$
Taking the trace of (34) at $x = 0$ and using the lemma with $\tilde{f} = \frac{1}{h_b^2} f_{NSW}$, we get that $q_0$ and $q_l$ satisfy the differential system

$$ egin{cases} 
\dot{q}_0 \mathfrak{s}'_{(b,0)}(0) + \dot{q}_\ell \mathfrak{s}'_{(b,\ell)}(0) - g(\partial_x R_b^0(h\partial_x b))_0 = \\ 
-\ddot{g}_0 - \frac{3}{h_b(0)^2} f_{NSW}(g_0, q_0) + 3(R_b^1(\frac{1}{h_b^2} f_{NSW}))_0 - \frac{3}{2} [\partial_x R_b^0(\frac{\partial_x b}{h_b} f_{NSW})]_0, \\ 
\dot{q}_\ell \mathfrak{s}'_{(b,0)}(\ell) + \dot{q}_\ell \mathfrak{s}'_{(b,\ell)}(\ell) - g(\partial_x R_b^0(h\partial_x b))_\ell = \\ 
-\ddot{g}_\ell - \frac{3}{h_b(\ell)^2} f_{NSW}(g_\ell, q_\ell) + 3(R_b^1(\frac{1}{h_b^2} f_{NSW}))_\ell - \frac{3}{2} [\partial_x R_b^0(\frac{\partial_x b}{h_b} f_{NSW})]_\ell, 
\end{cases} $$

which corresponds to (32). The end of the proof is as in Proposition 2.2. □

The well-posedness of the initial value problem given in the second point of the proposition is simply obtained as for Proposition 2.3 from Cauchy-Lipschitz theorem.

**Proposition 2.11.** Assume that $b$ is a smooth function satisfying (14) and (21). Let $g_0, g_\ell \in C^\infty(\mathbb{R}^+)$. Let also $n \in \mathbb{N}\setminus\{0\}$ and $(\zeta^{\text{in}}, q^{\text{in}}) \in H^n(0, \ell) \times H^{n+1}(0, \ell)$ be such that $\inf_{[0,\ell]}(h_0 - b + \zeta^{\text{in}}) > 0$ and assume that the compatibility conditions (27) hold. Then there exists a maximal existence time $T^* > 0$ and a unique solution $(\zeta, q, q_0, q_\ell) \in C^\infty([0, T^*]; H^n(0, \ell) \times H^{n+1}(0, \ell) \times \mathbb{R}^2)$ to (31) and (32) with initial conditions (26) and (33), and for all $t \in [0, T^*)$, one has $q_0(t) = q(t, 0)$ and $q_\ell(t) = q(t, \ell)$, and $\inf_{[0,\ell]}(h_0 - b + \zeta(t, \cdot)) > 0$.

*Proof.* The system (31)-(32) can be rewritten as

$$ \frac{d}{dt} \mathbb{U} = \mathbb{F}[t, \mathbb{U}] $$

with $\mathbb{U} = (\zeta, q, q_0, q_\ell)^\text{T}$ and

$$ \mathbb{F}[t, \mathbb{U}] = (-\partial_x q, -\partial_x(f_b[\zeta, q]) + g_b[\zeta, q] + \left(\mathfrak{s}^{(b,0)}_{\mathfrak{s}^{(b,\ell)}}\right) \cdot \mathfrak{h}_b[t, \mathbb{U}], \mathfrak{h}_b[t, \mathbb{U}]^\text{T})^\text{T}, $$

where

$$ \mathfrak{h}_b[t, \mathbb{U}] = (\mathfrak{S}'_b)^{-1} \left[ -V_{\text{bdry}}(g_0(t), q_0, g_\ell(t), q_\ell) + V_{\text{int}}[\zeta, q] - \begin{pmatrix} \ddot{g}_0(t) \\ \ddot{g}_\ell(t) \end{pmatrix} \right]. $$

Let $\mathbb{X} = H^n(0, \ell) \times H^{n+1}(0, \ell) \times \mathbb{R}^2$ and $\Omega \subset \mathbb{X}$ the open subset of $\mathbb{X}$ consisting of all $\mathbb{U} = (\zeta, q, q_0, q_\ell)^\text{T} \in \mathbb{X}$ such that $\inf_{(0,\ell)}(h_0 + \zeta - b) > 0$. Then it follows from standard product estimates, and because $R_b^0$ and $R_b^1$ map $H^n(0, \ell)$ to $H^{n+2}(0, \ell)$, that $\mathbb{F}: \mathbb{R}^+ \times \Omega \to \mathbb{X}$ is well-defined, continuous and locally-Lipschitz with respect to $\mathbb{U}$. By Cauchy-Lipschitz theorem, for all initial data $\mathbb{U}^{\text{in}} \in \Omega$ there exists therefore a maximal solution $\mathbb{U} \in C^1([0, T^*], \Omega)$, with $T^* > 0$, and moreover this solution belongs to $\mathbb{U} \in C^\infty([0, T^*], \Omega)$. The result then follows from Proposition 2.8. □
## 2.3 A reformulation adapted to well-balancedness

Since the steady state $(\zeta, q) = (0,0)$ (lake at rest) solves the system (31)-(32) with boundary data $g_0 = g_\ell = 0$, one gets in particular the following identity from the second equation of (31)

$$ \sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1(\frac{1}{h_b^2} f_{\text{NSW}}(0,0))) = R_b^0 B(0,0). $$

In the presence of topography, both terms of this identity are nonzero. Subtracting this identity to the second equation of (31), one obtains the following equivalent formulation

$$ (35) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1(\frac{1}{h_b^2} \tilde{f}_{\text{NSW}})) = R_b^0 \tilde{B} + \dot{q}_0 \mathfrak{s}^{(b,0)} + \dot{q}_\ell \mathfrak{s}^{(b,\ell)}, \end{cases} \quad \text{in } (0, \ell), $$

where, denoting $\tilde{f}_{\text{NSW}}(\zeta, q) := f_{\text{NSW}}(\zeta, q) - f_{\text{NSW}}(0,0)$, one has

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

Proceeding similarly with (32) we get that $q_0$ and $q_\ell$ solve the ODE

$$ (37) \qquad \mathfrak{S}'_b \left( \frac{\dot{q}_0}{\dot{q}_\ell} \right) + \tilde{V}_{\text{bdry}}(g_0, q_0, g_\ell, q_\ell) = \tilde{V}_{\text{int}}[\zeta, q] - \left( \begin{matrix} \ddot{g}_0 \\ \ddot{g}_\ell \end{matrix} \right), $$

with

$$ (38) \qquad \tilde{V}_{\text{bdry}}(\zeta_0, q_0, \zeta_\ell, q_\ell) = \left( \frac{\frac{3}{h_b(0)^2} \tilde{f}_{\text{NSW}}(\zeta_0, q_0)}{\frac{3}{h_b(\ell)^2} \tilde{f}_{\text{NSW}}(\zeta_\ell, q_\ell)} \right) $$

and

$$ (39) \qquad \tilde{V}_{\text{int}}[\zeta, q] = \left( \begin{matrix} 3(R_b^1(\frac{1}{h_b^2} \tilde{f}_{\text{NSW}}))|_{x=0} & - (\partial_x(R_b^0 \tilde{B}))|_{x=0} \\ 3(R_b^1(\frac{1}{h_b^2} \tilde{f}_{\text{NSW}}))|_{x=\ell} & - (\partial_x(R_b^0 \tilde{B}))|_{x=\ell} \end{matrix} \right), $$

In order to obtain well-balanced numerical schemes, we will use (35) and (37) rather than (31) and (32).

## 3 More general boundary conditions

In the previous sections, we considered the Boussinesq-Abbott system (15) on a finite interval $(0, \ell)$ and with imposed boundary conditions at $x = 0$ and $x = \ell$ on the surface elevation $\zeta$. In some situations, one may have to impose rather a boundary condition on $q$, or an a function of $\zeta$ and $q$ such as a Riemann invariant of the nonlinear shallow water system.
Denoting $(\zeta_0, q_0) = (\zeta, q)|_{x=0}$ and $(\zeta_\ell, q_\ell) = (\zeta, q)|_{x=\ell}$, we study in this section the possibility of imposing more general boundary conditions which take the form

$$ (40) \qquad \xi_0^+(\zeta_0, q_0) = g_0, \qquad \xi_\ell^-(\zeta_\ell, q_\ell) = g_\ell, $$

where the data $(g_0, g_\ell)$ are given functions of time, and where $\xi_0^+(\zeta, q)$ and $\xi_\ell^-(\zeta, q)$ are functions of $\zeta$ and $q$; we refer to $\xi_0^+(\zeta, q)$ and $\xi_\ell^-(\zeta, q)$ as the input functions.

In order to solve the initial boundary value problem associated with (40), an essential point is to compute the traces of $\zeta$ and $q$ at the boundaries $x = 0$ and $x = \ell$. This is done in two steps. We first need two output functions that can be found in terms of the input functions; we then reconstruct the traces of $\zeta$ and $q$ in terms of the input and output functions. These two steps are explained in §3.1. We then derive in §3.2 an ODE that can be used to compute the output functions. In §3.3, we discuss the issue of providing initial data to this ODE, as well as the issue of compatibility conditions between initial and boundary data. We can then state and prove in §3.4 our main well-posedness result. We finally introduce in §3.5 the issue of asymptotic stability which is an open theoretical problem that will be numerically investigated later in Section 5.

## 3.1 The output functions and the reconstruction mappings

We assume that there exist two other functions $\xi_0^-(\cdot, \cdot)$ and $\xi_\ell^+(\cdot, \cdot)$, referred to as the output functions, and such that $(\zeta_0, q_0)$ can be recovered from the knowledge of $g_0$ and $\xi_0^-(\zeta_0, q_0)$ and similarly, $(\zeta_\ell, q_\ell)$ can be recovered from the knowledge of $g_\ell$ and $\xi_0^+(\zeta_\ell, q_\ell)$. More precisely, we make the following assumption.

**Assumption 3.1.** There exists a non-empty connected open set $U \subset \mathbb{R}^2$ verifying $(0, 0) \in U$ such that the following conditions hold:

i. The input functions $\xi_0^+ : U \to \mathbb{R}$ and $\xi_\ell^- : U \to \mathbb{R}$ are well defined and smooth;

ii. There are two output functions $\xi_0^- : U \to \mathbb{R}$ and $\xi_\ell^+ : U \to \mathbb{R}$ and two open sets $V_0$ and $V_\ell$ that contain the range of the mappings $\xi_0 : U \ni (\zeta, q) \mapsto (\xi_0^+(\zeta, q), \xi_0^-(\zeta, q)) \in \mathbb{R}^2$ and $\xi_\ell : U \ni (\zeta, q) \mapsto (\xi_\ell^+(\zeta, q), \xi_\ell^-(\zeta, q)) \in \mathbb{R}^2$ respectively, as well as two smooth mappings $H_0 : V_0 \to \mathbb{R}^2$ and $H_\ell : V_\ell \to \mathbb{R}^2$ such that for all $(\zeta, q) \in U$, one has

$$ (41) \qquad H_0(\xi_0^+(\zeta, q), \xi_0^-(\zeta, q)) = (\zeta, q)^T \quad \text{and} \quad H_\ell(\xi_\ell^+(\zeta, q), \xi_\ell^-(\zeta, q)) = (\zeta, q)^T; $$

we refer to $H_0$ and $H_\ell$ as the reconstruction mappings.

**Example 3.2.** For instance, let us impose the boundary conditions $\zeta_0 = g_0$ and $R^-(\zeta_\ell, q_\ell) = g_\ell$, where $R^-(\zeta, q) = \frac{q}{h_0+\zeta-b} - 2\sqrt{g(h_0+\zeta-b)}$ is the left-going Riemann invariant associated with the nonlinear shallow water equations. With the above notations, this amounts to take

$$ \xi_0^+(\zeta, q) = \zeta \quad \text{and} \quad \xi_\ell^-(\zeta, q) = R^-(\zeta, q). $$
We can then choose

$$ \xi_0^-(\zeta, q) = q \quad \text{and} \quad \xi_\ell^+(\zeta, q) = R^+(\zeta, q), $$

where $R^+ = \frac{q}{h_0+\zeta-b} + 2\sqrt{g(h_0+\zeta-b)}$ is the right-going analogue of $R^-$. The corresponding maps $\mathcal{H}_0$ and $\mathcal{H}_\ell$ are then given by

$$ \mathcal{H}_0(\xi^+, \xi^-) = (\xi^+, \xi^-)^T, $$

$$ \mathcal{H}_\ell(\xi^+, \xi^-) = \left( \frac{1}{16g} (\xi^+ - \xi^-)^2 - h_0 + b, \frac{\xi^+ + \xi^-}{32g} (\xi^+ - \xi^-)^2 \right)^T. $$

## 3.2 Equations for the output functions

If Assumption 3.1 holds, the knowledge of the output functions $\xi_0^-(\zeta, q)$ and $\xi_\ell^+(\zeta, q)$, together with the boundary conditions set on the input functions $\xi_0^+(\zeta, q)$ and $\xi_\ell^-(\zeta, q)$ allow one to determine the traces of $\zeta$ and $q$ at $x = 0$ and $x = \ell$. The issue known is to be able to compute the output functions $\xi_0^-(\zeta, q)$ and $\xi_\ell^+(\zeta, q)$. In the (hyperbolic) case of the nonlinear shallow water equations, this can be done using the characteristic equations satisfied by the Riemann invariants, provided that the input functions $\xi_0^\pm$ and $\xi_\ell^\pm$ are correctly chosen (see [20] for a full analysis of 1D hyperbolic initial boundary value problems).

Due to the presence of dispersion, there are no Riemann invariants associated with the Boussinesq-Abbott model. However, we can notice that the relation (37) stems from a more general differential identity that relates $\zeta_0, q_0, \zeta_\ell$ and $q_\ell$, namely,

$$ (42) \qquad \mathcal{S}'_b \begin{pmatrix} \dot{q}_0 \\ \dot{q}_\ell \end{pmatrix} + \tilde{V}_{\text{bdry}}(\zeta_0, q_0, \zeta_\ell, q_\ell) = \tilde{V}_{\text{int}}[\zeta, q] - \begin{pmatrix} \ddot{\zeta}_0 \\ \ddot{\zeta}_\ell \end{pmatrix}; $$

in fact we see that (42) reduces to (37) when enforcing $(\zeta_0, \zeta_\ell) = (g_0, g_\ell)$. We now show how to use this differential identity to compute $\underline{\xi}_0^- := \xi_0^-(\zeta_0, q_0)$ and $\underline{\xi}_\ell^+ := \xi_\ell^+(\zeta_0, q_0)$ when enforcing the general boundary conditions (40). By definition of $\mathcal{H}_0$ and $\mathcal{H}_\ell$, and using the boundary conditions (40), one has

$$ \begin{pmatrix} \zeta_0 \\ q_0 \end{pmatrix} = \mathcal{H}_0(g_0, \underline{\xi}_0^-) \quad \text{and} \quad \begin{pmatrix} \zeta_\ell \\ q_\ell \end{pmatrix} = \mathcal{H}_\ell(\underline{\xi}_\ell^+, g_\ell); $$

denoting $\mathcal{H}_0 = (\mathcal{H}_{0,1}, \mathcal{H}_{0,2})^T$ and $\mathcal{H}_\ell = (\mathcal{H}_{\ell,1}, \mathcal{H}_{\ell,2})^T$, one deduces that

$$ (43) \qquad \dot{q}_0 = \nabla \mathcal{H}_{0,2}(g_0, \underline{\xi}_0^-) \cdot \begin{pmatrix} \dot{g}_0 \\ \dot{\underline{\xi}}_0^- \end{pmatrix} \quad \text{and} \quad \dot{q}_\ell = \nabla \mathcal{H}_{\ell,2}(\underline{\xi}_\ell^+, g_\ell) \cdot \begin{pmatrix} \dot{\underline{\xi}}_\ell^+ \\ \dot{g}_\ell^- \end{pmatrix} $$

as well as

$$ \ddot{\zeta}_0 = \nabla \mathcal{H}_{0,1}(g_0, \underline{\xi}_0^-) \cdot \begin{pmatrix} \ddot{g}_0 \\ \ddot{\underline{\xi}}_0^- \end{pmatrix} + \begin{pmatrix} \dot{g}_0 \\ \dot{\underline{\xi}}_0^- \end{pmatrix} \cdot \text{Hess}_{\mathcal{H}_{0,1}}(g_0, \underline{\xi}_0^-) \begin{pmatrix} \dot{g}_0 \\ \dot{\underline{\xi}}_0^- \end{pmatrix} $$
$$ \ddot{\zeta}_{\ell} = \nabla \mathcal{H}_{\ell,1}(\underline{\xi}_{\ell}^{+}, g_{\ell}) \cdot \left( \frac{\ddot{\underline{\xi}}_{\ell}^{+}}{\dot{g}_{\ell}} \right) + \left( \frac{\dot{\underline{\xi}}_{\ell}^{+}}{\dot{g}_{\ell}} \right) \cdot \mathrm{Hess}_{\mathcal{H}_{\ell,1}}(\underline{\xi}_{\ell}^{+}, g_{\ell}) \left( \frac{\dot{\underline{\xi}}_{\ell}^{+}}{\dot{g}_{\ell}} \right), $$

where $\mathrm{Hess}_{\mathcal{H}_{0,1}}$ and $\mathrm{Hess}_{\mathcal{H}_{\ell,1}}$ denote the $2 \times 2$ Hessian matrices of $\mathcal{H}_{0,1}$ and $\mathcal{H}_{\ell,1}$ respectively. Introducing the matrices

$$ D_j = \begin{pmatrix} \partial_1 \mathcal{H}_{0,j}(g_0, \underline{\xi}_0^-) & 0 \\ 0 & \partial_2 \mathcal{H}_{\ell,j}(\underline{\xi}_\ell^+, g_\ell) \end{pmatrix} $$

and

$$ \tilde{D}_j = \begin{pmatrix} \partial_2 \mathcal{H}_{0,j}(g_0, \underline{\xi}_0^-) & 0 \\ 0 & \partial_1 \mathcal{H}_{\ell,j}(\underline{\xi}_\ell^+, g_\ell) \end{pmatrix}, $$

and the quadratic forms defined on $\mathbb{R}^2$ by

$$ Q_1(u_0, v_0) = \begin{pmatrix} u_0 \\ v_0 \end{pmatrix} \cdot \mathrm{Hess}_{\mathcal{H}_{0,1}}(g_0, \underline{\xi}_0^-) \begin{pmatrix} u_0 \\ v_0 \end{pmatrix}, $$

$$ Q_2(u_\ell, v_\ell) = \begin{pmatrix} u_\ell \\ v_\ell \end{pmatrix} \cdot \mathrm{Hess}_{\mathcal{H}_{\ell,1}}(\underline{\xi}_\ell^+, g_\ell) \begin{pmatrix} u_\ell \\ v_\ell \end{pmatrix}, $$

we deduce from (42) the following systems of two scalar ODEs on $(\underline{\xi}_0^-, \underline{\xi}_\ell^+)$,

$$ \begin{align} \tilde{D}_1 \begin{pmatrix} \ddot{\underline{\xi}}_0^- \\ \dot{\underline{\xi}}_\ell^+ \end{pmatrix} + \mathcal{G}'_b \tilde{D}_2 \begin{pmatrix} \dot{\underline{\xi}}_0^- \\ \dot{\underline{\xi}}_\ell^+ \end{pmatrix} &+ \begin{pmatrix} Q_1(\dot{g}_0, \dot{\underline{\xi}}_0^-) \\ Q_2(\dot{\underline{\xi}}_\ell^+, \dot{g}_\ell) \end{pmatrix} + \tilde{V}_{\mathrm{bdry}}(\mathcal{H}_0(g_0, \underline{\xi}_0^-), \mathcal{H}_\ell(\underline{\xi}_\ell^+, g_\ell)) \tag{44} \\ & = \tilde{V}_{\mathrm{int}}[\zeta, q] - D_1 \begin{pmatrix} \ddot{g}_0 \\ \ddot{g}_\ell \end{pmatrix} - \mathcal{G}'_b D_2 \begin{pmatrix} \dot{g}_0 \\ \dot{g}_\ell \end{pmatrix}. \nonumber \end{align} $$

The order of this system of ODEs depends on whether the coefficients of $\tilde{D}_1$ vanish or not. We make the following assumption which ensures that (44) can be put as a system of explicit ODEs of first or second order on $\underline{\xi}_0^-$ and $\underline{\xi}_\ell^+$.

**Assumption 3.3.** Under Assumption 3.1 and with the same notations, we assume that:

i. Either $\partial_2 \mathcal{H}_{0,1} \equiv 0$ and $\partial_2 \mathcal{H}_{0,2}$ does not vanish on $V_0$, or $\partial_2 \mathcal{H}_{0,1}$ does not vanish on $V_0$;

ii. Either $\partial_1 \mathcal{H}_{\ell,1} \equiv 0$ and $\partial_1 \mathcal{H}_{\ell,2}$ does not vanish on $V_0$, or $\partial_1 \mathcal{H}_{\ell,1}$ does not vanish on $V_\ell$.

**Remark 3.4.** If $\partial_2 \mathcal{H}_{0,1} \equiv 0$ on $V_0$ then $Q_1(\dot{g}_0, \dot{\underline{\xi}}_0^-) = \partial_1^2 \mathcal{H}_{0,1}(g_0, \underline{\xi}_0^-) \dot{g}_0^2$ which is independent of $\underline{\xi}_0^-$. The system of ODEs (44) is therefore of order 1 in $\underline{\xi}_0^-$ and can be put in explicit form if $\partial_2 \mathcal{H}_{0,2}$ does not vanish. A similar comment can be made for the ODE on $\underline{\xi}_\ell^+$.

**Example 3.5.** Considering the same configuration as in Example 3.2, and with the same notations, one readily checks that $\partial_2 \mathcal{H}_{0,1} = 0$ and $\partial_2 \mathcal{H}_{0,2} = 1$ so that the first point of the assumption is satisfied. Moreover, in this case, $Q_1 = 0$. We also compute $\partial_1 \mathcal{H}_{\ell,1} = \frac{1}{8g}(\xi^+ - \xi^-)$. It is possible to choose the open set $V_\ell$
in such a way that $\xi^+ - \xi^- > 0$ for all $(\xi^+, \xi^-) \in V_\ell$ provided that we assume the total water height $h$ never vanishes. Assumption 3.3 is therefore satisfied. Recalling the definition (20) of $\mathfrak{S}'_b$, the system (44) can be put under the form

$$ (45) \qquad T \begin{pmatrix} \dot{\xi}_0^- \\ \ddot{\xi}_\ell^+ \end{pmatrix} + \mathfrak{S}'_b \tilde{D}_2 \begin{pmatrix} 0 \\ \dot{\xi}_\ell^+ \end{pmatrix} + \begin{pmatrix} 0 \\ Q_2(\dot{\xi}_\ell^+, \dot{g}_\ell) \end{pmatrix} + \tilde{V}_{\text{bdry}}(\mathcal{H}_0(g_0, \underline{\xi}_0^-), \mathcal{H}_\ell(\underline{\xi}_\ell^+, g_\ell)) \\ = \tilde{V}_{\text{int}}[\zeta, q] - D_1 \begin{pmatrix} \ddot{g}_0 \\ \ddot{g}_\ell \end{pmatrix} - \mathfrak{S}'_b D_2 \begin{pmatrix} \dot{g}_0 \\ \dot{g}_\ell \end{pmatrix}, $$

where $T$ is the triangular matrix

$$ T = \begin{pmatrix} (\mathfrak{s}^{b,0})'(0) & 0 \\ (\mathfrak{s}^{b,0})'(\ell) & \partial_1 \mathcal{H}_{\ell,1}(\underline{\xi}_\ell^+, g_\ell) \end{pmatrix}. $$

Under the assumption on $b$ made in the statement of Proposition 2.4, one has $(\mathfrak{s}^{b,0})'(0) \neq 0$ and $T$ is invertible, and the system (45) therefore furnishes an ODE in explicit form of first order in $\underline{\xi}_0^-$ and second order on $\underline{\xi}_\ell^+$. For all the applications considered in this paper, we will always have ODEs of second order, except when the boundary condition is imposed on the surface elevation (as in this example at $x = 0$).

We can notice that if $\xi_0^+(\zeta, q) = q$ (boundary condition on the discharge $q$ at $x = 0$), and if we choose $\xi_0^-(\zeta, q) = \zeta$ as output function, then $\mathcal{H}_0(\xi^+, \xi^-) = (\xi^-, \xi^+) = 0$ and in particular $\partial_1 \mathcal{H}_{0,1} = 0$. We make the following assumption to ensure that if we make another choice of input function, then $\partial_1 \mathcal{H}_{0,1}$ does not vanish.

**Assumption 3.6.** Under Assumption 3.1 and with the same notations, we assume that:

i. If $\partial_1 \mathcal{H}_{0,1}$ vanishes on $V_0$ then $\xi_0^+(\zeta, q) = \xi_0^+(q)$ on $U$.

ii. If $\partial_2 \mathcal{H}_{\ell,1}$ vanishes on $V_\ell$ then $\xi_\ell^-(\zeta, q) = \xi_\ell^-(q)$ on $U$.

### 3.3 Compatibility conditions and initial conditions for the output functions

If the initial condition $(\zeta, q)|_{t=0} = (\zeta^{\text{in}}, q^{\text{in}})$ is imposed for the Boussinesq-Abbott system, then a necessary compatibility condition with the boundary conditions (40) to allow solutions that are continuous at $(t = 0, x = 0, \ell)$ is that

$$ (46) \qquad \xi_0^+(\zeta^{\text{in}}(0), q^{\text{in}}(0)) = g_0(0) \quad \text{and} \quad \xi_\ell^-(\zeta^{\text{in}}(\ell), q^{\text{in}}(\ell)) = g_\ell(0); $$

these conditions generalize the first condition of (27). In order to solve (44), one also has to prescribe initial data on $\underline{\xi}_0^-$ and $\underline{\xi}_\ell^+$; one naturally takes

$$ (47) \qquad \underline{\xi}_0^-(0) = \underline{\xi}_0^-(\zeta^{\text{in}}(0), q^{\text{in}}(0)) \quad \text{and} \quad \underline{\xi}_\ell^-(0) = \underline{\xi}_\ell^+(\zeta^{\text{in}}(\ell), q^{\text{in}}(\ell)). $$

By definition of $\mathcal{H}_0$ and $\mathcal{H}_\ell$, we know that

$$ \dot{\zeta}_0 = \partial_1 \mathcal{H}_{0,1}(g_0, \underline{\xi}_0^-) \dot{g}_0 + \partial_2 \mathcal{H}_{0,1}(g_0, \underline{\xi}_0^-) \dot{\underline{\xi}}_0^-, $$
$$ \dot{\zeta}_{\ell} = \partial_{1}\mathcal{H}_{\ell,1}(\xi_{\ell}^{+}, g_{\ell})\dot{\xi}_{\ell}^{+} + \partial_{2}\mathcal{H}_{\ell,1}(\xi_{\ell}^{+}, g_{\ell})\dot{g}_{\ell}. $$

Using the first equation of the Boussinesq-Abbott system (15), we can replace $\dot{\zeta}_0(t)$ by $-\partial_x q(t, 0)$ (and proceed similarly at $x = \ell$) to obtain

$$ (48) \quad \left\{ \begin{array}{l} \partial_2 \mathcal{H}_{0,1}(g_0(0), \underline{\xi}_0^-(0)) \dot{\underline{\xi}}_0^- = -\partial_x q^{\text{in}}(0) - \partial_1 \mathcal{H}_{0,1}(g_0(0), \underline{\xi}_0^-(0)) \dot{g}_0(0), \\ \partial_1 \mathcal{H}_{\ell,1}(\underline{\xi}_\ell^+(0), g_\ell(0)) \dot{\underline{\xi}}_\ell^+ = -\partial_x q^{\text{in}}(\ell) - \partial_2 \mathcal{H}_{\ell,1}(\underline{\xi}_\ell^+(0), g_\ell(0)) \dot{g}_\ell(0), \end{array} \right. $$

where we assumed that the solution was regular enough to take the trace at $(t=0, x=0, \ell)$. Depending on the situation, these conditions can be either a compatibility condition in the boundary and initial data to allow the possibility of regular solutions, or an initial condition for $\dot{\underline{\xi}}_0^-(0)$ or $\dot{\underline{\xi}}_\ell^+(0)$. More precisely:

* If (44) is of second order on $\underline{\xi}_0^-$ (resp. on $\underline{\xi}_\ell^+$), then an initial condition is also needed on $\dot{\underline{\xi}}_0^-$ (resp. on $\dot{\underline{\xi}}_\ell^+$). Remarking further that under Assumption 3.3, $\partial_2 \mathcal{H}_{0,1}$ does not vanish if (44) is of second order in $\underline{\xi}_0^-$ (resp. $\partial_1 \mathcal{H}_{\ell,1}$ does not vanish if (44) is of second order in $\underline{\xi}_\ell^+$), (48) furnishes the needed initial data for $\dot{\underline{\xi}}_0^-$ (resp. on $\dot{\underline{\xi}}_\ell^+$).

* If (44) is of first order on $\underline{\xi}_0^-$ (resp. on $\underline{\xi}_\ell^+$) then by Assumption 3.3 one has $\partial_2 \mathcal{H}_{0,1} \equiv 0$ (resp. $\partial_1 \mathcal{H}_{\ell,1} \equiv 0$), and (48) is then a compatibility condition on the data that generalizes the second equations in the compatibility conditions (27) imposed in Proposition 2.11.

**Example 3.7.** For the configuration considered in Examples 3.2 and 3.5, and assuming that (46) holds, the ODE is of first order on $\underline{\xi}_0^-$ and of second order in $\underline{\xi}_\ell^+$. The initial conditions we need to impose are therefore (47) and also, at $x = \ell$, the second condition of (48). In order to expect $C^1$ solutions at the origin, we also need to impose at $x = 0$ a compatibility condition on the initial and boundary data, which is given by the first equation of (47) (and which coincides in this case with the compatibility condition (27), namely, $-\partial_x q^{\text{in}}(0) = \dot{g}_0$).

## 3.4 Well-posedness of the initial boundary value problem with general boundary conditions

We can now state the main result of this paper, which proves the well-posedness of the Boussinesq-Abbott system (15) with general boundary conditions (40). Proposition 2.11 is a particular case of the theorem, corresponding to $\xi_0^+(\zeta, q) = \xi_\ell^-(\zeta, q) = \zeta$.

**Theorem 3.8.** Assume that $b$ is a smooth function satisfying (14) and (21). Let also $\xi_0^+$ and $\xi_\ell^-$ be two input functions, $\xi_0^-$ and $\xi_\ell^+$ be two output functions, and $\mathcal{H}_0$ and $\mathcal{H}_\ell$ be reconstruction mappings that satisfy Assumptions 3.1, 3.3 and 3.6.

Let $g_0, g_\ell \in C^\infty(\mathbb{R}^+)$. Let also $n \in \mathbb{N}\setminus\{0\}$ and $(\zeta^{\text{in}}, q^{\text{in}}) \in H^n(0, \ell) \times H^{n+1}(0, \ell)$
be such that $\inf_{[0,\ell]}(h_0 - b + \zeta^{\text{in}}) > 0$ and assume that the compatibility condition (46) holds. If moreover $\partial_2 \mathcal{H}_{0,1} \equiv 0$ (resp. $\partial_1 \mathcal{H}_{\ell,1} \equiv 0$) then we also assume that the first (resp. the second) compatibility condition of (48) holds. Then there exists a maximal existence time $T^* > 0$ and a unique solution $(\zeta, q) \in C^\infty([0, T^*]; H^n(0, \ell) \times H^{n+1}(0, \ell))$ to the Boussinesq-Abbott system (15) with initial conditions (26) and boundary conditions (40), and moreover for all $t \in (0, T^*)$, one has $\inf_{[0,\ell]}(h_0 - b + \zeta(t, \cdot)) > 0$.

*Proof.* From the analysis of §3.2 and with the same notations, we know that if such a solution exists, then $(\zeta, q)$ solves the system (35). Substituting for $\dot{q}_0$ and $\dot{q}_\ell$ in the right-hand side of (35) using (43), we obtain that

$$ (49) \quad \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1(\frac{1}{h_b^2} \tilde{f}_{\text{NSW}})) = R_b^0 \tilde{B} \\ + \nabla \mathcal{H}_{0,2}(g_0, \underline{\xi}_0^-) \cdot \left( \begin{array}{c} \dot{g}_0 \\ \dot{\underline{\xi}}_0^- \end{array} \right) \mathfrak{s}^{(b,0)} + \nabla \mathcal{H}_{\ell,2}(\underline{\xi}_\ell^+, g_\ell) \cdot \left( \begin{array}{c} \dot{\underline{\xi}}_\ell^+ \\ \dot{g}_\ell \end{array} \right) \mathfrak{s}^{(b,\ell)}, \end{array} \right. $$

This system is complemented by the ODE (44) for $(\underline{\xi}_0^-, \underline{\xi}_\ell^+)$ which can be of first or second order in $\underline{\xi}_0^-$ or $\underline{\xi}_\ell^+$. We now distinguish several cases and, for the sake of clarity, focus our attention at $x = 0$, the adaptations to $x = \ell$ being straightforward.

* If $\partial_2 \mathcal{H}_{0,1} \equiv 0$, then by Assumption 3.3 we know that $\partial_2 \mathcal{H}_{0,2}$ does not vanish and the ODE (44) is explicit and of first order in $\underline{\xi}_0^-$. The needed initial condition on $\underline{\xi}_0^-$ is furnished by (47).

* If this is not the case then $\partial_2 \mathcal{H}_{0,1}$ does not vanish in virtue of Assumption 3.3, and the ODE (44) is explicit and of second order in $\underline{\xi}_0^-$. In addition to the initial condition on $\underline{\xi}_0^-$ furnished by (47), we need an initial condition on $\underline{\xi}_0^-$, which is furnished by (48).

In all cases, we can use as in the proof of Proposition 2.11 the Cauchy-Lipschitz theorem to construct a solution $(\zeta, q, \underline{\xi}_0^-, \underline{\xi}_\ell^+)$ to the initial value problem formed by (49), (44) and the initial conditions. We need to prove now that $(\zeta, q)$ solves the Boussinesq-Abbott system (15) with initial conditions (26) and boundary conditions (40).

Applying $(1+h_b T_b)$ to (49), we obtain that $(\zeta, q)$ solves (15); it remains to prove that the boundary conditions (40) are satisfied. Taking the trace of the second equation of (49) at $x = 0$ and $x = \ell$, we deduce from (18), (19) and (22) that

$$ \frac{d}{dt} q(t, 0) = \frac{d}{dt} \mathcal{H}_{0,2}(g_0, \underline{\xi}_0^-) \quad \text{and} \quad \frac{d}{dt} q(t, \ell) = \frac{d}{dt} \mathcal{H}_{\ell,2}(\underline{\xi}_\ell^+, g_\ell). $$

Since we also know from Assumption 3.1 and (47) that $\mathcal{H}_{0,2}(g_0, \underline{\xi}_0^-)|_{t=0} = q^{\text{in}}(0)$ and $\mathcal{H}_{\ell,2}(\underline{\xi}_\ell^+, g_\ell)|_{t=0} = q^{\text{in}}(\ell)$, it follows that for all time, one has $q(t, 0) = \mathcal{H}_{0,2}(g_0, \underline{\xi}_0^-)$ and $q(t, \ell) = \mathcal{H}_{\ell,2}(\underline{\xi}_\ell^+, g_\ell)$. We here again have to distinguish two cases, and as above we focus our attention at $x = 0$.
* If $\xi_0^+(\zeta, q) = q$ then one has $H_{0,2}(\xi^+, \xi^-) = \xi^+$ so that we deduce that $q(t, 0) = g_0$, which is exactly the boundary condition (40).

* If this is not the case, then we know by Assumption 3.6 that $\partial_1 H_{0,1}(\xi^+, \xi^-)$ does not vanish on $V_0$; therefore the ODE

$$ \tilde{D}_1 \begin{pmatrix} \ddot{\xi}_0 \\ \dot{\xi}_0 \\ \ddot{\xi}_\ell \\ \dot{\xi}_\ell \end{pmatrix} + \mathcal{S}'_b \tilde{D}_2 \begin{pmatrix} \dot{\xi}_0 \\ \xi_0 \\ \dot{\xi}_\ell \\ \xi_\ell \end{pmatrix} + \begin{pmatrix} Q_1(\dot{\xi}_0^+, \dot{\xi}_0^-) \\ Q_2(\dot{\xi}_\ell^+, \dot{\xi}_\ell^-) \end{pmatrix} + \tilde{V}_{\text{bdry}}(H_0(\xi_0^+, \xi_0^-), H_\ell(\xi_\ell^+, \xi_\ell^-)) \\ = \tilde{V}_{\text{int}}[\zeta, q] - D_1 \begin{pmatrix} \ddot{\xi}_0 \\ \dot{\xi}_0 \\ \ddot{\xi}_\ell \\ \dot{\xi}_\ell \end{pmatrix} - \mathcal{S}'_b D_2 \begin{pmatrix} \dot{\xi}_0 \\ \xi_0 \\ \dot{\xi}_\ell \\ \xi_\ell \end{pmatrix}. $$

which is obtained as (44), is of second order in $\xi_0^+ = \xi_0^+(\zeta_0, q_0)$. It follows that $(\xi_0^+, \xi_\ell^-)$ satisfy the same ODE as $(g_0, g_\ell)$. These two quantities are therefore equal (which implies that the boundary condition (40) is satisfied) if $(\xi_0^+(0), \xi_\ell^-(0)) = (g_0(0), g_\ell(0))$ and $(\xi_0^+(0), \dot{\xi}_\ell^-(0)) = (\dot{g}_0(0), \dot{g}_\ell(0))$. The first of these two conditions corresponds to (46). For the second one (focusing our attention on the case $x = 0$), we observe by time differentiating the relation $H_{0,1}(\xi_0^+, \xi_0^-) = \zeta(t, 0)$ and using the first equation of (15) that

$$ \partial_1 H_{0,1}(\xi_0^+(0), \xi_0^-(0)) \dot{\xi}_0^+(0) + \partial_2 H_{0,1}(\xi_0^+(0), \xi_0^-(0)) \dot{\xi}_0^-(0) = -\partial_x q^{\text{in}}(0). $$

Since $\xi_0^+(0) = g(0)$ by (46), we can use (48) (which holds by assumption if $\partial_2 H_{0,1} = 0$ and by the initial condition imposed on $\dot{\xi}_0^-$ otherwise) to obtain

$$ \partial_1 H_{0,1}(\xi_0^+(0), \xi_0^-(0)) (\dot{\xi}_0^+(0) - \dot{g}_0(0)) = 0, $$

and since $\partial_1 H_{0,1}$ does not vanish, this implies that $\dot{\xi}_0^+(0) = \dot{g}_0(0)$, so that the boundary condition (40) is satisfied.

The proof of the theorem is then complete.

□

## 3.5 Asymptotic stability

For some applications, we have some information on $\zeta$ and $q$ at $x = 0$ and $x = \ell$; this provides us with boundary conditions of the form (40) on some input functions $\xi_0^+$ and $\xi_\ell^-$. For instance, we may know the surface elevation $\zeta(t, 0)$ and $\zeta(t, \ell)$ through buoys located at $x = 0$ and $x = \ell$. By Theorem 3.8 we are able to compute the solution $(\zeta, q)$ of the Boussinesq-Abbott equations (15) on the full domain $(0, \ell)$ for some time interval $(0, T^*)$ with $T^* > 0$ provided that we know the initial data $(\zeta^{\text{in}}, q^{\text{in}})$. Unfortunately, for most applications, the initial data are not known (it is very complicated to measure the surface elevation and the horizontal discharge on the whole interval $(0, \ell)$). A question of high practical relevance is therefore the following: if we consider the solution $(\tilde{\zeta}, \tilde{q})$ of the initial boundary value problem with the same boundary conditions but with initial data $(\tilde{\zeta}^{\text{in}}, \tilde{q}^{\text{in}}) \neq (\zeta^{\text{in}}, q^{\text{in}})$ do we have $(\tilde{\zeta}, \tilde{q}) \sim (\zeta, q)$ for large times? If
this is the case, we shall say that the solution $(\zeta, q)$ of the initial boundary value problem formed by (15) with initial conditions (26) and boundary conditions (40) is asymptotically stable.

Asymptotic stability cannot be expected in general, even if we consider only the linear equation. Consider for instance the linearized Boussinesq equations with a flat topography

$$ (50) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ (1 - \frac{h_0^2}{3} \partial_x^2) \partial_t q + gh_0 \partial_x \zeta = 0, \end{cases} $$

with boundary conditions

$$ (51) \qquad \zeta(t, 0) = \zeta(t, \ell) = 0. $$

The rest state $(\zeta, q) = (0, 0)$ is obviously a solution to this problem, associated with a homogeneous initial condition. This following proposition shows that the rest state is not asymptotically stable. In the statement, we use the notation

$$ (52) \qquad E(t) = \frac{1}{2} g |\tilde{\zeta}(t)|_{L^2(0, \ell)}^2 + \frac{1}{2} \frac{1}{h_0} |\tilde{q}(t)|_{L^2(0, \ell)}^2 + \frac{h_0}{6} |\partial_x \tilde{q}(t)|_{L^2(0, \ell)}^2 $$

**Proposition 3.9.** If $(\tilde{\zeta}, \tilde{q})$ is a smooth solution to (50) and (51) then for all times, one has $E(t) = E(0)$. In particular, if $(\tilde{\zeta}, \tilde{q})|_{t=0} \neq (0, 0)$ then one cannot have $(\tilde{\zeta}, \tilde{q}) \to (0, 0)$ in $L^2(0, \ell) \times H^1(0, \ell)$.

*Proof.* Multiplying (50) by $(g\tilde{\zeta}, \frac{1}{h_0}\tilde{q})$ and integrating by parts and using the boundary conditions (51), we obtain

$$ \frac{d}{dt} E - \frac{h_0}{3} [\partial_x \partial_t \tilde{q}]_0^\ell = 0. $$

Using the first equation of (50), one has $\partial_x \partial_t \tilde{q} = -\partial_t^2 \tilde{\zeta}$ so that using (51) we deduce that $[\partial_x \partial_t \tilde{q}]_0^\ell = 0$, and the result follows easily. $\square$

The fact that we allow in (40) very general boundary conditions is important because it may allow us to use boundary conditions for which the solution furnished by Theorem 3.8 is asymptotically stable (provided that it is well defined globally in time of course). This question of asymptotic stability is a difficult open problem for the Boussinesq-Abbott equations. Even in the case of the much more widely studied nonlinear shallow water equations, very little is known. The most relevant result is the fact that the rest state is asymptotically stable if the input functions $\xi_0^+$ and $\xi_\ell^-$ at $x = 0$ and $x = \ell$ are respectively the right and left-going Riemann invariants $R^+$ and $R^-$ [16].

In §5.2, we investigate numerically the asymptotic stability of the solution of the initial boundary value problem associated with various types of boundary conditions.
## 4 Numerical schemes

In this section we detail the numerical approximation of the Boussinesq-Abbott system with topography (15), namely,

$$ \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0 \\ (1 + h_b T_b) \partial_t q + \partial_x f_{NSW} = -gh \partial_x b \quad \text{in } (0, \ell), \end{array} \right. $$

with the general boundary conditions (40), that is, with the notations of Section 3,

$$ \xi_0^+(\zeta_0, q_0) = g_0, \quad \xi_\ell^-(\zeta_\ell, q_\ell) = g_\ell, $$

and with initial condition

$$ (\zeta, q) = (\zeta^{\text{in}}, q^{\text{in}}) \quad \text{at} \quad t=0. $$

Under the assumptions of Theorem 3.8, we know that there is a unique solution to this initial boundary-value problem. We show in this section how to approximate numerically this solution. We actually solve the reformulation (49) of the equations that reads

$$ (53) \quad \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \sqrt{3} \frac{h_b}{\alpha_b} (h_b S)^* (R_b^1 (\frac{1}{h_b^2} \tilde{f}_{\text{NSW}})) = R_b^0 \tilde{B} \\ + \nabla \mathcal{H}_{0,2}(g_0, \underline{\xi}_0^-) \cdot \begin{pmatrix} \dot{g}_0 \\ \dot{\underline{\xi}}_0^- \end{pmatrix} \mathfrak{s}^{(b,0)} + \nabla \mathcal{H}_{\ell,2}(\underline{\xi}_\ell^+, g_\ell) \cdot \begin{pmatrix} \dot{\underline{\xi}}_\ell^+ \\ \dot{g}_\ell \end{pmatrix} \mathfrak{s}^{(b,\ell)}, \end{array} \right. $$

where $(\underline{\xi}_0^+, \underline{\xi}_\ell^-)$ solves the system of ODEs (44) that we rewrite here for the sake of clarity,

$$ (54) \quad \begin{aligned} & \tilde{D}_1 \begin{pmatrix} \ddot{\underline{\xi}}_0^- \\ \dot{\underline{\xi}}_\ell^+ \end{pmatrix} + \mathfrak{S}'_b \tilde{D}_2 \begin{pmatrix} \dot{\underline{\xi}}_0^- \\ \dot{\underline{\xi}}_\ell^+ \end{pmatrix} + \begin{pmatrix} Q_1(\dot{g}_0, \dot{\underline{\xi}}_0^-) \\ Q_2(\dot{\underline{\xi}}_\ell^+, \dot{g}_\ell) \end{pmatrix} + \tilde{V}_{\text{bdry}}(\mathcal{H}_0(g_0, \underline{\xi}_0^-), \mathcal{H}_\ell(\underline{\xi}_\ell^+, g_\ell)) \\ & = \tilde{V}_{\text{int}}[\zeta, q] - D_1 \begin{pmatrix} \ddot{g}_0 \\ \ddot{g}_\ell \end{pmatrix} - \mathfrak{S}'_b D_2 \begin{pmatrix} \dot{g}_0 \\ \dot{g}_\ell \end{pmatrix}. \end{aligned} $$

As explained in Section 2.3, this reformulation is adapted to well-balancedness in the sense that it enables to derive numerical schemes which naturally preserve the hydrostatic equilibrium $(\zeta, q) = (0, 0)$.

Our strategy consists in an hybrid approach mixing a finite volumes scheme for the interior equations (53), and a finite difference discretization of the nonlocal operators $R_b^0, R_b^1$ defined in (18)-(22) and of the system of ODEs (54) that relates the evolution of the output functions at the boundaries to the boundary data. We explain in Section 4.1 how to discretize the various operators involved in (53); we then propose in Section 4.2 a first order Lax-Friedrichs scheme, and a second order MacCormack scheme in Section 4.3.
**Notations.** We introduce a few notations used in the following lines. We consider the grid points $x_i = (i - 1)\Delta x$ for all $1 \le i \le N$ with $(N - 1)\Delta x = \ell$. A dual mesh is then obtained as the $N$ cells centered on the $(x_i)_{1 \le i \le N}$ and delimited by their interfaces $x_{i\pm 1/2} = x_i \pm \Delta x/2$. The use of this dual mesh is not mandatory, but it will allow us to deal with the boundary conditions (40) and the dispersive boundary layer directly in the border cells which are centered on $x_1 = 0$ and $x_N = \ell$. In particular, this rids us from the need to use interpolation formulas whenever a quantity has to be evaluated at the boundary of the domain. We denote the discrete times $t^n = n\Delta t$, and we consider $U_i^n = (\zeta_i^n, q_i^n)^T$ an approximation of $\frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} (\zeta, q)^T(t^n, s)ds$, the average of the solution in the cell $i$ at time $t^n$. A similar notation is adopted for the bathymetry $b$.

## 4.1 Discrete operators

In order to discretize $(h_bS)$, its adjoint $(h_bS)^*$ and the nonlocal operators $R_b^0, R_b^1$, we shall make use of the second order finite difference operator $\delta_x : \mathbb{R}^N \to \mathbb{R}^N$ defined such that for any $v \in \mathbb{R}^N$ and $1 \le i \le N$ we have

$$ (55) \qquad \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} $$

For any index $1 \le i \le N$, we define the quantities

$$ \alpha_{b,i} = 1 + \frac{1}{4}(\delta_x(b)_i)^2, \quad h_{b,i} = h_0 - b_i, $$

and whenever considering the discrete setting, $\alpha_b$ and $h_b$ will from now on refer to the vectors with corresponding coefficients $(\alpha_{b,i})_i$ and $(h_{b,i})_i$. It is then possible to approximate $(h_bS)$ by $(h_bS)$ as follows:

$$ (56) \qquad \forall v \in \mathbb{R}^N, \quad (h_bS)(v)_i = -\frac{h_{b,i}^2}{\sqrt{3}} \delta_x \left(\frac{v}{h_b}\right)_i + \frac{\sqrt{3}}{2} \delta_x(b)_i v_i, $$

which is consistent with the definition (17) of $S(\cdot)$. In a similar fashion we introduce $(h_bS)^*$ defined as

$$ (57) \qquad (h_bS)^*(f)_i = \frac{\alpha_{b,i}}{\sqrt{3}h_{b,i}} \delta_x \left(\frac{h_b^2 f}{\alpha_b}\right)_i + \left(\frac{h_{b,i}}{\sqrt{3}} \frac{\delta_x(\alpha_b)_i}{\alpha_{b,i}} + \frac{\sqrt{3}}{2} \delta_x b_i \right) f_i, $$

and that is shown to be consistent with $(h_bS)^*$ up to a $O(\Delta x^2)$ error.

Next we discretize the nonlocal operator $R_b^0$ by $\underline{R}_b^0$ defined implicitly as

$$ \underline{R}_b^0 : f \in \mathbb{R}^N \mapsto v = \underline{R}_b^0 f \in \mathbb{R}^N \text{ such that} $$

$$ (58) \qquad \left\{ \begin{array}{l} v_i - \delta_x \left( \frac{h_b^3}{3} \delta_x \left( \frac{v}{h_b} \right) \right)_i + \frac{h_{b,i}}{2\Delta x^2} (b_{i+1} - 2b_i + b_{i-1}) v_i = f_i, \quad 1 < i < N, \\ v_1 = 0, \quad v_N = 0. \end{array} \right. $$
and one readily checks that this is a consistent approximation of $(1 + h_b T_b)$ at order 2, with the operator $T_b$ given in (16). Likewise $R_b^1$ is discretized by $\underline{R}_b^1$ defined implicitly as

$$ (59) \qquad \underline{R}_b^1 : f \in \mathbb{R}^N \mapsto v = \underline{R}_b^1 f \in \mathbb{R}^N \text{ such that} \\ \left\{ \begin{array}{ll} v_i + \frac{1}{h_{b,i}} (\underline{h}_b S) \left( \frac{h_b}{\alpha_b} (\underline{h}_b S)^*(v) \right)_i = f_i, & 1 < i < N, \\ (\underline{h}_b S)^*(v)_1 = 0, & (\underline{h}_b S)^*(v)_N = 0, \end{array} \right. $$

which is consistent with (22). The definitions (58)-(59) are implicit since the practical computation of $\underline{R}_b^0$ and $\underline{R}_b^1$ applied to some vector will require the resolution of a linear system, and this will be the costliest step of the algorithm. However, these discrete operators do not evolve in time and it is thus possible to assemble them and perform their factorization in a prepossessing step, which greatly reduces the computational time.

**Remark 4.1.** Owing to the homogeneous conditions from (59) verified in the border cells, one has $(\underline{h}_b S)^*(f)_1 = (\underline{h}_b S)^*(f)_N = 0$ for any vector $f$ in $\text{Im}(\underline{R}_b^1)$. This is the discrete counterpart to the fact that $(\underline{h}_b S)^*(\underline{R}_b^1)_{|_{x=0,\ell}} = 0$ in $L^2(0, \ell)$, which enables to obtain the relation (54) — ensuing from Lemma 2.10 after applying $\partial_x$ to the second equation of (53) and taking the trace at $x = 0, \ell$ — and hereby to ensure that the boundary conditions (40) are satisfied. Therefore it is important to have this property at the discrete level to remain consistent with the continuous reformulation (53)-(54).

## 4.2 First order Lax-Friedrichs scheme

We begin by describing how the system (53) for the interior values is approximated. For the first equation, the free surface elevation is updated in cells $2 \le i \le N - 1$ using a finite volumes strategy

$$ (60) \qquad \frac{\zeta_i^{n+1} - \zeta_i^n}{\Delta t} + \frac{1}{\Delta x} (q_{i+1/2}^n - q_{i-1/2}^n) = 0, $$

where $q_{i\pm 1/2}^n$ is the Lax-Friedrichs numerical flux defined by

$$ (61) \qquad q_{i+1/2}^n = \frac{1}{2}(q_i^n + q_{i+1}^n) - \frac{\Delta x}{2\Delta t}(\zeta_{i+1}^n - \zeta_i^n). $$

The boundary values $\zeta_1^{n+1}, \zeta_N^{n+1}$ will be deduced upon approximating the ODE for the output functions. To approximate the second equation of (53) we first introduce the vectors $\tilde{f}_{\text{NSW}}^n, \tilde{B}^n \in \mathbb{R}^N$ with coefficients

$$ (62) \qquad \left\{ \begin{array}{ll} \tilde{f}_{\text{NSW},i}^n = \frac{(q_i^n)^2}{h_0 + \zeta_i^n - b_i} + \frac{g}{2}((\zeta_i^n)^2 + 2h_{b,i}\zeta_i^n), & \tilde{B}_i^n = -g\zeta_i^n\delta_x(b)_i + \frac{3}{2}\frac{\delta_x(b)_i}{h_{b,i}}\tilde{f}_{\text{NSW},i}^n, \end{array} \right. $$
which is consistent with the definition (36) of $\tilde{B}(\zeta, q)$. The update of the discharge in cells $1 \le i \le N$ makes use of the discrete operators introduced in Section 4.1 as follows:

$$ (63) \quad \frac{q_i^{n+1} - q_i^n}{\Delta t} + \sqrt{3} \frac{h_{b,i}}{\alpha_{b,i}} (h_b S)^* \left( \underline{R}_b^1 \left( \frac{1}{h_b^2} \tilde{f}_{NSW}^n \right) \right)_i = \\ \underline{R}_b^0 (\tilde{B}^n)_i + \delta_t q_1^n \cdot \underline{\mathfrak{s}}_i^{(b,0)} + \delta_t q_N^n \cdot \underline{\mathfrak{s}}_i^{(b,\ell)} + \frac{\Delta x^2}{2\Delta t} \frac{q_{i+1}^n - 2q_i^n + q_{i-1}^n}{\Delta x^2} \mathbf{1}_{i \notin \{1,N\}}, $$

where it remains to specify how to compute $\delta_t q_1^n$ and $\delta_t q_N^n$, and where the vectors $\underline{\mathfrak{s}}_i^{(b,0)}$ and $\underline{\mathfrak{s}}_i^{(b,\ell)}$ are taken as the solutions $v \in \mathbb{R}^N$ of

$$ (64) \quad v_i - \delta_x \left( \frac{h_b^3}{3} \delta_x \left( \frac{v}{h_b} \right) \right)_i + \frac{h_{b,i}}{2\Delta x^2} (b_{i+1} - 2b_i + b_{i-1}) v_i = 0, \quad 2 \le i \le N-1 $$

with respective conditions

$$ (65) \quad \begin{cases} \underline{\mathfrak{s}}_1^{(b,0)} = 1 \\ \underline{\mathfrak{s}}_N^{(b,0)} = 0 \end{cases} , \quad \begin{cases} \underline{\mathfrak{s}}_1^{(b,\ell)} = 0 \\ \underline{\mathfrak{s}}_N^{(b,\ell)} = 1 \end{cases} , $$

yielding a discrete counterpart of (19). Owing to the conditions (65), to the definition (58) of $\underline{R}_b^0$ and to Remark 4.1, the update (63) in border cells $i = 1, N$ reduces to

$$ (66) \quad \frac{q_1^{n+1} - q_1^n}{\Delta t} = \delta_t q_1^n, \qquad \frac{q_N^{n+1} - q_N^n}{\Delta t} = \delta_t q_N^n. $$

**Remark 4.2.** We also want to comment on the last term in the right hand side of (63), which is a first order numerical diffusion required for stability purposes. We can then show that the update (63) amounts to a finite volumes method with Lax-Friedrichs nonlocal numerical flux and with source terms accounting for the bathymetry and the boundary layer. In fact, defining the nonlocal flux vector as

$$ (67) \quad \mathfrak{f}^n = \left( \frac{h_{b,i}^2}{\alpha_{b,i}} \underline{R}_b^1 \left( \frac{1}{h_b^2} \tilde{f}_{NSW}^n \right) \right)_i, \quad 1 \le i \le N, $$

and the associated Lax-Friedrichs flux

$$ (68) \quad \mathfrak{f}_{i+1/2}^n = \frac{1}{2} (\mathfrak{f}_i^n + \mathfrak{f}_{i+1}^n) - \frac{\Delta x}{2\Delta t} (q_{i+1}^n - q_i^n), $$

we can show that for $2 \le i \le N-1$ the update (63) equivalently rewrites

$$ (69) \quad \frac{q_i^{n+1} - q_i^n}{\Delta t} + \frac{\mathfrak{f}_{i+1/2}^n - \mathfrak{f}_{i-1/2}^n}{\Delta x} = \\ \underline{R}_b^0 (\tilde{B}^n)_i - \left( \frac{\delta_x(\alpha_b)_i}{\alpha_{b,i}} + \frac{3}{2} \frac{\delta_x(b)_i}{h_{b,i}} \right) \mathfrak{f}_i^n + \delta_t q_1^n \cdot \underline{\mathfrak{s}}_i^{(b,0)} + \delta_t q_N^n \cdot \underline{\mathfrak{s}}_i^{(b,\ell)}. $$

When the bottom is flat, the vectors $\delta_x(b)$ and $\tilde{B}^n$ cancel, so that the discharge update (69) is conservative up to the boundary cells, which is coherent with the
continuous model. Note also that when considering the scheme for the Boussinesq-Abbott model with flat topography ($b \equiv 0$), we can directly use the expression (5) to compute $\mathfrak{s}^{(b,0)}$ and $\mathfrak{s}^{(b,\ell)}$ instead of approximating them through (64)-(65). However in practice this does not affect the numerical results in a noticeable way.

Finally we detail the handling of the border cells and the computation of $\delta_t q_1^n$, $\delta_t q_N^n$ in the framework of general boundary conditions discussed in Section 3. Assuming Assumption 3.1 holds, the knowledge of $U_1^{n+1}$, $U_N^{n+1}$ is deduced from the reconstruction formulas (41) reading as

$$U_1^{n+1} = \mathcal{H}_0(\xi_0^+(U_1^{n+1}), \xi_0^-(U_1^{n+1})) \quad \text{and} \quad U_N^{n+1} = \mathcal{H}_\ell(\xi_\ell^+(U_N^{n+1}), \xi_\ell^-(U_N^{n+1})).$$

In the above, we set the input functions according to (40) as $\xi_0^+(U_1^{n+1}) = g_0^{n+1}$ and $\xi_\ell^-(U_N^{n+1}) = g_\ell^{n+1}$. On the other hand, the output values $\xi_0^-(U_1^{n+1})$ and $\xi_\ell^+(U_N^{n+1})$ are obtained by discretizing the ODE (54) for the traces as follows. Given $X^n$ and $Y^n$ the respective approximations at time $t^n$ of the functions

$$(70) \qquad X : t \mapsto \begin{pmatrix} \xi_0^-(U(t, 0)) \\ \xi_\ell^+(U(t, \ell)) \end{pmatrix} \in \mathbb{R}^2, \qquad Y : t \mapsto \dot{X}(t),$$

we define the update $X^{n+1}, Y^{n+1}$ as

$$(71) \quad \left\{ \begin{aligned} \frac{X^{n+1} - X^n}{\Delta t} &= Y^{n+1} \\ \tilde{D}_1^n \frac{Y^{n+1} - Y^n}{\Delta t} + \mathfrak{S}_b' \tilde{D}_2^n Y^{n+1} &+ \left( \frac{\mathcal{Q}_1^n \left( \frac{g_0^{n+1} - g_0^{n-1}}{2\Delta t}, Y_1^n \right)}{\mathcal{Q}_2^n \left( Y_2^n, \frac{g_\ell^{n+1} - g_\ell^{n-1}}{2\Delta t} \right)} \right) + \tilde{V}_{\text{bdry}}^n \\ &= \tilde{V}_{\text{int}}^n - D_1^n \left( \frac{\frac{g_0^{n+1} - 2g_0^n + g_0^{n-1}}{\Delta t^2}}{\frac{g_\ell^{n+1} - 2g_\ell^n + g_\ell^{n-1}}{\Delta t^2}} \right) - \mathfrak{S}_b' D_2^n \left( \frac{\frac{g_0^{n+1} - g_0^{n-1}}{2\Delta t}}{\frac{g_\ell^{n+1} - g_\ell^{n-1}}{2\Delta t}} \right) \end{aligned} \right.$$

where the diagonal matrices $D_j^n, \tilde{D}_j^n$ are defined for $j \in \{1, 2\}$ as

$$D_j^n = \begin{pmatrix} \partial_1 \mathcal{H}_{0,j}(g_0^n, X_1^n) & 0 \\ 0 & \partial_2 \mathcal{H}_{\ell,j}(X_2^n, g_\ell^n) \end{pmatrix}, \\ \tilde{D}_j^n = \begin{pmatrix} \partial_2 \mathcal{H}_{0,j}(g_0^n, X_1^n) & 0 \\ 0 & \partial_1 \mathcal{H}_{\ell,j}(X_2^n, g_\ell^n) \end{pmatrix},$$

where we used the quadratic forms

$$\mathcal{Q}_1^n(u_0, v_0) = \begin{pmatrix} u_0 \\ v_0 \end{pmatrix} \cdot \text{Hess}_{\mathcal{H}_{0,1}}(g_0^n, X_1^n) \begin{pmatrix} u_0 \\ v_0 \end{pmatrix}, \\ \mathcal{Q}_2^n(u_\ell, v_\ell) = \begin{pmatrix} u_\ell \\ v_\ell \end{pmatrix} \cdot \text{Hess}_{\mathcal{H}_{\ell,1}}(X_2^n, g_\ell^n) \begin{pmatrix} u_\ell \\ v_\ell \end{pmatrix},$$

and where

$$\mathfrak{S}_b' = \begin{pmatrix} \delta_x(\mathfrak{s}^{(b,0)})_1 & \delta_x(\mathfrak{s}^{(b,\ell)})_1 \\ \delta_x(\mathfrak{s}^{(b,0)})_N & \delta_x(\mathfrak{s}^{(b,\ell)})_N \end{pmatrix}, \qquad \tilde{V}_{\text{bdry}}^n = \begin{pmatrix} \frac{3}{h_{b,1}^2} \tilde{f}_{\text{NSW}}(\mathcal{H}_0(g_0^n, X_1^n)) \\ \frac{3}{h_{b,N}^2} \tilde{f}_{\text{NSW}}(\mathcal{H}_\ell(X_2^n, g_\ell^n)) \end{pmatrix}$$
$$ \tilde{V}_{\text{int}}^{n} = \left( \begin{array}{c} 3(\underline{R}_{b}^{1}(\frac{1}{h_{b}^{2}}\tilde{f}_{\text{NSW}}^{n}))_{1} - (\delta_{x}(\underline{R}_{b}^{0}\tilde{B}^{n}))_{1} \\ 3(\underline{R}_{b}^{1}(\frac{1}{h_{b}^{2}}\tilde{f}_{\text{NSW}}^{n}))_{N} - (\delta_{x}(\underline{R}_{b}^{0}\tilde{B}^{n}))_{N} \end{array} \right). $$

The functions $g_0, g_\ell$ are defined over $t \in \mathbb{R}_+$, therefore the terms $g_0^{n-1}, g_\ell^{n-1}$ found in (71) are not known when $n=0$. We suggest to define them following the relations below:

$$ \frac{g_0(\Delta t) - g_0^{-1}}{2\Delta t} = \frac{-3g_0(0) + 4g_0(\Delta t) - g_0(2\Delta t)}{2\Delta t}, \\ \frac{g_\ell(\Delta t) - g_\ell^{-1}}{2\Delta t} = \frac{-3g_\ell(0) + 4g_\ell(\Delta t) - g_\ell(2\Delta t)}{2\Delta t}. $$

The sequence $(X^n)_n$ is initialized by taking $X^0 = (\xi_0^-(U_1^0), \xi_\ell^+(U_N^0))^T$. As in the continuous case (see the proof of Theorem 3.8), when $\partial_2 \mathcal{H}_{0,1}(g_0^0, X_1^0) \neq 0$ (respectively when $\partial_2 \mathcal{H}_{0,1}(X_2^0, g_\ell^0) \neq 0$), an initial value $Y_1^0$ (respectively $Y_2^0$) must also be provided in order to compute the second line of (71) for $n=0$. These quantities can be obtained by discretizing (48) in the following way

$$ (72) \quad \left\{ \begin{array}{l} \partial_2 \mathcal{H}_{0,1}(g_0(0), X_1^0) Y_1^0 = -\delta_x(q^0)_1 - \partial_1 \mathcal{H}_{0,1}(g_0(0), X_1^0) \frac{g_0(\Delta t) - g_0(0)}{\Delta t}, \\ \partial_1 \mathcal{H}_{\ell,1}(X_2^0, g_\ell(0)) Y_2^0 = -\delta_x(q^0)_N - \partial_2 \mathcal{H}_{\ell,1}(X_2^0, g_\ell(0)) \frac{g_\ell(\Delta t) - g_\ell(0)}{\Delta t}. \end{array} \right. $$

Remark that in (72), if for a given index $1 \le k \le 2$ the coefficient in factor of $Y_k^0$ cancels, then the quadratic function $Q_k^0$ does not depend on $Y_k^0$ anymore (this is the discrete counterpart of Remark 3.4); hence in the second equation of (71) the value $Y_k^0$ is not required to compute $Y_k^1$. The proposed discretization (71) thus automatically adapts to the order of the ODE for the output functions, which dictates as in the continuous case how many initial conditions are required on the output functions.

Once the approximation $X^{n+1}$ of the output functions is known, it is possible to compute the border values $U_1^{n+1} = \mathcal{H}(g_0^{n+1}, X_1^{n+1})$ and $U_N^{n+1} = \mathcal{H}(X_2^{n+1}, g_\ell^{n+1})$, and $\delta_t q_1^n, \delta_t q_N^n$ are deduced from (66) so that we can update the interior discharge through (69). The iteration is complete.

**Remark 4.3.** We need to check that the update (71) becomes well-defined for $\Delta t > 0$ small enough. Indeed, the second equation of this system admits a unique solution $Y^{n+1}$ if the matrix $\tilde{D}_1^n + \Delta t \underline{\mathcal{S}}'_b \tilde{D}_2^n$ is invertible. Its determinant is the second degree polynomial in $\Delta t$ given by

$$ \det(\tilde{D}_1^n) + \Delta t \left( (\underline{\mathcal{S}}'_b)_{1,1}(\tilde{D}_1^n)_{2,2}(\tilde{D}_2^n)_{1,1} + (\underline{\mathcal{S}}'_b)_{2,2}(\tilde{D}_1^n)_{1,1}(\tilde{D}_2^n)_{2,2} \right) \\ + \Delta t^2 \left( \det(\underline{\mathcal{S}}'_b) \det(\tilde{D}_2^n) \right); $$

as a consequence of Assumption 3.3 and Proposition 2.4 this determinant is non zero for $\Delta t$ small enough, so that the update (71) is well-defined.
When approximating the nonlinear Boussinesq-Abbott system (35)-(37) with the scheme described above, in practice we need to restrict the time-step with a CFL condition similar to that usually encountered when discretizing the hyperbolic nonlinear shallow water system, that is to say we take

$$ \frac{\Delta t}{\Delta x} \leq \frac{K}{\max_i |\lambda(U_i^n)|}, $$

with $|\lambda(U_i^n)| = |u_i^n| + \sqrt{gh_i^n}$ the maximum eigenvalue in absolute value of the Jacobian matrix $Df_{NSW}(U_i^n)$. In our numerical experiments such a constraint seems to be a requirement to achieve stable results.

## 4.3 Second order MacCormack scheme

Due to the dispersive nature of the Boussinesq-type models, it can be challenging to increase the order of their numerical approximations while keeping them stable and free of spurious oscillations in the presence of non homogeneous boundary conditions. A strategy introduced in [8] for the numerical simulation of waves interacting with floating objects is the MacCormack discretization method, which we describe in the lines below. It consists in alternating a first order prediction step with a correction step during which the predicted state is used. The final update is then obtained by averaging the prediction and correction steps. This procedure is explained in detail thereafter.

**PREDICTION STEP** — Given an approximation $(U_i^n)_{1 \le i \le N}$ of the solution to (35)-(37) at time $t^n$, we first compute $U_{P,i}^{n+1} = (\zeta_{P,i}^{n+1}, q_{P,i}^{n+1})^T$ the predicted state using a forward Euler method with a left upwinding for the numerical fluxes in every interior cell $2 \le i \le N-1$

$$ (73) \qquad \left\{ \begin{aligned} & \frac{\zeta_{P,i}^{n+1} - \zeta_i^n}{\Delta t} + \frac{1}{\Delta x} (q_i^n - q_{i-1}^n) = 0 \\ & \frac{q_{P,i}^{n+1} - q_i^n}{\Delta t} + \frac{f_i^n - f_{i-1}^n}{\Delta x} = \underline{R}_b^0 (\tilde{B}^n)_{i-1} + \delta_t q_{P,1}^n \cdot \underline{s}_{i-1}^{(b,0)} + \delta_t q_{P,N}^n \cdot \underline{s}_{i-1}^{(b,\ell)} , \\ & \quad - \frac{1}{2} \left( \frac{\alpha_{b,i+1} - \alpha_{b,i-1}}{\alpha_{b,i} \Delta x} + \frac{3}{2} \frac{b_{i+1} - b_{i-1}}{h_{b,i} \Delta x} \right) f_{i-1}^n \end{aligned} \right. $$

where we kept the same definitions for $\underline{R}_b^0, \underline{R}_b^1, f^n, \tilde{B}^n, \underline{s}^{(b,0)}, \underline{s}^{(b,\ell)}, \alpha_b$ as in the previous Section, which involve the second order centered operator $\delta_x$ given in (55). The boundary cells are treated similarly to the first order scheme described in Section 4.2; that is to say we enforce

$$ \xi_0^+(U_{P,1}^{n+1}) = g_0^{n+1}, \quad \xi_\ell^-(U_{P,N}^{n+1}) = g_\ell^{n+1}, $$

then compute the output functions according to (71) to get $(X, Y)_{P}^{n+1}$, and finally use the reconstruction formulas (41) to get $U_{P,1}^{n+1}, U_{P,N}^{n+1}$ and define

$$ \delta_t q_{P,1}^n = \frac{q_{P,1}^{n+1} - q_1^n}{\Delta t}, \quad \delta_t q_{P,N}^n = \frac{q_{P,N}^{n+1} - q_N^n}{\Delta t}. $$
Notice that, although in (73) we did not write the discharge update under a form involving an explicit discretization of $(h_b S)^*$, it is of course possible to rewrite it in a form similar to (63), where $(h_b S)^*$ defined in (57) has to be substituted with its left-upwinded counterpart, and without the numerical diffusion term. A similar comment will hold in the correction step.

CORRECTION STEP — Next the correction term $U_{C,i}^{n+1} = (\zeta_{C,i}^{n+1}, q_{C,i}^{n+1})^T$ is computed using a forward Euler update with right upwinding in the interior cells $2 \le i \le N-1$. It involves the predicted state obtained in the previous step in order to evaluate the flux and source terms:

$$ (74) \quad \left\{ \begin{aligned} & \frac{\zeta_{C,i}^{n+1} - \zeta_i^n}{\Delta t} + \frac{1}{\Delta x} (q_{P,i+1}^{n+1} - q_{P,i}^{n+1}) = 0 \\ & \frac{q_{C,i}^{n+1} - q_i^n}{\Delta t} + \frac{f_{P,i+1}^{n+1} - f_{P,i}^{n+1}}{\Delta x} = \underline{R}_b^0 (\tilde{B}_P^{n+1})_{i+1} + \delta_t q_{C,i}^n \cdot \underline{s}_{i+1}^{(b,0)} + \delta_t q_{C,N}^n \cdot \underline{s}_{i+1}^{(b,l)} \\ & \qquad - \frac{1}{2} \left( \frac{\alpha_{b,i+1} - \alpha_{b,i-1}}{\alpha_{b,i} \Delta x} + \frac{3}{2} \frac{b_{i+1} - b_{i-1}}{h_{b,i} \Delta x} \right) f_{P,i+1}^{n+1} \end{aligned} \right. $$

In the above, the fluxes $f_{P,i}^{n+1}$ and the source term $\tilde{B}_{P,i}^{n+1}$ are defined as in (67) and (62) but making use of the prediction state, that is to say

$$ \forall 1 \le i \le N, \quad \left\{ \begin{aligned} \tilde{B}_{P,i}^{n+1} &= -g \zeta_{P,i}^{n+1} \delta_x(b)_i + \frac{3}{2} \frac{\delta_x(b)_i}{h_{b,i}} \tilde{f}_{\text{NSW}}(U_{P,i}^{n+1}), \\ f_{P,i}^{n+1} &= \frac{h_{b,i}^2}{\alpha_{b,i}} \underline{R}_b^1 \left( \left\{ \frac{1}{h_{b,j}^2} \tilde{f}_{\text{NSW}}(U_{P,j}^{n+1}) \right\}_{1 \le j \le N} \right)_i. \end{aligned} \right. $$

The handling of the border cells is the same as before, except that in (71), for $k \in \{1, 2\}$ the matrices $D_k$, $\tilde{D}_k$ and the quadratic forms $Q_k^n$ need to be evaluated with the predicted values $(X, Y)_P^{n+1}$ obtained in the previous step instead of $(X, Y)^n$, and the vectors $\tilde{V}_{\text{bdry}}^n, \tilde{V}_{\text{int}}^n$ need to be evaluated with $U_P^{n+1}$ instead of $U^n$. Therefore the correction update for the output functions can be written

$$ (75) \quad \left\{ \begin{aligned} \frac{X_C^{n+1} - X^n}{\Delta t} &= Y_C^{n+1} \\ \tilde{D}_{P,1}^{n+1} \frac{Y_C^{n+1} - Y^n}{\Delta t} &+ \underline{\mathcal{S}}'_b \tilde{D}_{P,2}^{n+1} Y_C^{n+1} + \left( \begin{array}{l} Q_{P,1}^{n+1} \left( \frac{g_0^{n+1} - g_0^{n-1}}{2\Delta t}, Y_{P,1}^{n+1} \right) \\ Q_{P,2}^{n+1} \left( Y_{P,2}^{n+1}, \frac{g_\ell^{n+1} - g_\ell^{n-1}}{2\Delta t} \right) \end{array} \right) + \tilde{V}_{\text{bdry},P}^{n+1} \\ &= \tilde{V}_{\text{int},P}^{n+1} - D_{P,1}^{n+1} \left( \frac{\frac{g_0^{n+1} - 2g_0^n + g_0^{n-1}}{\Delta t^2}}{\frac{g_\ell^{n+1} - 2g_\ell^n + g_\ell^{n-1}}{\Delta t^2}} \right) - \underline{\mathcal{S}}'_b D_{P,2}^{n+1} \left( \frac{g_0^{n+1} - g_0^{n-1}}{2\Delta t} \right) \end{aligned} \right. , $$

and after computing $U_{C,1}^{n+1} = \mathcal{H}(g_0^{n+1}, X_{C,1}^{n+1})$ and $U_{C,N}^{n+1} = \mathcal{H}(X_{C,2}^{n+1}, g_\ell^{n+1})$ we set

$$ \delta_t q_{C,1}^n = \frac{q_{C,1}^{n+1} - q_1^n}{\Delta t}, \quad \delta_t q_{C,N}^n = \frac{q_{C,N}^{n+1} - q_N^n}{\Delta t}. $$
FINAL STEP — The final update is defined as the average between the prediction and correction states for all $2 \le i \le N - 1$, and likewise the output functions are averaged in the border cells $i = 1, N$:

(76)

$$ \left\{ \begin{array}{l} \zeta_i^{n+1} = \frac{1}{2}(\zeta_{P,i}^{n+1} + \zeta_{C,i}^{n+1}) \\ q_i^{n+1} = \frac{1}{2}(q_{P,i}^{n+1} + q_{C,i}^{n+1}) \end{array} \right. , \quad \left\{ \begin{array}{l} \xi_0^{-}(U_1^{n+1}) = \frac{1}{2}(\xi_0^{-}(U_{P,1}^{n+1}) + \xi_0^{-}(U_{C,1}^{n+1})) \\ \xi_{\ell}^{+}(U_N^{n+1}) = \frac{1}{2}(\xi_{\ell}^{+}(U_{P,N}^{n+1}) + \xi_{\ell}^{+}(U_{C,N}^{n+1})) \end{array} \right. $$

The final border elevations and discharges $\zeta_1^{n+1}, q_1^{n+1}, \zeta_N^{n+1}, q_N^{n+1}$ are recovered through the formula (41) involving the reconstruction maps $\mathcal{H}_0, \mathcal{H}_{\ell}$ evaluated respectively with $(g_0(t^{n+1}), \xi_0^{-}(U_1^{n+1}))$ and $(\xi_{\ell}^{+}(U_N^{n+1}), g_{\ell}(t^{n+1}))$.

**Remark 4.4.** Nothing prevents us from switching the upwinding direction in the first two steps, so that a right upwinding is performed in the prediction step together with a left upwinding in the correction step. This does not seem to affect the numerical results in practice. What is important is that the direction of upwinding alternates between these two stages, otherwise the method becomes first order in space.

**Remark 4.5.** Since we do not make use of ghost cells, the alternation of the upwinding direction characterizing the MacCormack scheme cannot be performed in the border cells, where we always need to upwind towards the interior of the domain. This does not prevent the scheme from reaching a second order of accuracy in space in the border cells. Indeed (71), (72) and (75) are already second order in space thanks to the use of the discrete operator $\delta_x$ given by (55) to discretize $S'_b, \tilde{V}_{\text{int}}$.

# 5 Numerical simulations

In this section we aim to validate the proposed first order Lax-Friedrichs and second order MacCormack schemes, and to assess the effect that different boundary conditions can have on the solutions; in particular, we numerically investigate the issue of asymptotic stability of the boundary conditions. To this end we shall consider three cases:

1. Input functions given by the elevation, output functions given by the discharge

$$ \xi_0^+(\zeta, q) = \xi_\ell^-(\zeta, q) = \zeta, \quad \xi_0^-(\zeta, q) = \xi_\ell^+(\zeta, q) = q; $$

2. Input functions given by the discharge, output functions given by the elevation

$$ \xi_0^+(\zeta, q) = \xi_\ell^-(\zeta, q) = q, \quad \xi_0^-(\zeta, q) = \xi_\ell^+(\zeta, q) = \zeta; $$
3. Input functions given by the incoming Riemann invariants, output functions given by the outgoing Riemann invariants

$$ \left\{ \begin{array}{l} \xi_0^+(\zeta, q) = u + 2\sqrt{gh} \\ \xi_\ell^-(\zeta, q) = u - 2\sqrt{gh} \end{array} \right. , \qquad \left\{ \begin{array}{l} \xi_0^-(\zeta, q) = u - 2\sqrt{gh} \\ \xi_\ell^+(\zeta, q) = u + 2\sqrt{gh} \end{array} \right. . $$

The gravitational acceleration g will be taken equal to 9.81, and unless specified otherwise the ratio $\Delta t/\Delta x$ will be taken equal to 0.8 for the Lax-Friedrichs scheme and to 0.45 for the MacCormack scheme.

We propose in §5.1 a numerical resolution of the initial boundary value problem with the different boundary conditions presented above. In the case of a flat bottom, considered in §5.1.1, the Boussinesq-Abbott system possesses solitary waves that can be used to study the convergence of the schemes. In the case of a non-flat topography investigated in §5.1.2, this is no longer true and we therefore use a solution computed in a wider domain and for a very refined mesh to study the convergence. The theoretical open problem of asymptotic stability is then numerically investigated in §5.2.

# 5.1 Numerical resolution of the initial boundary value problem

## 5.1.1 The Boussinesq-Abbott system with flat topography

We consider the case of a solitary wave travelling over a flat bottom, which consists in a solution of (1) of the form

$$ (77) \qquad \zeta(t, x) = \tilde{\zeta}(x - x_0 - ct), \qquad q(t, x) = \tilde{q}(x - x_0 - ct), $$

where c is the celerity at which the solitary wave propagates without deformation, and where $x_0 \in \mathbb{R}$ is the initial position of the maximum elevation $\zeta(0, x_0)$ denoted $\zeta_{max}$. Despite the apparent simplicity of this solution, it involves nonlinear and dispersive effects which make it an interesting test-case for numerical validation purposes. As explained in [28], the wave profile $\tilde{\zeta}$ verifies the following second order ODE

$$ (78) \qquad -c^2 \frac{h_0 \tilde{\zeta}}{h_0 + \tilde{\zeta}} + \frac{c^2 h_0^2}{3} \tilde{\zeta}'' + \frac{g}{2} (\tilde{\zeta}^2 + 2h_0 \tilde{\zeta}) = 0. $$

It is obtained by remarking that (77) implies $\tilde{q} = c\tilde{\zeta}$ from the first equation of (1). Injecting this in the second equation of (1) and upon integration in space with the hypothesis that $\tilde{\zeta}$ vanishes at infinity, one recovers (78). Furthermore, keeping in mind that $\tilde{\zeta}(0) = \zeta_{max}$, the celerity is shown to satisfy the relation

$$ c^2 = \frac{\frac{g}{6} \zeta_{\max}^3 + \frac{gh_0}{2} \zeta_{\max}^2}{h_0^2 (\zeta_{\max}/h_0 - \ln(1 + \zeta_{\max}/h_0))} $$
by integrating (78) against $\tilde{\zeta}'$, using again that $\tilde{\zeta}$ cancels at infinity and that $\tilde{\zeta}'(0) = 0$.

Once the solution of (78) has been approximated, it can be used to generate initial and boundary conditions in the domain $(0, \ell)$, and can also serve as a reference solution to measure the error of the proposed numerical schemes. The initial condition is taken as $(\zeta^{\text{in}}, q^{\text{in}})(x) = (\tilde{\zeta}, c\tilde{\zeta})(x - x_0)$ and the general boundary conditions write

$$ \begin{aligned} \xi_0^+(\zeta_0, q_0)(t) &= g_0(t) := \xi_0^+(\tilde{\zeta}(-x_0 - ct), c\tilde{\zeta}(-x_0 - ct)), \\ \xi_\ell^-(\zeta_\ell, q_\ell)(t) &= g_\ell(t) := \xi_\ell^-(\tilde{\zeta}(\ell - x_0 - ct), c\tilde{\zeta}(\ell - x_0 - ct)). \end{aligned} $$

In practice we take $\ell = 60$, $h_0 = 1$, $\zeta_{\text{max}} = h_0/5$, and $c > 0$. In order to validate the treatment of the boundary conditions we consider the two settings below.

**INCOMING SOLITARY WAVE** – The solitary wave is initially centered on $x_0 = -\ell/2$ outside of the computational domain $(0, \ell)$. It then enters the domain from the left boundary, and the simulation stops at time $t = \ell/c$ when the solitary wave is centered on $\ell/2$.

**OUTGOING SOLITARY WAVE** – The solitary wave is initially centered on $x_0 = \ell/2$ inside the computational domain, and the simulation is stopped once most of the wave left the domain through the right boundary at time $t = 3\ell/4c$.

Figure 1 compares the Lax-Friedrichs and MacCormack schemes at final time for both settings with incoming Riemann invariants enforced as boundary conditions. Unsurprisingly, the Lax-Friedrichs is much more diffusive than MacCormack, and also creates more reflections when trying to evacuate the wave through the right boundary.

Next we compare the three different choices for $\xi_0^+$, $\xi_\ell^-$ mentioned at the beginning of Section 5, and for each case we display the following $\ell^2$ error on the elevation

$$ (79) \qquad E_{\text{num}}^n = \sqrt{\frac{\Delta x}{\ell}} \left( \sum_{i=1}^N (\zeta_i^n - \tilde{\zeta}(x_i - x_0 - ct^n))^2 \right)^{1/2} $$

at final time $t = \ell/c$.

Tables 1 and 2 correspond respectively to the Lax-Friedrichs and MacCormack schemes for the incoming solitary wave test-case. Experimentally, both methods achieve the expected order of convergence, that is to say first order for Lax-Friedrichs and second order for MacCormack. The difference between the various boundary conditions in terms of the $\ell^2$ error is marginal, except when enforcing the discharge with the MacCormack scheme which leads to a substantially larger error than the two other choices.

Likewise, the outgoing solitary wave test-case is addressed in Table 3 for Lax-Friedrichs and Table 4 for MacCormack. The latter method shows good agreement with a second order convergence. For Lax-Friedrichs the situation is comparable to the incoming solitary wave test-case, except that a more refined
<figure><img src="image_4.png" /><figcaption>Figure 1: Free surface elevation for the incoming and outgoing solitary wave test-cases (respectively top and bottom); incoming Riemann invariants imposed at the boundaries</figcaption></figure>

<table><caption>Table 1: Lax-Friedrichs scheme for the incoming solitary wave with different boundary conditions. In the last column $R^{\pm}$ refers to the incoming Riemann invariant associated with the shallow water system.</caption><thead><tr><th rowspan="2">Δx</th><th colspan="2">ζ enforced</th><th colspan="2">q enforced</th><th colspan="2">R<sup>±</sup> enforced</th></tr><tr><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th></tr></thead><tbody><tr><td>8.82E-02</td><td>4.409E-03</td><td>—</td><td>4.728E-03</td><td>—</td><td>4.482E-03</td><td>—</td></tr><tr><td>6.25E-02</td><td>3.263E-03</td><td>0.87</td><td>3.484E-03</td><td>0.88</td><td>3.312E-03</td><td>0.87</td></tr><tr><td>4.42E-02</td><td>2.387E-03</td><td>0.90</td><td>2.539E-03</td><td>0.91</td><td>2.420E-03</td><td>0.91</td></tr><tr><td>3.12E-02</td><td>1.727E-03</td><td>0.93</td><td>1.833E-03</td><td>0.94</td><td>1.750E-03</td><td>0.93</td></tr><tr><td>2.21E-02</td><td>1.244E-03</td><td>0.95</td><td>1.318E-03</td><td>0.95</td><td>1.260E-03</td><td>0.95</td></tr><tr><td>1.56E-02</td><td>8.906E-04</td><td>0.97</td><td>9.411E-04</td><td>0.97</td><td>9.015E-04</td><td>0.97</td></tr></tbody></table>
<table><caption>Table 2: MacCormack scheme for the incoming solitary wave.</caption><thead><tr><th rowspan="2">Δx</th><th colspan="2">ζ enforced</th><th colspan="2">q enforced</th><th colspan="2">R<sup>±</sup> enforced</th></tr><tr><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th></tr></thead><tbody><tr><td>1.00E+00</td><td>7.846E-03</td><td>—</td><td>1.186E-02</td><td>—</td><td>9.424E-03</td><td>—</td></tr><tr><td>5.00E-01</td><td>2.473E-03</td><td>1.67</td><td>3.789E-03</td><td>1.65</td><td>2.783E-03</td><td>1.76</td></tr><tr><td>2.50E-01</td><td>6.670E-04</td><td>1.89</td><td>1.026E-03</td><td>1.89</td><td>7.199E-04</td><td>1.95</td></tr><tr><td>1.25E-01</td><td>1.696E-04</td><td>1.98</td><td>2.662E-04</td><td>1.95</td><td>1.832E-04</td><td>1.97</td></tr><tr><td>6.25E-02</td><td>4.312E-05</td><td>1.98</td><td>6.804E-05</td><td>1.97</td><td>4.624E-05</td><td>1.99</td></tr><tr><td>3.12E-02</td><td>1.107E-05</td><td>1.96</td><td>1.731E-05</td><td>1.97</td><td>1.162E-05</td><td>1.99</td></tr></tbody></table>

<table><caption>Table 3: Lax-Friedrichs scheme for the outgoing solitary wave.</caption><thead><tr><th rowspan="2">Δx</th><th colspan="2">ζ enforced</th><th colspan="2">q enforced</th><th colspan="2">R<sup>±</sup> enforced</th></tr><tr><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th></tr></thead><tbody><tr><td>6.25E-02</td><td>2.363E-03</td><td>—</td><td>2.666E-03</td><td>—</td><td>1.331E-04</td><td>—</td></tr><tr><td>4.42E-02</td><td>1.785E-03</td><td>0.81</td><td>2.031E-03</td><td>0.78</td><td>1.070E-04</td><td>0.63</td></tr><tr><td>3.12E-02</td><td>1.322E-03</td><td>0.87</td><td>1.514E-03</td><td>0.85</td><td>8.401E-05</td><td>0.70</td></tr><tr><td>2.21E-02</td><td>9.725E-04</td><td>0.89</td><td>1.119E-03</td><td>0.87</td><td>6.503E-05</td><td>0.74</td></tr><tr><td>1.56E-02</td><td>7.053E-04</td><td>0.93</td><td>8.145E-04</td><td>0.92</td><td>4.939E-05</td><td>0.79</td></tr><tr><td>1.11E-02</td><td>5.098E-04</td><td>0.94</td><td>5.903E-04</td><td>0.93</td><td>3.711E-05</td><td>0.83</td></tr></tbody></table>

<table><caption>Table 4: MacCormack scheme for the outgoing solitary wave.</caption><thead><tr><th rowspan="2">Δx</th><th colspan="2">ζ enforced</th><th colspan="2">q enforced</th><th colspan="2">R<sup>±</sup> enforced</th></tr><tr><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th></tr></thead><tbody><tr><td>1.00E+00</td><td>6.430E-03</td><td>—</td><td>1.031E-02</td><td>—</td><td>4.287E-03</td><td>—</td></tr><tr><td>5.00E-01</td><td>1.885E-03</td><td>1.77</td><td>4.437E-03</td><td>1.22</td><td>1.993E-03</td><td>1.11</td></tr><tr><td>2.50E-01</td><td>5.568E-04</td><td>1.76</td><td>1.257E-03</td><td>1.82</td><td>5.607E-04</td><td>1.83</td></tr><tr><td>1.25E-01</td><td>1.459E-04</td><td>1.93</td><td>3.252E-04</td><td>1.95</td><td>1.318E-04</td><td>2.09</td></tr><tr><td>6.25E-02</td><td>3.674E-05</td><td>1.99</td><td>8.223E-05</td><td>1.98</td><td>2.826E-05</td><td>2.22</td></tr><tr><td>3.12E-02</td><td>9.162E-06</td><td>2.00</td><td>2.067E-05</td><td>1.99</td><td>7.085E-06</td><td>2.00</td></tr></tbody></table>
mesh is required to start approaching a first order convergence. Note also that when enforcing the incoming Riemann invariants with the Lax-Friedrichs scheme, the order of convergence seems smaller than for the other boundary conditions, however we remark that in this case the errors differ by one order of magnitude in favor of the boundary conditions obtained by enforcing the incoming Riemann invariants.

### 5.1.2 The Boussinesq-Abbott system with varying bathymetry

Since no analytic expression is available for solutions of the Boussinesq-Abbott model (15) in the general case of unsteady flows over a non-flat bottom, we instead approximate a reference solution over a large domain $(-\ell, 2\ell)$ with periodic boundary conditions by the mean of a fine mesh; this approximated reference solution is then used to generate boundary conditions for the small domain $(0, \ell)$, in which the solution can be solved numerically with the proposed Lax-Friedrichs and MacCormack schemes.

The small domain length is set to $\ell = 25$, and we choose the following values for the characteristic depth, bathymetry and wave amplitudes:

$$h_0 = 1, \quad \beta = 1/4, \quad A = h_0/4.$$

The bathymetry profile is defined as follows

$$ (80) \qquad b(x) = \beta h_0 \cdot \begin{cases} 1 & \text{if } |x - \ell/2| < \frac{\ell}{4} \\ \sigma(1.5 - \frac{2}{\ell}|x - \ell/2|) & \text{if } \frac{\ell}{4} \le |x - \ell/2| \le \frac{3\ell}{4}, \\ 0 & \text{otherwise} \end{cases} $$

where $\sigma: [0, 1] \to [0, 1]$ is a smooth step function given here by $\sigma(x) = 6x^5 - 15x^4 + 10x^3$, which allows the bathymetry (80) to be $C^2$. The initial condition is taken as a Gaussian centered outside of the small domain:

$$ (81) \qquad \zeta^{\text{in}}(x) = A \cdot \exp(-(x-11)^2/3), \quad q^{\text{in}}(x) = 5 \cdot \zeta^{\text{in}}(x). $$

Recalling the shallowness parameter $\mu = (2\pi h_0)^2/\lambda^2$ with the characteristic wavelength defined here as $\lambda = \text{length}\{x, \zeta^{\text{in}}(x) > A/100\}$, one has $\mu \approx 0.6$ so that the flow is only moderately shallow. The nonlinearity parameter $A/h_0$ is $\epsilon = 0.25$.

The periodic reference solution is approximated over the large domain $(-\ell, 2\ell)$ with a spatial step $\Delta x \approx 8.14 \cdot 10^{-3}$, and the simulation stops at time $t = 15$; see Figure 2 for plots of the initial and final states. Note that the shape of the initial free surface elevation (81) is not preserved through time. More specifically, the Gaussian splits into right- and left-going signals; owing to the periodic condition, both of these signals will enter the small domain $(0, \ell)$ respectively from the left and right boundaries, moreover the right-going signal begins to exit the small domain before the simulation ends. Therefore this can be considered a rich test-case featuring wave generation, interaction and evacuation over a non-flat bottom $b$ whose derivative doesn't vanish at boundaries $x = 0, \ell$.
<figure><img src="image_5.png" /><figcaption>Figure 2: Left: initial condition. Right: reference solution, Lax-Friedrichs and MacCormack approximations at time $t = 15$, both obtained using the same CFL constant of 0.45.</figcaption></figure>

<table><caption>Table 5: Lax-Friedrichs scheme for the Gaussian over bump test-case.</caption><thead><tr><th rowspan="2">Δx</th><th colspan="2">ζ enforced</th><th colspan="2">q enforced</th><th colspan="2">R<sup>±</sup> enforced</th></tr><tr><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th></tr></thead><tbody><tr><td>3.26E-02</td><td>1.017E-02</td><td>—</td><td>1.300E-02</td><td>—</td><td>8.545E-03</td><td>—</td></tr><tr><td>2.30E-02</td><td>7.862E-03</td><td>0.74</td><td>1.016E-02</td><td>0.71</td><td>6.604E-03</td><td>0.74</td></tr><tr><td>1.63E-02</td><td>5.962E-03</td><td>0.80</td><td>7.759E-03</td><td>0.78</td><td>5.003E-03</td><td>0.80</td></tr><tr><td>1.15E-02</td><td>4.440E-03</td><td>0.85</td><td>5.802E-03</td><td>0.84</td><td>3.721E-03</td><td>0.85</td></tr><tr><td>8.14E-03</td><td>3.265E-03</td><td>0.89</td><td>4.278E-03</td><td>0.88</td><td>2.733E-03</td><td>0.89</td></tr><tr><td>5.75E-03</td><td>2.375E-03</td><td>0.92</td><td>3.118E-03</td><td>0.91</td><td>1.986E-03</td><td>0.92</td></tr></tbody></table>

<table><caption>Table 6: MacCormack scheme for the Gaussian over bump test-case.</caption><thead><tr><th rowspan="2">Δx</th><th colspan="2">ζ enforced</th><th colspan="2">q enforced</th><th colspan="2">R<sup>±</sup> enforced</th></tr><tr><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th><th>L<sup>2</sup>-error</th><th>Order</th></tr></thead><tbody><tr><td>2.60E-01</td><td>4.635E-03</td><td>—</td><td>5.708E-03</td><td>—</td><td>3.301E-03</td><td>—</td></tr><tr><td>1.84E-01</td><td>2.470E-03</td><td>1.81</td><td>3.076E-03</td><td>1.78</td><td>1.690E-03</td><td>1.92</td></tr><tr><td>1.30E-01</td><td>1.288E-03</td><td>1.89</td><td>1.629E-03</td><td>1.84</td><td>8.485E-04</td><td>2.00</td></tr><tr><td>9.19E-02</td><td>6.575E-04</td><td>1.93</td><td>8.402E-04</td><td>1.90</td><td>4.095E-04</td><td>2.09</td></tr><tr><td>6.51E-02</td><td>3.384E-04</td><td>1.93</td><td>4.316E-04</td><td>1.93</td><td>1.930E-04</td><td>2.18</td></tr><tr><td>4.60E-02</td><td>1.777E-04</td><td>1.85</td><td>2.211E-04</td><td>1.92</td><td>9.056E-05</td><td>2.17</td></tr></tbody></table>
Comparing the numerical solution obtained in the small domain with different mesh sizes to the reference solution, we can compute the $\ell^2$ error (79) and the associated experimental orders of convergence are shown in Tables 5 and 6. Both the Lax-Friedrichs and MacCormack schemes achieve the expected order for the different boundary conditions considered. The most advantageous choice seems to enforce the incoming Riemann invariants, as it leads to the smallest error.

## 5.2 Asymptotic stability

Finally we wish to assess the feasibility of achieving a reliable approximation of the solution to the Boussinesq-Abbott system (15) over a varying bathymetry when starting the simulation with a wrong initial condition. This is motivated by the fact that in real-life, measurements of the free surface elevation can only be performed at a limited number of points, preventing us from knowing the state of the flow everywhere in the domain of interest. Instead, one has to provide an initial guess, which might not be accurate at all.

To this end, we use a setup similar to the one from Section 5.1.2: a periodic reference solution computed on the large domain $(-\ell, 2\ell)$ is used to generate boundary conditions for the small domain $(0, \ell)$. We keep the same domain length $\ell = 25$ and the same definition (80) for the bottom, however the characteristic depth and wave amplitude are now scaled with respect to a parameter $K_\mu$ as follows:

$$h_0 = \sqrt{K_\mu}, \quad A = \sqrt{K_\mu}/4.$$

The initial data of the reference solution consists in a sine wave of period $\lambda = 3\ell/5$ and amplitude $A$ for the elevation, together with a nonzero discharge:

$$(82) \qquad \tilde{\zeta}^{\text{in}}(x) = A \cdot \sin(2\pi x/\lambda), \qquad \tilde{q}^{\text{in}}(x) = 2 \cdot \tilde{\zeta}^{\text{in}}(x).$$

The factor $K_\mu$ shall be taken equal to $(2\pi/\lambda)^{-2}\mu$, with $\mu$ the shallowness parameter valued in $\{0.01, 0.1, 1\}$. The nonlinearity parameter $\epsilon$ and bathymetry parameter $\beta$ are kept to 0.25. We also introduce the characteristic time $T_0$ given by $2\pi T_0 = \lambda/\sqrt{gh_0}$.

When approximating the reference solution over the small domain $(0, \ell)$ with the proposed MacCormack scheme, the simulation is initialized with a lake at rest $(\zeta^{\text{in}}, q^{\text{in}}) = (0, 0)$ which is different from the initial state of the reference solution. This initial setup is plotted in Figure 3 for $\mu = 1$. In order to comply with the system (48), which we recall provides compatibility conditions when enforcing $\zeta$, the reference solution $(\tilde{\zeta}, \tilde{q})$ is progressively enforced at the boundaries of the small domain by the mean of the smooth step function $\sigma$ defined in Section 5.1.2:

$$\xi_0^+(\zeta_0, q_0)(t) = g_0(t) := \sigma(5t/T_0) \cdot \xi_0^+(\tilde{\zeta}(t, 0), \tilde{q}(t, 0)),$$

$$\xi_\ell^-(\zeta_\ell, q_\ell)(t) = g_\ell(t) := \sigma(5t/T_0) \cdot \xi_\ell^-(\tilde{\zeta}(t, \ell), \tilde{q}(t, \ell)).$$

If on the contrary the boundary conditions on the elevation were enforced directly at $t = 0$, then system (48) would be violated, and the well-posedness result from Theorem 3.8 does not hold anymore in this situation.
<figure><img src="image_6.png" /><figcaption>Figure 3: Initial elevation for reference and small domain solutions ($\mu = 1$).</figcaption></figure>

We plot the results obtained at final time $t = 50 T_0$ in Figure 4 for $\mu = 1$ and $\mu = 0.1$, and in Figure 5 for $\mu = 0.01$. Out of the three boundary conditions tested, only the choice of enforcing the incoming Riemann invariants seems to allow the $\ell^2$ error (79) to decay to zero, which is consistent with the considerations of Section 3.5. On the other hand, this is not the case when enforcing the elevation or the discharge; our interpretation is that the initial information — which represents the perturbation of the reference solution — is prevented from leaving the domain due to reflections. Reducing the shallowness parameter to 0.01, one eventually observes a blow up of the numerical approximation for these two boundary conditions, and which is avoided with the incoming Riemann invariants. As a conclusion, the latter should clearly be preferred over the elevation or discharge in a situation where waves need to be evacuated from the computational domain.

**Remark 5.1.** In practice, enforcing the incoming Riemann invariant at only one of the two boundaries is enough to recover a good approximation of the reference solution. In this case the $\ell^2$ error on the elevation decays slightly slower towards zero compared to when the incoming Riemann invariants are enforced at both boundaries.
<figure><img src="image_7.png" /><figcaption>Figure 4: Flow over bar test-case for $\mu = 1$ (first column) and $\mu = 10^{-1}$ (second column). First and second rows correspond respectively to the MacCormack elevation and its deviation from the reference elevation $\tilde{\zeta}$ at final time $t = 50 T_0$ for various boundary conditions. The last row represents the $\ell^2$ error (79) over time. Out of the three boundary conditions investigated here, only enforcing the incoming Riemann invariants allows this error to decay towards zero.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 5: Flow over bar test-case for $\mu = 10^{-2}$. The first plot was obtained at final time $t = 50 T_0$. Enforcing $\zeta$ or $q$ as boundary conditions led to a blow up of the MacCormack approximation, therefore only the case of incoming Riemann invariants is displayed here. The $\ell^2$ error (79) decays over time but does not quite vanish; this is most likely due to the coarseness of the mesh which make it difficult to accurately capture the high frequencies that develop in the solution. A fix would be to refine the mesh further.</figcaption></figure>
# 6 Conclusion and perspectives

In this article, the initial boundary value problem linked to the Boussinesq-Abbott model with non flat bottom has been considered. When imposing the values of the surface elevation at the boundaries, we have seen that like in the case of a flat bottom over a half-line treated in [28], this problem can be equivalently reformulated as an initial value problem involving a nonlocal flux and an ODE on the trace of the discharge, even in presence of a varying topography and over a bounded domain. This formulation is however more complex because some commutation properties used in the case of a flat bottom are no longer true. We then proposed a new method allowing to prescribe a nonlinear function of the unknowns at the endpoints of the domain by a given function of time. This method is based on the notion of *input* functions, which represent the information that one wants to let in the domain, and that of *output* functions. We believe that the ability to enforce such general boundary conditions is quite relevant for complex applications such as the ones encountered in oceanography. Importantly, our reformulated model defines an ODE on some functional space; under regularity assumptions on the boundary data and on the initial condition, it is possible to obtain the well-posedness of this model in finite time by the Cauchy-Lipschitz theorem. We do not believe however that this approach can be transposed to the Serre-Green-Naghdi model, for which the infinite dimensional ODE structure of the Boussinesq-Abbott system is lost.

This theoretical study was then followed by a numerical discretization of the equations using a hybrid finite difference/finite volumes approach. New schemes of first and second orders were designed using respectively the Lax-Friedrichs numerical flux and the MacCormack strategy, similarly to what was proposed in [8] for wave-structure interactions over a flat bottom and for boundary conditions on the horizontal discharge. Owing to a simple adaptation of the equations, these schemes can be rendered well-balanced, that is to say that they preserve the lake at rest equilibrium. Next we performed numerical experiments validating the expected orders of convergence for various types of boundary conditions, in a complex setup featuring wave generation and evacuation over a varying topography. Finally we explored numerically the question of asymptotic stability, which is to know whether solutions arising from different initial data but with the same boundary conditions converge to each other after a transitory regime. We found that this is the case when enforcing the incoming Riemann invariants of the shallow water equations, but not when enforcing the elevation nor the discharge. Again this underlines the interest there is to be able to enforce general boundary conditions.

We consider several perspectives to the present work. The first one is related to a limitation inherent to dispersive models in general (including Boussinesq and Serre-Green-Naghdi), that is their inability to describe the breaking of waves; in fact for steep variations of the elevation it is well known that these models produce non physical spurious oscillations. Despite being less accurate than these dispersive models, the shallow water equations present the benefit of being
able to handle breaking waves in a physically relevant way by dissipating the energy in the form of shocks, see [13] for a physical justification and [32] for a general methodology on Riemann problems. Note also that the nonlinear shallow water equations are able to take into account the vanishing of the water height $h$ (see [27] for the mathematical analysis of this singularity and [6] for its numerical treatment); therefore this model is pertinent to describe the wave dynamics in the surf and swash areas. For these reasons, several authors propose to use a dispersive model and to switch the dispersive terms off in the vicinity of wave breaking [31, 14, 23], hereby working locally with the less precise but more robust nonlinear shallow water equations. The coupling between the nonlinear shallow water equations and dispersive models like the Boussinesq and Serre-Green-Naghdi equation is however not well understood essentially because the boundary conditions that must be imposed on the dispersive component is unclear. A natural perspective of this work is therefore to investigate this coupling.

As pointed in [19], the shoaling phenomenon is systematically underestimated (resp. overestimated) by Boussinesq-type models written in elevation-discharge form (resp. in elevation-velocity form). Preliminary tests show that our reformulated Boussinesq-Abbott model is no different; it would thus be interesting to try circumventing this issue. An idea that we would like to explore would consist to perform some form of averaging between the Boussinesq-Abbott model in $(\zeta, q)$ coordinates and the classical Boussinesq-Peregrine model in $(\zeta, v)$ coordinates.

To conclude, we mention a long term goal related to the statistical study of extreme waves. This would require to generate random waves at the boundaries of the domain so as to force a realistic wave field in the domain. We presume that the use of asymptotically stable boundary conditions is crucial in this setting.
References

[1] M. Abbott, H. Petersen, and O. Skovgaard. "Computations of short waves in shallow water". In: *Coastal Engineering Proceedings* 1.16 (Jan. 1978), p. 23.

[2] D. Antonopoulos and V. Dougalis. "Galerkin-finite element methods for the shallow water equations with characteristic boundary conditions". In: *IMA Journal of Numerical Analysis* 37.1 (2017), pp. 266-295.

[3] D. Antonopoulos and V. Dougalis. "Numerical solution of the 'classical'Boussinesq system". In: *Mathematics and Computers in Simulation* 82.6 (2012), pp. 984-1007.

[4] D. Antonopoulos, V. Dougalis, and D. Mitsotakis. "Initial-boundary-value problems for the Bona-Smith family of Boussinesq systems". In: *Advances in Differential Equations* 14.1/2 (2009), pp. 27-53.

[5] D. Antonopoulos, V. Dougalis, and D. Mitsotakis. "Numerical solution of Boussinesq systems of the Bona-Smith family". In: *Applied numerical mathematics* 60.4 (2010), pp. 314-336.

[6] E. Audusse and M.-O. Bristeau. "A well-balanced positivity preserving "second-order" scheme for shallow water flows on unstructured meshes". In: *Journal of Computational Physics* 206.1 (2005), pp. 311-333.

[7] G. Beck and D. Lannes. "Freely floating objects on a fluid governed by the Boussinesq equations". In: *Annales de l'Institut Henri Poincaré C* (2022).

[8] G. Beck, D. Lannes, and L. Weynans. "A numerical method for wave-structure interactions in the Boussinesq regime". In: *arXiv preprint* 2307.01749 (2023).

[9] S. Benzoni-Gavage and D. Serre. *Multi-dimensional hyperbolic partial differential equations: First-order Systems and Applications*. Oxford University Press on Demand, 2007.

[10] C. Besse, B. Mésognon-Gireau, and P. Noble. "Artificial boundary conditions for the linearized Benjamin-Bona-Mahony equation". In: *Numerische Mathematik* 139.2 (2018), pp. 281-314.

[11] C. Besse, P. Noble, and D. Sanchez. "Discrete transparent boundary conditions for the mixed KDV-BBM equation". In: *Journal of Computational Physics* 345 (2017), pp. 484-509.

[12] C. Besse et al. "Perfectly Matched Layers Methods for Mixed Hyperbolic-Dispersive Equations". In: *Water Waves* 4.3 (2022), pp. 313-343.

[13] P. Bonneton. "Modelling of periodic wave transformation in the inner surf zone". In: *Ocean Engineering* 34.10 (2007), pp. 1459-1471. ISSN: 0029-8018. DOI: https://doi.org/10.1016/j.oceaneng.2006.09.002. URL: https://www.sciencedirect.com/science/article/pii/S0029801806002459.
[14] P. Bonneton et al. “A splitting approach for the fully nonlinear and weakly dispersive Green–Naghdi model”. In: *J. Comput. Phys* 230 (2011), pp. 1479–1498.

[15] D. Bresch, D. Lannes, and G. Métivier. “Waves interacting with a partially immersed obstacle in the Boussinesq regime”. In: *Anal. PDE* 14.4 (2021), pp. 1085–1124.

[16] J.-M. Coron and G. Bastin. “Dissipative boundary conditions for one-dimensional quasi-linear hyperbolic systems: Lyapunov stability for the $C^1$-norm”. In: *SIAM Journal on Control and Optimization* 53.3 (2015), pp. 1464–1483.

[17] V. Dougalis, D. Mitsotakis, and J.-C. Saut. “Boussinesq systems of Bona-Smith type on plane domains: theory and numerical analysis”. In: *Journal of Scientific Computing* 44.2 (2010), pp. 109–135.

[18] V. Dougalis, D. Mitsotakis, and J.-C. Saut. “On initial-boundary value problems for a Boussinesq system of BBM-BBM type in a plane domain”. In: *Discrete Contin. Dyn. Syst* 23.4 (2009), pp. 1191–1204.

[19] 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.

[20] T. Iguchi and D. Lannes. “Hyperbolic free boundary problems and applications to wave-structure interactions”. In: *Indiana University Mathematics Journal* 70.1 (2021), pp. 353–464.

[21] M. Kazakova. “Dispersive models of ocean waves propagation: Numerical issues and modelling”. PhD thesis. Université Toulouse 3 Paul Sabatier, 2018.

[22] M. Kazakova and P. Noble. “Discrete transparent boundary conditions for the linearized Green–Naghdi system of equations”. In: *SIAM Journal on Numerical Analysis* 58.1 (2020), pp. 657–683.

[23] M. Kazolea and M. Ricchiuto. “On wave breaking for Boussinesq-type models”. In: *Ocean Modelling* 123 (2018), pp. 16–39.

[24] G. Kounadis, D. Antonopoulos, and V. Dougalis. “Galerkin finite element methods for the numerical solution of two classical-Boussinesq type systems over variable bottom topography”. In: *Wave Motion* 102 (2021), p. 102715.

[25] D. Lannes. “Modeling shallow water waves”. In: *Nonlinearity* 33.5 (2020), R1.

[26] D. Lannes. *The Water Waves Problem: Mathematical Analysis and Asymptotics*. Vol. 188. Mathematical Surveys and Monographs. AMS, 2013.

[27] D. Lannes and G. Métivier. “The shoreline problem for the one-dimensional shallow water and Green-Naghdi equations”. In: *J. Éc. polytech. Math.* 5 (2018), pp. 455–518.

[28] D. Lannes and L. Weynans. “Generating boundary conditions for a Boussinesq system”. In: *Nonlinearity* 33.12 (Oct. 2020), pp. 6868–6889.
[29] S. Noelle, M. Parisot, and T. Tscherpel. "A Class of Boundary Conditions for Time-Discrete Green-Naghdi Equations with Bathymetry". In: *SIAM Journal on Numerical Analysis* 60.5 (2022), pp. 2681-2712.

[30] J. Steinstraesser, G. Kemlin, and A. Rousseau. "A domain decomposition method for linearized Boussinesq-type equations". In: *Journal of Mathematical Study* 52.3 (2019), pp. 320-340.

[31] M. Tonelli and M. Petti. "Shock-capturing Boussinesq model for irregular wave propagation". In: *Coastal Engineering* 61 (2012), pp. 8-19.

[32] E.F. Toro. *Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction*. Springer Berlin Heidelberg, 2013.

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