Generating boundary conditions for a Boussinesq system

TL;DR

This paper introduces a novel numerical method for implementing boundary conditions in a 1D Boussinesq system, effectively handling dispersive boundary layers without increasing computational cost.

Contribution

A new reformulation-based approach for generating boundary conditions in Boussinesq systems that avoids additional computational overhead.

Findings

Method accurately handles dispersive boundary layers.

Validated through multiple numerical simulations.

No increase in computational time compared to periodic boundaries.

Abstract

We present a new method for the numerical implementation of generating boundary conditions for a one dimensional Boussinesq system. This method is based on a reformulation of the equations and a resolution of the dispersive boundary layer that is created at the boundary when the boundary conditions are non homogeneous. This method is implemented for a simple first order finite volume scheme and validated by several numerical simulations. Contrary to the other techniques that can be found in the literature, our approach does not cause any increase in computational time with respect to periodic boundary conditions.

Figures

Equations

$$\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\ \partial_{t}q+\partial_{x}\big{(}\frac{1}{2}gh^{2}+\frac{1}{h}q^{2}\big{)}=0,\end{cases}$$

$$\zeta (t,x=0)=f(t),\text{mbox}forallt\ge 0,$$

$$(\zeta ,q)(t=0,x)=({\zeta }^0,q^0)(x),\text{mbox}forallx\ge 0;$$

$$\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\ (1-\frac{H_{0}^{2}}{3}\partial_{x}^{2})\partial_{t}q+\partial_{x}\big{(}\frac{1}{2}gh^{2}+\frac{1}{h}q^{2}\big{)}=0,\qquad(h=H_{0}+\zeta)\end{cases}$$

$$\begin{cases}\displaystyle\partial_{t}\zeta+\partial_{x}q=0,\\ \displaystyle\partial_{t}q+\partial_{x}R_{1}\big{(}\frac{1}{2}g(h^{2}-H_{0}^{2})+\frac{1}{h}q^{2}\big{)}=\underline{{\mathcal{Q}}}(\underline{q},f,\ddot{f},\zeta,q)\exp\big{(}-\sqrt{3}\frac{x}{H_{0}}\big{)},\end{cases}$$

$$\underline{\mathcal{Q}}(\underline{q},f,\ddot{f},\zeta ,q)=$$

$$\displaystyle-\frac{\sqrt{3}}{H_{0}}\underline{R}_{1}\big{(}\frac{1}{2}g(h^{2}-H_{0}^{2})+\frac{1}{h}q^{2}\big{)},$$

$$(\zeta ,q)(t=0,x)=({\zeta }^0,q^0)(x),\underline{q}(t=0)=q^0(x=0),$$

$$\zeta (t,x=0)=f(t);$$

$$x_0=0,x_1={\delta }_x,\dots ,x_i=i{\delta }_x,\dots ,x_{n_x-1}=(n_x-1){\delta }_x,x_{n_x}=L,$$

$$\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\ \partial_{t}q+\partial_{x}\big{(}\frac{1}{2}gh^{2}+\frac{1}{h}q^{2}\big{)}=0,\end{cases}\qquad(h=H_{0}+\zeta).$$

$$(\zeta ,q)(t=0,x)=({\zeta }^0,q^0)(x),\zeta (t,x=0)=f(t).$$

$$u=\frac{q}{h}\text{mbox}withh=H_0+\zeta .$$

$$\{\begin{array}{l}{\partial }_t\zeta +{\partial }_x(hu)=0,\\ {\partial }_{t}u+g{\partial }_x\zeta +u{\partial }_{x}u=0\end{array}$$

$$\partial_{t}U+A(U)\partial_{x}U=0\quad\mbox{ with }\quad U=(\zeta,u)^{T}\quad\mbox{ and }\quad A(U)=\left(\begin{array}[]{cc}u&h\\ g&u\end{array}\right).$$

$$\lambda_{\pm}(U)=\pm u+\sqrt{gh}\quad\mbox{ and }\quad{\bf e}_{\pm}(U)=\big{(}\sqrt{\frac{g}{h}},\pm 1\big{)}^{T}.$$

$$\big{(}\sqrt{\frac{g}{h}}\partial_{t}h\pm\partial_{t}u\big{)}\pm(\pm u+\sqrt{gh})\big{(}\sqrt{\frac{g}{h}}\partial_{x}h\pm\partial_{x}u\big{)}=0.$$

$$R_{\pm}(U):=2\big{(}\sqrt{gh}-\sqrt{gH_{0}})\pm u,$$

$${\partial }_t{R}_{+}+{\lambda }_{+}(U){\partial }_x{R}_{+}=0,{\partial }_t{R}_{-}-{\lambda }_{-}(U){\partial }_x{R}_{-}=0.$$

$$\partial_{t}U+\partial_{x}\big{(}F(U)\big{)}=0\quad\mbox{ with }\quad F=\big{(}q,\frac{1}{2}g(h^{2}-H_{0}^{2})+\frac{1}{h}q^{2}\big{)}^{T},$$

$$\frac{U_i^{n+1}-U_i^n}{{\delta }_t}+\frac{1}{{\delta }_x}(F_{i+1/2}^n-F_{i-1/2}^n)=0,i\ge 1,$$

$$F_{i-1/2}^n=\frac{1}{2}(F_i^n+F_{i-1}^n)-\frac{{\delta }_x}{2{\delta }_t}(U_i^n-U_{i-1}^n),i\ge 1,$$

$$({\zeta }_i^0,q_i^0)=({\zeta }^0,q^0)(x_i)(i\ge 1),{\zeta }_0^n=f^n:=f(t^n),$$

$${\zeta }_0^n=f^n,$$

$$q=\frac{h}{2}(R_{+}-R_{-})\quad\mbox{ and }\quad R_{+}+R_{-}=4\big{(}\sqrt{gh}-\sqrt{gH_{0}})$$

$$q=h\big{(}2(\sqrt{gh}-\sqrt{gH_{0}})-R_{-}\big{)}.$$

$$\underline{q}=(H_{0}+f)\big{(}2\big{(}\sqrt{g(H_{0}+f)}-\sqrt{gH_{0}}\big{)}-{R_{-}}_{|_{x=0}}\big{)}$$

$$q^{n}_{0}=(H_{0}+f^{n})\big{(}2(\sqrt{g(H_{0}+f^{n})}-\sqrt{gH_{0}})-R_{-,0}^{n}\big{)}.$$

$$\frac{R_{-,0}^n-R_{-,0}^{n-1}}{{\delta }_t}-{\lambda }_{-}^{n-1}\frac{R_{-,1}^{n-1}-R_{-,0}^{n-1}}{{\delta }_x}=0;$$

$${\lambda }_{-}^{n-1}=(1-{\alpha }^{n-1}){\lambda }_{-}(U_0^{n-1})+{\alpha }^{n-1}{\lambda }_{-}(U_1^{n-1})$$

$$R_{-,0}^n=(1-{\lambda }_{-}^{n-1}\frac{{\delta }_t}{{\delta }_x})R_{-,0}^{n-1}+{\lambda }_{-}^{n-1}\frac{{\delta }_t}{{\delta }_x}{R}_{-,1}^{n-1},$$

Full text

D. L. is supported by the Fondation Del Duca de l’Académie des Sciences, the ANR grants ANR-17-CE40-0025 NABUCO and ANR-18-CE40-0027-01 Singflows, and the Conseil Régional d’Aquitaine.

1. Introduction

1.1. General setting

Among the many reduced models used to describe the evolution of waves at the surface of a fluid in shallow water, the nonlinear shallow water equations are certainly one of the most used for applications. They can be written in conservative form as

[TABLE]

where $\zeta$ is the surface elevation above the rest state and $q$ the horizontal discharge (equivalently, the vertical integral of the horizontal velocity), and where $h=H_{0}+\zeta$ is the total water depth ($H_{0}$ being the depth at rest) and $g$ the acceleration of gravity.

For many applications, the surface elevation is known at the entrance of the domain (through buoy measurements for instance, or it can be provided by offshore swell models) and is imposed as a boundary condition for the model

[TABLE]

as well as the initial values for $q$ and $\zeta$ in the domain,

[TABLE]

this type of boundary condition is often referred to as generating boundary condition (see for instance [13]). It is used a lot in coastal oceanography, where the offshore swell is imposed at the entrance of the domain of interest (see for instance the classical benchmark [12]).

The problem consisting in solving (1) together with (2) and (3) is a mixed initial-boundary value problem (IBVP); due to its hyperbolic nature, it can be solved theoretically (see for instance [23, 25], and more recently [16] for sharp well-posedness results). From the numerical viewpoint, solving this IBVP is also possible, using the decomposition of the solution into Riemann invariants (see for instance [24]).

The nonlinear shallow water equations provide a robust model used in many applications; it is known [2, 15, 18, 19] to provide an approximation of the full free surface Euler equations with a precision $O(\mu)$, where $\mu=H_{0}^{2}/L^{2}$ is the shallowness parameter ($L$ denotes here the typical horizontal scale of the waves). It omits however the dispersive effects that play an important role in coastal areas, in particular in the shoaling zone; in order to take them into account, one has to keep the $O(\mu^{2})$ terms that are neglected in the derivation of the nonlinear shallow water equations. The most simple models that reach such a precision and therefore take into account the dispersive effects while retaining nonlinear terms are the so-called Boussinesq models. There are actually many asymptotically equivalent Boussinesq models [9, 10, 8]; their simplicity is due to the fact that they are derived under the assumption that the waves are of small amplitude compared to the depth, which allows to neglect some of the nonlinear terms (without this assumption, one has to work with the much more complicated Serre-Green-Naghdi equations, see [18] for instance). We choose here to work with the so-called Boussinesq-Abbott model [1, 14] since its structure is obviously a dispersive perturbation of the nonlinear shallow water equations,

[TABLE]

(removing the dispersive term $-\frac{H_{0}^{2}}{3}\partial_{x}^{2}\partial_{t}q$, this model reduces to (1)). As above, we are interested in the initial-boundary value problem for this system, we therefore complement it with the boundary condition (2) and initial condition (3). Contrary to (1), this system is no longer hyperbolic, and there is no general theory to address the IBVP. Only some particular cases have been considered, such as in [28] with homogeneous boundary conditions, [7, 3] for a particular class of Boussinesq systems (the Bona-Smith family) where a regularizing dispersion is also present in the first equation, [22] for a higher order Boussinesq system or [21] for the shoreline problem (vanishing depth).

Due to its importance for numerical simulations in coastal oceanography, there has been a significant amount of work devoted to finding numerical answers to this issue in recent years. For the related problem of transparent boundary conditions in particular (i.e. which boundary conditions should be put at the boundary of the computational domain so that waves can pass through this artificial boundary without being affected by it), the linear problem has been considered for scalar equations (such as KdV or BBM) in [4, 5] as well as for the linearization of (4) around the rest state. For the nonlinear case, a different approach has been used recently, which consists in implementing a perfectly matched layer (PML) approach for a hyperbolic relaxation of the Green-Naghdi equations [17]. This approach can be used to deal with generating boundary conditions such as (2) but the size of the layer in which the PML approach is implemented is typically of two wavelength, which for applications to coastal oceanography can typically represent an increase of $100\%$ of the computational domain. Other methods such as the source function method [27] also require a significant increase of computational time.

Our goal in this note is to propose a new approach to the nonlinear Boussinesq system (4) with generating boundary condition (2) and initial condition (3), and which does not require any extension of the computational domain. It is based in a reformulation of the problem (4) and (2)-(3) into a non-homogeneous system of conservation laws for $\zeta$ and $q$ with a nonlocal flux, and with a source term accounting for the dispersive boundary layer,

[TABLE]

where $\underline{q}=q_{|_{x=0}}$ and

[TABLE]

with the initial condition

[TABLE]

and the boundary condition

[TABLE]

here, we denoted by $R_{1}$ the inverse of the operator $(1-\frac{H_{0}^{2}}{3}\partial_{x}^{2})$ on ${\mathbb{R}}_{+}$ with homogeneous Neumann boundary condition, and $\underline{R}_{1}f=(R_{1}f)_{|_{x=0}}$ (see Definition 1 below). The source term in the equation for $q$ is a dispersive boundary layer that appears because the time derivative of the trace $\underline{q}$ of $q$ at $x=0$ does not necessarily vanish.

We propose here a simple numerical scheme based on this new formulation of the problem (which we prove to be well-posed), very easy to implement and that does not require to work on an extended computational domain. The ability of this method to handle generating boundary conditions with a very good precision is illustrated by several computations where nonlinear and dispersive terms both play an important role.

1.2. Organization of the paper

We describe in this paper how to handle a generating boundary condition on the left border of the computational domain. For the sake of clarity, we consider the problem on the half line $(0,\infty)$ so that we do not have to deal with boundary conditions at the right boundary; for our numerical simulations, we either consider a domain that is large enough so that the influence of the right boundary condition is negligible, or take a wall boundary condition $q=0$ on the right boundary (of course, a generating boundary condition on the right-boundary can be handled by a straightforward adaptation of what is done at the left boundary).

In Section 2 we briefly recall the theory and numerical simulation of generating boundary conditions for the nonlinear shallow water equations in order to make clear that different mechanisms are at stake in the hyperbolic (shallow water) and dispersive (Boussinesq) cases. Note in particular that for the hyperbolic case considered in this section, the missing data at the boundary (i.e. the trace of the discharge $q$ at $x=0$) is deduced from the value at the boundary of the outgoing Riemann invariant which can itself be determined in terms of interior values by solving the characteristic equation.

In Section 3 we consider generating boundary conditions for the Boussinesq system (4). In order to make more apparent the structure of the dispersive boundary layer we shall construct, we first non-dimensionalize the equations in §3.1. The main step of the analysis is performed in §3.2 where the dispersive boundary layer is constructed and the reformulation (5) of the problem is derived. This reformulation is used in §3.3 to prove the local well-posedness of the initial-boundary value problem for the Boussinesq system (4) with generating boundary condition (2) which, to our knowledge, was not known so far. A discretization of the reformulation (5) is then proposed in §3.4; for the sake of clarity, it is based on the standard Lax-Friedrichs scheme. It is not possible to recover the missing boundary data using Riemann invariants as in the hyperbolic case, but the nonlocality of the operator $R_{1}$ allows us to express this missing information in terms of interior values. We insist on the fact that our numerical treatment of the generating boundary condition does not increase the computational time compared to simple boundary conditions (periodic, physical well, etc.), contrary to the previously used approaches mentioned in the introduction.

Finally, we provide in Section 4 several numerical computations showing the accuracy of our numerical scheme. Our validation method is first presented in §4.1; it consists in computing a reference solution in a large domain $[-L,L]$ with a very refined mesh, and to use the values of the water elevation at $x=0$ provided by this solution as a generating boundary condition for computations on the small domain $[0,L]$. The accuracy of this new solution is measured by comparing it with the reference solution. A first example is provided in §4.2 in a situation where both incoming and outgoing waves are present. In §4.3 we show that the Boussinesq system (4) admits solitary waves, and that we are able to generate them with good accuracy using the corresponding generating boundary condition. We then provide in §4.4 another example, relevant for applications to coastal oceanography [6, 26, 20], which consists in the generation of a sinusoidal wave train.

1.3. Notations

For the numerical computations, the computational domain $[0,L]$ is discretized using a uniform grid:

[TABLE]

with $\delta_{x}=\displaystyle{\frac{L}{n_{x}}}$. The time step is denoted $\delta_{t}$. The variables $\zeta_{i}^{n}$ and $q^{n}_{i}$ denote the values of the numerical solution for $\zeta$ and $q$ at the time $n\,\delta_{t}$ and at the location $x_{i}$. Generally speaking, the subscript $i$ and the superscript $n$ indicate respectively a discretization at the location $x_{i}$ and at the time $n\,\delta_{t}$.

2. The nonlinear shallow water equations

We recall that the nonlinear shallow water (or Saint-Venant) equations are a system of equations coupling the surface elevation $\zeta$ above the rest state to the horizontal discharge $q$,

[TABLE]

This system of equations is complemented by the initial and boundary conditions

[TABLE]

For the sake of completeness and as a basis for comparisons with the dispersive (Boussinesq) case, we briefly recall here how this problem can be handled numerically.

2.1. The Riemann invariants

The nonlinear shallow water equations (9) can be written under an equivalent quasilinear form by introducing the vertically averaged horizontal velocity $u$,

[TABLE]

The resulting system of equations on $(\zeta,u)$ is given by

[TABLE]

or, in more condensed form,

[TABLE]

The matrix $A(U)$ is diagonalizable with eigenvalues $\lambda_{+}(U)$ and $-\lambda_{-}(U)$ and associated left-eigenvectors ${\bf e}_{\pm}(U)$ given by

[TABLE]

Taking the scalar product of (12) and ${\bf e}_{\pm}(U)$, we obtain

[TABLE]

This leads us to introduce the Riemann invariants $R_{\pm}$ as

[TABLE]

which satisfy the transport equations

[TABLE]

These Riemann invariants play a central role in the numerical resolution of the IBVP (9)-(10) presented below.

2.2. The discrete equations

Writing $U=(\zeta,q)^{T}$, we first write (9) in the condensed form,

[TABLE]

for which a finite volume type discretization gives

[TABLE]

the choice of $F^{n}_{i+1/2}$ depending on the numerical scheme. Our focus here being on explaining how to handle the boundary condition (10), we consider here the most simple case of the Lax-Friedrichs scheme where the discrete flux is given by

[TABLE]

with $F^{n}_{i}=F(U^{n}_{i})$. For $i=1$, this equation involves $U^{n}_{0}$ that we need to express in terms of $U^{n}=(U^{n}_{i})_{1\leq i}$ and the initial-boundary condition (10), which, in discrete form, reads

[TABLE]

with $t^{n}:=n\delta_{t}$; this is done in the following section.

2.3. Data on the water depth on the left boundary

For $i=1$, the flux $F_{1/2}$ requires the knowledge of $U^{n}_{0}=(\zeta^{n}_{0},q^{n}_{0})$. From the initial-boundary condition (18), one takes

[TABLE]

but we need to determine $q_{0}^{n}$, which can be deduced from the knowledge of $R^{n}_{\pm,0}:=R_{\pm}(t^{n},0)$. From (13) one gets indeed

[TABLE]

and therefore

[TABLE]

Evaluating this relation at $x=0$ provides an expression for the trace $\underline{q}=q_{|_{x=0}}$ in terms of the boundary data $f=\zeta_{|_{x=0}}$ and of the trace of the outgoing Riemann invariant $R_{-}$,

[TABLE]

and at the discrete level, we get at $t=t^{n}$

[TABLE]

Therefore, we just need to determine $R_{-,0}^{n}$ in order to determine $q_{0}^{n}$. We use the characteristic equation (14) satisfied by $R_{-}$; after discretization, this gives

[TABLE]

as in [24], $\lambda^{n-1}_{-}$ is computed as a linear interpolation between $\lambda_{-}(U^{n-1}_{0})$ and $\lambda_{-}(U^{n-1}_{1})$,

[TABLE]

and $0\leq\alpha^{n-1}\leq 1$ computed such that $\lambda^{n-1}_{-}\,\delta_{t}=\alpha^{n-1}\,\delta_{x}$. Therefore

[TABLE]

which gives $R_{-,0}^{n}$ in terms of its values at the previous time step and in terms of interior points.

3. The Boussinesq equations

We consider here the following Boussinesq-Abbott system [1, 14], which includes the dispersive effects neglected by the nonlinear shallow water equations (1)

[TABLE]

complemented with the initial and boundary conditions

[TABLE]

The key step in our analysis is the reformulation of this IBVP into a system of two conservation laws with nonlocal flux and a source term accounting for the presence of a dispersive boundary layer, and whose coefficient is found through the resolution of a nonlinear ODE.

In order to make clearer the structure of the dispersive boundary layer, we work with a dimensionless version of (4). The non-dimensionalization is performed in §3.1. The reformulation of the equations is then derived in §3.2 and a numerical scheme based on this newly exhibited structure is proposed in §3.4.

3.1. Dimensionless equations

Denoting by $a$ the typical amplitude of the waves, by $L$ its typical horizontal scale, we introduce the following dimensionless quantities, denoted with a prime,

[TABLE]

Replacing in (22) (and omitting the primes for the sake of clarity), we obtain the dimensionless version of the Boussinesq equations (1)

[TABLE]

where $\varepsilon$ and $\mu$ are respectively called nonlinearity and shallowness parameters and defined as

[TABLE]

the Boussinesq equations are derived in the shallow water, weakly nonlinear regime characterized by the assumptions

[TABLE]

Under these smallness assumptions, the Boussinesq model (24) provides an approximation consistent with the full free surface Euler equations up to $O(\mu^{2})$ and the convergence error is of order $O(\mu^{2}t)$ for times of order $O(1/\varepsilon)$ [2, 18, 19].

3.2. Reformulation of the equations

Solving the equations (24) on the full line requires the inversion of the operator $(1-\frac{\mu}{3}\partial_{x}^{2})$, which does not raise any difficulty. The situation is different here since we need to invert this operator on the half-line $(0,\infty)$, and we therefore need a boundary condition on $\partial_{t}q$ which is not directly at our disposal. Our strategy is, as in [11] for the description of the interaction of a floating objects with waves governed by a Boussinesq model, to use the inverse of the operator $(1-\frac{\mu}{3}\partial_{x}^{2})$ with homogeneous Dirichlet boundary condition, and to construct the dispersive boundary layer due to the fact that the boundary value $\underline{q}$ of $q$ is not equal to zero in general; we shall denote

[TABLE]

and we also need to define the Dirichlet and Neumann inverses of the operator $(1-\frac{\mu}{3}\partial_{x}^{2})$.

Definition 1**.**

We denote by $R_{0}$ and $R_{1}$ the inverse of the operator $(1-\frac{\mu}{3}\partial_{x}^{2})$ with homogeneous Dirichlet and Neumann boundary conditions respectively,

[TABLE]

and

[TABLE]

We also introduce the boundary operator $\underline{R}_{1}$ as

[TABLE]

Recalling that the ODE

[TABLE]

admits a unique solution in $H^{2}({\mathbb{R}}_{+})$ given by

[TABLE]

the second equation in (24) can be written equivalently under the form

[TABLE]

The last step is therefore to express $\dot{\underline{q}}$ in terms of the data $f={\zeta}_{|_{x=0}}$ of the problem. This is done in the following proposition.

Proposition 1**.**

If $(\zeta,q)$ are a smooth enough solution of (24), then the boundary value $\underline{q}$ of $q$ are related to the boundary value $f=\zeta_{|_{x=0}}$ and to the interior value of $\zeta$ and $q$ by solving the ODE

[TABLE]

where $\underline{R}_{1}$ is the boundary operator introduced in Definition 1.

Proof.

Differentiating (26) with respect to $x$, one obtains

[TABLE]

Lemma 1**.**

For all $g\in L^{2}({\mathbb{R}}_{+})$, the following identity holds,

[TABLE]

Proof of the lemma.

Just remark that if $v=R_{1}g$, then one easily gets from the definition of $R_{1}$ that

[TABLE]

so that, by definition of $R_{0}$, one has $\partial_{x}v=R_{0}\partial_{x}g$ (note that by classical variational arguments, $R_{0}$ is well defined as a mapping $\partial_{x}L^{2}({\mathbb{R}}_{+})\to H^{1}({\mathbb{R}}_{+})$). ∎

Using the first equation of (24) to substitute $\partial_{t}\partial_{x}q=-\partial_{t}^{2}\zeta$ and the lemma, one then deduces from (27) that

[TABLE]

Remarking further that $-\partial_{x}^{2}=\frac{1}{\delta^{2}}(1-\frac{\mu}{3}\partial_{x}^{2})-\frac{1}{\delta^{2}}$ and recalling that $(1-\frac{\mu}{3}\partial_{x}^{2})R_{1}=\mbox{Id}$, we obtain that

[TABLE]

Taking the trace of this expression at $x=0$ then yields

[TABLE]

from which the result follows. ∎

Using once again the lemma to replace $R_{0}\partial_{x}$ by $\partial_{x}R_{1}$ in (26), it follows from the above that the dimensionless Boussinesq equations (24) with initial and boundary conditions (23) can be equivalently written under the form

[TABLE]

where $\underline{q}=q_{|_{x=0}}$ and ${\mathfrak{f}}(\zeta,q)$ is the flux in the momentum equation for the nonlinear shallow water equations (1) in dimensionless variables,

[TABLE]

and

[TABLE]

with the initial condition

[TABLE]

and the boundary condition

[TABLE]

Remark 1*.*

Recalling that by definition of $R_{1}$, the trace of $\partial_{x}R_{1}{\mathfrak{f}}$ vanishes at $x=0$, one can take the trace at $x=0$ in the second equation in (28) to obtain the following evolution equation on $\underline{q}=q_{|_{x=0}}$,

[TABLE]

This relation has to be compared to (19) in the hyperbolic case, where $\underline{q}$ is given in terms of $f=\zeta_{|_{x=0}}$ and the trace of the outgoing Riemann invariant $R_{-}$. The mechanisms that allow to express $\underline{q}$ in terms of $f$ and interior values of $\zeta$ and $q$ are therefore completely different in the hyperbolic and in the dispersive cases: in the former, the decomposition into Riemann invariants is used to propagate information from the interior domain, while in the latter, this is done by using the non local nature of the operator $R_{1}$.

3.3. Well-posedness of the initial boundary value problem

As said in the introduction, very few results exist regarding the local well-posedness result for Boussinesq systems, except in some special cases such as [3, 28]. To our knowledge, no result exist so far for the Boussinesq-Abbott system considered here. Our reformulation (28)-(30) of this IBVP allows an easy proof of local well-posedness since it forms a simple ODE on $(\zeta,q)$ (here again, this is in strong contrast with the hyperbolic case where, of course, the equations cannot be put recast as an ODE).

Theorem 1**.**

*Let $f\in C^{2}({\mathbb{R}}_{+})$, $n\in{\mathbb{N}}\backslash\{0\}$ and $(\zeta^{0},q^{0})\in H^{n}({\mathbb{R}}_{+})\times H^{n+1}({\mathbb{R}}_{+})$ be such that $\inf(1+\varepsilon\zeta^{0})>0$. Then there exist $T>0$ and a unique solution $(\zeta,q)\in C^{1}([0,T];H^{n}({\mathbb{R}}_{+})\times H^{n+1}({\mathbb{R}}_{+}))$ to (28)-(30).

If moreover ${\zeta^{0}}_{|_{x=0}}=f(0)$ and $-{\partial_{x}q^{0}}_{|_{x=0}}=\dot{f}(0)$, then the boundary condition (31) is also satisfied for all times.*

Remark 2*.*

The existence time furnished by the theorem depends on $\varepsilon$ and $\mu$. The relevant time scale for the existence of the solution is $O(1/\varepsilon)$ in the case of the full line [2, 18]. Proving that such a time scale is also achieved in our case would require much more effort and an in depth analysis of the dispersive boundary layer together with additional compatibility conditions. Such a study is performed in [11] in the related problem of waves interaction with a floating object in the Boussinesq regime.

Proof.

To prove the first part of the theorem, it is enough to prove that (28)-(30) is actually an ODE on $H^{n}({\mathbb{R}}^{+})\times H^{n+1}({\mathbb{R}}^{+})$ meeting the requirements of the Cauchy-Lipschitz theorem.

With $U=(\zeta,q)^{\rm T}$ we can write the equations under the form

[TABLE]

By standard product estimates, $(\zeta,q)\in H^{n}\times H^{n+1}\mapsto{\mathfrak{f}}(\zeta,q)\in H^{n}$ is regular in a neighborhood of $(\zeta^{0},q^{0})$; moreover, $\partial_{x}R_{1}$ maps $H^{n}$ into $H^{n+1}$ by definition of $R_{1}$. It follows easily that ${\mathcal{F}}(t,U)$ is continuous and locally Lipschitz with respect to the second variable, so that we can apply Cauchy-Lipschitz theorem.

We now need to check that $\zeta(t,0)=f(t)$ for all time. In order to do so, one computes from the first equation in (28) that $\partial_{t}^{2}\zeta=-\partial_{t}\partial_{x}q$. Using the second equation to compute $\partial_{t}\partial_{x}q$ and taking the trace at $x=0$ one gets (proceeding as in the proof of Proposition 1) that

[TABLE]

This can be seen as a second order non-autonomous ODE on $\zeta_{|_{x=0}}$ with a right-hand side that is locally Lipschitz with respect to $\zeta_{|_{x=0}}$. There is therefore a unique solution to this ODE satisfying the initial conditions $\zeta_{|_{x=0}}(0)=f(0)$ and $\frac{d}{dt}(\zeta_{|_{x=0}})(0)=-\partial_{x}q^{0}(0)=\dot{f}(0)$. This solution is obviously given by $\zeta_{|_{x=0}}=f$, so that the proof is complete. ∎

3.4. Discretization of the equations

The goal of this section is to derive a numerical scheme to solve the initial boundary value problem (28)-(31).

3.4.1. A discrete version of the operators $R_{1}$ and $\underline{R}_{1}$

We still denote by $R_{1}$ the discrete inverse of the operator $(1-\frac{\mu}{3}\partial_{x}^{2})$ with homogeneous Neumann condition at the boundary. We use here a standard centered second order finite difference approximation for the discretization of $\partial_{x}^{2}$. More precisely, if $F=(f_{i})_{i\geq 1}$, we denote by $R_{1}F$ the vector $R_{1}F=V$ where $V=(v_{i})_{i\geq 1}$ is given by the resolution of the equations

[TABLE]

while, for $i=1$ the Neumann boundary condition is taken into account as follows,

[TABLE]

Similarly, we still denote by $\underline{R}_{1}$ the discrete version of the boundary operator $\underline{R}_{1}$, naturally defined by the second order approximation

[TABLE]

3.4.2. A finite volume scheme with nonlocal flux

We first rewrite (28) in the condensed form

[TABLE]

with $U=(\zeta,q)^{T}$ and

[TABLE]

and where

[TABLE]

(${\mathfrak{f}}(U)$ is the flux in the momentum equation for the nonlinear shallow water equations (1) in dimensionless variables). The flux in (33) is therefore a nonlocal operator with respect to $U$. The source term $S$ in (33) is given by

[TABLE]

where we recall that $\underline{{\mathcal{Q}}}(\underline{q},f,\ddot{f},\zeta,q)$ is defined in (29).

Using a finite volume type discretization for the (33) and a standard Euler scheme for the ODE on $\underline{q}$, we obtain the following general discretization of the Boussinesq system (28),

[TABLE]

where $U^{n}=(\zeta^{n},q^{n})^{T}=(\zeta_{i}^{n},q_{i}^{n})_{i\geq 1}^{T}$ and the source term $S^{n}_{i}$ is being given by

[TABLE]

(note that the definition for the discretized version of $\underline{\mathcal{Q}}$ can straightforwardly be deduced from (29) along the lines of §3.4.1; see also Remark 3 below). The source term involves the quantity $q_{0}^{n}$ which cannot be computed by induction through (36) since in (36), one assumes that $i\geq 1$. However, a direct discretization of (32) yields

[TABLE]

It remains of course to explain how to compute the discrete fluxes ${\mathfrak{F}}_{\mu,i+1/2}$. As above for the nonlinear shallow water equations, we consider here the simplest case of the Lax-Friedrichs scheme where the discrete flux is given by

[TABLE]

where we write, when $i\geq 1$,

[TABLE]

the discrete operator $R_{1}$ being constructed as in §3.4.1.

When $i=0$, this definition is naturally adapted as follows,

[TABLE]

the discrete boundary operator $\underline{R}_{1}$ being constructed as in §3.4.1 while ${q}_{0}^{n}$ is provided by (38).

Remark 3*.*

The quantity $\underline{{\mathcal{Q}}}(q_{0}^{n},f^{n},\ddot{f}^{n},\zeta^{n},q^{n})$ that appears in the right-hand side of the momentum equation in (36) and in the discrete ODE (38) for $q_{0}^{n}$ can be written using the notation (41) as

[TABLE]

All these quantities are already known so that handling generating boundary condition can be done with no extra computational cost compared to, say, periodic boundary conditions.

4. Numerical validations

4.1. The validation method

Since the implementation of reflecting or periodic boundary conditions does not raise any problem for the Boussinesq equations (24) we compute first a solution $U^{L}=(\zeta^{L},q^{L})^{T}$ of the equations under consideration in a larger domain $[-L,L]$ until a final time $T_{f}$, with reflective or periodic boundary conditions at both extremities, and with a non trivial initial condition. We then define a reference solution as the restriction of $U^{L}$ on $[0,L]$, and a boundary data $f$ as

[TABLE]

We then use the scheme presented in §3.4 to compute the solution $U$ of the Boussinesq system (24) with initial data $U^{0}(x)=U^{\rm ref}(t=0,x)$ and boundary data $f$, and compare it with the reference solution $U^{\rm ref}$. We define in particular the errors ${\mathbf{e}^{\zeta}_{\delta_{x}}}(t)$ and ${\mathbf{e}^{q}_{\delta_{x}}}(t)$ as

[TABLE]

and we compute the overall errors $e^{\zeta}_{\delta_{x}}$ and $e^{q}_{\delta_{x}}$ on $[0,T_{f}]$ as

[TABLE]

The convergence order $p$ is computed with a least-squares linear regression, whose coefficient is plotted on the error curves.

4.2. Propagation of gaussian initial conditions

We recall that the Boussinesq equations (24) are derived under the smallness assumption (25) on $\varepsilon$ and $\mu$. We consider here the approximation error in different cases,

[TABLE]

the nonlinear and dispersive effect become more important when $\varepsilon$ and $\mu$ respectively become larger; in particular, the configuration (${\rm I}$) is quite stiff and in the limit of the range of validity of the Boussinesq equations (for strong nonlinearities, one should rather work with the more complicated Serre-Green-Naghdi equations [18, 19]).

The initial datum for $U^{L}$ in the larger domain is

[TABLE]

The reference solution is computed with the Lax-Friedrichs scheme with non local flux introduced above, on the domain $[-L,L]$, with a very refined mesh: $n_{x}=3600$, and a time step $\delta_{t}=0.9\delta_{x}$ in agreement with the CFL condition computed from the approximated velocities of the Riemann invariants. We take $L=5$. We compute the numerical solution with the nonlocal Lax-Friedrichs scheme in the domain $[0,L]$, on coarser meshes: $n_{x}=90,120,150,180,200,300,360,450$. The meshes are defined so that the points of the coarse meshes always coincide with the points of the finer mesh. The boundary conditions at $x=0$ are taken into account by imposing the reference solution and its second-order time derivative approximated with the classical centered second-order scheme.

As the initial data is zero near the boundaries of the large domain, no special effort is necessary for the computation with the coarse mesh at the right boundary $x=L$ if the final time of the simulation is not too large. We shall compare the solution over a time interval $t\in[0,2]$. The qualitative behavior of the solution is the following: each of the two gaussians decomposes into two waves roughly traveling at speed $1$ and $-1$ respectively. The gaussian located on the left being closer to the boundary $x=0$ of the small domain, this configuration is rich enough to contain the three main relevant cases,

1. (1)

The forcing $f$ corresponds to an essentially incoming wave. This is the situation that occurs for $t\sim 0.1$ (see Figure 1 left) 2. (2)

The forcing $f$ corresponds to the superposition of an outgoing and an incoming wave. This is the situation that occurs for $t\sim 1$ (see Figure 1 middle) 3. (3)

The forcing $f$ corresponds to an essentially outgoing wave. This is the situation that occurs for $t\sim 1.5$ (see Figure 1 right)

Numerical results for the initial condition (43)-(44) are presented on Figures 2, 3 and 4, in logarithmic scale with the slope obtained from a linear regression. On Figure 2, corresponding to the case $\mu=\epsilon=0.3$, the slope of the linear regression obtained with all error points is completed with the slope of the linear regression obtained with the four more refined error points. Globally, a first-order convergence in space is observed for both variables.

4.3. Soliton propagation

We test here our scheme on the propagation of a solitary wave, which involves both nonlinear and dispersive effects. The soliton for the non-linear Boussinesq system (24) is not explicit, but we compute it by solving numerically a second-order differential equation that we obtain as follows. We look for a solution of the non-linear Boussinesq equations such that $\zeta(x,t)=\tilde{\zeta}(x-ct)$ and $q(x,t)=\tilde{q}(x-ct)$. We inject this form in the first equation of (24) and find, omitting the tilde symbol for the sake of brevity

[TABLE]

Then we inject this relationship in the second equation of (24), integrate in space, and we can write (imposing that $\zeta$ vanishes at infinity),

[TABLE]

Multiplying this equation by $\zeta^{\prime}$ and using again that $\zeta$ tends to zero when $x$ tends to $\pm\infty$, we obtain

[TABLE]

Denoting $\zeta_{{\rm max}}$ the maximum value of $\zeta$ we can compute $c$ as a function of $\zeta_{{\rm max}}$ and $\varepsilon$.

[TABLE]

Once $c$ is computed, we solve the differential equation (45) with a high order numerical method in order to obtain our reference solution. We choose $\zeta_{{\rm max}}=1$ and $\mu=\varepsilon=0.3$ or $0.1$. We have checked that if we solve the Boussinesq system with this reference solution as an initial datum, with the nonlocal Lax-Friedrichs scheme and periodic boundary conditions, the numerical results show a first order convergence: see Figures 5 and 6. The space steps $\delta_{x}$ were computed as $\delta_{x}=L/n_{x}$, with $n_{x}=400,600,800,1000,1200,1400,1600,1800,2000$.

To test the imposition of the generating boundary condition we compute the numerical solution of the soliton on the small domain $[0,L]$ with L = 10. We use the nonlocal Lax-Friedrichs scheme and a constant time step $\delta_{t}=0.8\,\delta_{x}$ for $\mu=\varepsilon=0.3$, and $\delta_{t}=0.9\,\delta_{x}$ for $\mu=\varepsilon=0.1$, taking into account the values of the approximated eigenvalues. The space step is computed as $\delta_{x}=L/n_{x}$, with $n_{x}=400,600,800,1000,1200,1400,1600,1800,2000$. The maximum of the soliton is initially located on the left of the computational domain, at $x=-L/2$, so that the initial datum in the small domain is almost zero, and then the soliton propagates inside it. The boundary conditions on the left boundary of the small domain are taken into account by imposing the reference solution and its second-order time derivative approximated with the classical centered second-order scheme. As the initial datum is zero near the right boundary $x=L$, no special effort is necessary for the computation at this boundary if the final time of the simulation is not too large. The values of $\zeta$ in the small domain at the final time for the reference solution and the numerical solution are plotted on Figure 7. The numerical results are presented on Figures 8 and 9. On Figure 9, the slope of the linear regression obtained with all error points is completed with the slope of the linear regression obtained with the four more refined error points. Globally a first-order convergence is observed when the grid is sufficiently refined.

4.4. Sinusoidal boundary condition

We consider the cases

[TABLE]

Note that case (II) with different values of $\varepsilon$ and $\mu$ has been added here for its relevance for applications in coastal oceanography where a sinusoidal swell is imposed at the entrance of the domain in a region not so shallow (so that $\mu$ is not very small) but the waves are of small amplitude (they become bigger in the shoaling phase, nearer to the shore), so that $\varepsilon$ is small.

We first compute a numerical solution $U^{L}$ with a very refined mesh ($n_{x}$ =3600) on a larger domain $[-L,L]$, with $L=10$, with the Lax-Friedrichs scheme and a time step $\delta_{t}=0.9\delta_{x}$ The initial condition for $U^{L}$ in the larger domain is

[TABLE]

and we impose until the final time $T_{f}=15$ the generating boundary condition

[TABLE]

We define the reference solution on the slightly smaller domain $[-0.8\,L,L]$

[TABLE]

Then we compute a solution with coarse meshes on this smaller domain $[-0.8\,L,L]$. The boundary conditions at $x=-0.8\,L$ are taken into account by imposing the reference solution and its second-order time derivative approximated with the classical centered second-order scheme. No special effort is necessary for the computation with the coarse mesh at the right boundary $x=L$ until the final time $T_{f}=15$. The values of $\zeta$ in the small domain at the final time for the reference solution and the numerical solution are plotted on Figure 10. Because there is numerical dissipation, the numerical solution has a smaller amplitude than the reference solution after some time of propagation inside the small domain, but both solutions coincide well near the left boundary. The error between the reference solution and the solution on the coarse mesh is computed near the left boundary on the interval $[-0.8\,L,-0.6\,L]$ in order to measure the error due to the generating boundary condition rather than the dissipation error inherent to the Lax-Friedrichs scheme (which stronger in this numerical test than in the previous ones due to the fact the the reference solution involves higher frequencies).

The numerical results for the cases $({\rm I})$ and $({\rm II})$ are presented on Figures 11 and 12. The space steps $\delta_{x}$ were chosen as $\delta_{x}=2L/n_{x}$, with $n_{x}$ = 100, 120, 150, 200, 240, 300, 360, 400, 600, 720, 900. On both Figures the slope of the linear regression obtained with all error points is completed with the slope of the linear regression obtained with the three more refined error points. The numerical results for the cases $({\rm III})$ and $({\rm IV})$ are presented on 13 and 14. The space steps $\delta_{x}$ were chosen as $\delta_{x}=2L/n_{x}$, with $n_{x}=$ 90, 100, 120, 150, 200, 240, 300, 360, 400. Globally, a first-order convergence is observed when the grid is sufficiently refined.