# General boundary conditions for a Boussinesq model with varying bathymetry

David Lannes  Institut de Mathématiques de Bordeaux, Université de Bordeaux et
CNRS UMR 5251, 351 Cours de la libération 33405 Talence Cedex, France  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,

$\left\\{\begin{array}[]{l}\partial_{t}\zeta+\partial_{x}q=0\\\
(1+h_{b}{\mathcal{T}}_{b})\partial_{t}q+\partial_{x}(\frac{1}{h}q^{2})+{\mathtt{g}}h\partial_{x}\zeta=0\end{array}\right.\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; ${\mathtt{g}}$ is the acceleration
of gravity; and finally ${\mathcal{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
$\overline{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,

$\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
\partial_{t}q+\partial_{x}(\frac{1}{h}q^{2})+{\mathtt{g}}h\partial_{x}\zeta=0\end{cases}\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 $\zeta$). All these references considered the case of
the Boussinesq-Abbott model with a flat topography, and with boundary
conditions imposed on $\zeta$ 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=\ell$. 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 $\zeta$ 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 $\zeta$ 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 an 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 ${\mathcal{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

(3)
$\displaystyle\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\end{array}\right.\quad\text{in}\
(0,\ell)\ ,$

where we recall that

$f_{\text{NSW}}=f_{\rm
NSW}(\zeta,q)=\frac{1}{2}{\mathtt{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

(6) $\displaystyle R^{0}:f\in L^{2}(0,\ell)\longmapsto u\in
H^{2}(0,\ell)\quad\text{solving}\quad\left\\{\begin{array}[]{l}\big{[}1-\frac{h_{0}^{2}}{3}\partial_{x}^{2}\big{]}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,

(7)
$\left\\{\begin{array}[]{l}\big{[}1-\frac{h_{0}^{2}}{3}\partial_{x}^{2}\big{]}{\mathfrak{s}}^{(0)}=0\\\
{\mathfrak{s}}^{(0)}(0)=1,\qquad{\mathfrak{s}}^{(0)}(\ell)=0\end{array}\right.\quad\mbox{
and
}\quad\left\\{\begin{array}[]{l}\big{[}1-\frac{h_{0}^{2}}{3}\partial_{x}^{2}\big{]}{\mathfrak{s}}^{(\ell)}=0\\\
{\mathfrak{s}}^{(\ell)}(0)=0,\qquad{\mathfrak{s}}^{(\ell)}(\ell)=1\end{array}\right..$

We also denote by ${\mathfrak{S}}^{\prime}$ the matrix

(8)
${\mathfrak{S}}^{\prime}=\begin{pmatrix}({\mathfrak{s}}^{(0)})^{\prime}(0)&({\mathfrak{s}}^{(\ell)})^{\prime}(0)\\\
({\mathfrak{s}}^{(0)})^{\prime}(\ell)&({\mathfrak{s}}^{(\ell)})^{\prime}(\ell)\end{pmatrix}.$

######  Remark 2.1.

One can compute ${\mathfrak{s}}^{(0)}$ and ${\mathfrak{s}}^{(\ell)}$
explicitly, namely,

(9)
${\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}}^{\prime}$ 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$;

(12) $\displaystyle R^{1}:f\in L^{2}(0,\ell)\longmapsto u\in
H^{2}(0,\ell)\quad\text{solving}\quad\left\\{\begin{array}[]{l}\big{[}1-\frac{h_{0}^{2}}{3}\partial_{x}^{2}\big{]}u=f\\\
\partial_{x}u(0)=\partial_{x}u(\ell)=0.\end{array}\right.\ .$

An important observation is the following commutation property,

(13) $\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 (3) with
boundary conditions

(14) $\zeta(\cdot,0)=g_{0}\quad\mbox{ and }\quad\zeta(\cdot,\ell)=g_{\ell},$

and initial condition

(15) $(\zeta,q)(0,\cdot)=(\zeta^{\rm in},q^{\rm in})$

is well posed. Remarking that from the first equation of (23) one has
$\frac{{\rm d}}{{\rm d}t}(\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

(16) $\begin{cases}\zeta^{\rm in}(0)=g_{0}(0),\\\ -\partial_{x}q^{\rm
in}(0)=\dot{g}_{0}(0),\end{cases}\quad\mbox{ and }\quad\begin{cases}\zeta^{\rm
in}(\ell)=g_{\ell}(0),\\\ -\partial_{x}q^{\rm
in}(\ell)=\dot{g}_{\ell}(0).\end{cases}$

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^{\rm in},q^{\rm in})$ and
$(g_{0},g_{\ell})$ satisfy the compatibility condition (16). Then the two
following assertions are equivalent:
i. The couple $(\zeta,q)$ is a regular solution to (3) such that the depth $h$
never vanishes and with boundary conditions (14) and initial condition (15).
ii. The couple $(\zeta,q)$ is a regular solution such that the depth $h$ never
vanishes to

(17) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
\partial_{t}q+\partial_{x}(R_{1}f_{\rm
NSW})=\dot{q}_{0}{\mathfrak{s}}^{(0)}+\dot{q}_{\ell}{\mathfrak{s}}^{(\ell)},\end{cases}\quad\text{in}\
(0,\ell),$

with initial condition (15), and where $q_{0}$ and $q_{\ell}$ solve the ODE

(18) ${\mathfrak{S}}^{\prime}\begin{pmatrix}\dot{q}_{0}\\\
\dot{q}_{\ell}\end{pmatrix}+\frac{3}{h_{0}^{2}}\begin{pmatrix}f_{\rm
NSW}(g_{0},q_{0})\\\ f_{\rm
NSW}(g_{\ell},q_{\ell})\end{pmatrix}=-\begin{pmatrix}\ddot{g}_{0}\\\
\ddot{g}_{\ell}\end{pmatrix}+\frac{3}{h_{0}^{2}}\begin{pmatrix}(R_{1}f_{\rm
NSW})_{|_{x=0}}\\\ (R_{1}f_{\rm NSW})_{|_{x=\ell}}\end{pmatrix},$

with ${\mathfrak{S}}^{\prime}$ defined in (8), and with initial condition

(19) $q_{0}(0)=q^{\rm in}(0)\quad\mbox{ and }\quad q_{\ell}(0)=q^{\rm
in}(\ell).$

###### Proof.

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

$\partial_{t}q+R_{0}\partial_{x}f_{\rm
NSW}=\dot{q}_{0}{\mathfrak{s}}^{(0)}+\dot{q}_{\ell}{\mathfrak{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_{\rm
NSW}=\dot{q}_{0}{\mathfrak{s}}^{\prime(0)}+\dot{q}_{\ell}{\mathfrak{s}}^{\prime(\ell)}.$

Using the first equation one can replace
$\partial_{t}\partial_{x}q=-\partial^{2}_{t}\zeta$ while the identities
$R_{0}\partial_{x}=\partial_{x}R_{1}$ and
$(1-\frac{h_{0}^{2}}{3}\partial_{x}^{2})R_{1}={\rm Id}$ allow us to deduce
that

$-\partial_{t}^{2}\zeta-\frac{3}{h_{0}^{2}}(1-R_{1})f_{\rm
NSW}=\dot{q}_{0}{\mathfrak{s}}^{\prime}_{(0)}+\dot{q}_{\ell}{\mathfrak{s}}^{\prime}_{(\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}{\mathfrak{s}}^{\prime}_{(0)}(0)\dot{q}_{0}+{\mathfrak{s}}^{\prime}_{(\ell)}(0)\dot{q}_{\ell}&=-\frac{3}{h_{0}^{2}}f_{\rm
NSW}(g_{0},q_{0})+\frac{3}{h_{0}^{2}}(R_{1}f_{\rm
NSW})_{|_{x=0}}-\ddot{g}_{0}\\\
{\mathfrak{s}}^{\prime}_{(0)}(\ell)\dot{q}_{0}+{\mathfrak{s}}^{\prime}_{(\ell)}(\ell)\dot{q}_{\ell}&=-\frac{3}{h_{0}^{2}}f_{\rm
NSW}(g_{\ell},q_{\ell})+\frac{3}{h_{0}^{2}}(R_{1}f_{\rm
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_{l})$ solve the above system then $(\zeta,q)$ solves the Boussinesq
system (3). By taking the trace of the second equation of (17) at $x=0$ and
$x=\ell$ respectively, we also obtain that $\frac{{\rm d}}{{\rm
d}t}q(t,0)=\dot{q}_{0}(t)$ and $\frac{{\rm d}}{{\rm
d}t}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 (14) are not necessarily satisfied. However, as for the
derivation of (18), we obtain that
$(\zeta_{0}(t),\zeta_{\ell}(t)):=(\zeta(t,0),\zeta(t,\ell))$ satisfies

${\mathfrak{S}}^{\prime}\begin{pmatrix}\dot{q}_{0}\\\
\dot{q}_{\ell}\end{pmatrix}+\frac{3}{h_{0}^{2}}\begin{pmatrix}f_{\rm
NSW}(\zeta_{0},q_{0})\\\ f_{\rm
NSW}(\zeta_{\ell},q_{\ell})\end{pmatrix}=-\begin{pmatrix}\ddot{\zeta}_{0}\\\
\ddot{\zeta}_{\ell}\end{pmatrix}+\frac{3}{h_{0}^{2}}\begin{pmatrix}(R_{1}f_{\rm
NSW})_{|_{x=0}}\\\ (R_{1}f_{\rm NSW})_{|_{x=\ell}}\end{pmatrix},$

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\mbox{ 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 (16). It follows that the boundary conditions (14) are also
satisfied. This concludes the proof of the proposition. ∎

The initial value problem formed by (17) and (18) with initial conditions (15)
and (19) 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^{\rm in})>0$ ensures that the shallow water flux $f_{\rm
NSW}$ is well defined. By Proposition 2.2, this solution furnishes a solution
to the original initial boundary value problem (3), (14) and (15).

######  Proposition 2.3.

Let $g_{0},g_{\ell}\in C^{\infty}({\mathbb{R}}^{+})$. Let also
$n\in{\mathbb{N}}\backslash\\{0\\}$ and $(\zeta^{\rm in},q^{\rm in})\in
H^{n}(0,\ell)\times H^{n+1}(0,\ell)$ be such that $\inf(h_{0}+\zeta^{\rm
in})>0$ and satisfying the compatibility conditions (16).
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 (17) and (18) with initial
conditions (15) and (19), 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 (3). 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,

(20) $\inf_{(0,\ell)}(h_{0}-b)>0.$

Denoting $h_{b}=h_{0}-b$, the Boussinesq-Abbott model is given by

(23)
$\displaystyle\left\\{\begin{array}[]{l}\partial_{t}\zeta+\partial_{x}q=0\\\
(1+h_{b}{\mathcal{T}}_{b})\partial_{t}q+\partial_{x}f_{\text{NSW}}=-{\mathtt{g}}h\partial_{x}b\end{array}\right.\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_{\rm
NSW}(\zeta,q)=\frac{1}{2}{\mathtt{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 $\mathcal{T}_{b}$ is
given by

(24)
$\displaystyle{\mathcal{T}}_{b}(\cdot)=-\frac{1}{3h_{b}}\partial_{x}\Big{(}h_{b}^{3}\partial_{x}\frac{(\cdot)}{h_{b}}\Big{)}+\frac{1}{2}\partial_{x}^{2}b$

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

${\mathcal{T}}_{b}(\cdot)=S^{*}\big{(}h_{b}S(\cdot)\big{)}+\frac{1}{4}\frac{(\partial_{x}b)^{2}}{h_{b}}$

with

(25)
$S(\cdot)=-\frac{1}{\sqrt{3}}h_{b}\partial_{x}(\frac{1}{h_{b}}\cdot)+\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}{\mathcal{T}}_{b})$ with homogeneous Dirichlet boundary
conditions, that is to say

(28) $\displaystyle R^{0}_{b}:f\in L^{2}(0,\ell)\longmapsto u\in
H^{2}(0,\ell)\quad\text{solving}\quad\left\\{\begin{array}[]{l}\big{[}1+h_{b}{\mathcal{T}}_{b}\big{]}u=f\\\
u(0)=u(\ell)=0\end{array}\right.;$

here again, the operator $R^{0}_{b}$ 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,

(29)
$\left\\{\begin{array}[]{l}\big{[}1+h_{b}{\mathcal{T}}_{b}\big{]}{\mathfrak{s}}^{(b,0)}=0\\\
{\mathfrak{s}}^{(b,0)}(0)=1,\qquad{\mathfrak{s}}^{(b,0)}(\ell)=0\end{array}\right.\quad\mbox{
and
}\quad\left\\{\begin{array}[]{l}\big{[}1+h_{b}{\mathcal{T}}_{b}\big{]}{\mathfrak{s}}^{(b,\ell)}=0\\\
{\mathfrak{s}}^{(b,\ell)}(0)=0,\qquad{\mathfrak{s}}^{(b,\ell)}(\ell)=1\end{array}\right.$

and we also denote by ${\mathfrak{S}}_{b}^{\prime}$ the matrix

(30)
${\mathfrak{S}}_{b}^{\prime}=\begin{pmatrix}({\mathfrak{s}}^{(b,0)})^{\prime}(0)&({\mathfrak{s}}^{(b,\ell)})^{\prime}(0)\\\
({\mathfrak{s}}^{(b,0)})^{\prime}(\ell)&({\mathfrak{s}}^{(b,\ell)})^{\prime}(\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
${\mathfrak{S}}_{b}^{\prime}$ 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 (20) is satisfied and that in addition

(31)
$\min_{(0,\ell)}(1+\frac{1}{4}(\partial_{x}b)^{2}+\frac{1}{6}h_{b}\partial_{x}^{2}b)>0.$

Then the matrix ${\mathfrak{S}}_{b}^{\prime}$ is invertible and
$({\mathfrak{s}}^{(b,0)})^{\prime}(0)\neq 0$ and
$({\mathfrak{s}}^{(b,\ell)})^{\prime}(\ell)\neq 0$.

###### Proof.

We can remark that

$1+h_{b}{\mathcal{T}}_{b}=-\frac{1}{3}\partial_{x}(h_{b}^{2}\partial_{x}\cdot)-\frac{1}{3}h_{b}(\partial_{x}b)\partial_{x}+(1+\frac{1}{3}(\partial_{x}b)^{2}+\frac{1}{6}h_{b}\partial_{x}^{2}b).$

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}{\mathcal{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}(1+\frac{1}{3}(\partial_{x}b)^{2}+\frac{1}{6}h_{b}\partial_{x}^{2}b)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}{\mathcal{T}}_{b})u\geq\frac{1}{3}(1-\frac{\epsilon}{2})\int_{0}^{\ell}h_{b}^{2}(\partial_{x}u)^{2}+\int_{0}^{\ell}(1+\frac{1}{3}(1-\frac{1}{2\epsilon})(\partial_{x}b)^{2}+\frac{1}{6}h_{b}\partial_{x}^{2}b)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 ${\mathfrak{S}}_{b}^{\prime}$ 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}}^{\prime}(0)={\mathfrak{s}}^{\prime}(\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)})^{\prime}(0)\neq 0$, we can
modify the arguments above to get that if $u\in H^{1}(0,\ell)$ is such that
$u^{\prime}(0)=0$ and $u(\ell)=0$ then $u\equiv 0$. If we had
$({\mathfrak{s}}^{(b,0)})^{\prime}(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)})^{\prime}(\ell)\neq 0$, which
concludes the proof. ∎

Another important difference with the case of flat topography is that the
commutation property (13) is no longer true, that is, in general we have
$R^{0}_{b}\partial_{x}\neq R^{1}_{b}\partial_{x}$ if we define $R^{1}_{b}$ as
the inverse of $(1+h_{b}{\mathcal{T}}_{b})$ with homogeneous Neumann boundary
conditions. We can however define $R^{1}_{b}$ as follows

$\displaystyle R^{1}_{b}:f\in L^{2}(0,\ell)\longmapsto u\in H^{2}(0,\ell)$

where $u$ solves

(34)
$\displaystyle\left\\{\begin{array}[]{l}\big{[}1+\frac{1}{h_{b}}(h_{b}S)\big{(}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}\big{)}\big{]}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^{0}_{b}\partial_{x}f=\sqrt{3}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}(\frac{1}{h_{b}^{2}}f))-\frac{3}{2}R^{0}_{b}\big{(}\frac{\partial_{x}b}{h_{b}}f\big{)},$

where $\alpha_{b}=1+\frac{1}{4}(\partial_{x}b)^{2}$, and where we recall that
$S(\cdot)$ is defined in (25).

######  Remark 2.6.

When $b\equiv 0$ then $h_{b}=h_{0}$, $R^{0}_{b}=R^{0}$, $R^{1}_{b}=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

(35)
$1+h_{b}{\mathcal{T}}_{b}=\alpha_{b}+h_{b}S^{*}\big{(}h_{b}S(\cdot)\big{)};$

it follows that the equation $(1+h_{b}{\mathcal{T}}_{b})u=\partial_{x}f$ can
be equivalently written under the form

$\displaystyle\big{[}\frac{\alpha_{b}}{h_{b}}+S^{*}\big{(}h_{b}S(\cdot)\big{)}\big{]}u$
$\displaystyle=\frac{1}{h_{b}}\partial_{x}f$ (36)
$\displaystyle=S^{*}(\frac{\sqrt{3}}{h_{b}}f)-\frac{3}{2}\frac{1}{h_{b}^{2}}(\partial_{x}b)f.$

We now need the following lemma.

######  Lemma 2.7.

For all $\widetilde{f}\in L^{2}(0,\ell)$, the following identity holds,

$R^{0}_{b}\big{(}h_{b}S^{*}\widetilde{f}\big{)}=\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}(\frac{\widetilde{f}}{h_{b}})).$

###### Proof of the lemma.

Let us write $v=R^{1}_{b}(\frac{\widetilde{f}}{h_{b}})$; by definition, one
has

$\left\\{\begin{array}[]{l}\big{[}1+\frac{1}{h_{b}}(h_{b}S)\big{(}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(\cdot)\big{)}\big{]}v=\frac{1}{h_{b}}\widetilde{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

$\big{[}(h_{b}S)^{*}+S^{*}(h_{b}S)\big{(}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(\cdot)\big{)}\big{]}v=S^{*}\widetilde{f}$

or equivalently

$\big{[}\alpha_{b}+h_{b}S^{*}(h_{b}S)\big{]}\big{(}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(v)\big{)}=h_{b}S^{*}\widetilde{f};$

using (35), we deduce that

$\big{[}1+h_{b}{\mathcal{T}}_{b}\big{]}\big{(}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(v)\big{)}=h_{b}S^{*}\widetilde{f}$

Since moreover $\big{(}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(v)\big{)}$
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^{0}_{b}\big{(}h_{b}S^{*}\widetilde{f}\big{)},$

which proves the lemma. ∎

Owing to (35), one can write (36) under the form

$R^{0}_{b}\partial_{x}f=R^{0}_{b}\big{(}h_{b}S^{*}(\frac{\sqrt{3}}{h_{b}}f)\big{)}-\frac{3}{2}R^{0}_{b}\big{(}\frac{\partial_{x}b}{h_{b}}f\big{)},$

so that the result follows by applying the lemma with
$\widetilde{f}=\frac{\sqrt{3}}{h_{b}}f$. ∎

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

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

(37) $\zeta(\cdot,0)=g_{0}\quad\mbox{ and }\quad\zeta(\cdot,\ell)=g_{\ell},$

and initial condition

(38) $(\zeta,q)(0,\cdot)=(\zeta^{\rm in},q^{\rm 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

(39) $\begin{cases}\zeta^{\rm in}(0)=g_{0}(0),\\\ -\partial_{x}q^{\rm
in}(0)=\dot{g}_{0}(0),\end{cases}\quad\mbox{ and }\quad\begin{cases}\zeta^{\rm
in}(\ell)=g_{\ell}(0),\\\ -\partial_{x}q^{\rm
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

(40)
$B(\zeta,q)=-{\mathtt{g}}h\partial_{x}b+\frac{3}{2}\frac{\partial_{x}b}{h_{b}}f_{\rm
NSW}.$

We also introduce the two dimensional vectors

(41) $V_{\rm
bdry}(\zeta_{0},q_{0},\zeta_{\ell},q_{\ell})=\begin{pmatrix}\frac{3}{h_{b}(0)^{2}}f_{\rm
NSW}(\zeta_{0},q_{0})\\\ \frac{3}{h_{b}(\ell)^{2}}f_{\rm
NSW}(\zeta_{\ell},q_{\ell})\end{pmatrix}$

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

(42) $V_{\rm
int}[\zeta,q]=\begin{pmatrix}3(R_{b}^{1}(\frac{1}{h_{b}^{2}}f_{\rm
NSW}))_{|_{x=0}}-\big{(}\partial_{x}(R^{0}_{b}B)\big{)}_{|_{x=0}}\\\
3(R_{b}^{1}(\frac{1}{h_{b}^{2}}f_{\rm
NSW}))_{|_{x=\ell}}-\big{(}\partial_{x}(R^{0}_{b}B)\big{)}_{|_{x=\ell}}\end{pmatrix},$

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 (20) and (31). Assume that the
initial and boundary data $(\zeta^{\rm in},q^{\rm in})$ and $(g_{0},g_{\ell})$
satisfy the compatibility condition (39). Then the two following assertions
are equivalent:
i. The couple $(\zeta,q)$ is a regular solution to (23) such that the depth
$h$ never vanishes and with boundary conditions (37) and initial condition
(38).
ii. The couple $(\zeta,q)$ is a regular solution such that the depth $h$ never
vanishes to

(43) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
\partial_{t}q+\sqrt{3}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}(\frac{1}{h_{b}^{2}}f_{\rm
NSW}))=R^{0}_{b}B+\dot{q}_{0}{\mathfrak{s}}^{(b,0)}+\dot{q}_{\ell}{\mathfrak{s}}^{(b,\ell)},\end{cases}\quad\text{in}\
(0,\ell),$

where we recall that
$(h_{b}S)^{*}=\frac{1}{\sqrt{3}}\frac{1}{h_{b}}\partial_{x}(h_{b}^{2}\cdot)+\frac{\sqrt{3}}{2}\partial_{x}b$
and that $B(\zeta,q)$ is given by (40), and with initial condition (38), while
$q_{0}$ and $q_{\ell}$ solve the ODE

(44) $\displaystyle{\mathfrak{S}}_{b}^{\prime}\begin{pmatrix}\dot{q}_{0}\\\
\dot{q}_{\ell}\end{pmatrix}+V_{\rm bdry}(g_{0},q_{0},g_{\ell},q_{\ell})=V_{\rm
int}[\zeta,q]-\begin{pmatrix}\ddot{g}_{0}\\\ \ddot{g}_{\ell}\end{pmatrix},$

with ${\mathfrak{S}}^{\prime}_{b}$ defined in (30), and with initial condition

(45) $q_{0}(0)=q^{\rm in}(0)\quad\mbox{ and }\quad q_{\ell}(0)=q^{\rm
in}(\ell).$

######  Remark 2.9.

When $b=0$ (flat topography), one can check that
${\mathfrak{f}}_{b}=R^{1}f_{\rm NSW}$ and ${\mathfrak{g}}_{b}=0$, so that (43)
coincides as expected with (17).

###### Proof.

We proceed as in the proof of Proposition 2.2. Applying $R^{0}_{b}$ to the
equation in $q$, and using Proposition 2.5, we first obtain

$\displaystyle\partial_{t}q+\sqrt{3}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}(\frac{1}{h_{b}^{2}}f_{\rm
NSW}))$
$\displaystyle-\frac{3}{2}R^{0}_{b}\big{(}\frac{\partial_{x}b}{h_{b}}f_{\rm
NSW}\big{)}$
$\displaystyle=-{\mathtt{g}}R_{b}^{0}(h\partial_{x}b)+\dot{q}_{0}{\mathfrak{s}}^{(b,0)}+\dot{q}_{\ell}{\mathfrak{s}}^{(b,\ell)}.$

Regrouping the terms involving $R_{b}^{0}$, one deduces the second equation of
(43). Applying $\partial_{x}$ to this equation and using the equation on
$\zeta$, we then get

$\displaystyle-\partial_{t}^{2}\zeta+\sqrt{3}\partial_{x}\big{[}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}(\frac{1}{h_{b}^{2}}f_{\rm
NSW}))\big{]}-\frac{3}{2}\partial_{x}\big{[}R^{0}_{b}\big{(}\frac{\partial_{x}b}{h_{b}}f_{\rm
NSW}\big{)}\big{]}$ (46)
$\displaystyle=-{\mathtt{g}}\partial_{x}R_{b}^{0}(h\partial_{x}b)+\dot{q}_{0}({\mathfrak{s}}^{(b,0)})^{\prime}+\dot{q}_{\ell}({\mathfrak{s}}^{(b,\ell)})^{\prime}.$

We now need the following lemma.

######  Lemma 2.10.

For all $\widetilde{f}\in L^{2}(0,\ell)$, one has

$\sqrt{3}\partial_{x}\big{[}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}\widetilde{f})\big{]}_{|_{x=0,\ell}}=-3\big{[}(\widetilde{f}-R^{1}_{b}\widetilde{f})\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}\big{[}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}\widetilde{f})\big{]}=-\frac{3}{h_{b}}(h_{b}S)\big{[}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}\widetilde{f})\big{]}+\beta\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}\widetilde{f});$

by definition of $R_{b}^{1}$, one has
$(h_{b}S)\big{[}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}\widetilde{f})\big{]}=h_{b}\widetilde{f}-h_{b}R^{1}_{b}\widetilde{f}$
so that

$\sqrt{3}\partial_{x}\big{[}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}\widetilde{f})\big{]}=-3(\widetilde{f}-R^{1}_{b}\widetilde{f})+\beta\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}\widetilde{f}).$

Since by construction $(h_{b}S)^{*}(R^{1}_{b}\widetilde{f})$ vanishes at $x=0$
and $x=\ell$, the result follows upon taking the trace of the above identity
at the two boundary points. ∎

Taking the trace of (46) at $x=0$ ad $x=\ell$ and using the lemma with
$\widetilde{f}=\frac{1}{h_{b}^{2}}f_{\rm NSW}$, we get that $q_{0}$ and
$q_{l}$ satisfy the differential system

$\begin{cases}\dot{q}_{0}{\mathfrak{s}}^{\prime}_{(b,0)}(0)+\dot{q}_{\ell}{\mathfrak{s}}^{\prime}_{(b,\ell)}(0)-{\mathtt{g}}\big{(}\partial_{x}R_{b}^{0}(h\partial_{x}b)\big{)}_{0}=\\\
\qquad-\ddot{g}_{0}-\frac{3}{h_{b}(0)^{2}}f_{\rm
NSW}(g_{0},q_{0})+3\big{(}R^{1}_{b}(\frac{1}{h_{b}^{2}}f_{\rm
NSW})\big{)}_{0}-\frac{3}{2}\big{[}\partial_{x}R^{0}_{b}\big{(}\frac{\partial_{x}b}{h_{b}}f_{\rm
NSW}\big{)}\big{]}_{0},\\\
\dot{q}_{\ell}{\mathfrak{s}}^{\prime}_{(b,0)}(\ell)+\dot{q}_{\ell}{\mathfrak{s}}^{\prime}_{(b,\ell)}(\ell)-{\mathtt{g}}\big{(}\partial_{x}R_{b}^{0}(h\partial_{x}b)\big{)}_{\ell}=\\\
\qquad-\ddot{g}_{\ell}-\frac{3}{h_{b}(\ell)^{2}}f_{\rm
NSW}(g_{\ell},q_{\ell})+3\big{(}R^{1}_{b}(\frac{1}{h_{b}^{2}}f_{\rm
NSW})\big{)}_{\ell}-\frac{3}{2}\big{[}\partial_{x}R^{0}_{b}\big{(}\frac{\partial_{x}b}{h_{b}}f_{\rm
NSW}\big{)}\big{]}_{\ell},\end{cases}$

which corresponds to (44). 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 (20) and (31). Let
$g_{0},g_{\ell}\in C^{\infty}({\mathbb{R}}^{+})$. Let also
$n\in{\mathbb{N}}\backslash\\{0\\}$ and $(\zeta^{\rm in},q^{\rm in})\in
H^{n}(0,\ell)\times H^{n+1}(0,\ell)$ be such that
$\inf_{[0,\ell]}(h_{0}-b+\zeta^{\rm in})>0$ and assume that the compatibility
conditions (39) 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 (43) and (44) with initial
conditions (38) and (45), 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 (43)-(44) can be rewritten as

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

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

${\mathbb{F}}[t,{\mathbb{U}}]=(-\partial_{x}q,-\partial_{x}({\mathfrak{f}}_{b}[\zeta,q])+{\mathfrak{g}}_{b}[\zeta,q]+\begin{pmatrix}{\mathfrak{s}}^{(b,0)}\\\
{\mathfrak{s}}^{(b,\ell)}\end{pmatrix}\cdot\bm{\mathfrak{h}}_{b}[t,{\mathbb{U}}],\bm{\mathfrak{h}}_{b}[t,{\mathbb{U}}]^{\rm
T})^{\rm T},$

where

$\bm{\mathfrak{h}}_{b}[t,{\mathbb{U}}]=({\mathfrak{S}}_{b}^{\prime})^{-1}\big{[}-V_{\rm
bdry}(g_{0}(t),q_{0},g_{\ell}(t),q_{\ell})+V_{\rm
int}[\zeta,q]-\begin{pmatrix}\ddot{g}_{0}(t)\\\
\ddot{g}_{\ell}(t)\end{pmatrix}\big{]}.$

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})^{\rm T}\in{\mathbb{X}}$ such
that $\inf_{(0,\ell)}(h_{0}+\zeta-b)>0$. Then it follows from standard product
estimates, and because $R^{0}_{b}$ and $R^{1}_{b}$ 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}}^{\rm 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
(43)-(44) with boundary data $g_{0}=g_{\ell}=0$, one gets in particular the
following identity from the second equation of (43)

$\sqrt{3}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}(\frac{1}{h_{b}^{2}}f_{\rm
NSW}(0,0)))=R^{0}_{b}B(0,0).$

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

(47) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
\partial_{t}q+\sqrt{3}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}(\frac{1}{h_{b}^{2}}\widetilde{f}_{\rm
NSW}))=R^{0}_{b}\widetilde{B}+\dot{q}_{0}{\mathfrak{s}}^{(b,0)}+\dot{q}_{\ell}{\mathfrak{s}}^{(b,\ell)},\end{cases}\quad\text{in}\
(0,\ell),$

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

(48)
$\widetilde{B}(\zeta,q)=-{\mathtt{g}}\zeta\partial_{x}b+\frac{3}{2}\frac{\partial_{x}b}{h_{b}}\widetilde{f}_{\rm
NSW}.$

Proceeding similarly with (44) we get that $q_{0}$ and $q_{\ell}$ solve the
ODE

(49) $\displaystyle{\mathfrak{S}}_{b}^{\prime}\begin{pmatrix}\dot{q}_{0}\\\
\dot{q}_{\ell}\end{pmatrix}+\widetilde{V}_{\rm
bdry}(g_{0},q_{0},g_{\ell},q_{\ell})=\widetilde{V}_{\rm
int}[\zeta,q]-\begin{pmatrix}\ddot{g}_{0}\\\ \ddot{g}_{\ell}\end{pmatrix},$

with

(50) $\widetilde{V}_{\rm
bdry}(\zeta_{0},q_{0},\zeta_{\ell},q_{\ell})=\begin{pmatrix}\frac{3}{h_{b}(0)^{2}}\widetilde{f}_{\rm
NSW}(\zeta_{0},q_{0})\\\ \frac{3}{h_{b}(\ell)^{2}}\widetilde{f}_{\rm
NSW}(\zeta_{\ell},q_{\ell})\end{pmatrix}$

and

(51) $\widetilde{V}_{\rm
int}[\zeta,q]=\begin{pmatrix}3(R_{b}^{1}(\frac{1}{h_{b}^{2}}\widetilde{f}_{\rm
NSW}))_{|_{x=0}}-\big{(}\partial_{x}(R^{0}_{b}\widetilde{B})\big{)}_{|_{x=0}}\\\
3(R_{b}^{1}(\frac{1}{h_{b}^{2}}\widetilde{f}_{\rm
NSW}))_{|_{x=\ell}}-\big{(}\partial_{x}(R^{0}_{b}\widetilde{B})\big{)}_{|_{x=\ell}}\end{pmatrix},$

In order to obtain well-balanced numerical schemes, we will use (47) and (49)
rather than (43) and (44).

##  3 More general boundary conditions

In the previous sections, we considered the Boussinesq-Abbott system (23) 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

(52)
$\displaystyle\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 (52), 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
${\mathcal{U}}\subset{\mathbb{R}}^{2}$ verifying $(0,0)\in\mathcal{U}$ such
that the following conditions hold:
i. The input functions $\xi^{+}_{0}:{\mathcal{U}}\to{\mathbb{R}}$ and
$\xi^{-}_{\ell}:{\mathcal{U}}\to{\mathbb{R}}$ are well defined and smooth;
ii. There are two output functions $\xi^{-}_{0}:{\mathcal{U}}\to{\mathbb{R}}$
and $\xi^{+}_{\ell}:{\mathcal{U}}\to{\mathbb{R}}$ and two open sets
${\mathcal{V}}_{0}$ and ${\mathcal{V}}_{\ell}$ that contain the range of the
mappings
$\xi_{0}:{\mathcal{U}}\ni(\zeta,q)\mapsto(\xi_{0}^{+}(\zeta,q),\xi_{0}^{-}(\zeta,q))\in{\mathbb{R}}^{2}$
and
$\xi_{\ell}:{\mathcal{U}}\ni(\zeta,q)\mapsto(\xi_{\ell}^{+}(\zeta,q),\xi_{\ell}^{-}(\zeta,q))\in{\mathbb{R}}^{2}$
respectively, as well as two smooth mappings
${\mathcal{H}}_{0}:{\mathcal{V}}_{0}\to{\mathbb{R}}^{2}$ and
${\mathcal{H}}_{\ell}:{\mathcal{V}}_{\ell}\to{\mathbb{R}}^{2}$ such that for
all $(\zeta,q)\in{\mathcal{U}}$, one has

(53)
${\mathcal{H}}_{0}(\xi^{+}_{0}(\zeta,q),\xi^{-}_{0}(\zeta,q))=(\zeta,q)^{\rm
T}\quad\mbox{ and
}\quad{\mathcal{H}}_{\ell}(\xi^{+}_{\ell}(\zeta,q),\xi^{-}_{\ell}(\zeta,q))=(\zeta,q)^{\rm
T};$

we refer to ${\mathcal{H}}_{0}$ and ${\mathcal{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{{\mathtt{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\mbox{ and
}\quad\xi_{\ell}^{-}(\zeta,q)=R^{-}(\zeta,q).$

We can then choose

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

where $R^{+}=\frac{q}{h_{0}+\zeta-b}+2\sqrt{{\mathtt{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

$\displaystyle{\mathcal{H}}_{0}(\xi^{+},\xi^{-})$
$\displaystyle=(\xi^{+},\xi^{-})^{\rm T},$
$\displaystyle{\mathcal{H}}_{\ell}(\xi^{+},\xi^{-})$
$\displaystyle=\big{(}\frac{1}{16{\mathtt{g}}}(\xi^{+}-\xi^{-})^{2}-h_{0}+b,\frac{\xi^{+}+\xi^{-}}{32{\mathtt{g}}}(\xi^{+}-\xi^{-})^{2}\big{)}^{\rm
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 know 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^{\pm}_{0}$ and $\xi^{\pm}_{\ell}$ 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
(49) stems from a more general differential identity that relates $\zeta_{0}$,
$q_{0}$, $\zeta_{\ell}$ and $q_{\ell}$, namely,

(54) $\displaystyle{\mathfrak{S}}_{b}^{\prime}\begin{pmatrix}\dot{q}_{0}\\\
\dot{q}_{\ell}\end{pmatrix}+\widetilde{V}_{\rm
bdry}(\zeta_{0},q_{0},\zeta_{\ell},q_{\ell})=\widetilde{V}_{\rm
int}[\zeta,q]-\begin{pmatrix}\ddot{\zeta}_{0}\\\
\ddot{\zeta}_{\ell}\end{pmatrix};$

in fact we see that (54) reduces to (49) 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_{\ell},q_{\ell})$ when
enforcing the general boundary conditions (52). By definition of
${\mathcal{H}}_{0}$ and ${\mathcal{H}}_{\ell}$, and using the boundary
conditions (52), one has

$\begin{pmatrix}\zeta_{0}\\\
q_{0}\end{pmatrix}={\mathcal{H}}_{0}(g_{0},\underline{\xi}_{0}^{-})\quad\mbox{
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})^{\rm T}$
and ${\mathcal{H}}_{\ell}=({\mathcal{H}}_{\ell,1},{\mathcal{H}}_{\ell,2})^{\rm
T}$, one deduces that

(55)
$\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\mbox{ 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

$\displaystyle\ddot{\zeta}_{0}$
$\displaystyle=\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{\rm
Hess}_{{\mathcal{H}}_{0,1}}(g_{0},\underline{\xi}^{-}_{0})\begin{pmatrix}\dot{g}_{0}\\\
\dot{\underline{\xi}}_{0}^{-}\end{pmatrix}$ $\displaystyle\ddot{\zeta}_{\ell}$
$\displaystyle=\nabla{\mathcal{H}}_{\ell,1}(\underline{\xi}_{\ell}^{+},g_{\ell})\cdot\begin{pmatrix}\ddot{\underline{\xi}}_{\ell}^{+}\\\
\ddot{g_{\ell}}\end{pmatrix}+\begin{pmatrix}\dot{\underline{\xi}}_{\ell}^{+}\\\
\dot{g}_{\ell}\end{pmatrix}\cdot{\rm
Hess}_{{\mathcal{H}}_{\ell,1}}(\underline{\xi}_{\ell}^{+},g_{\ell})\begin{pmatrix}\dot{\underline{\xi}}_{\ell}^{+}\\\
\dot{g}_{\ell}\end{pmatrix},$

where ${\rm Hess}_{{\mathcal{H}}_{0,1}}$ and ${\rm
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

$\widetilde{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

$\displaystyle{\mathcal{Q}_{1}}(u_{0},v_{0})$
$\displaystyle=\begin{pmatrix}u_{0}\\\ v_{0}\end{pmatrix}\cdot{\rm
Hess}_{{\mathcal{H}}_{0,1}}(g_{0},\underline{\xi}^{-}_{0})\begin{pmatrix}u_{0}\\\
v_{0}\end{pmatrix},$ $\displaystyle{\mathcal{Q}_{2}}(u_{\ell},v_{\ell})$
$\displaystyle=\begin{pmatrix}u_{\ell}\\\ v_{\ell}\end{pmatrix}\cdot{\rm
Hess}_{{\mathcal{H}}_{\ell,1}}(\underline{\xi}_{\ell}^{+},g_{\ell})\begin{pmatrix}u_{\ell}\\\
v_{\ell}\end{pmatrix},$

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

$\displaystyle\widetilde{D}_{1}\begin{pmatrix}\ddot{\underline{\xi}}^{-}_{0}\\\
\ddot{\underline{\xi}}^{+}_{\ell}\end{pmatrix}+{\mathfrak{S}}^{\prime}_{b}\widetilde{D}_{2}\begin{pmatrix}\dot{\underline{\xi}}^{-}_{0}\\\
\dot{\underline{\xi}}^{+}_{\ell}\end{pmatrix}$
$\displaystyle+\begin{pmatrix}{\mathcal{Q}}_{1}(\dot{g}_{0},\dot{\underline{\xi}}_{0}^{-})\\\
{\mathcal{Q}}_{2}(\dot{\underline{\xi}}_{\ell}^{+},\dot{g}_{\ell})\end{pmatrix}+\widetilde{V}_{\rm
bdry}\big{(}{\mathcal{H}}_{0}(g_{0},\underline{\xi}_{0}^{-}),{\mathcal{H}}_{\ell}(\underline{\xi}_{\ell}^{+},g_{\ell})\big{)}$
(56) $\displaystyle=\widetilde{V}_{\rm
int}[\zeta,q]-{D}_{1}\begin{pmatrix}\ddot{g}_{0}\\\
\ddot{g}_{\ell}\end{pmatrix}-{\mathfrak{S}}^{\prime}_{b}{D}_{2}\begin{pmatrix}\dot{g}_{0}\\\
\dot{g}_{\ell}\end{pmatrix}.$

The order of this system of ODEs depends on whether the coefficients of
$\widetilde{D}_{1}$ vanish or not. We make the following assumption which
ensures that (56) 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 ${\mathcal{V}}_{0}$, or
$\partial_{2}{\mathcal{H}}_{0,1}$ does not vanish on ${\mathcal{V}}_{0}$;
ii. Either $\partial_{1}{\mathcal{H}}_{\ell,1}\equiv 0$ and
$\partial_{1}{\mathcal{H}}_{\ell,2}$ does not vanish on ${\mathcal{V}}_{0}$,
or $\partial_{1}{\mathcal{H}}_{\ell,1}$ does not vanish on
${\mathcal{V}}_{\ell}$

######  Remark 3.4.

If $\partial_{2}{\mathcal{H}}_{0,1}\equiv 0$ on ${\mathcal{V}}_{0}$ then
${\mathcal{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 (56) 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, ${\mathcal{Q}}_{1}=0$. We also compute
$\partial_{1}{\mathcal{H}}_{\ell,1}=\frac{1}{8{\mathtt{g}}}(\xi^{+}-\xi^{-})$.
It is possible to choose the open set ${\mathcal{V}}_{\ell}$ in such a way
that $\xi^{+}-\xi^{-}>0$ for all $(\xi^{+},\xi^{-})\in{\mathcal{V}}_{\ell}$
provided that we assume the total water height $h$ never vanishes. Assumption
3.3 is therefore satisfied. Recalling the definition (30) of
${\mathfrak{S}}^{\prime}_{b}$, the system (56) can be put under the form

$\displaystyle T\begin{pmatrix}\dot{\underline{\xi}}_{0}^{-}\vspace{1mm}\\\
\ddot{\underline{\xi}}_{\ell}^{+}\end{pmatrix}+{\mathfrak{S}}^{\prime}_{b}\widetilde{D}_{2}\begin{pmatrix}0\vspace{1mm}\\\
\dot{\xi}^{+}_{\ell}\end{pmatrix}$
$\displaystyle+\begin{pmatrix}0\vspace{1mm}\\\
{\mathcal{Q}}_{2}(\dot{\underline{\xi}}_{\ell}^{+},\dot{g}_{\ell})\end{pmatrix}+\widetilde{V}_{\rm
bdry}\big{(}{\mathcal{H}}_{0}(g_{0},\underline{\xi}_{0}^{-}),{\mathcal{H}}_{\ell}(\underline{\xi}_{\ell}^{+},g_{\ell})\big{)}$
(57) $\displaystyle=\widetilde{V}_{\rm
int}[\zeta,q]-{D}_{1}\begin{pmatrix}\ddot{g}_{0}\\\
\ddot{g}_{\ell}\end{pmatrix}-{\mathfrak{S}}^{\prime}_{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})^{\prime}(0)&0\\\
({\mathfrak{s}}^{b,0})^{\prime}(\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})^{\prime}(0)\neq 0$ and $T$ is invertible, and the
system (57) 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^{+})$
and in particular $\partial_{1}{\mathcal{H}}_{0,1}\equiv 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 ${\mathcal{V}}_{0}$ then
$\xi_{0}^{+}(\zeta,q)=\xi_{0}^{+}(q)$ on ${\mathcal{U}}$.
ii. If $\partial_{2}{\mathcal{H}}_{\ell,1}$ vanishes on ${\mathcal{V}}_{\ell}$
then $\xi_{\ell}^{-}(\zeta,q)=\xi_{\ell}^{-}(q)$ on ${\mathcal{U}}$.

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

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

(58) $\xi_{0}^{+}(\zeta^{\rm in}(0),q^{\rm in}(0))=g_{0}(0)\quad\mbox{ and
}\quad\xi_{\ell}^{-}(\zeta^{\rm in}(\ell),q^{\rm in}(\ell))=g_{\ell}(0);$

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

(59) $\underline{\xi}^{-}_{0}(0)=\xi^{-}_{0}(\zeta^{\rm in}(0),q^{\rm
in}(0))\quad\mbox{ and
}\quad\underline{\xi}^{-}_{\ell}(0)=\xi^{+}_{\ell}(\zeta^{\rm in}(\ell),q^{\rm
in}(\ell)).$

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

$\displaystyle\dot{\zeta}_{0}$
$\displaystyle=\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}^{-},$
$\displaystyle\dot{\zeta}_{\ell}$
$\displaystyle=\partial_{1}{\mathcal{H}}_{\ell,1}(\underline{\xi}_{\ell}^{+},g_{\ell})\dot{\underline{\xi}}_{\ell}^{+}+\partial_{2}{\mathcal{H}}_{\ell,1}(\underline{\xi}_{\ell}^{+},g_{\ell})\dot{g}_{\ell}.$

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

(60)
$\begin{cases}\partial_{2}{\mathcal{H}}_{0,1}(g_{0}(0),\underline{\xi}_{0}^{-}(0))\,\dot{\underline{\xi}}_{0}^{-}&=-\partial_{x}q^{\rm
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^{\rm
in}(\ell)-\partial_{2}{\mathcal{H}}_{\ell,1}(\underline{\xi}_{\ell}^{+}(0),g_{\ell}(0))\dot{g}_{\ell}(0),\end{cases}$

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 (56) 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 (56) is of second order in $\underline{\xi}_{0}^{-}$ (resp.
$\partial_{1}{\mathcal{H}}_{\ell,1}$ does not vanish if (56) is of second
order in $\underline{\xi}_{\ell}^{+}$), (60) furnishes the needed initial data
for $\dot{\underline{\xi}}^{-}_{0}$ (resp. on
$\dot{\underline{\xi}}_{\ell}^{+}$).

  * •

If (56) 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 (60) is then a
compatibility condition on the data that generalizes the second equations in
the compatibility conditions (39) imposed in Proposition 2.11.

######  Example 3.7.

For the configuration considered in Examples 3.2 and 3.5, and assuming that
(58) 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 (59) and also, at $x=\ell$, the second condition of
(60). 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 (59) (and which coincides in this case
with the compatibility condition (39), namely, $-\partial_{x}q^{\rm
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 (23) with general boundary
conditions (52). 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 (20) and (31). 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}}\backslash\\{0\\}$ and $(\zeta^{\rm in},q^{\rm in})\in
H^{n}(0,\ell)\times H^{n+1}(0,\ell)$ be such that
$\inf_{[0,\ell]}(h_{0}-b+\zeta^{\rm in})>0$ and assume that the compatibility
condition (58) 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 (60) 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 (23) with initial conditions (38) and boundary
conditions (52), 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 (47). Substituting for
$\dot{q}_{0}$ and $\dot{q}_{\ell}$ in the right-hand side of (47) using (55),
we obtain that

(61) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
\partial_{t}q+\sqrt{3}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}(\frac{1}{h_{b}^{2}}\widetilde{f}_{\rm
NSW}))=R^{0}_{b}\widetilde{B}\\\ \hskip
42.67912pt+\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{cases}$

This system is complemented by the ODE (56) 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 (56) is
explicit and of first order in $\underline{\xi}_{0}^{-}$. The needed initial
condition on $\underline{\xi}_{0}^{-}$ is furnished by (59).

  * •

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 (56) is explicit and of second order
in $\underline{\xi}_{0}^{-}$. In addition to the initial condition on
$\underline{\xi}_{0}^{-}$ furnished by (59), we need an initial condition on
$\dot{\underline{\xi}}_{0}^{-}$, which is furnished by (60).

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 (61), (56) and the initial conditions. We need to
prove now that $(\zeta,q)$ solves the Boussinesq-Abbott system (23) with
initial conditions (38) and boundary conditions (52).
Applying $(1+h_{b}{\mathcal{T}}_{b})$ to (61), we obtain that $(\zeta,q)$
solves (23); it remains to prove that the boundary conditions (52) are
satisfied. Taking the trace of the second equation of (61) at $x=0$ and
$x=\ell$, we deduce from (28), (29) and (34) that

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

Since we also know from Assumption 3.1 and (59) that
${\mathcal{H}}_{0,2}(g_{0},\underline{\xi}^{-}_{0})_{|_{t=0}}=q^{\rm in}(0)$
and
${\mathcal{H}}_{\ell,2}(\underline{\xi}^{+}_{\ell},g_{\ell})_{|_{t=0}}=q^{\rm
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
${\mathcal{H}}_{0,2}(\xi^{+},\xi^{-})=\xi^{+}$ so that we deduce that
$q(t,0)=g_{0}$, which is exactly the boundary condition (52).

  * •

If this is not the case, then we know by Assumption 3.6 that
$\partial_{1}{\mathcal{H}}_{0,1}(\xi^{+},\xi^{-})$ does not vanish on
${\mathcal{V}}_{0}$; therefore the ODE

$\displaystyle\widetilde{D}_{1}\begin{pmatrix}\ddot{\underline{\xi}}^{-}_{0}\\\
\ddot{\underline{\xi}}^{+}_{\ell}\end{pmatrix}+{\mathfrak{S}}^{\prime}_{b}\widetilde{D}_{2}\begin{pmatrix}\dot{\underline{\xi}}^{-}_{0}\\\
\dot{\underline{\xi}}^{+}_{\ell}\end{pmatrix}$
$\displaystyle+\begin{pmatrix}{\mathcal{Q}}_{1}(\dot{\underline{\xi}}^{+}_{0},\dot{\underline{\xi}}_{0}^{-})\\\
{\mathcal{Q}}_{2}(\dot{\underline{\xi}}_{\ell}^{+},\dot{\underline{\xi}}^{-}_{\ell})\end{pmatrix}+\widetilde{V}_{\rm
bdry}({\mathcal{H}}_{0}(\underline{\xi}^{+}_{0},\underline{\xi}_{0}^{-}),{\mathcal{H}}_{\ell}(\underline{\xi}_{\ell}^{+},\underline{\xi}^{+}_{\ell}))$
$\displaystyle=\widetilde{V}_{\rm
int}[\zeta,q]-{D}_{1}\begin{pmatrix}\ddot{\underline{\xi}}^{+}_{0}\\\
\ddot{\underline{\xi}}^{-}_{\ell}\end{pmatrix}-{\mathfrak{S}}^{\prime}_{b}{D}_{2}\begin{pmatrix}\dot{\underline{\xi}}^{+}_{0}\\\
\dot{\underline{\xi}}^{-}_{\ell}\end{pmatrix}.$

which is obtained as (56), is of second order in
$\underline{\xi}_{0}^{+}=\xi^{+}_{0}(\zeta_{0},q_{0})$. It follows that
$(\underline{\xi}_{0}^{+},\underline{\xi}_{\ell}^{-})$ satisfy the same ODE as
$(g_{0},g_{\ell})$. These two quantities are therefore equal (which implies
that the boundary condition (52) is satisfied) if
$(\underline{\xi}_{0}^{+}(0),\underline{\xi}_{\ell}^{-}(0))=(g_{0}(0),g_{\ell}(0))$
and
$(\dot{\underline{\xi}}_{0}^{+}(0),\dot{\underline{\xi}}_{\ell}^{-}(0))=(\dot{g}_{0}(0),\dot{g}_{\ell}(0))$.
The first of these two conditions corresponds to (58). For the second one
(focusing our attention on the case $x=0$), we observe by time differentiating
the relation
${\mathcal{H}}_{0,1}(\underline{\xi}_{0}^{+},\underline{\xi}_{0}^{-})=\zeta(t,0)$
and using the first equation of (23) that

$\partial_{1}{\mathcal{H}}_{0,1}(\underline{\xi}_{0}^{+}(0),\underline{\xi}_{0}^{-}(0))\dot{\underline{\xi}}_{0}^{+}(0)+\partial_{2}{\mathcal{H}}_{0,1}(\underline{\xi}_{0}^{+}(0),\underline{\xi}_{0}^{-}(0))\dot{\underline{\xi}}_{0}^{-}(0)=-\partial_{x}q^{\rm
in}(0).$

Since $\underline{\xi}_{0}^{+}(0)=g(0)$ by (58), we can use (60) (which holds
by assumption if $\partial_{2}{\mathcal{H}}_{0,1}\equiv 0$ and by the initial
condition imposed on $\dot{\underline{\xi}}_{0}^{-}$ otherwise) to obtain

$\partial_{1}{\mathcal{H}}_{0,1}(\underline{\xi}_{0}^{+}(0),\underline{\xi}_{0}^{-}(0))\big{(}\dot{\underline{\xi}}_{0}^{+}(0)-\dot{g}_{0}(0)\big{)}=0,$

and since $\partial_{1}{\mathcal{H}}_{0,1}$ does not vanish, this implies that
$\dot{\underline{\xi}}_{0}^{+}(0)=\dot{g}_{0}(0)$, so that the boundary
condition (52) 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 (52) 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 (23) on the full
domain $(0,\ell)$ for some time interval $(0,T^{*})$ with $T^{*}>0$ provided
that we know the initial data $(\zeta^{\rm in},q^{\rm 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 $(\widetilde{\zeta},\widetilde{q})$ of
the initial boundary value problem with the same boundary conditions but with
initial data $(\widetilde{\zeta}^{\rm in},\widetilde{q}^{\rm
in})\neq(\zeta^{\rm in},q^{\rm in})$ do we have
$(\widetilde{\zeta},\widetilde{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 (23) with initial conditions (38) and boundary
conditions (52) 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

(62) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
(1-\frac{h_{0}^{2}}{3}\partial_{x}^{2})\partial_{t}q+{\mathtt{g}}h_{0}\partial_{x}\zeta=0,\end{cases}$

with boundary conditions

(63) $\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

(64)
$E(t)=\frac{1}{2}{\mathtt{g}}|\widetilde{\zeta}(t)|^{2}_{L^{2}(0,\ell)}+\frac{1}{2}\frac{1}{h_{0}}|\widetilde{q}(t)|^{2}_{L^{2}(0,\ell)}+\frac{h_{0}}{6}|\partial_{x}\widetilde{q}(t)|^{2}_{L^{2}(0,\ell)}$

######  Proposition 3.9.

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

###### Proof.

Multiplying (62) by
$({\mathtt{g}}\widetilde{\zeta},\frac{1}{h_{0}}\widetilde{q})$ and integrating
by parts and using the boundary conditions (63), we obtain

$\frac{{\rm d}}{{\rm
d}t}E-\frac{h_{0}}{3}[\partial_{x}\partial_{t}\widetilde{q}\widetilde{q}]_{0}^{\ell}=0.$

Using the first equation of (62), one has
$\partial_{x}\partial_{t}\widetilde{q}=-\partial^{2}_{t}\widetilde{\zeta}$ so
that using (63) we deduce that
$[\partial_{x}\partial_{t}\widetilde{q}\widetilde{q}]_{0}^{\ell}=0$, and the
result follows easily. ∎

The fact that we allow in (52) 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 (23), namely,

$\displaystyle\left\\{\begin{array}[]{l}\partial_{t}\zeta+\partial_{x}q=0\\\
(1+h_{b}{\mathcal{T}}_{b})\partial_{t}q+\partial_{x}f_{\text{NSW}}=-{\mathtt{g}}h\partial_{x}b\end{array}\right.\quad\text{in}\
(0,\ell)\ ,$

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

$\displaystyle\xi^{+}_{0}(\zeta_{0},q_{0})=g_{0},\qquad\xi^{-}_{\ell}(\zeta_{\ell},q_{\ell})=g_{\ell},$

and with initial condition

$(\zeta,q)=(\zeta^{\rm in},q^{\rm in})\quad\mbox{ 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
(61) of the equations that reads

(65) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
\partial_{t}q+\sqrt{3}\frac{h_{b}}{\alpha_{b}}(h_{b}S)^{*}(R^{1}_{b}(\frac{1}{h_{b}^{2}}\widetilde{f}_{\rm
NSW}))=R^{0}_{b}\widetilde{B}\\\ \hskip
42.67912pt+\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{cases}$

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

$\displaystyle\widetilde{D}_{1}\begin{pmatrix}\ddot{\underline{\xi}}^{-}_{0}\\\
\ddot{\underline{\xi}}^{+}_{\ell}\end{pmatrix}+{\mathfrak{S}}^{\prime}_{b}\widetilde{D}_{2}\begin{pmatrix}\dot{\underline{\xi}}^{-}_{0}\\\
\dot{\underline{\xi}}^{+}_{\ell}\end{pmatrix}$
$\displaystyle+\begin{pmatrix}{\mathcal{Q}}_{1}(\dot{g}_{0},\dot{\underline{\xi}}_{0}^{-})\\\
{\mathcal{Q}}_{2}(\dot{\underline{\xi}}_{\ell}^{+},\dot{g}_{\ell})\end{pmatrix}+\widetilde{V}_{\rm
bdry}\big{(}{\mathcal{H}}_{0}(g_{0},\underline{\xi}_{0}^{-}),{\mathcal{H}}_{\ell}(\underline{\xi}_{\ell}^{+},g_{\ell})\big{)}$
(66) $\displaystyle=\widetilde{V}_{\rm
int}[\zeta,q]-{D}_{1}\begin{pmatrix}\ddot{g}_{0}\\\
\ddot{g}_{\ell}\end{pmatrix}-{\mathfrak{S}}^{\prime}_{b}{D}_{2}\begin{pmatrix}\dot{g}_{0}\\\
\dot{g}_{\ell}\end{pmatrix}.$

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 (65), and a finite difference discretization of the
nonlocal operators $R_{b}^{0},R_{b}^{1}$ defined in (28)-(34) and of the
system of ODEs (66) 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 (65); 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\leq i\leq N$ with
$(N-1)\Delta x=\ell$. A dual mesh is then obtained as the $N$ cells centered
on the $(x_{i})_{1\leq i\leq 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 (52) 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)\mathrm{d}s$, 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_{b}S)$, its adjoint $(h_{b}S)^{*}$ 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}\rightarrow\mathbb{R}^{N}$
defined such that for any $v\in\mathbb{R}^{N}$ and $1\leq i\leq N$ we have

(67) $\delta_{x}(v)_{i}=\left\\{\begin{array}[]{ll}\frac{1}{2\Delta
x}(v_{i+1}-v_{i-1})&\text{if}\ 2\leq i\leq N-1,\\\\[5.0pt] \frac{1}{2\Delta
x}(-3v_{1}+4v_{2}-v_{3})&\text{if}\ i=1,\\\\[5.0pt] \frac{1}{2\Delta
x}(v_{N-2}-4v_{N-1}+3v_{N})&\text{if}\ i=N.\end{array}\right.$

For any index $1\leq i\leq 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_{b}S)$ by $(\underline{h_{b}S})$ as follows:

(68) $\forall
v\in\mathbb{R}^{N},\quad(\underline{h_{b}S})(v)_{i}=-\frac{h_{b,i}^{2}}{\sqrt{3}}\delta_{x}\Big{(}\frac{v}{h_{b}}\Big{)}_{i}+\frac{\sqrt{3}}{2}\delta_{x}(b)_{i}v_{i},$

which is consistent with the definition $\eqref{defS}$ of $S(\cdot)$. In a
similar fashion we introduce $(\underline{h_{b}S})^{*}$ defined as

(69) $\displaystyle(\underline{h_{b}S})^{*}(f)_{i}={}$
$\displaystyle\frac{\alpha_{b,i}}{\sqrt{3}h_{b,i}}\delta_{x}\Big{(}\frac{h_{b}^{2}f}{\alpha_{b}}\Big{)}_{i}+\Big{(}\frac{h_{b,i}}{\sqrt{3}}\frac{\delta_{x}(\alpha_{b})_{i}}{\alpha_{b,i}}+\frac{\sqrt{3}}{2}\delta_{x}b_{i}\Big{)}f_{i},$

and that is shown to be consistent with $(h_{b}S)^{*}$ 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

$\displaystyle\underline{R}_{b}^{0}:f$
$\displaystyle\in\mathbb{R}^{N}\longmapsto
v=\underline{R}_{b}^{0}f\in\mathbb{R}^{N}\text{ such that}$ (70)
$\displaystyle\begin{cases}\displaystyle
v_{i}-\delta_{x}\Big{(}\frac{h_{b}^{3}}{3}\delta_{x}\Big{(}\frac{v}{h_{b}}\Big{)}\Big{)}_{i}+\frac{h_{b,i}}{2\Delta
x^{2}}(b_{i+1}-2b_{i}+b_{i-1})v_{i}=f_{i},&1<i<N,\\\ v_{1}=0,\quad
v_{N}=0.\end{cases}$

and one readily checks that this is a consistent approximation of
$(1+h_{b}\mathcal{T}_{b})$ at order $2$, with the operator $\mathcal{T}_{b}$
given in (24). Likewise $R_{b}^{1}$ is discretized by $\underline{R}_{b}^{1}$
defined implicitly as

$\displaystyle\underline{R}_{b}^{1}:f$
$\displaystyle\in\mathbb{R}^{N}\longmapsto
v=\underline{R}_{b}^{1}f\in\mathbb{R}^{N}\text{ such that}$ (71)
$\displaystyle\begin{cases}\displaystyle
v_{i}+\frac{1}{h_{b,i}}(\underline{h_{b}S})\Big{(}\frac{h_{b}}{\alpha_{b}}(\underline{h_{b}S})^{*}(v)\Big{)}_{i}=f_{i},&1<i<N,\\\
(\underline{h_{b}S})^{*}(v)_{1}=0,\quad(\underline{h_{b}S})^{*}(v)_{N}=0,\end{cases}$

which is consistent with (34). The definitions (70)-(71) 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 (71) 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 $\mathrm{Im}(\underline{R}_{b}^{1})$. This is the
discrete counterpart to the fact that
$(h_{b}S)^{*}(R_{b}^{1}\cdot)_{|_{x=0,\ell}}\equiv 0$ in $L^{2}(0,\ell)$,
which enables to obtain the relation (66) — ensuing from Lemma 2.10 after
applying $\partial_{x}$ to the second equation of (65) and taking the trace at
$x=0,\ell$ — and hereby to ensure that the boundary conditions (52) are
satisfied. Therefore it is important to have this property at the discrete
level to remain consistent with the continuous reformulation (65)-(66).

###  4.2 First order Lax-Friedrichs scheme

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

(72) $\frac{\zeta_{i}^{n+1}-\zeta_{i}^{n}}{\Delta t}+\frac{1}{\Delta
x}\Big{(}q_{i+1/2}^{n}-q_{i-1/2}^{n}\Big{)}=0,$

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

(73) $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 (65) we first introduce the vectors
$\widetilde{f}_{\mathrm{NSW}}^{n},\widetilde{B}^{n}\in\mathbb{R}^{N}$ with
coefficients

(74)
$\displaystyle\begin{cases}\displaystyle\widetilde{f}_{\mathrm{NSW},i}^{n}=\frac{(q_{i}^{n})^{2}}{h_{0}+\zeta_{i}^{n}-b_{i}}+\frac{\mathtt{g}}{2}((\zeta_{i}^{n})^{2}+2h_{b,i}\zeta_{i}^{n}),\\\
\displaystyle\widetilde{B}_{i}^{n}=-{\mathtt{g}}\zeta_{i}^{n}\delta_{x}(b)_{i}+\frac{3}{2}\frac{\delta_{x}(b)_{i}}{h_{b,i}}\widetilde{f}_{\mathrm{NSW},i}^{n},\end{cases}$

which is consistent with the definition (48) of $\widetilde{B}(\zeta,q)$. The
update of the discharge in cells $1\leq i\leq N$ makes use of the discrete
operators introduced in Section 4.1 as follows:

(75) $\displaystyle\frac{q_{i}^{n+1}-q_{i}^{n}}{\Delta
t}+\sqrt{3}\frac{h_{b,i}}{\alpha_{b,i}}(\underline{h_{b}S})^{*}\Big{(}\underline{R}_{b}^{1}\Big{(}\frac{1}{h_{b}^{2}}\widetilde{f}_{\mathrm{NSW}}^{n}\Big{)}\Big{)}_{i}={}$
$\displaystyle\quad\underline{R}_{b}^{0}(\widetilde{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}}\mathds{1}_{i\not\in\\{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}}^{(b,0)}$ and $\underline{\mathfrak{s}}^{(b,\ell)}$
are taken as the solutions $v\in\mathbb{R}^{N}$ of

(76) $\displaystyle
v_{i}-\delta_{x}\Big{(}\frac{h_{b}^{3}}{3}\delta_{x}\Big{(}\frac{v}{h_{b}}\Big{)}\Big{)}_{i}+\frac{h_{b,i}}{2\Delta
x^{2}}(b_{i+1}-2b_{i}+b_{i-1})v_{i}=0,\quad 2\leq i\leq N-1$

with respective conditions

(77) $\displaystyle\begin{cases}\underline{\mathfrak{s}}_{1}^{(b,0)}=1\\\
\underline{\mathfrak{s}}_{N}^{(b,0)}=0\end{cases},\qquad\begin{cases}\underline{\mathfrak{s}}_{1}^{(b,\ell)}=0\\\
\underline{\mathfrak{s}}_{N}^{(b,\ell)}=1\end{cases},$

yielding a discrete counterpart of (29). Owing to the conditions (77), to the
definition (70) of $\underline{R}_{b}^{0}$ and to Remark 4.1, the update (75)
in border cells $i=1,N$ reduces to

(78) $\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 (75), which
is a first order numerical diffusion required for stability purposes. We can
then show that the update (75) 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

(79)
$\mathfrak{f}^{n}=\Big{(}\frac{h_{b,i}^{2}}{\alpha_{b,i}}\underline{R}_{b}^{1}\Big{(}\frac{1}{h_{b}^{2}}\widetilde{f}_{\mathrm{NSW}}^{n}\Big{)}_{i}\Big{)}_{1\leq
i\leq N},$

and the associated Lax-Friedrichs flux

(80)
$\displaystyle\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\leq i\leq N-1$ the update (75) equivalently rewrites

(81) $\displaystyle\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}={}$
$\displaystyle\quad\underline{R}_{b}^{0}(\widetilde{B}^{n})_{i}-\Big{(}\frac{\delta_{x}(\alpha_{b})_{i}}{\alpha_{b,i}}+\frac{3}{2}\frac{\delta_{x}(b)_{i}}{h_{b,i}}\Big{)}\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 $\widetilde{B}^{n}$
cancel, so that the discharge update (81) 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 (9) to compute
$\underline{\mathfrak{s}}^{(b,0)}$ and $\underline{\mathfrak{s}}^{(b,\ell)}$
instead of approximating them through (76)-(77). 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 (53) reading as

$U_{1}^{n+1}={\mathcal{H}}_{0}(\xi^{+}_{0}(U_{1}^{n+1}),\xi^{-}_{0}(U_{1}^{n+1}))\quad\mbox{
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 (52) 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 (66) for the traces as follows. Given $X^{n}$
and $Y^{n}$ the respective approximations at time $t^{n}$ of the functions

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

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

(83) $\displaystyle\begin{cases}\cfrac{X^{n+1}-X^{n}}{\Delta t}=Y^{n+1}\\\
\widetilde{D}_{1}^{n}\cfrac{Y^{n+1}-Y^{n}}{\Delta
t}+\underline{\mathfrak{S}}^{\prime}_{b}\widetilde{D}_{2}^{n}Y^{n+1}+\begin{pmatrix}{\mathcal{Q}}_{1}^{n}(\frac{g_{0}^{n+1}-g_{0}^{n-1}}{2\Delta
t},Y_{1}^{n})\\\
{\mathcal{Q}}_{2}^{n}(Y_{2}^{n},\frac{g_{\ell}^{n+1}-g_{\ell}^{n-1}}{2\Delta
t})\end{pmatrix}+\widetilde{V}_{\rm bdry}^{n}\\\ \hskip
60.00009pt=\widetilde{V}_{\rm
int}^{n}-{D}_{1}^{n}\begin{pmatrix}\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}}\end{pmatrix}-\underline{\mathfrak{S}}^{\prime}_{b}{D}_{2}^{n}\begin{pmatrix}\frac{g_{0}^{n+1}-g_{0}^{n-1}}{2\Delta
t}\\\ \frac{g_{\ell}^{n+1}-g_{\ell}^{n-1}}{2\Delta t}\end{pmatrix}\end{cases}$

where the diagonal matrices $D_{j}^{n},\widetilde{D}_{j}^{n}$ are defined for
$j\in\\{1,2\\}$ as

$\displaystyle
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},$
$\displaystyle\widetilde{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

$\displaystyle{\mathcal{Q}_{1}^{n}}(u_{0},v_{0})$
$\displaystyle=\begin{pmatrix}u_{0}\\\ v_{0}\end{pmatrix}\cdot{\rm
Hess}_{{\mathcal{H}}_{0,1}}(g_{0}^{n},X_{1}^{n})\begin{pmatrix}u_{0}\\\
v_{0}\end{pmatrix},$ $\displaystyle{\mathcal{Q}_{2}^{n}}(u_{\ell},v_{\ell})$
$\displaystyle=\begin{pmatrix}u_{\ell}\\\ v_{\ell}\end{pmatrix}\cdot{\rm
Hess}_{{\mathcal{H}}_{\ell,1}}(X_{2}^{n},g_{\ell}^{n})\begin{pmatrix}u_{\ell}\\\
v_{\ell}\end{pmatrix},$

and where

$\displaystyle\underline{\mathfrak{S}}_{b}^{\prime}=\begin{pmatrix}\delta_{x}(\underline{\mathfrak{s}}^{(b,0)})_{1}&\delta_{x}(\underline{\mathfrak{s}}^{(b,\ell)})_{1}\\\
\delta_{x}(\underline{\mathfrak{s}}^{(b,0)})_{N}&\delta_{x}(\underline{\mathfrak{s}}^{(b,\ell)})_{N}\end{pmatrix},\qquad\widetilde{V}_{\rm
bdry}^{n}=\begin{pmatrix}\frac{3}{h_{b,1}^{2}}\widetilde{f}_{\rm
NSW}({\mathcal{H}}_{0}(g_{0}^{n},X_{1}^{n}))\\\
\frac{3}{h_{b,N}^{2}}\widetilde{f}_{\rm
NSW}({\mathcal{H}}_{\ell}(X_{2}^{n},g_{\ell}^{n}))\end{pmatrix}$
$\displaystyle\widetilde{V}_{\rm
int}^{n}=\begin{pmatrix}3(\underline{R}_{b}^{1}(\frac{1}{h_{b}^{2}}\widetilde{f}_{\rm
NSW}^{n}))_{1}-\big{(}\delta_{x}(\underline{R}^{0}_{b}\widetilde{B}^{n})\big{)}_{1}\\\
3(\underline{R}_{b}^{1}(\frac{1}{h_{b}^{2}}\widetilde{f}_{\rm
NSW}^{n}))_{N}-\big{(}\delta_{x}(\underline{R}^{0}_{b}\widetilde{B}^{n})\big{)}_{N}\end{pmatrix}.$

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 (83) are not known
when $n=0$. We suggest to define them following the relations below:

$\displaystyle\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},$
$\displaystyle\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 (83) for $n=0$. These quantities can be obtained by
discretizing (60) in the following way

(84)
$\begin{cases}\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})\cfrac{g_{0}(\Delta
t)-g_{0}(0)}{\Delta t},\\\\[2.5pt]
\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))\cfrac{g_{\ell}(\Delta
t)-g_{\ell}(0)}{\Delta t}.\end{cases}$

Remark that in (84), if for a given index $1\leq k\leq 2$ the coefficient in
factor of $Y_{k}^{0}$ cancels, then the quadratic function
$\mathcal{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 (83) the
value $Y_{k}^{0}$ is not required to compute $Y_{k}^{1}$. The proposed
discretization (83) 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 (78) so that we can
update the interior discharge through (81). The iteration is complete.

######  Remark 4.3.

We need to check that the update (83) 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 $\widetilde{D}_{1}^{n}+\Delta
t\underline{\mathfrak{S}}_{b}^{\prime}\widetilde{D}_{2}^{n}$ is invertible.
Its determinant is the second degree polynomial in $\Delta t$ given by

$\displaystyle\det(\widetilde{D}_{1}^{n})+\Delta
t\Big{(}(\underline{\mathfrak{S}}_{b}^{\prime})_{1,1}(\widetilde{D}_{1}^{n})_{2,2}(\widetilde{D}_{2}^{n})_{1,1}+(\underline{\mathfrak{S}}_{b}^{\prime})_{2,2}(\widetilde{D}_{1}^{n})_{1,1}(\widetilde{D}_{2}^{n})_{2,2}\Big{)}$
$\displaystyle{}+\Delta
t^{2}\Big{(}\det(\underline{\mathfrak{S}}_{b}^{\prime})\det(\widetilde{D}_{2}^{n})\Big{)};$

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 (83) is well-defined.

When approximating the nonlinear Boussinesq-Abbott system (47)-(49) 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}\lvert\lambda(U_{i}^{n})\lvert},$

with $\lvert\lambda(U_{i}^{n})\rvert=\lvert
u_{i}^{n}\rvert+\sqrt{\mathtt{g}h_{i}^{n}}$ the maximum eigenvalue in absolute
value of the Jacobian matrix $Df_{\rm 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\leq i\leq N}$ of the
solution to (47)-(49) at time $t^{n}$, we first compute
$U_{\mathrm{P},i}^{n+1}=(\zeta_{\mathrm{P},i}^{n+1},q_{\mathrm{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\leq i\leq N-1$

(85)
$\displaystyle\begin{cases}\displaystyle\frac{\zeta_{\mathrm{P},i}^{n+1}-\zeta_{i}^{n}}{\Delta
t}+\frac{1}{\Delta x}\Big{(}q_{i}^{n}-q_{i-1}^{n}\Big{)}=0\\\\[7.5pt]
\displaystyle\frac{q_{\mathrm{P},i}^{n+1}-q_{i}^{n}}{\Delta
t}+\frac{\mathfrak{f}_{i}^{n}-\mathfrak{f}_{i-1}^{n}}{\Delta
x}=\underline{R}_{b}^{0}(\widetilde{B}^{n})_{i-1}+\delta_{t}q_{\mathrm{P},1}^{n}\cdot\underline{\mathfrak{s}}_{i-1}^{(b,0)}+\delta_{t}q_{\mathrm{P},N}^{n}\cdot\underline{\mathfrak{s}}_{i-1}^{(b,\ell)}\\\\[7.5pt]
\displaystyle\hskip
112.50017pt{}-\frac{1}{2}\Big{(}\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}\Big{)}\mathfrak{f}_{i-1}^{n}\end{cases},$

where we kept the same definitions for
$\underline{R}_{b}^{0},\underline{R}_{b}^{1},\mathfrak{f}^{n},\widetilde{B}^{n},\underline{\mathfrak{s}}^{(b,0)},\underline{\mathfrak{s}}^{(b,\ell)},\alpha_{b}$
as in the previous Section, which involve the second order centered operator
$\delta_{x}$ given in (67). 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_{\mathrm{P},1}^{n+1})=g_{0}^{n+1},\qquad\xi_{\ell}^{-}(U_{\mathrm{P},N}^{n+1})=g_{\ell}^{n+1},$

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

$\delta_{t}q_{\mathrm{P},1}^{n}=\frac{q_{\mathrm{P},1}^{n+1}-q_{1}^{n}}{\Delta
t},\qquad\delta_{t}q_{\mathrm{P},N}^{n}=\frac{q_{\mathrm{P},N}^{n+1}-q_{N}^{n}}{\Delta
t}.$

Notice that, although in (85) 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 (75), where
$(\underline{h_{b}S})^{*}$ defined in (69) 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_{\mathrm{C},i}^{n+1}=(\zeta_{\mathrm{C},i}^{n+1},q_{\mathrm{C},i}^{n+1})^{T}$
is computed using a forward Euler update with right upwinding in the interior
cells $2\leq i\leq N-1$. It involves the predicted state obtained in the
previous step in order to evaluate the flux and source terms:

(86)
$\displaystyle\begin{cases}\displaystyle\frac{\zeta_{\mathrm{C},i}^{n+1}-\zeta_{i}^{n}}{\Delta
t}+\frac{1}{\Delta
x}\Big{(}q_{\mathrm{P},i+1}^{n+1}-q_{\mathrm{P},i}^{n+1}\Big{)}=0\\\\[7.5pt]
\displaystyle\frac{q_{\mathrm{C},i}^{n+1}-q_{i}^{n}}{\Delta
t}+\frac{\mathfrak{f}_{\mathrm{P},i+1}^{n+1}-\mathfrak{f}_{\mathrm{P},i}^{n+1}}{\Delta
x}=\underline{R}_{b}^{0}(\widetilde{B}_{\mathrm{P}}^{n+1})_{i+1}+\delta_{t}q_{\mathrm{C},1}^{n}\cdot\underline{\mathfrak{s}}_{i+1}^{(b,0)}+\delta_{t}q_{\mathrm{C},N}^{n}\cdot\underline{\mathfrak{s}}_{i+1}^{(b,\ell)}\\\\[7.5pt]
\displaystyle\hskip
130.0002pt{}-\frac{1}{2}\Big{(}\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}\Big{)}\mathfrak{f}_{\mathrm{P},i+1}^{n+1}\end{cases}\\!.$

In the above, the fluxes $\mathfrak{f}_{\mathrm{P},i}^{n+1}$ and the source
term $\widetilde{B}_{\mathrm{P},i}^{n+1}$ are defined as in (79) and (74) but
making use of the prediction state, that is to say

$\displaystyle\forall 1\leq i\leq
N,\qquad\begin{cases}\displaystyle\widetilde{B}_{\mathrm{P},i}^{n+1}&=\displaystyle-{\mathtt{g}}\zeta_{\mathrm{P},i}^{n+1}\delta_{x}(b)_{i}+\frac{3}{2}\frac{\delta_{x}(b)_{i}}{h_{b,i}}\widetilde{f}_{\mathrm{NSW}}(U_{\mathrm{P},i}^{n+1}),\\\
\displaystyle\mathfrak{f}_{\mathrm{P},i}^{n+1}&=\displaystyle\frac{h_{b,i}^{2}}{\alpha_{b,i}}\underline{R}_{b}^{1}\Big{(}\Big{\\{}\frac{1}{h_{b,j}^{2}}\widetilde{f}_{\mathrm{NSW}}(U_{\mathrm{P},j}^{n+1})\Big{\\}}_{1\leq
j\leq N}\Big{)}_{i}.\end{cases}$

The handling of the border cells is the same as before, except that in (83),
for $k\in\\{1,2\\}$ the matrices $D_{k},\widetilde{D}_{k}$ and the quadratic
forms $\mathcal{Q}_{k}^{n}$ need to be evaluated with the predicted values
$(X,Y)_{\mathrm{P}}^{n+1}$ obtained in the previous step instead of
$(X,Y)^{n}$, and the vectors $\widetilde{V}_{\rm bdry}^{n},\widetilde{V}_{\rm
int}^{n}$ need to be evaluated with $U_{\mathrm{P}}^{n+1}$ instead of $U^{n}$.
Therefore the correction update for the output functions can be written

(87) $\displaystyle\begin{cases}\cfrac{X_{\mathrm{C}}^{n+1}-X^{n}}{\Delta
t}=Y_{\mathrm{C}}^{n+1}\\\
\widetilde{D}_{\mathrm{P},1}^{n+1}\cfrac{Y_{\mathrm{C}}^{n+1}-Y^{n}}{\Delta
t}+\underline{\mathfrak{S}}^{\prime}_{b}\widetilde{D}_{\mathrm{P},2}^{n+1}Y_{\mathrm{C}}^{n+1}+\begin{pmatrix}{\mathcal{Q}}_{\mathrm{P},1}^{n+1}(\frac{g_{0}^{n+1}-g_{0}^{n-1}}{2\Delta
t},Y_{\mathrm{P},1}^{n+1})\\\
{\mathcal{Q}}_{\mathrm{P},2}^{n+1}(Y_{\mathrm{P},2}^{n+1},\frac{g_{\ell}^{n+1}-g_{\ell}^{n-1}}{2\Delta
t})\end{pmatrix}+\widetilde{V}_{\rm bdry,P}^{n+1}\\\ \hskip
60.00009pt=\widetilde{V}_{\rm
int,P}^{n+1}-{D}_{\mathrm{P},1}^{n+1}\begin{pmatrix}\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}}\end{pmatrix}-\underline{\mathfrak{S}}^{\prime}_{b}{D}_{\mathrm{P},2}^{n+1}\begin{pmatrix}\frac{g_{0}^{n+1}-g_{0}^{n-1}}{2\Delta
t}\\\ \frac{g_{\ell}^{n+1}-g_{\ell}^{n-1}}{2\Delta
t}\end{pmatrix}\end{cases}\\!,$

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

$\delta_{t}q_{\mathrm{C},1}^{n}=\frac{q_{\mathrm{C},1}^{n+1}-q_{1}^{n}}{\Delta
t},\qquad\delta_{t}q_{\mathrm{C},N}^{n}=\frac{q_{\mathrm{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\leq i\leq N-1$, and likewise the output
functions are averaged in the border cells $i=1,N$:

(88)
$\begin{cases}\zeta_{i}^{n+1}=\cfrac{1}{2}(\zeta_{\mathrm{P},i}^{n+1}+\zeta_{\mathrm{C},i}^{n+1})\\\
q_{i}^{n+1}=\cfrac{1}{2}(q_{\mathrm{P},i}^{n+1}+q_{\mathrm{C},i}^{n+1})\end{cases},\quad\begin{cases}\xi_{0}^{-}(U_{1}^{n+1})=\cfrac{1}{2}(\xi_{0}^{-}(U_{\mathrm{P},1}^{n+1})+\xi_{0}^{-}(U_{\mathrm{C},1}^{n+1}))\\\
\xi_{\ell}^{+}(U_{N}^{n+1})=\cfrac{1}{2}(\xi_{\ell}^{+}(U_{\mathrm{P},N}^{n+1})+\xi_{\ell}^{+}(U_{\mathrm{C},N}^{n+1}))\end{cases}$

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 (53) 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 (83), (84) and (87) are already
second order in space thanks to the use of the discrete operator $\delta_{x}$
given by (67) to discretize $\mathfrak{S}_{b}^{\prime},\widetilde{V}_{\rm
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. 1.

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

$\xi_{0}^{+}(\zeta,q)=\xi_{\ell}^{-}(\zeta,q)=\zeta,\qquad\xi_{0}^{-}(\zeta,q)=\xi_{\ell}^{+}(\zeta,q)=q;$

  2. 2.

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

$\xi_{0}^{+}(\zeta,q)=\xi_{\ell}^{-}(\zeta,q)=q,\qquad\xi_{0}^{-}(\zeta,q)=\xi_{\ell}^{+}(\zeta,q)=\zeta;$

  3. 3.

Input functions given by the incoming Riemann invariants, output functions
given by the outgoing Riemann invariants

$\begin{cases}\xi_{0}^{+}(\zeta,q)=u+2\sqrt{\mathtt{g}h}\\\
\xi_{\ell}^{-}(\zeta,q)=u-2\sqrt{\mathtt{g}h}\end{cases},\qquad\begin{cases}\xi_{0}^{-}(\zeta,q)=u-2\sqrt{\mathtt{g}h}\\\
\xi_{\ell}^{+}(\zeta,q)=u+2\sqrt{\mathtt{g}h}\end{cases}.$

The gravitational acceleration $\mathtt{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 (3) of the form

(89) $\zeta(t,x)=\widetilde{\zeta}(x-x_{0}-ct),\qquad
q(t,x)=\widetilde{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 $\widetilde{\zeta}$ verifies
the following second order ODE

(90)
$-c^{2}\frac{h_{0}\widetilde{\zeta}}{h_{0}+\widetilde{\zeta}}+\frac{c^{2}h_{0}^{2}}{3}\widetilde{\zeta}^{\prime\prime}+\frac{\mathtt{g}}{2}(\widetilde{\zeta}^{2}+2h_{0}\widetilde{\zeta})=0.$

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

$c^{2}=\frac{\frac{\mathtt{g}}{6}\zeta_{\max}^{3}+\frac{\mathtt{g}h_{0}}{2}\zeta_{\max}^{2}}{h_{0}^{2}(\zeta_{\max}/h_{0}-\ln(1+\zeta_{\max}/h_{0}))}$

by integrating (90) against $\widetilde{\zeta}^{\prime}$, using again that
$\widetilde{\zeta}$ cancels at infinity and that
$\widetilde{\zeta}^{\prime}(0)=0$.

Once the solution of (90) 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^{\rm in},q^{\rm
in})(x)=(\widetilde{\zeta},c\,\widetilde{\zeta})(x-x_{0})$ and the general
boundary conditions write

$\displaystyle\xi_{0}^{+}(\zeta_{0},q_{0})(t)$
$\displaystyle=g_{0}(t):=\xi_{0}^{+}(\widetilde{\zeta}(-x_{0}-ct),c\,\widetilde{\zeta}(-x_{0}-ct)),$
$\displaystyle\xi_{\ell}^{-}(\zeta_{\ell},q_{\ell})(t)$
$\displaystyle=g_{\ell}(t):=\xi_{\ell}^{-}(\widetilde{\zeta}(\ell-
x_{0}-ct),c\,\widetilde{\zeta}(\ell-x_{0}-ct)).$

In practice we take $\ell=60$, $h_{0}=1$, $\zeta_{\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: Free surface elevation for the incoming and outgoing solitary wave
test-cases (respectively top and bottom); incoming Riemann invariants imposed
at the boundaries

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

(91) $E_{\rm num}^{n}=\sqrt{\frac{\Delta
x}{\ell}}\bigg{(}\sum_{i=1}^{N}\big{(}\zeta_{i}^{n}-\widetilde{\zeta}(x_{i}-x_{0}-ct^{n})\big{)}^{2}\bigg{)}^{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.

$\Delta x$ |  $\zeta$ enforced |  $q$ enforced |  $R^{\pm}$ enforced
---|---|---|---
|  $L^{2}$-error | Order |  $L^{2}$-error | Order |  $L^{2}$-error | Order
8.82E-02 | 4.409E-03 | – | 4.728E-03 | – | 4.482E-03 | –
6.25E-02 | 3.263E-03 | 0.87 | 3.484E-03 | 0.88 | 3.312E-03 | 0.87
4.42E-02 | 2.387E-03 | 0.90 | 2.539E-03 | 0.91 | 2.420E-03 | 0.91
3.12E-02 | 1.727E-03 | 0.93 | 1.833E-03 | 0.94 | 1.750E-03 | 0.93
2.21E-02 | 1.244E-03 | 0.95 | 1.318E-03 | 0.95 | 1.260E-03 | 0.95
1.56E-02 | 8.906E-04 | 0.97 | 9.411E-04 | 0.97 | 9.015E-04 | 0.97
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. $\Delta x$ |  $\zeta$ enforced |  $q$ enforced |  $R^{\pm}$ enforced
---|---|---|---
|  $L^{2}$-error | Order |  $L^{2}$-error | Order |  $L^{2}$-error | Order
1.00E+00 | 7.846E-03 | – | 1.186E-02 | – | 9.424E-03 | –
5.00E-01 | 2.473E-03 | 1.67 | 3.789E-03 | 1.65 | 2.783E-03 | 1.76
2.50E-01 | 6.670E-04 | 1.89 | 1.026E-03 | 1.89 | 7.199E-04 | 1.95
1.25E-01 | 1.696E-04 | 1.98 | 2.662E-04 | 1.95 | 1.832E-04 | 1.97
6.25E-02 | 4.312E-05 | 1.98 | 6.804E-05 | 1.97 | 4.624E-05 | 1.99
3.12E-02 | 1.107E-05 | 1.96 | 1.731E-05 | 1.97 | 1.162E-05 | 1.99
Table 2: MacCormack scheme for the incoming solitary wave. $\Delta x$ |  $\zeta$ enforced |  $q$ enforced |  $R^{\pm}$ enforced
---|---|---|---
|  $L^{2}$-error | Order |  $L^{2}$-error | Order |  $L^{2}$-error | Order
6.25E-02 | 2.363E-03 | – | 2.666E-03 | – | 1.331E-04 | –
4.42E-02 | 1.785E-03 | 0.81 | 2.031E-03 | 0.78 | 1.070E-04 | 0.63
3.12E-02 | 1.322E-03 | 0.87 | 1.514E-03 | 0.85 | 8.401E-05 | 0.70
2.21E-02 | 9.725E-04 | 0.89 | 1.119E-03 | 0.87 | 6.503E-05 | 0.74
1.56E-02 | 7.053E-04 | 0.93 | 8.145E-04 | 0.92 | 4.939E-05 | 0.79
1.11E-02 | 5.098E-04 | 0.94 | 5.903E-04 | 0.93 | 3.711E-05 | 0.83
Table 3: Lax-Friedrichs scheme for the outgoing solitary wave. $\Delta x$ |  $\zeta$ enforced |  $q$ enforced |  $R^{\pm}$ enforced
---|---|---|---
|  $L^{2}$-error | Order |  $L^{2}$-error | Order |  $L^{2}$-error | Order
1.00E+00 | 6.430E-03 | – | 1.031E-02 | – | 4.287E-03 | –
5.00E-01 | 1.885E-03 | 1.77 | 4.437E-03 | 1.22 | 1.993E-03 | 1.11
2.50E-01 | 5.568E-04 | 1.76 | 1.257E-03 | 1.82 | 5.607E-04 | 1.83
1.25E-01 | 1.459E-04 | 1.93 | 3.252E-04 | 1.95 | 1.318E-04 | 2.09
6.25E-02 | 3.674E-05 | 1.99 | 8.223E-05 | 1.98 | 2.826E-05 | 2.22
3.12E-02 | 9.162E-06 | 2.00 | 2.067E-05 | 1.99 | 7.085E-06 | 2.00
Table 4: MacCormack scheme for the outgoing solitary wave.

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
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 (23) 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:

$\displaystyle h_{0}=1,\quad\beta=1/4,\quad A=h_{0}/4.$

The bathymetry profile is defined as follows

(92) $\displaystyle b(x)=\beta h_{0}\cdot\begin{cases}1&\text{if }\lvert
x-\ell/2\rvert<\frac{\ell}{4}\\\ \sigma(1.5-\frac{2}{\ell}\lvert
x-\ell/2\rvert)&\text{if }\frac{\ell}{4}\leq\lvert
x-\ell/2\rvert\leq\frac{3\ell}{4}\\\ 0&\text{otherwise}\end{cases},$

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

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

Recalling the shallowness parameter $\mu=(2\pi h_{0})^{2}/\lambda^{2}$ with
the characteristic wavelength defined here as $\lambda=\mathrm{length}\,\\{x,\
\zeta^{\rm 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$.

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.

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 (93)
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$.

$\Delta x$ |  $\zeta$ enforced |  $q$ enforced |  $R^{\pm}$ enforced
---|---|---|---
|  $L^{2}$-error | Order |  $L^{2}$-error | Order |  $L^{2}$-error | Order
3.26E-02 | 1.017E-02 | – | 1.300E-02 | – | 8.545E-03 | –
2.30E-02 | 7.862E-03 | 0.74 | 1.016E-02 | 0.71 | 6.604E-03 | 0.74
1.63E-02 | 5.962E-03 | 0.80 | 7.759E-03 | 0.78 | 5.003E-03 | 0.80
1.15E-02 | 4.440E-03 | 0.85 | 5.802E-03 | 0.84 | 3.721E-03 | 0.85
8.14E-03 | 3.265E-03 | 0.89 | 4.278E-03 | 0.88 | 2.733E-03 | 0.89
5.75E-03 | 2.375E-03 | 0.92 | 3.118E-03 | 0.91 | 1.986E-03 | 0.92
Table 5: Lax-Friedrichs scheme for the Gaussian over bump test-case. $\Delta x$ |  $\zeta$ enforced |  $q$ enforced |  $R^{\pm}$ enforced
---|---|---|---
|  $L^{2}$-error | Order |  $L^{2}$-error | Order |  $L^{2}$-error | Order
2.60E-01 | 4.635E-03 | – | 5.708E-03 | – | 3.301E-03 | –
1.84E-01 | 2.470E-03 | 1.81 | 3.076E-03 | 1.78 | 1.690E-03 | 1.92
1.30E-01 | 1.288E-03 | 1.89 | 1.629E-03 | 1.84 | 8.485E-04 | 2.00
9.19E-02 | 6.575E-04 | 1.93 | 8.402E-04 | 1.90 | 4.095E-04 | 2.09
6.51E-02 | 3.384E-04 | 1.93 | 4.316E-04 | 1.93 | 1.930E-04 | 2.18
4.60E-02 | 1.777E-04 | 1.85 | 2.211E-04 | 1.92 | 9.056E-05 | 2.17
Table 6: MacCormack scheme for the Gaussian over bump test-case.

Comparing the numerical solution obtained in the small domain with different
mesh sizes to the reference solution, we can compute the $\ell^{2}$ error (91)
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 (23) 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 (92) for the bottom, however
the characteristic depth and wave amplitude are now scaled with respect to a
parameter $K_{\mu}$ as follows:

$\displaystyle 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:

(94) $\widetilde{\zeta}^{\rm in}(x)=A\cdot\sin(2\pi
x/\lambda),\qquad\widetilde{q}^{\rm in}(x)=2\cdot\widetilde{\zeta}^{\rm
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{\mathtt{g}h_{0}}$.

Figure 3: Initial elevation for reference and small domain solutions
($\mu=1$).

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^{\rm in},q^{\rm 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 (60), which we recall
provides compatibility conditions when enforcing $\zeta$, the reference
solution $(\widetilde{\zeta},\widetilde{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:

$\displaystyle\xi_{0}^{+}(\zeta_{0},q_{0})(t)=g_{0}(t):=\sigma(5t/T_{0})\cdot\xi_{0}^{+}(\widetilde{\zeta}(t,0),\widetilde{q}(t,0)),$
$\displaystyle\xi_{\ell}^{-}(\zeta_{\ell},q_{\ell})(t)=g_{\ell}(t):=\sigma(5t/T_{0})\cdot\xi_{\ell}^{-}(\widetilde{\zeta}(t,\ell),\widetilde{q}(t,\ell)).$

If on the contrary the boundary conditions on the elevation were enforced
directly at $t=0$, then system (60) would be violated, and the well-posedness
result from Theorem 3.8 does not hold anymore in this situation.

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 (91) 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 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
$\widetilde{\zeta}$ at final time $t=50\,T_{0}$ for various boundary
conditions. The last row represents the $\ell^{2}$ error (91) over time. Out
of the three boundary conditions investigated here, only enforcing the
incoming Riemann invariants allows this error to decay towards zero.

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 (91) 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.

##  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