# Finite volume methods for unidirectional dispersive wave models

Denys Dutykh, Theodoros Katsaounis, Dimitrios Mitsotakis

► To cite this version:

Denys Dutykh, Theodoros Katsaounis, Dimitrios Mitsotakis. Finite volume methods for unidirectional dispersive wave models. International Journal for Numerical Methods in Fluids, 2012, 71 (6), pp.717-736. 10.1002/fld.3681 . hal-00538043v3

HAL Id: hal-00538043

https://hal.science/hal-00538043v3

Submitted on 3 Feb 2012

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# FINITE VOLUME METHODS FOR UNIDIRECTIONAL DISPERSIVE WAVE MODELS

DENYS DUTYKH, THEODOROS D. KATSAOUNIS*, AND DIMITRIOS MITSOTAKIS

ABSTRACT. We extend the framework of the finite volume method to dispersive unidirectional water wave propagation in one space dimension. In particular we consider a KdV-BBM type equation. Explicit and IMEX Runge-Kutta type methods are used for time discretizations. The fully discrete schemes are validated by direct comparisons to analytic solutions. Invariants conservation properties are also studied. Main applications include important nonlinear phenomena such as dispersive shock wave formation, solitary waves and their various interactions.

## 1. INTRODUCTION

Water wave modeling is a complicated process and usually leads to models which are hard to analyze mathematically as well as to solve numerically. Under certain simplifying assumptions approximate models are obtained, e.g. the KdV equation [28], the BBM equation [4] and Boussinesq systems [11, 40, 8]. All these models assume the wave to be weakly nonlinear and weakly dispersive, propagating mainly in one space direction. These approximate models consider mainly unidirectional or bidirectional wave propagation on flat or complex bathymetries.

In this paper we study the application of some finite volume schemes to a scalar non-linear dispersive partial differential equation modeling unidirectional wave propagation. Specifically, we consider the KdV-BBM equation in its general form:

$$u_t + \alpha u_x + \beta uu_x - \gamma u_{xxt} + \delta u_{xxx} = 0, \quad (1.1)$$

for $x \in \mathbb{R}$, $t > 0$, where $\alpha, \beta, \gamma, \delta$ are positive real numbers, [4]. The finite volume method is well known for its accuracy, efficiency, robustness and excellent local conservative properties. Most often this method is employed to approximate solutions to hyperbolic conservation laws. The system of Nonlinear Shallow Water Equations (NSWE) is a classical example of the successful application of modern finite volume schemes to water wave problems.

A wide range of numerical methods have been employed to compute approximate solutions to dispersive wave equations of KdV-BBM type: finite difference schemes [10, 50], finite element methods [9, 33, 2] and spectral methods [38, 39, 15, 35]. Recently discontinuous Galerkin schemes have also been employed to dispersive wave equations [51, 30, 18], (the list is far from being exhaustive). However, the application of finite volume or hybrid

<footnote>Key words and phrases. finite volume method; nonlinear dispersive waves; unidirectional propagation; solitary waves; water waves.</footnote>

<footnote>* Corresponding author.</footnote>
FV/FD methods remain most infrequent for this type of problems. To our knowledge, only a few recent works are in this direction [3, 17, 6, 47, 42, 12].

In order to apply the finite volume method to the KdV-BBM equation (1.1), we rewrite it in a conservative form, including a nontrivial evolution operator, an advective and a dispersive flux functions. In the finite volume literature there exist several ways to approximate these fluxes. For the advective part we test three different numerical fluxes each one representing a particular family of finite volume method:

• *average flux* (m-scheme),

• *central flux*, (KT-scheme) as a representative of central schemes, [34, 29],

• *characteristic flux* (CF-scheme), as a representative of upwind schemes and linearized Riemann solvers, [20, 22].

The dispersive term is discretized using simply the *average flux*, while high order approximations are used for the BBM term ($\gamma u_{xxt}$). The *central flux* and the *characteristic flux* are widely used in the case of conservation laws. On the other hand the *average flux*, known to be unstable for conservation laws, performs equally well.

The evaluation of the numerical flux functions require approximate values of the solution at the cell interfaces. The order of the approximation determines the space accuracy of the underlying finite volume scheme. We consider first order, taking simply piecewise constant approximations, as well as high order schemes. The high order accuracy is achieved through application of various reconstruction techniques such as TVD [46], UNO [26] and WENO [31].

The time discretization of (1.1) is based on Runge-Kutta methods. The stability of the resulting system of ode's depends on the interplay between the BBM term ($\gamma u_{xxt}$) and the KdV type dispersive term ($\delta u_{xxx}$). An explicit discretization of the ode system is sufficient when these terms are of the same order. Thus, Strong Stability Preserving Runge-Kutta (SSP-RK) methods, which preserve the TVD property of the finite volume scheme, [44, 24] are used for the explicit discretization.

However, when $\gamma \ll \delta$ the resulting semidiscrete system of ode's is highly stiff and therefore implicit methods with strong stability characteristics are preferable. To balance the high computational cost of fully implicit methods and stability considerations we rely on Implicit-Explicit Runge-Kutta (IMEX) methods, [1]. Indeed IMEX RK methods turned out to be well suited for the time discretization of the KdV-BBM equation (1.1) exhibiting excellent stability behavior.

The validated numerical method is applied to study the KdV-BBM equation (1.1) in a systematic way through a series of numerical experiments. In particular, we focus on the following issues:

• accuracy of the finite volume method for solitary wave propagation and invariants conservation

• dispersive shock formation (we underline that the finite element as well as spectral methods break down for this experiment while the finite volume method provides robust and accurate results)

• interactions of solitary waves (overtaking collisions)
The paper is organized as follows. In Section 2 the governing equation (1.1) is presented briefly along with its basic properties. In Section 3 the finite volume discretization as well as fully discrete schemes are presented in details. In Section 4 we validate the discretization procedure by comparisons with analytical solution. Several important test cases are also presented.

## 2. DISPERSIVE WATER WAVE MODEL EQUATION

We present briefly the mathematical model under consideration and some of its basic properties. The KdV-BBM equation takes the following general form:

$$u_t + \alpha u_x + \beta uu_x - \gamma u_{xxt} + \delta u_{xxx} = 0, \quad (2.1)$$

where $x \in \mathbb{R}$, $t > 0$, $u$ denotes the free surface elevation above the still water level $u = 0$ and $\alpha, \beta, \gamma, \delta$ are positive real numbers. Equation (2.1) incorporates nonlinear and dispersive effects and has been suggested as a model for surface water waves in a uniform channel with flat bottom, cf. ([4, 19]).

When $\delta = 0$, (2.1) reduces to the BBM equation [4], while taking $\gamma = 0$ leads the celebrated KdV equation [28]. The KdV-BBM model (2.1) has been studied thoroughly in the past and the Cauchy problem is known to be well-posed in appropriate Sobolev spaces, at least locally in time. Also the well-posedness of some initial-boundary value problems, including the initial-periodic boundary value problem, can be proved, cf. e.g. [4, 7, 19] and the references therein.

One may easily check that equation (2.1) admits exact solitary wave solutions of the form:

$$u(x,t) = 3 \frac{c_s - \alpha}{\beta} \operatorname{sech}^2 \left( \frac{1}{2} \sqrt{\frac{c_s - \alpha}{\gamma c_s + \delta}} (x - c_s t) \right), \quad (2.2)$$

that travel rightwards with a given speed $c_s$. We are going to exploit this solution below in order to validate our discretization procedure and measure the order of convergence of proposed numerical schemes. Further it is well known that (2.1) possesses two quantities invariant under its evolution dynamics. Assuming either the solution has compact support or $u \to 0$ $x \to \pm\infty$, one can easily check that quantities

$$I_1(t) = \int_{\mathbb{R}} u(x,t) dx, \quad I_2(t) = \int_{\mathbb{R}} (u^2(x,t) + \gamma u_x^2(x,t)) dx, \quad (2.3)$$

are conserved in time, i.e. $I_1(t) = I_1(0)$, $I_2(t) = I_2(0)$, $\forall t > 0$. The invariant $I_1$ reflects the physical property of the mass conservation, while invariant $I_2$ can be assimilated to the generalized kinetic energy. Invariants conservation is a fundamental property important not only for theoretical investigations but also for numerics since it allows to validate numerical schemes and to quantify the accuracy of the obtained results.

For more realistic situation one has to consider bidirectional models with uniform or variable bathymetry cf. e.g. [8, 40]. For a systematic numerical study of such Boussinesq type systems using finite volume methods analogous to those presented in this paper, including the runup algorithm we refer to [16].
### 3. FINITE VOLUME DISCRETIZATION

We proceed to the discretization of (2.1) by a finite volume method. Our motivation stems from the observation that the KdV-BBM equation can be seen as a dispersive perturbation<fnref n="1" /> of the following inviscid Burgers equation:

$$u_t + (\alpha u + \frac{\beta}{2} u^2)_x = 0.$$

Consequently, the proposed finite volume schemes are based on the corresponding schemes for scalar conservation laws. A special treatment is introduced for the discretization of dispersive terms.

Let $\mathcal{T} = \{x_i\}, i \in \mathbb{Z}$ be a partition of $\mathbb{R}$ into cells $C_i = (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}})$, where $x_i = (x_{i+\frac{1}{2}} + x_{i-\frac{1}{2}})/2$ denotes the midpoint of the cell $C_i$. Let $\Delta x_i = x_{i+\frac{1}{2}} - x_{i-\frac{1}{2}}$ denote the length of the cell $C_i$ and let $\Delta x_{i+\frac{1}{2}} = x_{i+1} - x_i$. Herein, we assume the partition $\mathcal{T}$ to be uniform, i.e. $\Delta x_i = \Delta x_{i+\frac{1}{2}} = \Delta x$, $i \in \mathbb{Z}$. For a scalar function $w(x,t)$ let $w_i$ denotes its cell average on $C_i$:

$$w_i(t) = \frac{1}{\Delta x} \int_{C_i} w(x, t) dx.$$

We rewrite (2.1) in a conservative-like form:

$$(I - \gamma \partial_x^2) u_t + [F(u)]_x + [G(u_{xx})]_x = 0, \quad (3.1)$$

where the advective flux is $F(u) = \alpha u + \frac{\beta}{2} u^2$ and the dispersive flux is $G(v) = \delta v$. We underline that $F$ is a convex flux function. A simple integration of (3.1) over a cell $C_i$ yields:

$$\begin{aligned} & \frac{d}{dt} \left[ u_i(t) - \frac{\gamma}{\Delta x} \left( u_x(x_{i+\frac{1}{2}}, t) - u_x(x_{i-\frac{1}{2}}, t) \right) \right] \\ & + \frac{1}{\Delta x} \left[ F(u(x_{i+\frac{1}{2}}, t)) - F(u(x_{i-\frac{1}{2}}, t)) \right] + \frac{1}{\Delta x} \left[ G(u_{xx}(x_{i+\frac{1}{2}}, t)) - G(u_{xx}(x_{i-\frac{1}{2}}, t)) \right] = 0, \end{aligned} \quad (3.2)$$

where the values of the advective and dispersive fluxes on the cell interfaces have to be properly defined.

### 3.1. Semidiscrete scheme.
We proceed to the construction of the semidiscrete finite volume approximation. Let $\chi_{C_i}$ be the characteristic function of the cell $C_i$. We define a piecewise constant function $u_h(x, t) = \sum_{i \in \mathbb{Z}} U_i(t) \chi_{C_i}(x)$, where $U_i(t)$ are solutions of the following system of ordinary differential equations:

$$\frac{d}{dt} \left[ U_i - \frac{\gamma}{\Delta x} \left( \frac{U_{i+1} - 2U_i + U_{i-1}}{\Delta x} \right) \right] + \frac{1}{\Delta x} \left( \mathcal{F}_{i+\frac{1}{2}} - \mathcal{F}_{i-\frac{1}{2}} \right) + \frac{1}{\Delta x} \left( \mathcal{G}_{i+\frac{1}{2}} - \mathcal{G}_{i-\frac{1}{2}} \right) = 0, \quad (3.3)$$

with initial conditions defined as a projection onto the space of piecewise constant functions on $\mathcal{T}$:

$$U_i(0) = \frac{1}{\Delta x} \int_{C_i} u(x, 0) dx, \quad i \in \mathbb{Z}.$$

<footnote n="1">¹Since the wave is assumed to be weakly nonlinear and weakly dispersive.</footnote>
In (3.3) $\mathcal{F}$ and $\mathcal{G}$ denote the advective and the (KdV-type) dispersive numerical fluxes respectively. More specifically, $\mathcal{F}_{i+\frac{1}{2}} = \mathcal{F}(U_{i+\frac{1}{2}}^L, U_{i+\frac{1}{2}}^R)$ and $\mathcal{G}_{i+\frac{1}{2}} = \mathcal{G}(W_{i+\frac{1}{2}}^L, W_{i+\frac{1}{2}}^R)$ are approximations of $F(u(x_{i+\frac{1}{2}}, t))$ and $G(u_{xx}(x_{i+\frac{1}{2}}, t))$ respectively at cell interfaces. Values $U_{i+\frac{1}{2}}^L, U_{i+\frac{1}{2}}^R$ are approximations to the point value $u(x_{i+\frac{1}{2}}, t)$ from cells $C_i, C_{i+1}$ respectively, while $W_{i+\frac{1}{2}}^L$ and $W_{i+\frac{1}{2}}^R$ are corresponding approximations to the point value of the second derivative $u_{xx}(x_{i+\frac{1}{2}}, t)$. All quantities $U_{i+\frac{1}{2}}^L, U_{i+\frac{1}{2}}^R$ as well as $W_{i+\frac{1}{2}}^L, W_{i+\frac{1}{2}}^R$ are computed by a reconstruction process described below (see Section 3.1.2).

### 3.1.1. Advective and dispersive numerical fluxes.
Over the last twenty years numerous numerical fluxes $\mathcal{F}$ have been proposed to discretize advective operators [41, 25, 37, 23, 5]. We select three quite different flux functions. Namely, we consider a simple *average flux* $\mathcal{F}^m$, a *central type flux* $\mathcal{F}^{KT}$, [29, 34] and a *characteristic flux* $\mathcal{F}^{CF}$,[20, 21, 22] :

$$ \mathcal{F}^m(U, V) = F\left(\frac{U+V}{2}\right), \qquad (3.4) $$

$$ \mathcal{F}^{KT}(U, V) = \frac{1}{2}\{[F(U) + F(V)] - \mathcal{A}(U, V)[V - U]\}, \qquad (3.5) $$

$$ \mathcal{F}^{CF}(U, V) = \frac{1}{2}\{[F(U) + F(V)] - \mathcal{A}(U, V)[F(V) - F(U)]\}. \qquad (3.6) $$

The *average* flux is perhaps the simplest one and is known to be unconditionally unstable for nonlinear conservation laws. However, this flux shows very good performance for dispersive waves (see Section 4).

The *central* flux is of Lax-Friedrichs type and is a representative of the family of central schemes. The operator $\mathcal{A}$ in the KT-scheme is related to characteristic speeds of the flow and is given by this expression:

$$ \mathcal{A}(U, V) = \max [|F'(U)|, |F'(V)|]. \qquad (3.7) $$

The *characteristic* flux function is somehow similar to the Roe scheme [41] and the operator $\mathcal{A}$ in this case is defined as:

$$ \mathcal{A}(U, V) = \text{sign}\left(F'\left(\frac{U+V}{2}\right)\right) = \text{sign}\left(\alpha + \beta \frac{U+V}{2}\right). \qquad (3.8) $$

For the dispersive numerical flux $\mathcal{G}$ we choose to work with the *average flux function* (3.4):

$$ \mathcal{G}(W, R) = \delta \frac{W + R}{2}, \qquad (3.9) $$

where $W$ and $R$ are standard central approximations of the second derivative from each side. The numerical flux $\mathcal{G}$ can be evaluated either using simple cell averages, denoted by $\mathcal{G}^m$, or higher order approximation based on a reconstruction procedure, denoted by $\mathcal{G}^{lm}$.
3.1.2. Reconstruction process. The values $U_{i+\frac{1}{2}}^L, U_{i+\frac{1}{2}}^R$ are approximations to $u(x_{i+\frac{1}{2}}, t)$ from cells $C_i$ and $C_{i+1}$ respectively. The simplest choice is to take the piecewise constant approximation in each cell:

$$U_{i+\frac{1}{2}}^L = U_i, \quad U_{i+\frac{1}{2}}^R = U_{i+1}. \qquad (3.10)$$

The resulting semidiscrete finite volume scheme is formally first order accurate in space. To achieve a higher order accuracy in space, we have to adopt more elaborated reconstruction process. The main idea is to use the cell averages $U_i$ to reconstruct more accurate approximation to the solution at cell interfaces $u(x_{i+\frac{1}{2}}, t)$. For this purpose we consider three different reconstruction methods: the classical MUSCL type (TVD2) piecewise linear reconstruction [27, 48], the UNO2 reconstruction [26] and WENO type reconstructions, [31].

* The classical TVD2 scheme uses a linear reconstruction :

$$U_{i+\frac{1}{2}}^L = U_i + \frac{1}{2}\phi(r_i)(U_{i+1} - U_i), \quad U_{i+\frac{1}{2}}^R = U_{i+1} - \frac{1}{2}\phi(r_{i+1})(U_{i+2} - U_{i+1}), \qquad (3.11)$$

where $r_i = \frac{U_i - U_{i-1}}{U_{i+1} - U_i}$, and $\phi$ is an appropriate slope limiter function, [46]. There exist many possible choices of the slope limiter. Some of the usual choices are

- MinMod (MM) limiter : $\phi(\theta) = \max(0, \min(1, \theta))$,

- VanLeer (VL) limiter : $\phi(\theta) = \frac{\theta+|\theta|}{1+|\theta|}$,

- Monotonized Central (MC) limiter : $\phi(\theta) = \max(0, \min((1+\theta)/2, 2, 2\theta))$,

- Van Albada (VA) limiter : $\phi(\theta) = \frac{\theta+\theta^2}{1+\theta^2}$.

The last three limiters have been shown to produce sharper resolution of discontinuities, and in our case less dissipative numerical results. The TVD2 reconstruction is formally second order accurate except at local extrema where it reduces to the first order. Reconstructions considered below were proposed to remove this shortcoming.

* The UNO2, like the TVD2, is also a linear reconstruction process which is second order accurate even at local extrema. The values $U_{i+\frac{1}{2}}^L, U_{i+\frac{1}{2}}^R$ are defined as

$$U_{i+\frac{1}{2}}^L = U_i + \frac{1}{2}S_i, \quad U_{i+\frac{1}{2}}^R = U_{i+1} - \frac{1}{2}S_{i+1}, \qquad (3.12)$$

where

$$S_i = m(S_i^+, S_i^-), \quad S_i^\pm = d_{i\pm\frac{1}{2}}U \mp \frac{1}{2}D_{i\pm\frac{1}{2}}U,$$

$$d_{i+\frac{1}{2}}U = U_{i+1} - U_i, \quad D_{i+\frac{1}{2}}U = m(D_iU, D_{i+1}U),$$

$$D_iU = U_{i+1} - 2U_i + U_{i-1}, \quad m(x, y) = \frac{1}{2}(\mathrm{sign}(x) + \mathrm{sign}(y)) \min(|x|, |y|)$$

The UNO2 reconstruction is formally second accurate even at local extrema.

* We also consider WENO type reconstructions [31, 43]. Namely, we implement the 3rd and 5th order accurate WENO methods, hereafter referred to as WENO3 and WENO5 respectively. For the sake of clarity, we present here only WENO3 scheme.
First of all we compute the 3rd order reconstructed values:

$$U_{i+\frac{1}{2}}^{(0)} = \frac{1}{2}(U_i + U_{i+1}), \quad U_{i+\frac{1}{2}}^{(1)} = \frac{1}{2}(-U_{i-1} + 3U_i),$$

$$U_{i-\frac{1}{2}}^{(0)} = \frac{1}{2}(3U_i - U_{i+1}), \quad U_{i-\frac{1}{2}}^{(1)} = \frac{1}{2}(U_{i-1} + U_i).$$

Then, we define the smoothness indicators:

$$\beta_0 = (U_{i+1} - U_i)^2, \quad \beta_1 = (U_i - U_{i-1})^2,$$

and constants $d_0 = \frac{2}{3}$, $d_1 = \frac{1}{3}$, $\tilde{d}_0 = d_1$, $\tilde{d}_1 = d_0$. The weights are defined as:

$$\omega_0 = \frac{\alpha_0}{\alpha_0 + \alpha_1}, \quad \omega_1 = \frac{\alpha_0}{\alpha_0 + \alpha_1}, \quad \tilde{\omega}_0 = \frac{\tilde{\alpha}_0}{\tilde{\alpha}_0 + \tilde{\alpha}_1}, \quad \tilde{\omega}_1 = \frac{\tilde{\alpha}_1}{\tilde{\alpha}_0 + \tilde{\alpha}_1},$$

where $\alpha_i = \frac{d_i}{\epsilon + \beta_i}$, $\tilde{\alpha}_i = \frac{\tilde{d}_i}{\epsilon + \beta_i}$ and $\epsilon$ is a small, positive number (in our computations we set $\epsilon = 10^{-15}$).

Finally, the reconstructed values are given by formulas:

$$U_{i+\frac{1}{2}}^{L} = \sum_{r=0}^{1} \omega_r U_{i+\frac{1}{2}}^{(r)}, \quad U_{i-\frac{1}{2}}^{R} = \sum_{r=0}^{1} \tilde{\omega}_r U_{i-\frac{1}{2}}^{(r)}. \qquad (3.13)$$

**Remark 1.** The elliptic operator approximation in (3.3) is only second order accurate. In the case where a high order WENO reconstruction is used, we need to increase also the elliptic solver accuracy. For example, the following semidiscrete scheme:

$$\frac{d}{dt} \left[ \frac{U_{i-1} + 10U_i + U_{i+1}}{12} - \gamma \frac{U_{i+1} - 2U_i + U_{i-1}}{\Delta x^2} \right] + \frac{\mathcal{H}_{i-1} + 10\mathcal{H}_i + \mathcal{H}_{i+1}}{12} = 0 \quad (3.14)$$

where $\mathcal{H}_i = \frac{1}{\Delta x}(\mathcal{F}_{i+\frac{1}{2}} - \mathcal{F}_{i-\frac{1}{2}}) + \frac{1}{\Delta x}(\mathcal{G}_{i+\frac{1}{2}} - \mathcal{G}_{i-\frac{1}{2}})$ is a fourth order approximation. Thus in the WENO3 case a global third order accuracy is observed, while for WENO5 interpolation, we profit only locally by the 5th order accuracy of the reconstruction, cf. Section 4.1.

**Remark 2.** In computation of the dispersive flux we distinguish between the simple averaging of cell centered values in $\mathcal{G}^m$ and of $\mathcal{G}^{lm}$, where higher order reconstructions of the second order derivatives are used.

**3.2. Fully discrete schemes.** We consider now fully discrete schemes for the ode system (3.3). The time discretization is based on Runge-Kutta type methods. Explicit schemes based on TVD preserving RK-methods are presented. In certain cases where stiffness becomes dominant, we use an implicit-explicit strategy based on IMEX type RK-methods.

**3.2.1. Explicit schemes.** The initial value problem (3.3) can be discretized by various methods. When the parameter $\gamma$ is of the same order as $\delta$ the system of ode's appeared to be non-stiff and therefore can be integrated numerically by any explicit time-stepping method. We use a special class of Runge-Kutta methods that preserve the TVD property of the finite volume scheme, [44, 24, 45].
Let $\Delta t$ be the temporal stepsize and let $t^{n+1} = t^n + \Delta t$, $n \ge 0$ be discrete time levels, then (3.3) is an initial value problem of the form

$$ \mathbf{T}\mathbf{U}' = L(\mathbf{U}), \qquad (3.15) $$

where $\mathbf{U} = \{U_i\}$, $i \in \mathbb{Z}$, $\mathbf{T} = \mathbf{I} + [-\gamma, 2\gamma, -\gamma]/\Delta x^2$ is a tridiagonal matrix and $L$ is a nonlinear operator incorporating the contribution of the numerical fluxes $\mathcal{F}$, $\mathcal{G}$. Assuming at time $t^n$, $\mathbf{U}^n$ is known then $\mathbf{U}^{n+1}$ is defined by

$$ \begin{aligned} \mathbf{U}^{n+1} &= \mathbf{U}^n - \frac{\Delta t}{\Delta x} \sum_{j=1}^{s} b_j \mathbf{T}^{-1} L(\mathbf{U}^{n,j}), \\ \mathbf{U}^{n,j} &= \mathbf{U}^n - \frac{\Delta t}{\Delta x} \sum_{\ell=1}^{s-1} a_{j\ell} \mathbf{T}^{-1} L(\mathbf{U}^{n,\ell}), \end{aligned} \qquad (3.16) $$

where the set of constants $A = (a_{j\ell})$, $b = (b_1, \dots, b_s)$ define a $s$-stage Runge-Kutta method. The following tableau are examples of explicit TVD RK-methods which are of 2nd and 3rd order respectively

$$ \begin{array}{c|c} \begin{array}{cc|c} 0 & 0 & 0 \\ 1 & 0 & 1 \\ \hline \frac{1}{2} & \frac{1}{2} & \end{array} & \begin{array}{cc|c} 0 & 0 & 0 \\ 1 & 0 & 0 \\ \hline \frac{1}{4} & \frac{1}{4} & 0 \\ \hline \frac{1}{6} & \frac{1}{6} & \frac{2}{3} \end{array} & \begin{array}{c} 0 \\ 1 \\ \frac{1}{2} \\ \hline \end{array} \end{array} \qquad (3.17) $$

In our computations we mainly use the 3-stage third order method.

### 3.2.2. Implicit-Explicit schemes.
As the parameter $\gamma$ decreases to zero the semidiscretization of the KdV-BBM equation leads to a stiff system of ode's. To solve efficiently this system we apply an IMEX type RK-method, [1]. The linear dispersive terms are treated in an implicit way while the rest of the terms are treated explicitly. Numerical evidence shows that IMEX methods exhibit excellent stability and handle stiffness in an efficient and robust way even in the limiting case $\gamma = 0$.

We consider an $s$-stage Diagonally Implicit Runge-Kutta (DIRK) method, properly chosen, that is given by the tableau

$$ \frac{A}{b} \rvert_{\tau} = \begin{array}{c} \begin{array}{ccccc} a_{11} & 0 & \cdots & 0 & \tau_1 \\ a_{21} & a_{22} & \cdots & 0 & \tau_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{s1} & a_{s2} & \cdots & a_{ss} & \tau_s \end{array} \\ \hline \begin{array}{ccccc} b_1 & b_2 & \cdots & b_s & \end{array} \end{array} \qquad (3.18) $$
and an $s + 1$ explicit Runge-Kutta method

$$ \begin{gathered} \begin{array}{c|c} \hat{A} & \hat{\tau} \\ \hline \hat{b} & \end{array} = \begin{array}{cccccc} 0 & 0 & \cdots & 0 & 0 & 0 \\ \hat{a}_{11} & 0 & \cdots & 0 & 0 & \hat{\tau}_1 \\ \hat{a}_{21} & \hat{a}_{22} & \cdots & 0 & 0 & \hat{\tau}_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ \hat{a}_{s1} & \hat{a}_{s2} & \cdots & \hat{a}_{ss} & 0 & \hat{\tau}_s \\ \hline \hat{b}_1 & \hat{b}_2 & \cdots & \hat{b}_s & 0 & \end{array} \end{gathered} \tag{3.19} $$

We rewrite system (3.15) in the form

$$ \mathbf{T}\mathbf{U}' = \mathcal{F}(\mathbf{U}) + \mathbf{D}\mathbf{U}, \tag{3.20} $$

where $\mathbf{D}$ is the five-diagonal matrix $\delta[-1/2, 1, 0, -1, 1/2]/\Delta x^3$ coming from the discretization of the KdV term when we use the numerical flux function $\mathcal{G}^m$. Then the fully discrete scheme can be written in the form

$$ (\mathbf{T} + \Delta t a_{ii} \mathbf{D}) \mathbf{U}^{(i)} = \mathbf{T} \mathbf{U}^n - \Delta t \sum_{j=1}^{i} \hat{a}_{ij} \mathcal{F}(\mathbf{U}^{(j)}) - \Delta t \sum_{j=1}^{i-1} a_{ij} \mathbf{D} \mathbf{U}^{(j)}, \quad i = 1, \dots, s, \tag{3.21} $$

$$ \mathbf{T}\mathbf{U}^{n+1} = \mathbf{T}\mathbf{U}^n - \Delta t \sum_{j=1}^{s} \hat{b}_j \mathcal{F}(\mathbf{U}^{(j)}) - \Delta t \sum_{j=1}^{s} b_j \mathbf{D}\mathbf{U}^{(j)}. \tag{3.22} $$

We employ four IMEX RK-methods of different number of stages, orders of accuracy and stability properties. In particular we consider the following pairs, [1]

* A two stage third order DIRK method and a corresponding three stage, third order accurate ERK method with $\gamma = (3 + \sqrt{3})/6$. The resulting IMEX method is third order accurate.

$$ \begin{array}{c|c} \begin{array}{cc} \gamma & 0 \\ 1 - 2\gamma & \gamma \end{array} & \begin{array}{cc} \gamma & 0 \\ 1 - \gamma & 0 \end{array} \\ \hline \frac{1}{2} & \frac{1}{2} \end{array}, \quad \begin{array}{c|c} \begin{array}{ccc} 0 & 0 & 0 \\ \gamma & 0 & 0 \\ 1 - \gamma & 2(1 - \gamma) & 0 \end{array} & \begin{array}{c} 0 \\ \gamma \\ 1 - \gamma \end{array} \\ \hline \frac{1}{2} & \frac{1}{2} \end{array} \tag{3.23} $$

* A two stage second order DIRK method which is stiffly accurate, with $\gamma = (2 - \sqrt{2})/2$. The corresponding ERK is a three stage second order accurate method with $\delta = -2\sqrt{2}/3$. The resulting IMEX combination is second order accurate.

$$ \begin{array}{c|c} \begin{array}{cc} \gamma & 0 \\ 1 - \gamma & \gamma \end{array} & \begin{array}{cc} \gamma & 0 \\ 1 & \gamma \end{array} \\ \hline \frac{1}{2} & \frac{1}{2} \end{array}, \quad \begin{array}{c|c} \begin{array}{ccc} 0 & 0 & 0 \\ \gamma & 0 & 0 \\ \delta & 1 - \delta & 0 \end{array} & \begin{array}{c} 0 \\ \gamma \\ 1 \end{array} \\ \hline \frac{1}{2} & \frac{1}{2} \end{array}, \tag{3.24} $$

* A three stage third order DIRK stiffly accurate method with larger dissipative region than (3.24). The corresponding ERK is a three stage third order method.
<table><caption>TABLE 1. Stability of IMEX for the KdV equation ($α = β = δ = 1$, $γ = 0$)</caption><thead><tr><th>Method</th><th>Δt/Δx ≤</th></tr></thead><tbody><tr><td>(3.23)</td><td>1/4</td></tr><tr><td>(3.24)</td><td>1/5</td></tr><tr><td>(3.25)</td><td>1</td></tr><tr><td>(3.26)</td><td>1</td></tr></tbody></table>

The resulting IMEX pair is third order accurate.

$$ \begin{array}{c c c c c c} \begin{array}{c} 0.4358665215 \\ 0.2820667392 \\ 1.208496649 \end{array} & \begin{array}{c} 0 \\ 0.4358665215 \\ -0.644363171 \end{array} & \begin{array}{c} 0 \\ 0 \\ 0.4358665215 \end{array} & \begin{array}{c} 0.4358665215 \\ 0.7179332608 \\ 1 \end{array} & & \\ \hline \begin{array}{c} 0.4358665215 \\ 0.3212788860 \\ -0.105858296 \end{array} & \begin{array}{c} 0 \\ 0.3966543747 \\ 0.5529291479 \end{array} & \begin{array}{c} 0 \\ 0 \\ 0.5529291479 \end{array} & \begin{array}{c} 0 \\ 0 \\ 0 \end{array} & & \\ \hline \begin{array}{c} 0 \\ 1.208496649 \end{array} & \begin{array}{c} 0 \\ -0.644363171 \end{array} & \begin{array}{c} 0.4358665215 \end{array} & \begin{array}{c} 0 \\ 0.4358665215 \\ 0.7179332608 \\ 1 \end{array} & & \end{array} \quad (3.25) $$

• A four stage, L-stable DIRK method with rational coefficients. The corresponding ERK is a five stage third order method. The resulting IMEX method is third order.

$$ \begin{array}{c c c c c c} \begin{array}{c} \frac{1}{2} \\ \frac{1}{6} \\ -\frac{1}{2} \\ \frac{3}{2} \\ \frac{3}{2} \end{array} & \begin{array}{c} 0 \\ \frac{1}{2} \\ \frac{1}{2} \\ -\frac{3}{2} \\ -\frac{3}{2} \end{array} & \begin{array}{c} 0 \\ 0 \\ \frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \end{array} & \begin{array}{c} 0 \\ \frac{1}{2} \\ \frac{3}{2} \\ \frac{1}{2} \\ 1 \end{array} & \begin{array}{c} 0 \\ \frac{1}{2} \\ \frac{11}{18} \\ \frac{5}{6} \\ \frac{1}{4} \end{array} & \begin{array}{c} 0 \\ 0 \\ \frac{1}{18} \\ -\frac{5}{6} \\ \frac{7}{4} \end{array} & \begin{array}{c} 0 \\ 0 \\ 0 \\ \frac{1}{2} \\ \frac{3}{4} \end{array} & \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{array} & \begin{array}{c} 0 \\ \frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \\ 1 \end{array} \end{array} \quad (3.26) $$

We tested these IMEX methods in the case of the KdV equation with $α = β = δ = 1$, $γ = 0$. In Table 1, we summarize the constraints for the timestep $Δt$, purely in term of $Δx$, to obtain a stable solution. IMEX methods (3.25) and (3.26) exhibit excellent stability behavior.

## 4. NUMERICAL RESULTS

In this section we present a series of numerical results aiming to show the performance and robustness of discretization procedures described above. There are many possible combinations of numerical fluxes, types of reconstruction and slope limiter functions. We begin by examining the accuracy of the methods by measuring the convergence rates in Section 4.1 and the preservation of the invariants in Section 4.2. The ability of the schemes to capture a solitary wave solution is demonstrated in Section 4.3. Solitary wave collisions are studied in Section 4.4. Finally, a dispersive shock wave formation is investigated in Section 4.5.
<table><caption>TABLE 2. Rates of convergence : CF-flux</caption><thead><tr><th colspan="3">(a) UNO2 MinMod</th><th colspan="3">(b) WENO3</th></tr><tr><th>Δx</th><th>Rate(E<sub>h</sub><sup>2</sup>)</th><th>Rate(E<sub>h</sub><sup>∞</sup>)</th><th>Δx</th><th>Rate(E<sub>h</sub><sup>2</sup>)</th><th>Rate(E<sub>h</sub><sup>∞</sup>)</th></tr></thead><tbody><tr><td>0.5</td><td>2.000</td><td>2.015</td><td>0.5</td><td>2.604</td><td>2.561</td></tr><tr><td>0.25</td><td>2.001</td><td>2.014</td><td>0.25</td><td>2.790</td><td>2.810</td></tr><tr><td>0.125</td><td>2.001</td><td>2.012</td><td>0.125</td><td>2.905</td><td>2.913</td></tr><tr><td>0.0625</td><td>2.001</td><td>2.010</td><td>0.0625</td><td>2.974</td><td>2.981</td></tr><tr><td>0.03125</td><td>2.001</td><td>2.008</td><td>0.03125</td><td>2.968</td><td>2.995</td></tr></tbody></table>

**Remark 3.** The solution of the linear system involved in (3.15) and (3.21) is obtained by a variation of Gauss elimination for tridiagonal systems with computational complexity $O(d)$, $d$-being the dimension of the system.

**4.1. Rates of convergences, accuracy test.** We consider an initial value problem for (2.1) with periodic boundary conditions in $[-100, 100]$. We take for simplicity $\alpha = \beta = \gamma = \delta = 1$ and consider a solitary wave solution of the form (2.2) with $c_s = 1.1$. We take a uniform mesh $h = \Delta x = 200/N$ and compute the solution up to $T = 100$ using the three stage third order explicit SSP-RK method (3.17) with time step $\Delta t = T/M$. The errors are measured using the discrete scaled norms $E_h^2$ and $E_h^\infty$, [30]

$$ 
\begin{aligned}
E_h^2(k) &= \|U^k\|_h / \|U^0\|_h, \quad \|U^k\|_h = \left( \sum_{i=1}^N \Delta x |U_i^k|^2 \right)^{1/2}, \\
E_h^\infty(k) &= \|U^k\|_{h,\infty} / \|U^0\|_{h,\infty}, \quad \|U^k\|_{h,\infty} = \max_{i=1,...,N} |U_i^k|,
\end{aligned}
$$

where $U^k = \{U_i^k\}_{i=1}^N$ denotes the solution of the fully-discrete scheme (3.16) at the time $t^k = k \Delta t$. The numerical rate of convergence is defined by

$$ \text{Rate} = \frac{\log(E_{h_1}/E_{h_2})}{\log(h_1/h_2)}, $$

for two different mesh sizes $h_1, h_2$.

We perform several tests using the TVD2, UNO2 and WENO3 reconstructions. Numerical solutions are computed with CF, KT or average fluxes. Table 2 shows the rates of convergence for the CF-scheme along with UNO2 and WENO3 reconstructions. We observe the theoretical 2nd order convergence for the average, TVD2 (not reported) and UNO2 schemes. The WENO3 reconstruction in conjunction with improved elliptic inversion scheme (3.14) gives us the expected 3rd order convergence. Rates in Table 2 are obtained with the most dissipative MinMod limiter function, while other limiters yield slightly sharper results. Moreover, the convergence results for the average $m$-flux and the KT numerical flux are qualitatively identical to those of CF. Analogous convergence rates were obtained using the IMEX methods.
4.2. **Invariants preservation.** As already mentioned in Section 2, (2.1) admits at least two quantities (2.3) which remain constant under the equation dynamics. We investigate the conservation of these quantities by computing their discrete counterparts:

$$I_1^h = \Delta x \sum_i U_i, \quad I_2^h = \Delta x \sum_i \left( U_i^2 + \gamma \left[ \frac{U_{i+1} - U_i}{\Delta x} \right]^2 \right). \quad (4.1)$$

The observation of invariants during numerical computations (3.16) may also give an idea on the overall discretization accuracy.

The initial value problem for (2.1) with periodic boundary conditions is considered. We set $\alpha = \beta = \gamma = \delta = 1$ and consider a solitary wave solution with celerity $c_s = 1.5$. We compute its evolution up to $T = 200$ using $\Delta x = 0.1$ and $\Delta t = \Delta x/2$.

The first observation is that the mass of the solitary wave $I_1^h = 13.41640786499$ is preserved in all computations independently from the choice either of the numerical flux, reconstruction method, or the slope limiter function.

The behavior of $I_2^h$ is quite different. Figure 1 shows the evolution of the solitary wave amplitude and of the invariant $I_2^h$. The numerical solution is obtained using $\mathcal{F}^m$, $\mathcal{F}^{CF}$ and $\mathcal{F}^{KT}$ numerical fluxes along with TVD2 and UNO2 reconstructions. The limiter MinMod is used and the dispersive flux is computed with $\mathcal{G}^{lm}$ flux function. The behavior of CF and KT schemes is almost identical. Perhaps, the CF-scheme is slightly less dissipative than the KT-scheme. However, the m-scheme appears to be the least dissipative.

For both KT and CF fluxes, the TVD2 reconstruction preserves neither the invariant $I_2^h$ nor the amplitude of the solitary wave. In the same time UNO2 reconstruction shows excellent behavior. Despite its simplicity, the m-scheme, using $\mathcal{F}^m$ and $\mathcal{G}^m$, performs very well too in preserving $I_2^h$ and the solitary wave amplitude.

In Figure 2 we show the influence of the dispersive flux $\mathcal{G}^m$, $\mathcal{G}^{lm}$ choice. One observes that $\mathcal{G}^{lm}$ flux shows better behavior than the simpler $\mathcal{G}^m$ flux. A comparable performance is achieved with CF-scheme using WENO3 and WENO5 reconstructions.

Finally, in Figure 3 we show a comparison between the various slope limiter functions (Minmod, Van Albada, Van Leer and MC) tested with CF-scheme. MinMod limiter exhibits a small dissipative effect, while other limiters we tested show comparable behavior. The choice of the time-stepping method do not induce any difference.

4.3. **Propagation of solitary waves.** We continue the presentation of numerical results by the classical test-case of a solitary wave propagation. This class of solutions (2.2) plays a very important role in the nonlinear physics and any practical numerical scheme should be able to compute with good accuracy this type of solutions. For simplicity, we will set to unity all coefficients $\alpha = \beta = \gamma = \delta = 1$ in (2.1).

A large-amplitude solitary wave travels rightwards with the speed $c_s = 1.5$. Its propagation is computed up to $T = 100$ with discretization parameters $\Delta x = \Delta t = 0.1$ using KT and CF numerical fluxes and TVD2 reconstruction. In both cases we use the Van Albada limiter. In Figure 4 we compare the analytical solution with the numerical one. Figure 4(b) is a magnification of the solitary pulse showing that the solitary wave shape is
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Evolution of amplitude and $I_2^h$ with $G^{lm}$ flux and Minmod limiter. '$∇$: CF-TVD2, '$◇$: CF-UNO2, '$□$: KT-TVD2, '$×$: KT-UNO2, '$○$: m-scheme</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Evolution of amplitude and $I_2^h$, UNO2 reconstruction with Minmod limiter: '$∇$: $\mathcal{F}^{CF}-G^{lm}$, '$◇$: $\mathcal{F}^{CF}-G^m$, '$□$: $\mathcal{F}^{KT}-G^{lm}$, '$×$: $\mathcal{F}^{KT}-G^m$, '$▷$: $\mathcal{F}^{CF}$-WENO3, '$▷$: $\mathcal{F}^{CF}$-WENO5, '$○$: $F^m-G^m$. (Notice the scale difference on the vertical axis with respect to Figure 1).</figcaption></figure>

perfectly retained. Also we note that up to the graphical resolution, all curves are undistinguishable. In order to observe the differences between these solutions we present in 4(c) the error $E_\ell = \log_{10}|u_{exact}(x, 100) - U(x, 100)|$. This shows that the difference between the numerical and the exact solution is analogous in all the cases and very small.

The behavior of the numerical solutions can be better understood by analyzing the so-called *effective* equation, that is the p.d.e that the numerical scheme satisfies up to the order of the method. Obtaining an *effective* equation is not always feasible. In the case
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Evolution of amplitude and $I_2^h$, $\mathcal{F}^{CF}-\mathcal{G}^{lm}$ fluxes and UNO2 reconstruction: '∇': Minmod, '◊': MC, '□': Van Albada, '○': Van Leer. (Notice the scale difference on the vertical axis with respect to Figure 1).</figcaption></figure>

of the m-scheme for the KdV-BBM equation (2.1), the numerical solution $u_h$ satisfies the following effective equation:

$$ u_{h,t} + \alpha u_{h,x} + \beta u_h u_{h,x} - \gamma u_{h,xxt} + \delta u_{h,xxx} + \Delta x^2 \left( \frac{\alpha}{6} u_{h,xxx} + \frac{\beta}{6} u_h u_{h,xxx} + \frac{\beta}{4} u_{h,x} u_{h,xx} + \frac{\delta}{4} u_{h,xxxx} - \frac{\gamma}{12} u_{h,xxxxt} \right) = 0. \quad (4.2) $$

On Figure 5 we illustrate some artifacts of the numerical discretization for the pure BBM equation ($\delta = 0$). In Figure 5(a) one can observe a small dispersive tail coming mainly from nonlinear terms discretization. The amplitude of the tail is related to the order of the method. Taking $\Delta x$ ten times smaller leads the reduction of the amplitude by two orders of magnitude, as it can be observed on Figure 5(b). The explanation of these phenomena is contained in the straightforward analysis of the effective equation (4.2).

We underline that the smallest tail is produced by the m-scheme and the largest by the KT-scheme. This shortcoming can be further reduced by UNO2 or WENO3 reconstruction procedures. We conclude that a detailed study of solitary wave interactions would require a combination of a higher order method with a finer grid resolution.

### 4.4. Solitary wave overtaking collisions.
The solitary wave solutions (also known as solitons) of the celebrated KdV equation ($\alpha = \beta = \delta = 1$, $\gamma = 0$) have a well-known property to interact in an elastic way during an overtaking collision. In other words, the solitary waves retain their initial shape after the interaction, cf. [14]. Contrary to the KdV equation, the overtaking collision of two solitary waves of the BBM model and in general of the KdV-BBM equation is not elastic. Interacting solitary waves change in shape and also a small dispersive tail appears after the process. However, a nonlinear phase shift can be still observed even in the KdV-BBM equation.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Comparison between the analytical and numerical solutions: ...: analytical solution, -: CF-TVD2, - -: KT-TVD2, - -: m-scheme.</figcaption></figure>

Here we study the overtaking collision of two solitary waves of the KdV-BBM equation with $\alpha = \beta = \gamma = \delta = 1$. Solitary waves are located initially at $X_1 = -50$ and $X_2 = 50$ with speeds $c_s = 1.5$ and $c_s = 1.1$ respectively. At $t = 0$ we have two well separated pulses and the wave behind (left) propagates faster. Space and time variables are discretized with $\Delta x = \Delta t = 0.01$ to capture this process accurately. The solution is computed using the CF-scheme and three types of reconstruction: TVD2 with Van Albada limiter, UNO2 reconstruction with MinMod limiter and WENO3 method, and with the third order explicit SSP-RK method.

The invariant $I_1^h = 18.915498698$ is conserved with the digits shown in all cases. With the invariant $I_2^h$ the situation is slightly different: UNO2 and WENO3 schemes preserved the value $I_2^h = 15.0633$, while the more dissipative TVD2 reconstruction yields $I_2^h = 15.063$.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Dispersive artifacts of the equivalent equation: ...: Analytical solution, -: FCVF-TVD2, --:KT-TVD2, .-.:FCVF-WENO3</figcaption></figure>

Figure 6 shows the interaction process at several time instances in the left column, while the right column shows the corresponding magnification of the dispersive tail. Essentially no difference can be observed among various numerical solutions even in the magnified region, up to the graphical resolution. Additional snapshots aiming to illustrate the interaction process are shown on Figure 7. We observe that the solitary waves propagate connected as a single pulse with a single maximum for a small time interval contrary to bidirectional models [16] and to Euler equations (cf. [13]).

Figure 8 shows the “elastic” collision of two solitons of the KdV equation ($α = β = δ = 1$, $γ = 0$) up to $t = 600$. In this experiment we took $Δx = Δt = 0.01$ and 0.005 using IMEX method (3.25). Contrary to the analogous collision in the case of the BBM equation, we do not observe any new dispersive tails. Further magnification of the images show small artifacts of the order $O(10^{-6})$. The invariants are $I_1^h = 12.280014566440$ and $I_2^h = 9.244$ for all the computations with $Δx = 0.01$. When a finer grid is considered, $Δx = 0.005$ we do not observe any improvement in the conservation of the invariant $I_1^h$ while $I_2^h$ was 9.2442. Analogous conservation properties observed when we studied the collision for the KdV-BBM equation with the IMEX method (3.25) we observed that $I_1^h = 18.915498698945$ but no other improvement in the invariant $I_2^h = 15.0633$.

### 4.5. Dispersive shock formation.
It was proven that smooth solutions to the KdV equation tend to become highly oscillatory as the parameter $δ$ tends to zero, cf. [49]. These oscillatory solutions are sometimes referred to in the literature as dispersive shock waves. In this section we study numerically this special class of solutions. Recently, a discontinuous Galerkin method was employed to study the same problem [51] in the classical setting of the KdV equation.

Namely we consider the KdV-BBM equation with $α = β = 1$, $γ = 10^{-5}$ and $δ = 0$. A solitary wave solution (2.2) is taken as an initial condition with parameters $α = β = γ = 1$, $δ = 0$ and $c_s = 1.3$. We underline that this initial condition is not an exact solution to
<figure><img src="image_7.png" /><figcaption>(a) $t = 0$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) $t = 0$ (magnification)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) $t = 200$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) $t = 200$ (magnification)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(e) $t = 350$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(f) $t = 350$ (magnification)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(g) $t = 600$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(h) $t = 600$ (magnification)</figcaption></figure>

<caption>FIGURE 6. Inelastic overtaking collision of two solitary waves for the KdV-BBM equation</caption>
<figure><img src="image_15.png" /><figcaption>FIGURE 7. Inelastic overtaking collision of two solitary waves for the KdV-BBM equation (detailed view)</figcaption></figure>

the BBM equation under consideration, since the coefficient $\gamma$ is different. A fine grid with $\Delta x = 0.001$ is required to observe this phenomenon. We note that even much more accurate schemes [51] require almost the same resolution. Figure 9 shows the formation of a dispersive shock wave. The numerical solution is computed with four different methods: the m-scheme and CF-scheme with TVD2, UNO2 and WENO5 reconstructions. The KT flux was also tested, producing almost identical to that of the CF-scheme. In all the cases we took $\Delta t = \Delta x/10$ except in the case of the WENO5 reconstruction where $\Delta t = \Delta x/2$.

The invariant $I_1^h = 7.493997530$ conserving the digits shown during all simulations for all numerical schemes we tested. The behavior of $I_2^h$ is considerably different. Figure 10 (left) shows that from the time the dispersive shock was formed, all numerical schemes, except the m-scheme, lose the conservation of the invariant $I_2^h$. As for the m-scheme the $I_2^h$ invariant was conserved to one decimal digit, during the whole simulation, see Figure 10 (b).

On the other hand, when a solitary wave solution evolves for longer time intervals, using for example the m-scheme, we observe that solitary-wave-like structures are formed, cf. Figure 11, while retaining the conservation of the invariant $I_2^h$ up to one digit. Analogous
<figure><img src="image_16.png" /><figcaption>(a) $t = 0$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) $t = 0$ (magnification)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) $t = 200$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(d) $t = 200$ (magnification)</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(e) $t = 350$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(f) $t = 350$ (magnification)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(g) $t = 600$</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(h) $t = 600$ (magnification)</figcaption></figure>

<caption>FIGURE 8. Elastic overtaking collision of two solitary waves computed with the KdV equation</caption>
<figure><img src="image_24.png" /><figcaption>FIGURE 9. Near the zero dispersion limit, BBM equation</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>FIGURE 10. Evolution of $I_2^h$</figcaption></figure>

behavior is observed for the KdV equation where general initial conditions evolved into
series of solitary waves, cf. [14].

In Figure 12 we present the same experiment for the KdV equation ($\alpha = \beta = 1$, $\gamma = 0$, $\delta = 10^{-5}$) where the time integration is performed with the IMEX method (3.25) up to
<figure><img src="image_26.png" /><figcaption>FIGURE 11. Near the zero dispersion limit</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>FIGURE 12. Near the zero dispersion limit, KdV equation</figcaption></figure>
$T = 20$ with discretization parameters $\Delta x = 0.001$ and $\Delta t = \Delta x/2$. We observe that the invariant $I_2^h$ is conserved with slightly, less accuracy while the $I_1^h = 6.572670686045$. When we use the IMEX method (3.25) and the m-scheme in the case of the BBM equation we observe that the invariant $I_2^h$ conserves 2 digits $I_2^h = 4.49$, while $I_1^h = 7.493997530374$ conserving the digits shown. Thus we conclude that in this experiment (as also observed in all previous ones) the use of the IMEX method might improve the conservation of mass.

## 5. CONCLUSIONS

The main scope of the present article is to extend the framework of finite volume methods to scalar unidirectional dispersive models. We chose the celebrated BBM-KdV equation (2.1) as an important representative model arising in the water wave theory and having all main features of dispersive wave equations.

The BBM-KdV equation can be also viewed as a dispersive perturbation of the inviscid Burgers equation. Consequently, our method relies on classical finite volume schemes which discretize the advection operator. Then, a special treatment was proposed for the KdV-dispersion term, while the BBM-dispersion required an elliptic operator inversion per each time step, hence, providing a physical regularization to numerical solutions. We propose and implement also several methods to obtain high order accurate schemes.

The proposed discretization procedure is validated by comparisons with an analytical solitary wave solution. The order of convergence is measured as well as invariant preservation is studied extensively. The numerical method is applied to several important test cases such as a solitary wave propagation and a dispersive shock formation. We make also use of proposed higher order extensions to study the overtaking solitary waves collision for the KdV-BBM equation.

The extension to more realistic bi-directional wave propagation models such as Boussinesq type equations [40, 36, 32, 16].

## ACKNOWLEDGEMENTS

D. Dutykh acknowledges the support from French Agence Nationale de la Recherche, project MathOcean (Grant ANR-08-BLAN-0301-01), Ulysses Program of the French Ministry of Foreign Affairs under the project 23725ZA and CNRS PICS project No. 5607. The work of Th. Katsaounis was partially supported by European Union FP7 program Capacities(Regpot 2009-1), through ACMAC (http://acmac.tem.uoc.gr).

## REFERENCES

[1] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri. Implicit-Explicit Runge-Kutta methods for time-dependent partial differential equations. *Applied Numerical Mathematics*, 25:151–167, 1997. 2, 8, 9

[2] P. Avilez-Valente and F. J. Seabra-Santos. A high-order Petrov-Galerkin finite element method for the classical Boussinesq wave model. *Int. J. Numer. Meth. Fluids*, 59:969–1010, 2009. 1

[3] G. Bellotti and M. Brocchini. On the shoreline boundary conditions for Boussinesq-type models. *Int. J. Num. Meth. in Fluids*, 37(4):479–500, 2001. 2
[4] T. B. Benjamin, J. L. Bona, and J. J. Mahony. Model equations for long waves in nonlinear dispersive systems. *Philos. Trans. Royal Soc. London Ser. A*, 272:47–78, 1972. 1, 3

[5] F. Benkhaldoun and L. Quivy. A Non Homogeneous Riemann Solver for shallow water and two phase flows. *Flow, Turbulence and Combustion*, 76:391–402, 2006. 5

[6] F. Benkhaldoun and M. Seaid. New finite-volume relaxation methods for the third-order differential equations. *Commun. Comput. Phys.*, 4:820–837, 2008. 2

[7] J. Bona and V. Dougalis. An initial- and boundary value problem for a model equation for propagation of long waves. *J. Math. Anal. and Applic.*, 75:503–522, 1980. 3

[8] J. L. Bona, M. Chen, and J.-C. Saut. Boussinesq equations and other systems for small-amplitude long waves in nonlinear dispersive media. I: Derivation and linear theory. *Journal of Nonlinear Science*, 12:283–318, 2002. 1, 3

[9] J. L. Bona, V. A. Dougalis, and D. E. Mitsotakis. Numerical solution of KdV-KdV systems of Boussinesq equations: I. The numerical scheme and generalized solitary waves. *Mat. Comp. Simul.*, 74:214–228, 2007. 1

[10] J. L. Bona, W. G. Pritchard, and L. R. Scott. Numerical schemes for a model for nonlinear dispersive waves. *J. Comput. Phys.*, 60:167–186, 1985. 1

[11] J. Boussinesq. Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond. *J. Math. Pures Appl.*, 17:55–108, 1872. 1

[12] Y. Brenier and D. Levy. Dissipative behavior of some fully non-linear KdV-type of equations. *Physica D*, 137(3-4):277-294, 2000. 2

[13] W. Craig, P. Guyenne, J. Hammack, D. Henderson, and C. Sulem. Solitary water wave interactions. *Phys. Fluids*, 18(5):57106, 2006. 16

[14] P. G. Drazin and R. S. Johnson. *Solitons: An introduction*. Cambridge, UK: Univ. Pr. (1989) 226 p, 1989. 14, 20

[15] D. Dutykh and F. Dias. Dissipative Boussinesq equations. *C. R. Mecanique*, 335:559-583, 2007. 1

[16] D. Dutykh, Th. Katsaounis, and D. Mitsotakis. Finite volume schemes for dispersive wave propagation and runup. *J. Comput. Phys*, 230:3035–3061, 2011. 3, 16, 22

[17] K. S. Erduran, S. Ilic, and V. Kutija. Hybrid finite-volume finite-difference scheme for the solution of Boussinesq equations. *Int. J. Numer. Meth. Fluids*, 49:1213–1232, 2005. 2

[18] C. Eskilsson and S. J. Sherwin. Spectral/hp discontinuous Galerkin methods for modelling 2D Boussinesq equations. *J. Comput. Phys*, 212(2):566–589, 2006. 1

[19] R. Fetecau and D. Levy. Aproximate model equations for water waves. *Comm. Math. Sci.*, 3:159–170, 2005. 3

[20] J.-M. Ghidaglia. Flux schemes for solving nonlinear systems of conservation laws. In J J Chattot and M Hafez, editors, *Proceedings of the meeting in honor of P.L. Roe*, Arcachon, July 1998. 2, 5

[21] J.-M. Ghidaglia. *Innovative Methods for Numerical Solution of Partial Differential Equations*, chapter Flux schem. World Scientific, Singapore, 2001. 5

[22] J.-M. Ghidaglia, A. Kumbaro, and G. Le Coq. Une méthode volumes-finis à flux caractéristiques pour la résolution numérique des systèmes hyperboliques de lois de conservation. *C. R. Acad. Sci. I*, 322:981–988, 1996. 2, 5

[23] J.-M. Ghidaglia, A. Kumbaro, and G. Le Coq. On the numerical solution to two fluid models via cell centered finite volume method. *Eur. J. Mech. B/Fluids*, 20:841–867, 2001. 5

[24] S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong Stability-Preserving High-Order Time Discretization Methods. *SIAM Review*, 43:89–112, 2001. 2, 7

[25] A. Harten, P. D. Lax, and B. van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. *SIAM Review*, 25:35–61, 1983. 5

[26] A. Harten and S. Osher. Uniformly high-order accurate nonscillatory schemes, I. *SIAM J. Numer. Anal.*, 24:279–309, 1987. 2, 6
[27] N. E. Kolgan. Application of the minimum-derivative principle in the construction of finite-difference schemes for numerical analysis of discontinuous solutions in gas dynamics. *Uchenye Zapiski TsaGI* [Sci. Notes Central Inst. Aerodyn], 3(6):68–77, 1972. 6

[28] D. J. Korteweg and G. de Vries. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves. *Phil. Mag.*, 39(5):422–443, 1895. 1, 3

[29] A. Kurganov and E. Tadmor. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations. *J. Comput. Phys.*, 160(1):241–282, 2000. 2, 5

[30] D. Levy, C.-W. Shu, and J. Yan. Local discontinuous Galerkin methods for nonlinear dispersive equations. *J. Comput. Phys.*, 196(2):751–772, 2004. 1, 11

[31] X.-D. Liu, S. Osher, and T. Chan. Weighted essentially non-oscillatory schemes. *J. Comp. Phys.*, 115:200–212, 1994. 2, 6

[32] P. A. Madsen, H. B. Bingham, and H. A. Schaffer. Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis. *Proc. R. Soc. Lond. A*, 459:1075–1104, 2003. 22

[33] D. E. Mitsotakis. Boussinesq systems in two space dimensions over a variable bottom for the generation and propagation of tsunami waves. *Math. Comp. Simul.*, 80:860–873, 2009. 1

[34] H. Nessyahu and E. Tadmor. Nonoscillatory central differencing for hyperbolic conservation laws. *J. Comp. Phys.*, 87(2):408–463, 1990. 2, 5

[35] H. Y. Nguyen and F. Dias. A Boussinesq system for two-way propagation of interfacial waves. *Physica D*, 237(18):2365–2389, 2008. 1

[36] O. Nwogu. Alternative form of Boussinesq equations for nearshore wave propagation. *J. Waterway, Port, Coastal and Ocean Engineering*, 119:618–638, 1993. 22

[37] S. Osher. Riemann solvers, the entropy condition, and difference approximations. *SIAM J. Numer. Anal.*, 21(2):217–235, 1984. 5

[38] H. T. Ozkan-Haller and J. T. Kirby. A Fourier-Chebyshev collocation method for the shallow water equations including shoreline runup. *Applied Ocean Research*, 19:21–34, 1997. 1

[39] B. Pelloni and V. Dougalis. Numerical modelling of two-way propagation of nonlinear dispersive waves. *Mat. Comp. Simul.*, 55:595–606, 2001. 1

[40] D. H. Peregrine. Long waves on a beach. *J. Fluid Mech.*, 27:815–827, 1967. 1, 3, 22

[41] P. L. Roe. Approximate Riemann solvers, parameter vectors and difference schemes. *J. Comput. Phys.*, 43:357–372, 1981. 5

[42] J. B. Shiach and C. G. Mingham. A temporally second-order accurate Godunov-type scheme for solving the extended Boussinesq equations. *Coastal Engineering*, 56:32–45, 2009. 2

[43] C.-W. Shu. Essentially non-oscillatory and weighted Essentially non-oscillatory schemes for hyperbolic conservation laws. In *Advanced Numerical Approximation of Nonlinear Hyperbolic Equations*, pages 325–432. Springer Berlin / Heidelberg, 1997. 6

[44] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes. *J. Comput. Phys.*, 77:439–471, 1988. 2, 7

[45] R. J. Spiteri and S. J. Ruuth. A New Class of Optimal High-Order Strong-Stability-Preserving Time Discretization Methods. *SIAM Journal on Numerical Analysis*, 40:469–491, 2002. 7

[46] P. K. Sweby. High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws. *SIAM J. Numer. Anal.*, 21(5):995–1011, 1984. 2, 6

[47] M. Tonelli and M. Petti. Hybrid finite-volume finite-difference scheme for 2DH improved Boussinesq equations. *Coastal Engineering*, 56:609–620, 2009. 2

[48] B. van Leer. Towards the ultimate conservative difference scheme V: a second order sequel to Godunov' method. *J. Comput. Phys.*, 32:101–136, 1979. 6

[49] S. Venakides. The zero dispersion limit of the Korteweg-de Vries equation with periodic initial data. *AMS Transactions*, 301:189–226, 1987. 16
[50] G. Wei, J. T. Kirby, S. T. Grilli, and R. Subramanya. A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves. *J. Fluid Mech.*, 294:71–92, 1995. 1

[51] J. Yan and C.-W. Shu. A local discontinuous Galerkin method for KdV type equations. *SIAM J. Num. Anal.*, 40:769–791, 2002. 1, 16, 18

LAMA, UMR 5127 CNRS, UNIVERSITÉ DE SAVOIE, CAMPUS SCIENTIFIQUE, 73376 LE BOURGET-DU-LAC CEDEX, FRANCE

*E-mail address:* Denys.Dutykh@univ-savoie.fr

URL: http://www.lama.univ-savoie.fr/~dutykh/

DEPARTMENT OF APPLIED MATHEMATICS, UNIVERSITY OF CRETE, HERAKLION, 71409 GREECE, INST. OF APP. AND COMP. MATH. (IACM), FORTH, HERAKLION, 71110, GREECE

*E-mail address:* thodoros@tem.uoc.gr

URL: http://www.tem.uoc.gr/~thodoros/

IMA, UNIVERSITY OF MINNESOTA, MINNEAPOLIS MN 55455, USA

*E-mail address:* dmitsot@gmail.com

URL: http://sites.google.com/site/dmitsot/