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

HAL
open science

# Finite volume schemes for dispersive wave propagation and runup

Denys Dutykh, Theodoros Katsaounis, Dimitrios Mitsotakis

► To cite this version:

Denys Dutykh, Theodoros Katsaounis, Dimitrios Mitsotakis. Finite volume schemes for dispersive wave propagation and runup. Journal of Computational Physics, 2011, 230 (8), pp.3035-3061. 10.1016/j.jcp.2011.01.003. hal-00472431v5

HAL Id: hal-00472431

https://hal.science/hal-00472431v5

Submitted on 9 Jan 2011

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# FINITE VOLUME SCHEMES FOR DISPERSIVE WAVE PROPAGATION AND RUNUP

DENYS DUTYKH*, THEODOROS KATSAOUNIS, AND DIMITRIOS MITSOTAKIS

ABSTRACT. Finite volume schemes are commonly used to construct approximate solutions to conservation laws. In this study we extend the framework of the finite volume methods to dispersive water wave models, in particular to Boussinesq type systems. We focus mainly on the application of the method to bidirectional nonlinear, dispersive wave propagation in one space dimension. Special emphasis is given to important nonlinear phenomena such as solitary waves interactions, dispersive shock wave formation and the runup of breaking and non-breaking long waves.

## 1. INTRODUCTION

The simulation of water waves in realistic and complex environments is a very challenging problem. Most of the applications arise from the areas of coastal and naval engineering, but also from natural hazards assessment. These applications may require the computation of the wave generation [DD07, KDD07], propagation [TG97], interaction with solid bodies, the computation of long wave runup [TS94, TS98] and even the extraction of the wave energy [Sim81]. Issues like wave breaking, robustness of the numerical algorithm in wet-dry processes along with the validity of the mathematical models in the near-shore zone are some basic problems in this direction [HP79]. During past years, the classical shallow water equations have been employed to solve some of these problems [AC99, DPD10]:

$$ 
\begin{aligned}
& H_t + (Hu)_x = 0, \\
& (Hu)_t + (Hu^2 + \frac{g}{2}H^2)_x = gHD_x,
\end{aligned} 
\quad (1.1) 
$$

where $H(x, t) := \eta(x, t) + D(x)$ is the total water depth, $D(x)$ describes the depth below the mean sea level while $\eta(x, t)$ is the free surface elevation, $u(x, t)$ denotes the depth-averaged fluid velocity and $g$ is the gravity acceleration constant. Mathematically, equations (1.1) represent a system of conservation laws describing the propagation of infinitely long waves with a hydrostatic pressure assumption. The wave breaking phenomenon is commonly assimilated to the formation of shock waves (or hydraulic jumps) which is a common feature of hyperbolic p.d.e's. Consequently, the finite volume(FV) method has become the method of choice for these problems due to its excellent intrinsic conservative and shock-capturing properties [AC99, DK03, DKK08, DPD10]. Furthermore the shallow water equations have been proven in practice to predict accurately the maximum runup of long waves [HH70, Syn87, TS94, TS98, KS98, CFVCP04].

*Key words and phrases.* finite volume method; dispersive waves; solitary waves; runup; water waves.

* Corresponding author.
On the other hand, various studies have shown that the inclusion of dispersive effects is beneficial for the description of long wave propagation and runup processes [Zel91, WKGS95, LR02]. Moreover, J.A. Zelt [Zel91] reported a divergence in the prediction of the rundown and in the prediction of the reflected wave-train after the wave climbing on the shore when a dispersionless model is employed. According to J.A. Zelt, the results of the nonlinear dispersive model considered in [Zel91] showed better performance compared to (1.1). During the last fifty years numerous dispersive models have been proposed for the simulation of long waves [Ser53, Per67, BS76, Nwo93, KCKD00, MBS03, Mit09].

In this work we will study numerically bidirectional water wave models. Specifically, we consider the following family of Boussinesq type systems of water wave theory, introduced in [BCS02], written in nondimensional, unscaled variables

$$ 
\begin{aligned}
\eta_t + u_x + (\eta u)_x + a u_{xxx} - b \eta_{xxt} &= 0, \\
u_t + \eta_x + uu_x + c \eta_{xxx} - d u_{xxt} &= 0,
\end{aligned} 
\qquad (1.2) 
$$

where $a, b, c, d \in \mathbb{R}$, $\eta = \eta(x, t)$, $u = u(x, t)$ are real functions defined for $x \in \mathbb{R}$ and $t \ge 0$.

For more realistic situations we introduce a new Boussinesq type system with variable bottom topography based on Peregrine's system, [Per67]. The new system incorporates a very important property — the invariance under vertical translations, thus more appropriate for practical applications such as wave runup on non-uniform shores. In dimensional variables the model reads

$$ 
\begin{aligned}
H_t + Q_x &= 0, \\
Q_t + \left(\frac{Q^2}{H} + \frac{g}{2}H^2\right)_x - \left[\frac{H^2}{3}Q_{xxt} - \left(\frac{1}{3}H_x^2 - \frac{1}{6}HH_{xx}\right)Q_t + \frac{1}{3}HH_xQ_{xt}\right] &= gHD_x,
\end{aligned} 
\qquad (1.3) 
$$

where $H(x, t) = \eta(x, t) + D(x)$, $Q(x, t) = H(x, t)u(x, t)$.

There is a wide range of numerical methods in the literature for computing approximate solutions to these models. Finite difference (FD) schemes [KWC+98, JFBM06, HMK+09], finite element methods [BDM07, Mit09, DM08, AVSS09] and spectral methods [PD01, Ngu08] have been proposed. More contemporary discontinuous Galerkin (DG) schemes have also been adapted with some success to dispersive wave equations [YS02, LSY04, ES06, ESB06] while the application of Finite Volumes (FV) or hybrid FV/FD methods remain most infrequent for this type of problems. To our knowledge, only a few very recent works are in this direction [BB01, EIK05, BS08, TP09, SM09].

Finite volume method is well known for its accuracy, efficiency and robustness for approximating solutions to conservation laws and in particular to nonlinear shallow water equations (1.1). The aforementioned bidirectional models (1.2) and (1.3) are rewritten in a conservative form and discretization by the finite volume method follows. Three different numerical fluxes are employed

* a simple *average flux* (m-scheme),

* a *central flux*, (KT-scheme) [NT90, KT00], as a representative of central schemes,

* a *characteristic flux* (CF-scheme), as a representative of the linearized Riemann solvers, [GKC96],
along with TVD, UNO and WENO reconstruction techniques, [Swe84, HO87, LOC94]. Time discretization is based on Runge-Kutta (RK) methods which preserve the total variation diminishing (TVD) property of the finite volume scheme, [SO88, GST01, SR02]. We use explicit RK methods since we work with BBM type systems (1.2) and not with KdV-type systems which is well known to be notoriously stiff, [Mit09]. These methods have been studied thoroughly in the case of nonlinear conservation laws. The average flux although is known to be unstable for conservation laws is proved to be very accurate for nonlinear dispersive waves. On the other hand finite volume methods based on the central flux as well as on characteristic flux work equally well for the numerical simulation of waves even in realistic environments.

The performance of the finite volume method applied to models (1.2) and to the new system (1.3) is studied in a systematic way through a series of numerical experiments. Our main focus is the evaluation, in terms of accuracy, efficiency and robustness of second order finite volume methods compared to high order schemes. In particular, in this study we take up on the following points

* accuracy of the finite volume method in the propagation of solitary waves with very satisfactory results.

* conservation of various invariant quantities during the formation of dispersive shocks is studied numerically. The finite element as well as spectral methods break down for these experiments. The finite volume method provides very accurate results.

* interactions of solitary waves are computed with high accuracy. It is shown numerically that Boussinesq type systems describe better overtaking collisions of solitary waves than unidirectional models like KdV-BBM. We compare our results, whenever possible, with experimental measurements with very good agreement.

* finite volume method allows to use appropriate techniques to treat the transition from wet to dry regions and vice versa. These techniques are applied successfully to systems with dispersive terms modeling runup of long waves. On the other hand, when the model fails due to wave breaking, the method allows to use locally the nonlinear shallow water system, thus enabling us to resolve a wide spectrum of hydrodynamic phenomena using a single computational framework.

* it is shown numerically the advantage of using dispersive models over standard nonlinear shallow water equations in computing the wave runup and, in particular, in capturing the reflected wave. It's also illustrated by an example the importance of the system being invariant under vertical translations.

The paper is organized as follows. In Section 2 Boussinesq type systems are presented along with some of their basic properties. A new system with uneven bottom and invariant under vertical translations is derived. In Section 3 the finite volume method is presented for a general framework incorporating all models.

Section 4 presents a series of numerical experiments for the Boussinesq systems (1.2). In this mathematical setting we validate the finite volume method and measure its accuracy. We study the propagation as well as the interaction of solitary waves: we consider in particular head-on and overtaking collisions, but also we present results concerning the small
dispersion effect. Finally, in Section 5 the new system with variable bottom, (1.3) is studied. Numerical simulations of non-breaking and breaking long wave runup are presented and compared with experimental data.

## 2. MATHEMATICAL MODELS

We present briefly the mathematical models being considered and some of their main properties.

### 2.1. Dispersive models with flat bottom.
We consider the following family of Boussinesq type systems of water wave theory, introduced in [BCS02], which may be written in nondimensional, unscaled variables

$$ 
\begin{aligned}
\eta_t + u_x + (\eta u)_x + au_{xxx} - b\eta_{xxt} &= 0, \\
u_t + \eta_x + uu_x + c\eta_{xxx} - du_{xxt} &= 0,
\end{aligned} 
\qquad (2.1) $$

where $\eta = \eta(x, t)$, $u = u(x, t)$ are real functions defined for $x \in \mathbb{R}$ and $t \ge 0$. Coefficients $a, b, c$ and $d$ are defined as

$$ a = \frac{1}{2}(\theta^2 - \frac{1}{3})\nu, \quad b = \frac{1}{2}(\theta^2 - \frac{1}{3})(1 - \nu), \quad c = \frac{1}{2}(1 - \theta^2)\mu, \quad d = \frac{1}{2}(1 - \theta^2)(1 - \mu), \qquad (2.2) $$

where $0 \le \theta \le 1$ and $\mu, \nu \in \mathbb{R}$. The variables in (2.1) are non-dimensional and unscaled: $x$ and $t$ are proportional to position along the channel and time, respectively, while $\eta(x, t)$ and $u(x, t)$ are proportional to the deviation of the free surface above an undisturbed level and to the horizontal velocity of the fluid at a height $y = -1 + \theta(1 + \eta(x, t))$, respectively. In terms of these variables the channel bottom is located at $y = -1$, ($\theta = 0$), while the free surface corresponds to $\theta = 1$. Boussinesq systems (2.1) with $b = d$ conserve the energy functional:

$$ I_1(t) = \int_{\mathbb{R}} (\eta^2(x, t) + (1 + \eta(x, t))u^2(x, t) - c\eta_x^2(x, t) - a u_x^2(x, t)) dx, \qquad (2.3) $$

i.e. $I_1(t) = I_1(0)$ for $t \ge 0$. System (2.1) is derived under the following assumptions:

$$ \varepsilon := A/h_0 \ll 1, \quad \sigma := h_0/\lambda \ll 1 \quad S := \frac{A\lambda^2}{h_0^3} = O(1), $$

where $S$ is the Stokes (or Ursell) number, $A$ is a typical wave amplitude of fluid of depth $h_0$ and $\lambda$ is a characteristic wavelength. If one takes $S = 1$ and switches to scaled, dimensionless variables, one may derive from Euler equations a scaled version of (2.1) by appropriate asymptotic expansion in powers of $\varepsilon$, cf. [BCS02]:

$$ 
\begin{aligned}
\eta_t + u_x + \varepsilon(\eta u)_x + \varepsilon[au_{xxx} - b\eta_{xxt}] &= O(\varepsilon^2), \\
u_t + \eta_x + \varepsilon uu_x + \varepsilon[c\eta_{xxx} - du_{xxt}] &= O(\varepsilon^2),
\end{aligned} 
\qquad (2.4) $$

from which we obtain (2.1) by unscaling and neglecting higher order terms $O(\varepsilon^2)$.

We list several examples of particular Boussinesq systems of the form (2.1) that we will refer to in the sequel. The initial-value problem for all these systems has been shown to be at least nonlinearly well-posed locally in time, cf. [BCS04].
(i) The 'classical' Boussinesq system ($\mu = 0, \theta^2 = 1/3$, i.e. $a = b = c = 0, d = 1/3$ in (2.1)), whose initial-value problem is globally well-posed, [Ami84, Sch81],

$$ 
\begin{aligned}
\eta_t + u_x + (\eta u)_x &= 0, \\
u_t + \eta_x + uu_x - \frac{1}{3}u_{xxt} &= 0.
\end{aligned} 
\qquad (2.5) $$

(ii) The BBM-BBM system ($\nu = \mu = 0, \theta^2 = 2/3$, i.e. $a = c = 0, b = d = 1/6$ in (2.1)), whose initial-value problem is locally well-posed, [BC98],

$$ 
\begin{aligned}
\eta_t + u_x + (\eta u)_x - \frac{1}{6}\eta_{xxt} &= 0, \\
u_t + \eta_x + uu_x - \frac{1}{6}u_{xxt} &= 0.
\end{aligned} 
\qquad (2.6) $$

(iii) The Bona-Smith system ($\nu = 0, \mu = (4 - 6\theta^2)/3(1 - \theta^2)$, i.e. $a = 0, b = d = (3\theta^2 - 1)/6, c = (2 - 3\theta^2)/3, 2/3 < \theta^2 < 1$ in (2.1)), whose initial-value problem is globally well-posed, cf. [BS76]. The limiting form of this system as $\theta \to 1$, corresponding to $a = 0, b = d = 1/3, c = -1/3$, is the system actually studied by Bona and Smith, [BS76]. These systems are given by

$$ 
\begin{aligned}
\eta_t + u_x + (\eta u)_x - \frac{3\theta^2 - 1}{6}\eta_{xxt} &= 0, \\
u_t + \eta_x + uu_x + \frac{2-3\theta^2}{3}\eta_{xxx} - \frac{3\theta^2 - 1}{6}u_{xxt} &= 0.
\end{aligned} 
\qquad (2.7) $$

The existence of solitary wave solutions to the above systems, in some cases the uniqueness also, has been proved in [Che00, Che98, DM04] and in the case of the Bona-Smith type systems (2.7), for each $\theta^2 \in (7/9, 1)$, there exists one solitary wave in closed form, [Che98]

$$ 
\begin{aligned}
\eta(\xi) &= \eta_0 \operatorname{sech}^2(\lambda\xi), \\
u(\xi) &= B \eta(\xi),
\end{aligned} 
\qquad (2.8) $$

with

$$ 
\begin{aligned}
\eta_0 &= \frac{9}{2} \cdot \frac{\theta^2 - 7/9}{1 - \theta^2}, & c_s &= \frac{4(\theta^2 - 2/3)}{\sqrt{2(1 - \theta^2)(\theta^2 - 1/3)}}, \\
\lambda &= \frac{1}{2} \sqrt{\frac{3(\theta^2 - 7/9)}{(\theta^2 - 1/3)(\theta^2 - 2/3)}}, & B &= \sqrt{\frac{2(1 - \theta^2)}{\theta^2 - 1/3}}.
\end{aligned} 
\qquad (2.9) $$

## 2.2. Dispersive models with variable bottom.
For more realistic applications one should consider Boussinesq systems with variable bottom. In his pioneering work Peregrine, [Per67], derived the following Boussinesq type system

$$ 
\begin{aligned}
\eta_t + [(D + \eta)u]_x &= 0, \\
u_t + g\eta_x + uu_x - \frac{D}{2}(Du)_{xxt} + \frac{D^2}{6}u_{xxt} &= 0,
\end{aligned} 
\qquad (2.10) $$

where $\eta(x, t)$ and $u(x, t)$ are defined as before and $D(x)$ describes the water depth below its rest position. Many other systems have been derived also, including systems with improved dispersion characteristics [Nwo93], high-order Boussinesq systems [MBS03] and other generalizations of (2.1), cf. [Mit09]. Most of these systems break Galilean invariance and the invariance under vertical translations. This is a restrictive property especially in the studies of realistic problems like the water wave runup on non-uniform shores. We note also that the complete water wave problem possesses these symmetries, [BO82].
To overcome this deficiency we develop a new system, analogous to the original Peregrine's system, [Per67], which is invariant under vertical translations. To derive the system we begin with (2.10) written in dimensionless scaled variables (in analogy with (2.4))

$$ 
\begin{aligned}
\eta_t + [(D + \varepsilon\eta)u]_x &= 0, \\
u_t + \eta_x + \varepsilon u u_x - \sigma^2 \left[ \frac{D}{2}(Du)_{xxt} - \frac{D^2}{6}u_{xxt} \right] &= O(\varepsilon^2, \varepsilon\sigma^2).
\end{aligned} 
\quad (2.11) $$

Then by setting $H = D + \varepsilon\eta$, we obtain

$$ 
\begin{aligned}
H_t + \varepsilon(Hu)_x &= 0, \\
(Hu)_t + (\varepsilon Hu^2 + \frac{1}{2}\varepsilon H^2)_x \\
&\quad -\sigma^2 \left[ \frac{HD}{2}(Du)_{xxt} - \frac{HD^2}{6}(\frac{Du}{D})_{xxt} \right] = \frac{1}{\varepsilon}HD_x + O(\varepsilon^2, \varepsilon\sigma^2).
\end{aligned} 
\quad (2.12) $$

Observing that $(\frac{Du}{D})_{xx} = [2\frac{D_x^2}{D^3} - \frac{D_{xx}}{D^2}](Du) - 2\frac{D_x}{D^2}(Du)_x + \frac{1}{D}(Du)_{xx}$ and that $H = D + O(\varepsilon)$ we have that

$$ 
\begin{aligned}
H_t + \varepsilon(Hu)_x &= 0, \\
(Hu)_t + (\varepsilon Hu^2 + \frac{1}{2}\varepsilon H^2)_x \\
&\quad -\sigma^2 \left[ \frac{D^2}{3}(Du)_{xxt} - \left(\frac{1}{3}D_x^2 - \frac{1}{6}DD_{xx}\right)Du_t + \frac{1}{3}DD_x(Du)_{xt} \right] \\
&= \frac{1}{\varepsilon}HD_x + O(\varepsilon^2, \varepsilon\sigma^2).
\end{aligned} 
\quad (2.13) $$

By setting $Q = Hu$, and using again the relation $H = D + O(\varepsilon)$ we have

$$ 
\begin{aligned}
H_t + \varepsilon Q_x &= 0, \\
Q_t + (\varepsilon Q^2/H + \frac{1}{2}\varepsilon H^2)_x \\
&\quad -\sigma^2 \left[ \frac{H^2}{3}Q_{xxt} - \left(\frac{1}{3}H_x^2 - \frac{1}{6}HH_{xx}\right)Q_t + \frac{1}{3}HH_xQ_{xt} \right] = \frac{1}{\varepsilon}HD_x + O(\varepsilon^2, \varepsilon\sigma^2).
\end{aligned} 
\quad (2.14) $$

In dimensional variables, neglecting the higher order terms at the right-hand side we obtain

$$ 
\begin{aligned}
H_t + Q_x &= 0, \\
Q_t + (Q^2/H + \frac{g}{2}H^2)_x - P(H, Q) &= gHD_x. \\
P(H, Q) &= \frac{H^2}{3}Q_{xxt} - \left(\frac{1}{3}H_x^2 - \frac{1}{6}HH_{xx}\right)Q_t + \frac{1}{3}HH_xQ_{xt}
\end{aligned} 
\quad (2.15) $$

where $H(x, t) = \eta(x, t) + D(x)$, $Q(x, t) = H(x, t)u(x, t)$. We underline that system (2.15) is invariant under vertical translations and therefore more appropriate for studying long wave runup. Moreover, the linearization of the system (2.15) coincides with the original Peregrine's system (2.11) and therefore, inherits all its linear dispersive characteristics. On the other hand system (2.15) cannot be regarded as a correct asymptotic model to the Euler equations since it contains terms of the order $O(\varepsilon\sigma^2)$. However, such terms considered in the correct (small amplitude and long wave) regime are negligible, hence their contribution will be negligible. Finally we note that ignoring the dispersive terms $P(H, Q)$ of system (2.15) we obtain the shallow water equations (1.1).

We also note that even though Boussinesq systems are not valid in the near-shore region, in practice they appear to predict well the behavior of small amplitude waves from moderately deep to shallower waters and for smooth flows, cf. [Zel91]. Of course, more accurate systems in the near-shore zone have been derived such as the Sérre equations (sometimes referred also as Green-Naghdi equations), cf. [Ser53, LB09, DM10]. These
systems appeared in practice to model better the breaking phenomena in the near shore zone but recent numerical studies of the Sérre system showed that unphysical oscillations might appear in analogy with the Boussinesq equations during the wave breaking and the runup process, [CBB06, CBB07].

2.3. **Source terms.** Nonlinear shallow water model (1.1) and Boussinesq system (1.3) may be completed to take into account some dissipative or friction effects which are very beneficial in describing the wave breaking phenomena. Usually this is accomplished by including appropriate source or dissipative terms into momentum conservation equations (1.1) or (1.3). Possible choices are the following :

$$
\text{Friction:} \qquad F(u, H) = -c_m g \frac{u|u|}{H^{1/3}}, \qquad (2.16)
$$

$$
\text{Viscosity:} \qquad V(u, H) = \mu \frac{\partial}{\partial x} \left( H \frac{\partial u}{\partial x} \right), \qquad (2.17)
$$

where $c_m$ is the Manning roughness coefficient and $\mu$ denotes the kinematic viscosity of the fluid. The particular form of the source terms is suggested by empirical laws, which are generally obtained for steady state flows. Similar models have been derived from then Navier-Stokes system for incompressible flows with a free surface. More complex friction laws can be also formulated to model bottom rugosity effects, etc.

3. NUMERICAL SCHEMES

In the present section we generalize the finite volume method to systems (1.2) and (1.3)
of dispersive PDEs. In our work we rely on corresponding schemes for conservation laws.
Next we present briefly the finite volume framework for conservation laws. Based on this
framework we introduce finite volume schemes for the dispersive models.

3.1. Finite volume method for conservation laws. We consider the initial value problem

$$
\begin{align}
w_t + F(w)_x &= S(w), \quad x \in \mathbb{R}, t > 0 \tag{3.1} \\
w(x, 0) &= w_0(x), \nonumber
\end{align}
$$

where $w(x,t)$ is the state variable, $F$ denotes the flux and $S$ is the source term. Let $\mathcal{T} = \{x_i\}$, $i \in \mathbb{Z}$ denotes 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 $C_i$. Let $\Delta x_i = x_{i+\frac{1}{2}} - x_{i-\frac{1}{2}}$ be the length of the cell $C_i$, $\Delta x_{i+\frac{1}{2}} = x_{i+1} - x_i$. Without loss of generality we assume a uniform partition $\mathcal{T}$, that is $\Delta x_i = \Delta x_{i+\frac{1}{2}} = \Delta x$, $i \in \mathbb{Z}$. Let $w_i$ denotes the cell average of $w$ on $C_i$ i.e $w_i(t) = \frac{1}{\Delta x} \int_{C_i} w(x,t) dx$. Then a simple integration of (3.1) over a cell $C_i$ yields

$$
\frac{d}{dt} w_i(t) + \frac{1}{\Delta x} \left( F(w(x_{i+\frac{1}{2}}, t)) - F(w(x_{i-\frac{1}{2}}, t)) \right) = \frac{1}{\Delta x} \int_{C_i} S(w(x,t)) \, dx. \quad (3.2)
$$
### 3.1.1. Semidiscrete schemes.
We now define the semidiscrete finite volume approximation of $w(x,t)$. Let $\chi_{C_i}$ denotes the characteristic function of the cell $C_i$, we seek a piecewise constant function $w_h(x,t) = \sum_{i \in \mathbb{Z}} W_i(t)\chi_{C_i}(x)$ with

$$ \begin{aligned} \frac{d}{dt} W_i(t) + \frac{1}{\Delta x} (\mathcal{F}_{i+\frac{1}{2}} - \mathcal{F}_{i-\frac{1}{2}}) &= S_i, \quad i \in \mathbb{Z}, \\ W_i(0) &= \frac{1}{\Delta x} \int_{C_i} w(x, 0) dx, \quad i \in \mathbb{Z}, \end{aligned} \quad (3.3) $$

where $\mathcal{F}_{i+\frac{1}{2}} = \mathcal{F}(W_{i+\frac{1}{2}}^L, W_{i+\frac{1}{2}}^R)$ is an approximation to $F(w(x_{i+\frac{1}{2}}, t))$ while $S_i$ approximates the source term $S_i = S_i(W_i) \approx \frac{1}{\Delta x} \int_{C_i} S(w(x, t)) dx$. The values $W_{i+\frac{1}{2}}^L, W_{i+\frac{1}{2}}^R$ are approximations to the point value $w(x_{i+\frac{1}{2}}, t)$ from cells $C_i, C_{i+1}$ respectively and $\mathcal{F}$ is a numerical flux function which is consistent and monotone. The values $W_{i+\frac{1}{2}}^L, W_{i+\frac{1}{2}}^R$ are computed by a reconstruction process described below (see Section 3.1.3).

### 3.1.2. The numerical fluxes.
There are many possible choices for the numerical flux function $\mathcal{F}$. In the present study we choose to work with three following fluxes

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

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

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

The *average* flux (3.4) is the simplest one. It is well known that although this flux is unstable for nonlinear conservation laws, it is proven very stable and accurate for nonlinear dispersive models.

The *central* flux (3.5) is a Lax-Friedrichs type flux and is a representative of central schemes [KT00, NT90]. The operator $\mathcal{A}$ is related to the characteristic speeds of the flow and is defined as

$$ \mathcal{A}(W, V) = \max[\rho(DF(W)), \rho(DF(V))], \quad (3.7) $$

where $DF$ denotes the Jacobian matrix and $\rho(A)$ is the spectral radius of $A$.

The *characteristic* flux function (3.6), [GKC96, GKC01], is similar to the upwind flux and the operator $\mathcal{A}$ in this case is defined by

$$ \mathcal{A}(W, V) = \text{sign} \left( DF \left( \frac{W+V}{2} \right) \right). \quad (3.8) $$

### 3.1.3. The reconstruction process.
The values $W_{i+\frac{1}{2}}^L, W_{i+\frac{1}{2}}^R$ are approximations to $w(x_{i+\frac{1}{2}}, t)$ from cells $C_i$ and $C_{i+1}$ respectively. The simplest possible choice is to take the piecewise constant approximation in each cell,

$$ W_{i+\frac{1}{2}}^L = W_i, \quad W_{i+\frac{1}{2}}^R = W_{i+1}. \quad (3.9) $$

The resulting semidiscrete finite volume scheme is formally first order accurate in space.
To construct a higher order scheme in space, the piecewise constant data is replaced by a piecewise polynomial representation. The main idea here is to construct higher order approximations to $w(x_{i+1/2}, t)$ using the computed cell averages $W_i$. For this purpose the classical MUSCL type (TVD2) linear reconstruction [Kol75, vL79] as well as UNO2, [HO87] or WENO type reconstructions, [LOC94], have been developed.

The classical TVD2 linear reconstruction is given by the following formulas:

$$W_{i+1/2}^{L} = W_{i} + \frac{1}{2}\phi(r_{i})(W_{i+1} - W_{i}), \quad W_{i+1/2}^{R} = W_{i+1} - \frac{1}{2}\phi(r_{i+1})(W_{i+2} - W_{i+1}), \quad (3.10)$$

where $r_i = \frac{W_i - W_{i-1}}{W_{i+1} - W_i}$, and $\phi$ is an appropriate slope limiter, [Swe84]. There are many options for a limiter function. Some of the most 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 exhibit sharper resolution of discontinuities since they do not reduce the slope as severely as (MM) near a discontinuity. The TVD2 reconstruction is second order accurate except at the local extrema where it reduces to the first order. A remedy is to consider the UNO2 type reconstruction.

The UNO2 reconstruction is a linear interpolation which is second order accurate even at local extrema, [HO87]. The values $W_{i+1/2}^L, W_{i+1/2}^R$ are defined as

$$W_{i+1/2}^{L} = W_{i} + \frac{1}{2}S_{i}, \quad W_{i+1/2}^{R} = W_{i+1} - \frac{1}{2}S_{i+1}, \quad (3.11)$$

where

$$S_i = m(S_i^+, S_i^-), \quad S_i^{\pm} = d_{i\pm\frac{1}{2}}W \mp \frac{1}{2}D_{i\pm\frac{1}{2}}W, \\ d_{i+\frac{1}{2}}W = W_{i+1} - W_i, \quad D_{i+\frac{1}{2}}W = m(D_i W, D_{i+1} W), \\ D_i W = W_{i+1} - 2W_i + W_{i-1}, \quad m(x, y) = \frac{1}{2}(\mathrm{sign}(x) + \mathrm{sign}(y)) \min(|x|, |y|).$$

Using either (TVD2) or (UNO2) reconstructions the semidiscrete finite volume scheme (3.3) is formally second order accurate.

In order to achieve higher order accuracy we also employ WENO type reconstructions for the values $W_{i\pm\frac{1}{2}}^R, W_{i\pm\frac{1}{2}}^L$. We implemented 3rd and 5th order accurate WENO methods (also referred to as WENO3 and WENO5, respectively) as they are described in [LOC94]. For the sake of simplicity we only present the WENO3 case. In order to compute the approximations $W_{i+1/2}^L$ and $W_{i-1/2}^R$, we first compute the 3rd order reconstructed values

$$W_{i+1/2}^{(0)} = \frac{1}{2}(W_i + W_{i+1}), \quad W_{i+1/2}^{(1)} = \frac{1}{2}(-W_{i-1} + 3W_i), \\ W_{i-1/2}^{(0)} = \frac{1}{2}(3W_i - W_{i+1}), \quad W_{i-1/2}^{(1)} = \frac{1}{2}(W_{i-1} + W_i).$$
We define the smoothness parameters

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

and the parameters $d_0 = \frac{2}{3}$, $d_1 = \frac{1}{3}$ and $\tilde{d}_0 = d_1$, $\tilde{d}_1 = d_0$, along with the weights

$$ \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+\tilde{\beta}_i}$ and $\epsilon$ to be a small, positive number (in our computations we set $\epsilon = 10^{-15}$). Then the reconstructed values are given by the following formulas

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

### 3.1.4. Discretization of source terms.
The finite volume discretization of the source term $S(w)$ in (3.1) depends on the particular choice. On the other hand the resulting approximation should preserve the upwind nature and the overall scheme should be well balanced. One possible discretization of the source term $S(w)$ is given by:

$$ \frac{1}{\Delta x} \int_{C_i} S(w) dx \approx \frac{S_{i-\frac{1}{2}} + S_{i+\frac{1}{2}}}{2}, \quad S_{i+\frac{1}{2}} = S\left(\frac{W_{i+\frac{1}{2}}^{L} + W_{i+\frac{1}{2}}^{R}}{2}\right). \qquad (3.13) $$

### 3.1.5. Fully discrete schemes.
Equation (3.3) is an initial value problem and can be discretized by various methods. In our case we use a special class of Runge-Kutta methods which ensure the TVD property of the finite volume scheme, [SO88, GST01, SR02].

Let $\Delta t$ be the time step and let $t^{n+1} = t^n + \Delta t$, $n \ge 0$ be discrete time levels. Assuming at $t^n$ the approximations $\{W_i^n\}$, $i \in \mathbb{Z}$ are known then $W_i^{n+1}$ are defined by

$$ W_i^{n+1} = W_i^n - \frac{\Delta t}{\Delta x} \sum_{j=1}^{s} b_j (\mathcal{F}_{i+\frac{1}{2}}^{n,j} - \mathcal{F}_{i-\frac{1}{2}}^{n,j}) + \Delta t \sum_{j=1}^{s} b_j S_i^{n,j}, \qquad (3.14) $$

$$ W_i^{n,j} = W_i^n - \frac{\Delta t}{\Delta x} \sum_{\ell=1}^{s} a_{j\ell} (\mathcal{F}_{i+\frac{1}{2}}^{n,\ell} - \mathcal{F}_{i-\frac{1}{2}}^{n,\ell}) + \Delta t \sum_{\ell=1}^{s} a_{j\ell} S_i^{n,\ell}, \qquad (3.14) $$

where $\mathcal{F}_{i+\frac{1}{2}}^{n,j} = \mathcal{F}(W_i^{n,j}, W_{i+1}^{n,j})$, $S_i^{n,j} = S(W_i^{n,j})$. The set of constants $A = (a_{j\ell})$, $b = (b_1, \dots, b_s)$ define an 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|cc|cc|cc} & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline & 1 & 0 & 1 & \frac{1}{4} & \frac{1}{4} & \frac{1}{2} \\ \hline & \frac{1}{2} & \frac{1}{2} & & \frac{1}{6} & \frac{1}{6} & \frac{2}{3} \end{array} \qquad (3.15) $$

In our computations we mainly use the three stage 3rd order method.
3.2. Finite volume schemes for dispersive models. To construct the finite volume schemes for the dispersive models the main idea is to rewrite the governing equations or systems in a conservative like form and discretize the resulting conservation laws using the aforementioned framework. One can use any of the numerical fluxes, $\mathcal{F}^m$, $\mathcal{F}^{KT}$, $\mathcal{F}^{CF}$ and reconstruction techniques TVD2, UNO2 or WENO. Temporal discretization is based on the TVD-Runge-Kutta methods, (3.15).

3.2.1. Boussinesq systems with flat bottom. Boussinesq systems (1.2) can be rewritten in a conservative like form as follows:

$$ (\mathbf{I} - \mathbf{D})\mathbf{v}_t + [\mathbf{F}(\mathbf{v})]_x + [\mathbf{G}(\mathbf{v})]_x = 0, \quad (3.16) $$

where $\mathbf{v} = (\eta, u)^T$, $\mathbf{F}(\mathbf{v}) = ((1 + \eta)u, \eta + \frac{1}{2}u^2)^T$, $\mathbf{G}(\mathbf{v}) = (a u_{xx}, c \eta_{xx})^T$, and $\mathbf{D} = \text{diag}(b \partial_x^2, d \partial_x^2)$. The simplest discretization is based on the average fluxes $\mathcal{F}^m$ for $\mathbf{F}$ and $\mathcal{G}^m$ for $\mathbf{G}$. For the other two choices of the numerical flux $\mathcal{F}$ the evaluation of a Jacobian is needed. Let $A$ denotes the Jacobian of $\mathbf{F}$, then

$$ A = \begin{pmatrix} u & 1 + \eta \\ 1 & u \end{pmatrix}, $$

with eigenvalues $\lambda_i = u \pm \sqrt{1 + \eta}$, $i = 1, 2$. It is readily seen, since $\mathbf{F}$ is a hyperbolic flux, that $A$ can be decomposed as $A = L\Lambda R$ thus for the characteristic flux $\mathcal{F}^{CF}$ we have with $\mu = \frac{W+V}{2}$, $s_i = \text{sign}(\lambda_i)$, $i = 1, 2$

$$ \mathcal{A}(W, V) = \begin{pmatrix} \frac{1}{2}(s_1 + s_2) & \frac{1}{2}\sqrt{1 + \mu_1}(s_1 - s_2) \\ \frac{s_1 - s_2}{2\sqrt{1 + \mu_1}} & \frac{1}{2}(s_1 + s_2) \end{pmatrix}. $$

For evaluating the numerical fluxes $\mathcal{F}$, $\mathcal{G}$ simple cell averages or higher order approximations such as UNO2 (3.11) or WENO (3.12) can be used.

**Remark 1.** The discretization of the elliptic operator $\mathbf{D}$ is based on the standard centered difference. This is a second order accurate approximation and it is compatible with the TVD2 and UNO2 reconstructions. For higher order interpolation we modify the elliptic and flux discretization. Indeed, the finite volume scheme is modified as

$$ \frac{d}{dt} \left[ \frac{\mathbf{V}_{i-1} + 10\mathbf{V}_i + \mathbf{V}_{i+1}}{12} - (b, d) \frac{\mathbf{V}_{i+1} - 2\mathbf{V}_i + \mathbf{V}_{i-1}}{\Delta x^2} \right] + \frac{\mathcal{H}_{i-1} + 10\mathcal{H}_i + \mathcal{H}_{i+1}}{12} = 0, $$

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}})$, resulting in a high order accurate 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 the sequel for the discretization of the dispersive term $\mathbf{G}$ we use mainly the average numerical flux $\mathcal{G}^m$ defined as $\mathcal{G}_{i+\frac{1}{2}}^m = (a, c)\frac{1}{2}(\mathbf{Y}_i + \mathbf{Y}_{i+1})$, where $\mathbf{Y}_i = \frac{1}{\Delta x^2}(\mathbf{V}_{i+1} - 2\mathbf{V}_i + \mathbf{V}_{i-1})$. In case of higher order WENO reconstructions we use the average numerical
flux based on the reconstructed values of $\mathbf{Y}_i$ i.e. the flux $G_{i+\frac{1}{2}}^{lm} = (a, c)\frac{1}{2}(\mathbf{Y}_{i+\frac{1}{2}}^L + \mathbf{Y}_{i+\frac{1}{2}}^R)$, where $\mathbf{Y}_{i+\frac{1}{2}}^L$ and $\mathbf{Y}_{i+\frac{1}{2}}^R$ are reconstructed values of $\mathbf{Y}_i$.

3.2.2. Boussinesq system with variable bottom. We write system (2.15) in terms of dependent variables $\mathbf{v} := (H, Q)^T$ in the following conservative form

$$[\mathbf{D}(\mathbf{v}_t)] + [\mathbf{F}(\mathbf{v})]_x = \mathbf{S}(\mathbf{v}), \qquad (3.17)$$

where

$$\mathbf{D}(\mathbf{v}_t) = \left( \left(1 + \frac{1}{3}H_x^2 - \frac{1}{6}HH_{xx}\right)Q_t - \frac{1}{3}HH_xQ_{xt} - \frac{H^2}{3}Q_{xxt} \right), \qquad (3.18)$$

$$\mathbf{F}(\mathbf{v}) = \begin{pmatrix} Q \\ \frac{Q^2}{H} + \frac{g}{2}H^2 \end{pmatrix}, \qquad \mathbf{S}(\mathbf{v}) = \begin{pmatrix} 0 \\ gHD_x \end{pmatrix}. \qquad (3.19)$$

We consider a uniform mesh and we denote by $H_i$, $U_i$ and $D_i$ the corresponding cell averages. To discretize the dispersive terms in (3.18) we consider the following approximations:

$$\frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \left[ 1 + \frac{1}{3}(H_x)^2 - \frac{1}{6}HH_{xx} \right] Q \, dx \approx \left( 1 + \frac{1}{3} \left( \frac{H_{i+1} - H_{i-1}}{2\Delta x} \right)^2 - \frac{1}{6}H_i \frac{H_{i+1} - 2H_i + H_{i-1}}{\Delta x^2} \right) Q_i, \\ \frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \frac{1}{3}HH_x Q_x \, dx \approx \frac{1}{3}H_i \frac{H_{i+1} - H_{i-1}}{2\Delta x} \frac{Q_{i+1} - Q_{i-1}}{2\Delta x}, \qquad (3.20)$$

$$\frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \frac{1}{3} H^2 Q_{xx} \, dx \approx \frac{1}{3} H_i^2 \frac{Q_{i+1} - 2Q_i + Q_{i-1}}{\Delta x^2}. \qquad (3.21)$$

The aforementioned discretizations lead to a linear system with tridiagonal matrix denoted by $\mathbf{L}$ that can be inverted efficiently by a variation of Gauss elimination for tridiagonal systems with computational complexity $O(d)$, d-being the dimension of the system. We note that on the dry cells the matrix becomes diagonal since $H_i$ is zero on dry cells. For the time integration the explicit third-order TVD-RK method, (3.15) is used. In the numerical experiments we observed that the fully discrete scheme is stable and preserves the positivity of $H$ during the runup under a mild restriction on the time step $\Delta t$.

Therefore, the semidiscrete problem of (3.18) - (3.19) is written as a system of o.d.e's in the form

$$\mathbf{L}_i \mathbf{v}_{it} + \frac{1}{\Delta x} (\mathcal{F}_{i+\frac{1}{2}} - \mathcal{F}_{i-\frac{1}{2}}) = \frac{1}{\Delta x} \mathbf{S}_i, \qquad (3.22)$$

where $\mathbf{L}_i$ is the i-th row of matrix $\mathbf{L}$ and $\mathcal{F}_{i+\frac{1}{2}}$ can be chosen as one of the numerical flux functions mentioned in the previous sections. In the sequel we will use the KT and the CF
numerical fluxes. In this case the Jacobian of $\mathbf{F}$ is given by the matrix

$$A = \begin{pmatrix} 0 & 1 \\ gH - (Q/H)^2 & 2Q/H \end{pmatrix},$$

and the eigenvalues are $\lambda_{1,2} = Q/H \pm \sqrt{gH}$. Therefore, the CF numerical flux takes the form

$$\mathcal{F}_{i+\frac{1}{2}} = \frac{\mathbf{F}(\mathbf{V}_{i+\frac{1}{2}}^L) + \mathbf{F}(\mathbf{V}_{i+\frac{1}{2}}^R)}{2} - \mathbf{U}(\boldsymbol{\mu}) \frac{\mathbf{F}(\mathbf{V}_{i+\frac{1}{2}}^R) - \mathbf{F}(\mathbf{V}_{i+\frac{1}{2}}^L)}{2}, \quad (3.23)$$

where $\boldsymbol{\mu} = (\mu_1, \mu_2)^T$ are the Roe average values,

$$\mu_1 = \frac{H_{i+\frac{1}{2}}^L + H_{i+\frac{1}{2}}^R}{2}, \quad \mu_2 = \frac{\sqrt{H_{i+\frac{1}{2}}^L U_{i+\frac{1}{2}}^L + \sqrt{H_{i+\frac{1}{2}}^R U_{i+\frac{1}{2}}^R}}{\sqrt{H_{i+\frac{1}{2}}^L} + \sqrt{H_{i+\frac{1}{2}}^R}}$$

and

$$\mathbf{U}(\boldsymbol{\mu}) = \begin{pmatrix} \frac{s_2(\mu_2+c)-s_1(\mu_2-c)}{2c} & \frac{s_1-s_2}{2c} \\ \frac{(s_2-s_1)(\mu_2^2-c^2)}{2c} & \frac{s_1(\mu_2+c)-s_2(\mu_2-c)}{2c} \end{pmatrix}, \quad c = \sqrt{g\mu_1}, \quad s_i = \mathrm{sign}(\lambda_i). \quad (3.24)$$

In order to guarantee the positivity of the reconstructed values $H_{i+\frac{1}{2}}$ on the cell interfaces we employ the well balanced hydrostatic reconstruction algorithm, [ABB$^+$04]. Here we briefly recall the great lines of this reconstruction algorithm.

In the cell $C_i$ we compute first the reconstructions $\mathbf{V}_{i,r}$ and $\mathbf{V}_{i,l}$ at $(i+\frac{1}{2})^-$ and $(i-\frac{1}{2})^+$, respectively using either TVD2 or UNO2 with MinMod limiter. Moreover, we compute in the same way the values $\eta_{i,l}$ and $\eta_{i,r}$ of the free surface elevation $\eta_i = H_i - D_i$. Now we can deduce the values $D_{i,l} = H_{i,l} - \eta_{i,l}$ and $D_{i,r} = H_{i,r} - \eta_{i,r}$. Letting $D_{i+\frac{1}{2}} = \min(D_{i,r}, D_{i,l})$ we compute

$$H_{i+\frac{1}{2}}^R = \max(0, H_{i,r} + D_{i,r} - D_{i+\frac{1}{2}}), \quad H_{i+\frac{1}{2}}^L = \max(0, H_{i+1,l} + D_{i+1,l} - D_{i+\frac{1}{2}}), \quad (3.25)$$

and we deduce conservative reconstructed variables

$$\mathbf{V}_{i+\frac{1}{2}}^L = \begin{pmatrix} H_{i+\frac{1}{2}}^L \\ H_{i+\frac{1}{2}}^L u_{i,r} \end{pmatrix}, \quad \mathbf{V}_{i+\frac{1}{2}}^R = \begin{pmatrix} H_{i+\frac{1}{2}}^R \\ H_{i+\frac{1}{2}}^R u_{i+1,l} \end{pmatrix}. \quad (3.26)$$

Then the term $\mathbf{S}_i$ can be written as $\mathbf{S}_i = \mathbf{S}_{i+\frac{1}{2}}^L + \mathbf{S}_{i+\frac{1}{2}}^R + \mathbf{S}_{ci}$, where

$$\mathbf{S}_{i+\frac{1}{2}}^L = \begin{pmatrix} 0 \\ \frac{g}{2} \left[ (H_{i+\frac{1}{2}}^L)^2 - (H_{i,r})^2 \right] \end{pmatrix}, \quad \mathbf{S}_{i+\frac{1}{2}}^R = \begin{pmatrix} 0 \\ \frac{g}{2} \left[ (H_{i,l})^2 - (H_{i+\frac{1}{2}}^R)^2 \right] \end{pmatrix}$$

and

$$\mathbf{S}_{ci} = \begin{pmatrix} 0 \\ g \frac{H_{i,l}+H_{i,r}}{2} (D_{i,l} - D_{i,r}) \end{pmatrix}.$$

Numerical experiments show that the resulting scheme is well-balanced even for Boussinesq system of equations.
3.2.3. Boundary conditions. In the case of Bona-Smith type systems with flat bottom we consider herein only the initial-periodic boundary value problem which is known to be well-posed [ADM09].

In case of the modified Peregrine's system with an uneven bottom we use reflective boundary conditions. We note that for the classical Boussinesq system posed in a bounded domain $I = [b_1, b_2]$, one needs to impose boundary conditions only in one of the two dependent variables, cf. [FP05]. In the case of reflective boundary conditions it is sufficient to take $u(b_1, t) = u(b_2, t) = 0$ cf. [AD]. In [AD] it was also observed that during solitary waves reflection the derivatives $\eta_x(b_1, t) = \eta_x(b_2, t) \to 0$, while for other wave types these derivatives remained very small.

In our case we consider analogous reflective boundary conditions taking the cell averages of $u$ on the first and the last cell to be $u_0 = u_{N+1} = 0$. We don't impose explicitly boundary conditions on $H$. The reconstructed values on the first and the last cell are computed using neighboring ghost cells and taking odd and even extrapolation for $u$ and $H$ respectively. These specific boundary conditions appeared to reflect incident waves on the boundaries while conserving the mass.

4. INTERACTIONS OF SOLITARY WAVES

For the Boussinesq system (1.2) we present initially results demonstrating the accuracy of the finite volume scheme. We study the propagation as well as the interaction of solitary waves. In particular we consider head-on and overtaking collisions.

4.1. Accuracy test, validation. We consider the initial value problem with periodic boundary conditions for the Bona-Smith systems (2.7) with known solitary wave solutions (2.8) – (2.9) to study the accuracy of the finite volume method. We fix $\theta^2 = 8/10$ in the system and an analytic solitary wave of amplitude $\eta_0 = 1/2$ is used as the exact solution in $[-50, 50]$ computed up to $T = 200$. The error is measured with respect to discrete $L^2$ and $L^\infty$ scaled norms $E_h^2, E_h^\infty$, namely

$$ 
\begin{aligned}
E_h^2(k) &= \|u - U^k\|_h / \|U^0\|_h, \quad \|u - U^k\|_h = \left( \sum_i \Delta x |u(x_i) - U_i^k|^2 \right)^{1/2}, \\
E_h^\infty(k) &= \|u - U^k\|_{h,\infty} / \|U^0\|_{h,\infty}, \quad \|u - U^k\|_{h,\infty} = \max_i |u(x_i) - U_i^k|,
\end{aligned}
$$

where $U^k = \{U_i^k\}_i$ denotes the solution of the fully-discrete scheme at the time $t^k = k \Delta t$. The expected theoretical order of convergence was confirmed for all finite volume methods presented above. Three indicative cases, demonstrating the order of convergence, are reported in Table 1: a) for the average flux, b) for the KT flux with TVD2 reconstruction using the minmod limiter and c) for CF flux with WENO3 reconstruction. The order of convergence for the WENO5 method cannot be obtained since a 4th order discretization is used for the elliptic operator.
<table><caption>TABLE 1. Rates of convergence.</caption><thead><tr><th colspan="3">(a) Average Flux</th><th colspan="3">(b) KT-TVD2(MinMod)</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>1.910</td><td>1.978</td><td>0.5</td><td>2.042</td><td>2.032</td></tr><tr><td>0.25</td><td>1.910</td><td>1.954</td><td>0.25</td><td>2.033</td><td>2.029</td></tr><tr><td>0.125</td><td>1.923</td><td>1.937</td><td>0.125</td><td>2.026</td><td>2.023</td></tr><tr><td>0.0625</td><td>1.936</td><td>1.941</td><td>0.0625</td><td>2.021</td><td>2.019</td></tr><tr><td>0.03125</td><td>1.946</td><td>1.948</td><td>0.03125</td><td>2.017</td><td>2.016</td></tr></tbody></table><table><thead><tr><th colspan="3">(c) CF-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></tr></thead><tbody><tr><td>0.5</td><td>2.976</td><td>2.975</td></tr><tr><td>0.25</td><td>3.017</td><td>3.022</td></tr><tr><td>0.125</td><td>3.031</td><td>3.044</td></tr><tr><td>0.0625</td><td>3.042</td><td>3.059</td></tr><tr><td>0.03125</td><td>3.051</td><td>3.073</td></tr></tbody></table>

We also check the preservation of the invariant (2.3) by computing its discrete counterpart:

$$I_1^h = \sum_i \Delta x \left( \eta_i^2 + [(1+\eta_i)u_i]^2 - c \left[ \frac{\eta_{i+1} - \eta_i}{\Delta x} \right]^2 - a \left[ \frac{u_{i+1} - u_i}{\Delta x} \right]^2 \right), \quad (4.1)$$

as well as the discrete mass $I_0^h = \Delta x \sum_i \eta_i$. Figure 1 shows the evolution of the amplitude and the invariant $I_1^h$ of the solitary wave up to $T = 200$. The comparison of various methods is performed. We observe that the UNO2 reconstruction is more accurate while KT and the CF schemes show comparable performance. (We note that the invariant $I_0^h = 1.932183566158$ conserved the digits shown for all numerical schemes. In this experiment we took $\Delta x = 0.1$, $\Delta t = \Delta x/2$.) Figure 1 (a) and (c) show the evolution of the amplitude of the analytical solitary wave of the Bona-Smith system ($\theta^2 = 8/10$) and of the solitary wave produced by the solution of the analogous ordinary differential equations system of the classical Boussinesq system respectively. In the case of the classical Boussinesq system we took $c_s = 1.2$ and we used the method described in [DM08].

The well balanced property of the finite volume schemes is also verified numerically. We consider a uniform shore, cf. Section 5, including a wet and dry region. The bottom is also modified by adding a small parabolic type hump located at $x = 40$. We tested the steady state preservation of all fluxes and possible reconstruction techniques. The results are similar in all cases. In Figure 2 we present the case of FC flux along with UNO2.

4.2. **Head-on collisions.** The head-on collision of two counter-propagating solitary waves is characterized by the change of the shape along with a small phase-shift of the waves
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Preservation of the solitary wave amplitude and conservation of the invariant $I_1^h$: $G^m$ flux with Minmod limiter</figcaption></figure>

as a consequence of the nonlinearity and dispersion. These effects have been studied extensively before by numerical means using high order numerical methods such as finite differences, [BC98], spectral and finite element methods [ADM10, DDLMM07, PD01] and experimentally in [CGH+06]. In Figure 3 we present the numerical solutions of the BBM-BBM system (2.6) and the Bona-Smith system (2.7) with $\theta^2 = 9/11$ (in dimensional and unscaled variables) along with the experimental data from [CGH+06]. The spatial variable is expressed in centimeters while the time in seconds. The solutions were obtained using the CF-scheme with UNO2 and WENO3 reconstruction using $\Delta x = 0.05$ cm and $\Delta t = 0.01$ s. For this experiment we constructed solitary waves for Boussinesq systems by solving the respective o.d.e's system in the spirit of [BDM07] such that they fit to experimentally
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Steady state preservation</figcaption></figure>

generated solitary waves before the collision. The speeds of the right and left-traveling
solitary waves are $c_{r,s} = 0.854$ m/s and $c_{l,s} = 0.752$ m/s respectively.

We observe that Boussinesq models converge to the same numerical solution with all nu-
merical schemes we tested. A very good agreement with the experimental data is observed.
The maximum height predicted by the numerical solution during the collision process is
slightly higher in the case of the BBM-BBM system but the difference is negligible within
the specific experimental scale. Furthermore, we observe similar underestimation of the
maximum amplitude of the colliding waves compared to the experimental data, [CGH+06].
This discrepancy might be explained by a possible "splash" phenomenon during the colli-
sion reported also earlier by T. Maxworthy, [Max76]. After the collision we observe that
the phase shift of the solitary waves is the same in both numerical and experimental data,
while the shape of the experimental solitary waves were not stabilized due to interactions
with other small amplitude dispersive waves. We note that after the head-on collision of
the waves small amplitude dispersive tails were developed, [BC98, ADM10, DDLMM07].

The discrete mass for the Bona-Smith system is $I_0^h = 0.0059904310418$ and for the
BBM-BMM system is $I_0^h = 0.0059199389479$ for all fluxes and reconstructions used. The
variances in $I_1^h$ are mainly due to different types of reconstruction and not to the choice of
numerical fluxes. In Table 2 these values are reported.

4.3. **Overtaking collisions.** The overtaking collision of two solitary waves similarly to the head-on collision incorporates nonlinear and dispersive effects. Overtaking collision has been studied recently in the case of bidirectional models in [ADM10]. The interaction is similar to that of the unidirectional models but it was found that a new N-shape wavelet is generated during the interaction. This wavelet is of small amplitude and travels in the opposite direction to solitary waves and its shape depends on the Boussinesq system in use. Furthermore, as it was observed numerically and experimentally in [CGH+06],
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Head-on collision of two solitary waves: —: BBM-BBM, --: Bona-Smith ($\theta^2$ = 9/11), •: experimental data of [CGH$^+$06]</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 3. (Cont'd) Head-on collision of two solitary waves. —: BBM-BBM, --: Bona-Smith ($\theta^2 = 9/11$), •: experimental data of [CGH$^+$06]</figcaption></figure>

<table><caption>TABLE 2. Preservation of the invariant $I_1^h$.</caption><thead><tr><th></th><th>(a) Bona-Smith<br>$I_1^h$</th><th>(b) BBM-BBM<br>$I_1^h$</th></tr></thead><tbody><tr><td>m-flux</td><td>0.000944236</td><td>0.00092793</td></tr><tr><td>UNO2</td><td>0.00094423</td><td>0.00092793</td></tr><tr><td>TVD2</td><td>0.00094</td><td>0.00092</td></tr><tr><td>WENO3</td><td>0.00094423</td><td>0.00092793</td></tr></tbody></table>

the interaction of two solitary waves during an overtaking collision is characterized by a mass exchange and not by a simple superposition of the solitary pulses. These pulses remain separate retaining two different maxima contrary to unidirectional models where they merge into a single pulse momentarily.

To study this interaction we solve numerically the Bona-Smith system (2.7) with $\theta^2 = 9/11$. Following the same process as before two solitary waves were generated numerically with speeds $c_{1,s} = 1.2$ and $c_{2,s} = 1.4$. We solved the system using all fluxes using UNO2 and WENO3 reconstructions with discretization parameters $\Delta x = 0.01$, $\Delta t = 0.005$ up to $T = 600$. During simulations we were able to observe the generation and propagation of a small N-shape wavelet. In all computations the invariants were $I_0^h = 4.6098804880$, $I_1^h = 5.116$ conserving the digits shown for all methods.

In Figure 4 we present the interaction of two solitary waves. Figure 5 shows a magnification on the generation of a small wavelet along with the generation of dispersive tails as an effect of the inelastic interaction of two waves. In Figure 6 we observe that the overtaking collision is accompanied by an exchange of mass between pulses while both peaks are permanently present. The situation is different for unidirectional models where two pulses
<figure><img src="image_7.png" /><figcaption>FIGURE 4. Overtaking collision of two solitary waves of the Bona-Smith system with $\theta^2 = 9/11$.</figcaption></figure>

merge during a few time-steps to travel as a single pulse. Up to the graphic resolution we could not observe any difference in numerical solutions between UNO2 and WENO3 reconstructions.

4.4. **Small dispersion effect.** In this section we study the small dispersion effects on solitary waves of the classical Boussinesq system. The motivation for this study is the lack of theory supporting the breaking phenomena in Boussinesq systems contrary to the KdV equation. For this reason we employ the Boussinesq system with $a = b = c = 0$, $d = 10^{-5}$ and we take the solitary wave of the Boussinesq systems (2.5) as an initial condition. In Figure 7 we present numerical results obtained with CF-UNO2 and CF-WENO3 schemes. In these experiments we take $\Delta x = 0.001$ and $\Delta t = \Delta x/2$. The invariant $I_0^h$ is 1.629096452537 preserving the digits shown during all simulations. The invariant $I_1^h$ is not preserved by this model since the coefficient $b$ is not equal to $d$. The oscillations generated in the case of the WENO3 reconstruction were larger compared to those generated by the UNO2 reconstruction. Moreover, a new W-shaped wavelet is
<figure><img src="image_8.png" /><figcaption>FIGURE 5. Generation of a wavelet during the overtaking collision of two solitary waves of the Bona-Smith system with $\theta^2 = 9/11$.</figcaption></figure>

generated traveling to the left. This small wavelet finishes by producing a secondary breaking very similar to that of the initial solitary wave.

## 5. BOUSSINESQ SYSTEM WITH VARIABLE BOTTOM: RUNUP OF LONG WAVES

The shallow water equations are routinely used to predict a tsunami wave runup and, subsequently, constitute inundation maps for tsunami hazard areas. One of the main questions we address in this study is whether the inclusion of dispersive effects is beneficial for the description of the wave/shore interaction. In this section we perform a comparison of numerical solutions to Boussinesq equations (2.15), shallow water equations (1.1) (solved by the same numerical method) and experimental measurements made by C.E. Synolakis [Syn87] and J.A. Zelt [Zel91]. In these experiments we consider a bottom of the form, Figure 8,

$$ -D(x) = \begin{cases} -x \tan \beta, & x \le \cot \beta, \\ -1, & x > \cot \beta, \end{cases} $$
<figure><img src="image_9.png" /><figcaption>FIGURE 6. Overtaking collision of two solitary waves of the Bona-Smith system with $\theta^2 = 9/11$: mass exchange process.</figcaption></figure>

In all experiments over a flat bottom, $D(x) = D_0$, we use an approximate solitary wave solution of the following form:

$$ \eta_0(x) = A_s \operatorname{sech}^2(\lambda(x - X_0)), \quad u_0(x) = -c_s \frac{\eta_0(x)}{D_0 + \eta_0(x)}, \qquad (5.1) $$

$$ \lambda = \sqrt{\frac{3A_s}{4(1+A_s)}}, \quad c_s = \sqrt{g} \frac{\sqrt{6(1+A_s)}}{\sqrt{3+2A_s}} \cdot \frac{\sqrt{(1+A_s)\log(1+A_s)-A_s}}{A_s}, \qquad (5.2) $$

where $A_s$ denotes the amplitude, $c_s$ is the correct speed of the solitary wave propagation for classical Boussinesq equations.

The first three experiments we tested are described in [Syn87] and deal with the runup of solitary waves on a shore with a mild slope of 1 : 19.85. The first is a non-breaking solitary wave with dimensionless and scaled amplitude $A_s/D_0 = 0.0185$, the second one is
<figure><img src="image_10.png" /><figcaption>FIGURE 7. The small dispersion effect onto classical Boussinesq equations solutions.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 8. Sketch of the problem setup.</figcaption></figure>
a nearly breaking solitary wave with $A_s/D_0 = 0.04$, while the third experimental setup is a breaking solitary wave with $A_s/D_0 = 0.28$.

System (2.15) has some advantages over other asymptotically equivalent models with variable bottom. Namely, it shows excellent stability properties even for nearly breaking waves on the shore. However, for the simulation of strong breaking events, it is beneficial to include friction or dissipative terms taking into account turbulence generation.

We also considered two experiments from [Zel91] concerning the runup of solitary waves on a shore with steep slope 1 : 2.74. These experiments shed some light on the differences between dispersive and non-dispersive models.

Finally we consider a non-uniform sloping shore that contains a small pond demonstrating the capability of the modified Peregrine’s system to handle simultaneously and correctly dispersive effects in two basins with different mean sea levels.

In the sequel $t$ denotes the dimensionless time scaled by the quantity $\sqrt{g/D_0}$. Furthermore, we denote by $R$ the height of the last dry cell at a specific time instance. In our computations a cell is considered as dry if the total water depth $H_i$ inside is less than $5 \cdot 10^{-14}$. The quantity $R$ will also be referred to as *runup*. The maximum runup will be denoted by $R_\infty$. In all experiments the discretization parameters were taken to be equal $\Delta x = 0.05$, $\Delta t = \Delta x/10$, unless otherwise mentioned. Further, we compute in all cases the discrete mass $I_0^h$ and show the preserved digits. We use the KT and CF schemes combined with the TVD2 and UNO2 reconstructions. The CF-scheme appeared to be less dissipative and we emphasize the results of this method.

### 5.1. Runup of a solitary wave on a gradual slope $\beta = 2.88^\circ$ with $A_s/D_0 = 0.0185$
We consider first the simplest case — the runup of a non-breaking solitary wave. In this experiment we take an initial solitary wave with the amplitude $A_s = 0.0185$, $D_0 = 1$ and $X_0 = 19.85$ in $I = [-10, 70]$ and a mildly sloping shore 1 : 19.85. This specific solitary wave does not break [Syn87] and the solution remains smooth during the runup and the rundown processes. In Figure 9 we show several profiles of numerical solutions to Boussinesq and shallow water equations along with the experimental data of [Syn87]. We observe that both models converge to the same solution. The runup as well as the rundown in this experiment is predicted very well. The runup value $R$ for both models is almost the same. The maximum runup is $R_\infty \approx 0.085$ for the Boussinesq system, while for NSWE is $R_\infty \approx 0.088$. The experimental value reported in [Syn87] is equal to $R_\infty \approx 0.078$. In Figure 10 the runup $R$ as a function of time is represented. The discrete mass is preserved $I_0^h = 60.3667671231$ conserving the digits shown for both models.

### 5.2. Runup of a solitary wave on a gradual slope $\beta = 2.88^\circ$ with $A_s/D_0 = 0.04$
We consider the same sloping shore as before. We study the runup of a solitary wave with amplitude $A_s = 0.04$, placed initially at $X_0 = 19.85$ in $I = [-10, 70]$. The solitary wave does not break during the runup phase. Breaking occurs during the rundown process as in experimental observations [Syn87]. Results of the numerical simulations are presented in Figure 11. In Figure 12 the evolution of the runup value is shown. The maximum runup
<figure><img src="image_12.png" /><figcaption>FIGURE 9. Solitary wave runup on a sloping shore: $A_s = 0.0185$ case.</figcaption></figure>

for the Boussinesq system is $R_\infty \approx 0.20$ and $R_\infty \approx 0.21$ for shallow water system. The experimental value reported in [Syn87] is $R_\infty \approx 0.156$.

In Figure 13 we perform a comparison with tide gauge data (free surface elevation measured in [Syn87]) collected at 32.1 m from the still shoreline position. We observe again a good agreement between the dispersive and nondispersive models. The discrete mass is preserved, $I_0^h = 60.5210181987$ conserving the digits shown.

5.3. **Runup of a solitary wave on a gradual slope ($\beta = 2.88^\circ$) with $A_s = 0.28$.**
Finally we present the stiffest case of a solitary wave with amplitude $A_s = 0.28$, placed initially at $X_0 = 19.85$ in $I = [-10, 60]$. This specific initial condition is characterized by the wave breaking phenomenon before even reaching the shoreline. Strictly speaking, in this case Boussinesq model is not valid unless a wave breaking mechanism is considered, cf. [Zel91]. In this case the approximate solitary wave given by formulas (5.1)-(5.2) with $A_s = 0.28$ it is outside the range of validity of the specific system. We proceed by constructing numerically a more accurate solitary wave following the *cleaning procedure* of
<figure><img src="image_13.png" /><figcaption>FIGURE 10. Runup value $R$ as a function of time: $A_s = 0.0185$ case.</figcaption></figure>

[BC98]. (Note that since the linearized system (2.15) coincides with the linearized classical Boussinesq system we conclude that there exist classical solitary wave solutions of at least small amplitude, cf. [DM08]). We constructed a solitary wave with $A_s \approx 0.28$ with $\Delta x = 0.1$ and $\Delta x = 0.05$ while $\Delta t = \Delta x/10$.

Specifically, we consider the interval $[-200, 200]$, and the initial condition $\eta_0(x) = 0.5 \exp(-(x-150)^2/25)$, $u_0 = 0$ (with reflective boundary conditions). The initial condition is resolved into two solitary waves traveling in opposite directions. We observed that the left-traveling solitary wave at $t = 88$ was separated enough from the rest of the solution. This solitary wave is isolated by keeping the numerical solution in the interval $[-200, -151]$ and setting it equal to zero to the rest of the interval. The solution is then translated to the right at $X_0 = 19.85$ and we let it propagate. We observe that the clean solitary wave propagates like the analytical solution of the Bona-Smith system with analogous behavior between the TVD2 and UNO2 reconstruction. In Figure 14 we present the results for the KT scheme since the results of the CF scheme are completely analogous and no difference can be observed. In Figure 14 (a) the rightward traveling waves were reflected by the right boundary.

In order to ensure the stability of the simulation and to study the runup, instead of smoothing, filtering or adding extra dissipative terms, we simply exclude the contribution of the term $Q_{xxt}$ in the vicinity of the shoreline (where $D_i < 0.3$). Wave transition between
<figure><img src="image_14.png" /><figcaption>(a) $t = 20$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) $t = 26$</figcaption></figure>

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

<figure><img src="image_17.png" /><figcaption>(d) $t = 38$</figcaption></figure>

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

<figure><img src="image_19.png" /><figcaption>(f) $t = 50$</figcaption></figure>

<caption>FIGURE 11. Solitary wave runup on a sloping shore: $A_s = 0.04$ case.</caption>
<figure><img src="image_20.png" /><figcaption>FIGURE 11. (Cont'd) Solitary wave runup on a sloping shore: $A_s = 0.04$ case.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIGURE 12. Runup value $R$ as a function of time: $A_s = 0.04$ case.</figcaption></figure>

these two regions appeared to be smooth as one may witness on Figure 15. After this slight modification, the algorithm became more robust for large amplitude breaking waves without creating any unphysical oscillations.

In this experiment friction appeared to play a significant role during the runup process, contrary to previous cases. The maximum runup computed without taking into account the friction of the bottom was far away from the experimentally measured values. For this reason, and only in this specific test case we included the empirical friction term (2.16) into the momentum conservation equation (2.15), with coefficient $c_m = 2 \cdot 10^{-4}$. The friction term is discretized according to (3.13). This discretization preserves the positivity
<figure><img src="image_22.png" /><figcaption>FIGURE 13. Free surface elevation at $x = 32.1$ m.</figcaption></figure>

of all numerical schemes we tested. Mass conservation in this experiment was perfect $I_0^h = 51.7504637472$ preserving the digits shown.

In Figure 15 we show the propagation of a breaking wave including its runup and run-
down. We observe a significant difference between shallow water system and the dispersive
model during the wave propagation. Discrepancies are present in the amplitude and in
the phase speed simultaneously. However the dispersive model solution approximates bet-
ter the measurements of J.A. Zelt [Zel91]. Nevertheless, we have to underline that the
runup and rundown are fairly well described by both models. The maximum runup value
according to the dispersive and nondispersive models is $R_\infty \approx 0.47$ which is in the range
of [0.42, 0.53] of the theoretical prediction of C.E. Synolakis, [Syn87]. There is no single
experimental value reported for the maximum runup in [Syn87] due to practical difficulties
in generating a solitary wave of such large amplitude. Finally, we mention that the specific
technique for handling the breaking wave leads to more accurate results for the rundown
process than one presented in [Zel91].

5.4. **Solitary wave runup on a steep slope (β = 20°).** Now we present two experiments pointing out some further differences in solutions to dispersive and nondispersive models. These experiments were performed by J.A. Zelt, [Zel91]. We consider two waves in $I = [-10, 30]$ with amplitudes $A_s = 0.12$ and $A_s = 0.2$ initially located at $X_0 = 8.85$ and $X_0 = 10.62$ respectively. These waves propagate onto a steep sloping shore 1 : 2.74. A very fine grid of $\Delta x = 0.01$, $\Delta t = \Delta x/100$ is used to guarantee the accuracy and stability of simulations.

As it was observed in [Zel91] both waves do not break during the runup but the second one generates a strong breaking event during the rundown. Friction does not play an important role in these experiments. Consequently, no friction term is included into the models.
<figure><img src="image_23.png" /><figcaption>FIGURE 14. Generation of a solitary wave with $A_s \approx 0.28$ applying cleaning</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>FIGURE 15. Solitary wave runup on a sloping shore: $A_s = 0.28$ case.</figcaption></figure>

Figure 16 shows the runup value $R$ as a function of time. We observe for both models that there is a phase lag compared to the experimental data. We believe that this discrepancy can be removed by changing the definition of the last dry cell. We also observe that shallow water equations over-predict the maximum runup and the minimum rundown while the Boussinesq model predicts correctly the extrema in both cases.

Figure 17 shows the rundown of the solitary wave of amplitude $A_s = 0.2$ during breaking. One may observe that the experimental data consist of two curves due to the difficulty in measuring the surface elevation of the breaking wave due to 3D effects which become important. On Figure 18 the free surface elevation at a gauge is presented. The gauge is located 8.85 meters away from the still shoreline position. The reflected wave appears in both cases to be highly dispersive thus, the Boussinesq model provides a much better approximation. In the case $A_s = 0.2$, the mass remained equal to $I_0^h = 29.770808175$, while in the case $A_s = 0.12$, $I_0^h = 29.4861671693$ conserving the digits shown.
<figure><img src="image_25.png" /><figcaption>FIGURE 15. (Cont'd) Solitary wave runup on a sloping shore: A_s = 0.28 case.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>FIGURE 16. Runup value $R$ as a function of time.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>FIGURE 17. Rundown of the wave with amplitude $A_s = 0.2$.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>FIGURE 18. The amplitude at the wave gauge A</figcaption></figure>

5.5. **Solitary wave runup on a gradual slope ($\beta = 2.88^{\circ}$) with a pond.** We repeat the experiment of Section 5.2 in $I = [-10, 50]$ with solitary wave amplitude equal to $A_s = 0.05$. However, we modify the bottom by adding a small pond over the shoreline described by the exponential function $0.1 e^{-(x+4)^2}$. The Boussinesq system preserves the correct dispersion characteristics for the waves reaching the pond. In Figure 19 we present the overall process. It is worth noting that after the pond was filled a breaking wave was reflected back. As the wave slides down, a small hydraulic jump appears. In the case of shallow water system this jump propagates as a shock wave due to hyperbolic character of equations. On the other hand, the Boussinesq system develops into an Airy type wave according to its dispersive characteristics. In Figure 20 we show the solution at two wave gauges located at $x = -3.4$ and $x = 8$ for both the dispersive and nondispersive models. The mass during the simulations is constantly equal to $I_0^h = 40.5198087147$.

## 6. CONCLUSIONS

In the present study we extend the finite volume framework, developed for hyperbolic conservations laws, to approximate solutions of dispersive wave equations. This type of equations arises naturally in many physical problems. In the water wave theory dispersive equations have been well known since the pioneering work of J. Boussinesq [Bou71] and Korteweg-de Vries [KdV95]. Currently, the so-called Boussinesq-type models become more and more popular as an operational model for coastal hydrodynamics and other fields of engineering.

We extend the finite volume framework to dispersive models. We tested several choices of numerical fluxes, various reconstruction methods ranging from classical MUSCL type to modern approaches such as WENO. Various choices of limiters have been also tested out. Advantages of specific methods are discussed and some suggestions are outlined.

For operational modeling of the wave runup we derived a new system which has some advantages over its classical counterpart. The new system together with proposed novel
<figure><img src="image_29.png" /><figcaption>(a) $t = 37.58$</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(b) $t = 48.31$</figcaption></figure>

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

<figure><img src="image_32.png" /><figcaption>(d) $t = 58.67$</figcaption></figure>

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

<figure><img src="image_34.png" /><figcaption>(f) $t = 72.47$</figcaption></figure>

<caption>FIGURE 19. Long wave runup on a shore with a pond.</caption>
<figure><img src="image_35.png" /><figcaption>FIGURE 19. (Cont'd): Long wave runup on a shore with a pond.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>FIGURE 20. Evolution of the free surface elevation at two wave gauges.</figcaption></figure>
discretization procedure are validated by extensive comparisons with experimental data of C.E. Synolakis [Syn87] and J.A. Zelt [Zel91].

We paid a special attention to the comparison of dispersive (Boussinesq) and nondispersive (shallow water) models. Nowadays shallow water equations have become the model of choice for operational tsunami modeling including the inundation zone estimation [TG97, SBT+08]. The question of dispersive effects importance arises recurrently in the tsunami wave modeling community [KML05, DT07]. Our results show that shallow water equations are sufficient to predict maximum runup values. However, the dispersive effects can be beneficial for more accurate description of long wave propagation, runup and rundown.

## ACKNOWLEDGMENT

D. Dutykh acknowledges the support from French Agence Nationale de la Recherche, project MathOcean (Grant ANR-08-BLAN-0301-01) and Ulysses Program of the French Ministry of Foreign Affairs under the project 23725ZA. The work of Th. Katsaounis was partially supported by European Union FP7 program Capacities (Regpot 2009-1), through ACMAC (http://acmac.tem.uoc.gr). The work of D. Mitsotakis was supported by Marie Curie Fellowship No. PIEF-GA-2008-219399 of the European Commission. We would like to thank also professors Diane Henderson and Costas Synolakis for providing us their experimental data and professors Jerry Bona and Vassilios Dougalis for very helpful discussions.

## REFERENCES

[ABB+04] E. Audusse, F. Bouchut, O.M. Bristeau, R. Klein, and B. Perthame. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. *SIAM J. of Sc. Comp.*, 25:2050–2065, 2004. 13

[AC99] K. Anastasiou and C. T. Chan. Solution of the 2d shallow water equations using the finite volume method on unstructured triangular meshes. *International Journal for Numerical Methods in Fluids*, 24:1225–1245, 1999. 1

[AD] D. C. Antonopoulos and V. D. Dougalis. Numerical solution of the ‘classical’ boussinesq system. *to appear*. 14

[ADM09] D. C. Antonopoulos, V. A. Dougalis, and D. E. Mitsotakis. Initial-boundary-value problems for the Bona-Smith family of Boussinesq systems. *Advances in Differential Equations*, 14:27–53, 2009. 14

[ADM10] D. C. Antonopoulos, V. A. Dougalis, and D. E. Mitsotakis. Numerical solution of Boussinesq systems of the Bona-Smith family. *Appl. Numer. Math.*, 30:314–336, 2010. 16, 17

[Ami84] C. J. Amick. Regularity and uniqueness of solutions to the Boussinesq system of equations. *J. Diff. Eqns.*, 54:231–247, 1984. 5

[AVSS09] 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. 2

[BB01] 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

[BC98] J. L. Bona and M. Chen. A Boussinesq system for two-way propagation of nonlinear dispersive waves. *Physica D*, 116:191–224, 1998. 5, 16, 17, 26
[BCS02] 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. 2, 4

[BCS04] J.L. Bona, M. Chen, and J.-C. Saut. Boussinesq equations and other systems for small-amplitude long waves in nonlinear dispersive media: II. the nonlinear theory. *Nonlinearity*, 17:925-952, 2004. 4

[BDM07] 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. 2, 16

[BO82] T.B. Benjamin and P. Olver. Hamiltonian structure, symmetries and conservation laws for water waves. *J. Fluid Mech*, 125:137-185, 1982. 5

[Bou71] J. Boussinesq. Théorie de l'intumescence liquide appelée onde solitaire ou de translation se propageant dans un canal rectangulaire. *C.R. Acad. Sci. Paris Sér. A-B*, 72:755-759, 1871. 34

[BS76] J. L. Bona and R. Smith. A model for the two-way propagation of water waves in a channel. *Math. Proc. Camb. Phil. Soc.*, 79:167-182, 1976. 2, 5

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

[CBB06] R. Cienfuegos, E. Barthelemy, and P. Bonneton. A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part I: Model development and analysis. *Int. J. Numer. Meth. Fluids*, 51:1217-1253, 2006. 7

[CBB07] R. Cienfuegos, E. Barthelemy, and P. Bonneton. A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: Boundary conditions and model validation. *Int. J. Numer. Meth. Fluids*, 53:1423-1455, 2007. 7

[CFVCP04] L. Cea, A. Ferreiro, M. E. Vazquez-Cendon, and J. Puertas. Experimental and numerical analysis of solitary waves generated by bed and boundary movement. *Int. J. Numer. Meth. Fluids*, 46:793-813, 2004. 1

[CGH+06] W. Craig, P. Guyenne, J. Hammack, D. Henderson, and C. Sulem. Solitary water wave interactions. *Phys. Fluids*, 18:57-106, 2006. 16, 17, 18, 19

[Che98] M. Chen. Exact traveling-wave solutions to bidirectional wave equations. *International Journal of Theoretical Physics*, 37:1547-1567, 1998. 5

[Che00] M. Chen. Solitary-wave and multi pulsed traveling-wave solutions of Boussinesq systems. *Applic. Analysis.*, 75:213-240, 2000. 5

[DD07] D. Dutykh and F. Dias. Water waves generated by a moving bottom. In Anjan Kundu, editor, *Tsunami and Nonlinear waves*. Springer Verlag (Geo Sc.), 2007. 1

[DDLMM07] V. Dougalis, A. Duran, M. A. Lopez-Marcos, and D. E. Mitsotakis. A numerical study of the stability of solitary waves of Bona-Smith family of Boussinesq systems. *J. Nonlinear Sci.*, 17:595-607, 2007. 16, 17

[DK03] A. I. Delis and Th. Katsaounis. Relaxation schemes for the shallow water equations. *Int. J. Numer. Meth. Fluids*, 41:695-719, 2003. 1

[DKK08] A. I. Delis, M. Kazolea, and N. A. Kampanis. A robust high-resolution finite volume scheme for the simulation of long waves over complex domains. *Int. J. Numer. Meth. Fluids*, 56:419-452, 2008. 1

[DM04] V.A. Dougalis and D.E. Mitsotakis. *Advances in scattering theory and biomedical engineering*, chapter Solitary waves of the Bona-Smith system, pages 286-294. World Scientific, New Jersey, 2004. 5

[DM08] V. A. Dougalis and D. E. Mitsotakis. Theory and numerical analysis of Boussinesq systems: A review. In N. A. Kampanis, V. A. Dougalis, and J. A. Ekaterinaris, editors, *Effective Computational Methods in Wave Propagation*, pages 63-110. CRC Press, 2008. 2, 15, 26
[DM10] F. Dias and P. Milewski. On the fully-nonlinear shallow-water generalized Serre equations. *Physics Letters A*, 374(8):1049-1053, 2010. 6

[DPD10] D. Dutykh, R. Poncet, and F. Dias. Complete numerical modelling of tsunami waves: generation, propagation and inundation. *Submitted*, http://arxiv.org/abs/1002.4553, 2010. 1

[DT07] M.H. Dao and P. Tkalich. Tsunami propagation modelling - a sensitivity study. *Nat. Hazards Earth Syst. Sci.*, 7:741-754, 2007. 37

[EIK05] 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

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

[ESB06] C. Eskilsson, S.J. Sherwin, and L. Bergdahl. An unstructured spectral/hp element model for enhanced Boussinesq-type equations. *Coastal Engineering*, 53(11):947-963, 2006. 2

[FP05] A. S. Fokas and B. Pelloni. Boundary value problems for Boussinesq type systems. *Math. Phys. Anal. Geom.*, 8:59-96, 2005. 14

[GKC96] 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, 8

[GKC01] 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. 8

[GST01] S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong stability-preserving high-order time discretization methods. *SIAM Review*, 43:89-112, 2001. 3, 10

[HH70] K.L. Heitner and G.W. Housner. Numerical model for tsunami runup. *J. Waterway, Port, Coastal and Ocean Engineering*, 96:701-719, 1970. 1

[HMK+09] M. Hamidou, B. Molin, M. Kadri, O. Kimmoun, and A. Tahakourt. A coupling method between extended Boussinesq equations and the integral equation method with application to a two-dimensional numerical wave-tank. *Ocean Engineering*, 36:1377-1385, 2009. 2

[HO87] A. Harten and S. Osher. Uniformly high-order accurate nonscillatory schemes, I. *SIAM J. Numer. Anal.*, 24:279-309, 1987. 3, 9

[HP79] S. Hibberd and D.H. Peregrine. Surf and run-up on a beach: a uniform bore. *J. Fluid Mech.*, 95:323-345, 1979. 1

[JFBM06] E. Jamois, D.R. Fuhrman, H.B. Bingham, and B. Molin. A numerical study of nonlinear wave run-up on a vertical plate. *Coastal Engineering*, 53:929-945, 2006. 2

[KCKD00] A.B. Kennedy, Q. Chen, J.T. Kirby, and R.A. Dalrymple. Boussinesq modelling of wave transformation, breaking, and runup. *J. Waterway, Port, Coastal and Ocean Engineering*, 126:39-47, 2000. 2

[KDD07] Y. Kervella, D. Dutykh, and F. Dias. Comparison between three-dimensional linear and nonlinear tsunami generation models. *Theor. Comput. Fluid Dyn.*, 21:245-269, 2007. 1

[KdV95] 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. 34

[KML05] E.A. Kulikov, P.P. Medvedev, and S.S. Lappo. Satellite recording of the Indian Ocean tsunami on December 26, 2004. *Doklady Earth Sciences A*, 401:444-448, 2005. 37

[Kol75] N.E. Kolgan. Finite-difference schemes for computation of three dimensional solutions of gas dynamics and calculation of a flow over a body under an angle of attack. *Uchenye Zapiski TsaGI [Sci. Notes Central Inst. Aerodyn]*, 6(2):1-6, 1975. 9

[KS98] U. Kanoglu and C.E. Synolakis. Long wave runup on piecewise linear topographies. *J. Fluid Mech.*, 374:1-28, 1998. 1

[KT00] 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, 8
[KWC$^{+98}$] J. T. Kirby, G. Wei, Q. Chen, A. B. Kennedy, and R. A. Dalrymple. Funwave 1.0, fully nonlinear boussinesq wave model documentation and user's manual. Research Report No. CACR-98-06, 1998. 2

[LB09] D. Lannes and P. Bonneton. Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation. *Phys. Fluids*, 21:016601, 2009. 6

[LOC94] X.-D. Liu, S. Osher, and T. Chan. Weighted essentially non-oscillatory schemes. *J. Comp. Phys.*, 115:200-212, 1994. 3, 9

[LR02] Y. Li and F. Raichlen. Non-breaking and breaking solitary wave run-up. *J. Fluid Mech.*, 456:295-318, 2002. 2

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

[Max76] T. Maxworthy. Experiments on collisions between solitary waves. *J Fluid Mech*, 76:177-185, 1976. 17

[MBS03] 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. 2, 5

[Mit09] 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. 2, 3, 5

[Ngu08] H.Y. Nguyen. Modèles pour les ondes interfaciales et leur intégration numérique. PhD thesis, Ecole Normale Supérieure de Cachan, 2008. 2

[NT90] H. Nessyahu and E. Tadmor. Nonoscillatory central differencing for hyperbolic conservation laws. *J. Computational Physics*, 87(2):408-463, 1990. 2, 8

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

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

[Per67] D. H. Peregrine. Long waves on a beach. *J. Fluid Mech.*, 27:815-827, 1967. 2, 5, 6

[SBT$^{+08}$] C.E. Synolakis, E.N. Bernard, V.V. Titov, U. Kânoglu, and F.I. González. Validation and verification of tsunami numerical models. *Pure Appl. Geophys.*, 165:2197-2228, 2008. 37

[Sch81] M. E. Schonbek. Existence of solutions for the Boussinesq system of equations. *J. Diff. Eqns.*, 42:325-352, 1981. 5

[Ser53] F. Serre. Contribution à l'étude des écoulements permanents et variables dans les canaux. *La Houille blanche*, 8:374-388 & 830-872, 1953. 2, 6

[Sim81] M.J. Simon. Wave-energy extraction by a submerged cylindrical resonant duct. *Journal of Fluid Mechanics*, 104:159-187, 1981. 1

[SM09] 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

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

[SR02] 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. 3, 10

[Swe84] P.K. Sweby. High resolution schemes using flux limiters for hyperbolic conservation laws. *SIAM J. Numer. Anal.*, 21(5):995-1011, 1984. 3, 9

[Syn87] C. Synolakis. The runup of solitary waves. *J. Fluid Mech.*, 185:523-545, 1987. 1, 21, 22, 24, 25, 29, 37

[TG97] V.V. Titov and F.I. González. Implementation and testing of the method of splitting tsunami (MOST) model. Technical Report ERL PMEL-112, Pacific Marine Environmental Laboratory, NOAA, 1997. 1, 37
[TP09] M. Tonelli and M. Petti. Hybrid finite-volume finite-difference scheme for 2DH improved Boussinesq equations. *Coastal Engineering*, 56:609-620, 2009. 2

[TS94] S. Tadepalli and C. E. Synolakis. The run-up of N-waves on sloping beaches. *Proc. R. Soc. Lond. A*, 445:99-112, 1994. 1

[TS98] V. V. Titov and C. E. Synolakis. Numerical modeling of tidal wave runup. *J. Waterway, Port, Coastal, and Ocean Engineering*, 124:157-171, 1998. 1

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

[WKGS95] 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. 2

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

[Zel91] J.A. Zelt. The run-up of nonbreaking and breaking solitary waves. *Coastal Engineering*, 15:205-246, 1991. 2, 6, 21, 24, 25, 29, 37

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/

UMR DE MATHÉMATIQUES, UNIVERSITÉ DE PARIS-SUD, BÂTIMENT 425, P.O. BOX, 91405 ORSAY, FRANCE

E-mail address: Dimitrios.Mitsotakis@math.u-psud.fr

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