# A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part I: Model development and analysis

Rodrigo Cienfuegos, Eric Barthélemy, Philippe Bonneton

► To cite this version:

Rodrigo Cienfuegos, Eric Barthélemy, Philippe Bonneton. A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part I: Model development and analysis. International Journal for Numerical Methods in Fluids, 2006, 51 (11), pp.1217-1253. 10.1002/fld.1141. hal-00182841

HAL Id: hal-00182841

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

Submitted on 23 Jan 2020

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
# A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations.
## Part I: Model development and analysis

R. Cienfuegos¹,³,*,†, E. Barthélemy¹,‡ and P. Bonneton²,§

¹Laboratoire des Écoulements Géophysiques et Industriels, Institut National Polytechnique de Grenoble,
BP53 38041 Grenoble Cedex 9, France

²Département de Géologie et d'Océanographie, Université Bordeaux I, Av. des Facultés 33405 Talence, France

³Departamento de Ingeniería Hidráulica y Ambiental, Pontificia Universidad Católica de Chile,
Vicuña Mackenna 4860, casilla 306, correo 221, Santiago, Chile

## SUMMARY

A high-order finite volume scheme is developed to numerically integrate a fully nonlinear and weakly dispersive set of Boussinesq-type equations (the so-called Serre equations) (J. Fluid Mech. 1987; 176:117-134; Surveys Geophys. 2004; 25(3-4):315-337). The choice of this discretization strategy is motivated by the fact that this particular set of equations is recasted in a convenient quasi-conservative form. Cell face values are reconstructed using implicit compact schemes (J. Comput. Phys. 1999; 156:137-180; J. Comput. Phys. 2004; 198:535-566) and time integration is performed with the help of a four-stage Runge-Kutta method. Numerical properties of the proposed scheme are investigated both, analytically using linear spectral analysis, and numerically for highly nonlinear cases. The numerical analysis indicates that the newly developed scheme has wider stability regions and better spectral resolution than most of the previously published numerical methods used to handle equivalent set of equations. Moreover, it was also noticed that the use of mixed-order strategies to discretize convective and dispersive terms may result in an important overall reduction of the spectral resolution of the scheme. Additionally, there is some numerical evidence, which seems to indicate that the incorporation of a high-order dispersion correction term as given by Madsen et al. (Coastal Eng. 1991; 15:371-388) may introduce instability in the system.

KEY WORDS: Boussinesq-type equations; Serre equations; finite volume method; compact schemes

<footnote>*Correspondence to: R. Cienfuegos, Departamento de Ingeniería Hidráulica y Ambiental, Pontificia Universidad Católica de Chile, Vicuña Mackenna 4860, casilla 306, correo 221, Santiago, Chile.</footnote>

<footnote>†E-mail: racienfu@ing.puc.cl</footnote>

<footnote>‡E-mail: eric.barthelemy@hmg.inpg.fr</footnote>

<footnote>§E-mail: p.bonneton@epoc.u-bordeax1.fr</footnote>

<footnote>Contract/grant sponsor: French Foreign Office</footnote>

<footnote>Contract/grant sponsor: Chilean Research Council (CONICYT)</footnote>
# 1. INTRODUCTION

In recent years, important research efforts have provided a solid theoretical background for a new generation of the so-called Boussinesq-type wave equations. This new class of equations, which takes advantage of pioneering contributions made by Boussinesq [1], Korteweg and de Vries [2], Serre [3] and Peregrine [4], have proven to adequately describe most of the water-wave phenomena taking place in the near-shore zone. Recent developments concern the extension of these shallow water sets of equations to deeper water propagation problems [5–10], and to describe surf-zone dynamics by adding extra wave-breaking terms [11–14] and removing the originally embedded weakly nonlinear hypothesis [15–17]. Even though those important improvements have provided a powerful set of mathematical tools to describe the hydrodynamics of coastal zones, the question of constructing efficient, stable and accurate numerical schemes to solve them have received considerably less attention.

Numerical solutions of Boussinesq-type equations have been mostly tackled with finite difference methods [18–24], and to some extent, using finite element techniques [25–28]. On the contrary, the application of the finite volume method to discretize them has been only reported very recently [29–32]. Even though finite volume techniques have been successfully applied to the strictly hyperbolic dispersionless set of nonlinear shallow water equations (NSWE) using shock-capturing solvers (see Reference [33] for a comprehensive review), parabolic terms responsible for the dispersive nature of Boussinesq-type equations make the extrapolation of such approaches to coastal applications rather difficult. However, shock-capturing methods may provide an efficient and elegant way to numerically integrate Boussinesq equations but novel strategies are to be developed in order to adequately handle the mixed hyperbolic–parabolic mathematical nature of the system. Several authors are working in that direction, developing in particular *hybrid* or *operator splitting* techniques where convective hyperbolic-type terms are treated using shock-capturing approaches while high-order dispersive terms are handled with help of standard finite difference approximations (e.g. References [31, 34]). Such strategies have only been reported in applications concerning weakly nonlinear forms of Boussinesq-type equations. On the other hand, achieving high-order accuracy in the latter is quite demanding from a computational point of view and there is still no clear knowledge of what type of flux reconstruction strategy would be the most appropriate in the framework of dispersive wave propagation problems (see, for instance, Reference [31]).

In the present work, a finite volume technique is applied to the fully nonlinear and weakly dispersive set of Serre equations first derived by Seabra-Santos et al. [15]. The inclusion of a dispersion correction term following the lines given by Madsen and Schäffer [35] allows for an extension of the range of application of this set into moderately deep waters. The choice of integrating equations with a finite volume technique appeared naturally because it is possible to recast continuity and momentum equations in an attractive weak conservative form using an auxiliary variable as given in Reference [36]. Here, *weak conservative* is used in order to make a distinction between a formulation written in terms of primitive variables, i.e. the water depth and the depth-averaged velocity, and a *strongly conservative* formulation, which would require to express the momentum conservation equation in terms of the momentum flux. We recall that recent applications of finite volume methods to Boussinesq-type equations mostly concern weakly dispersive and weakly nonlinear versions of these sets (e.g. a KdV-like system in References [29, 34] and a low-order Boussinesq model in References [31, 32]). Similarly, most of discretization strategies implemented in previous finite volume approaches appears
to be 2nd-order accurate in space. On the contrary, our aim here is to develop an efficient high-order finite volume scheme applicable even for strongly nonlinear situations. However, no attempt is made to take advantage of shock-capturing methods mainly because: (i) this convenient conservative set of equations is not really in a strong conservative mathematical form, (ii) high-order accuracy in space and time was intended at low computational cost, and (iii) shocks do not naturally emerge in Boussinesq-type equations because there is a balance between nonlinear and dispersive terms [37]. Nevertheless, our numerical approach borrows some recently developed ideas from modern computational fluid dynamics (CFD) which make high-order accuracy and high spectral resolution possible at an affordable computational cost. Such techniques are based on a compact reconstruction of cell face values at a high-order by solving implicit algebraical equations [38, 39]. Indeed, 4th-order accuracy is achieved by solving an implicit tridiagonal system. Similarly, spatial derivatives are estimated with a 4th-order truncation error using compact differencing [40]. Finally, numerical solution is advanced in time with help of the four-stage Runge-Kutta method.

This paper is organized as follows: in Section 2 governing equations are derived by vertically integrating Euler equations, while Section 3 deals with the numerical implementation of the newly developed compact finite volume scheme. A comprehensive linear numerical analysis of the method is given in Section 4 and additional information on the performance of the finite volume scheme is obtained with help of numerical experiments in Section 5. In particular, the influence that strong nonlinearities may have on stability and accuracy is investigated.

## 2. GOVERNING EQUATIONS

Serre equations describe the 2D, irrotational and shallow water flow of an incompressible and inviscid fluid over uneven bottom bathymetries. This set of equations rules the total water depth and the depth-averaged horizontal velocity. In the following we expand on the derivation of Serre equations for horizontal bathymetries as given by Barthélemy in Reference [41].

Generally speaking, Boussinesq-like equations are valid only for long waves. Long wave dynamics are characterized by a small parameter $\sigma = h_0/L \ll 1$ where $h_0$ and $L$ are, respectively, a vertical length scale (say a typical water depth) and a horizontal length scale of the flow. The continuity equation implies the following order of magnitude relationship:

$$ \frac{w}{u} \sim \sigma \qquad (1) $$

where $u$ and $w$ are, respectively, the horizontal and vertical fluid velocities ($\sim$ stands for order of magnitude). The bottom boundary condition suggests that bottom slopes $\xi_x$ should be of order $\sigma$, since $w/u = \xi_x$ at the bottom. Obviously, bottom horizontal scales are not disconnected from the scales of the waves. In other words, the flow length scales adapt to the bottom ones. Hence, bottom slopes for long wave propagation should be mild, which is what is assumed here. Even though this assumption is implicitly adopted in most of the available long wave propagation models, thus introducing practical limitations in the latter, Serre equations have been successfully used to describe solitary wave propagation over a submerged step in Reference [15]. This settled, the scaling detailed in Reference [41] leads
to the following dimensionless forms of continuity, irrotationality, momentum equations and boundary conditions:

$$ \frac{\partial u^*}{\partial x^*} + \frac{\partial w^*}{\partial z^*} = 0 \qquad (2) $$

$$ \frac{\partial u^*}{\partial z^*} - \sigma^2 \frac{\partial w^*}{\partial x^*} = 0 \qquad (3) $$

$$ \varepsilon \frac{\partial u^*}{\partial t^*} + \varepsilon^2 \frac{\partial}{\partial x^*} (u^{*2}) + \varepsilon^2 \frac{\partial}{\partial z^*} (u^* w^*) = -\frac{\partial p^*}{\partial x^*} \qquad (4) $$

$$ \varepsilon \sigma^2 \frac{\partial w^*}{\partial t^*} + \varepsilon^2 \sigma^2 u^* \frac{\partial w^*}{\partial x^*} + \varepsilon^2 \sigma^2 w^* \frac{\partial w^*}{\partial z^*} = -\frac{\partial p^*}{\partial z^*} - 1 \qquad (5) $$

$$ w^* = u^* \frac{\partial \xi^*}{\partial x^*} \quad \text{at the bottom } z^* = \xi^* \qquad (6) $$

$$ w^* = \frac{\partial \eta^*}{\partial t^*} + \varepsilon u^* \frac{\partial \eta^*}{\partial x^*} \quad \text{at the free surface } z^* = \varepsilon \eta^* \qquad (7) $$

$$ p^* = 0 \quad \text{at the free surface } z^* = \varepsilon \eta^* \qquad (8) $$

where * indicates dimensionless variables. Here, $p^*$ represents the pressure, $\eta^*$ is the free surface elevation and $\varepsilon = a/h_0$ is the nonlinear parameter with a being the local wave amplitude relative to the chosen reference level (see Figure 1). For simplicity, we introduce the following function:

$$ \Gamma^*(x^*, z^*, t^*) = \frac{\partial w^*}{\partial t^*} + \varepsilon \left( u^* \frac{\partial w^*}{\partial x^*} + w^* \frac{\partial w^*}{\partial z^*} \right) \qquad (9) $$

which is related to the vertical acceleration of fluid particles.

Recall, we seek a set of equations for the depth-averaged horizontal velocity $\overline{u}^*$ and the water depth $h^* = \varepsilon\eta^* - \xi^*$. Depth averaging a function $f^*$ results from applying the following

<figure><img src="image_2.png" /><figcaption>Figure 1. Definition sketch for wave propagation over uneven bathymetries. $z = \eta$ is the free surface elevation, $z = \xi$ is the bottom location, and $h = \eta - \xi$ is the local water depth (here $\xi$ has a negative value).</figcaption></figure>
operator:

$$ \overline{f^*} = \frac{1}{h^*} \int_{\xi^*}^{\epsilon\eta^*} f^* \, dz^* \qquad (10) $$

Depth-averaging continuity equation yields

$$ \frac{\partial \eta^*}{\partial t^*} + \frac{\partial}{\partial x^*} [h^* \overline{u^*}] = 0 \qquad (11) $$

Depth averaging the x-momentum equation (4) and following the derivation along similar lines as in Reference [41] gives

$$ \begin{aligned} & h^* \frac{\partial \overline{u^*}}{\partial t^*} + \epsilon h^* \overline{u^*} \frac{\partial \overline{u^*}}{\partial x^*} + h^* \frac{\partial \eta^*}{\partial x^*} + \sigma^2 \frac{\partial}{\partial x^*} \int_{\xi^*}^{\epsilon\eta^*} z^* \Gamma^* dz^* + \sigma^2 \frac{\partial \xi^*}{\partial x^*} \int_{\xi^*}^{\epsilon\eta^*} \Gamma^* dz^* \\ & = - \epsilon \frac{\partial}{\partial x^*} \int_{\xi^*}^{\epsilon\eta^*} (u^{*2} - (\overline{u^*})^2) dz^* \end{aligned} \qquad (12) $$

So far no approximation has been invoked, but the integrals in (12) cannot be reduced to partial differential expressions without approximating the integrands. We aim at obtaining a set of equations comprising terms at most of $O(\sigma^2)$. To be consistent, the two integrals on the left-hand side should be expanded to order 0. Further inspection is required for the one on the right-hand side, which also reads

$$ \int_{\xi^*}^{\epsilon\eta^*} (u^{*2} - (\overline{u^*})^2) dz^* = h^*(\overline{u^{*2}} - (\overline{u^*})^2) \qquad (13) $$

This term has given rise in the past to discussions as to whether it is of $O(\sigma^4)$ or not for uneven bathymetries. Dingemans [42, volume II, p. 613] argues it is not. It has been then taken advantage of this argument to mistrust Serre equations for long wave propagation over uneven topographies. Nevertheless, it can be readily shown that (13) is of order $\sigma^4$. Using series expansions for $u^*$ and $w^*$ given in Appendix A we have,

$$ u^*(x^*, z^*, t^*) = u_0^*(x^*, t^*) + \sigma^2(z^* - \xi^*)u_2^*(x^*, t^*) + O(\sigma^4) \qquad (14) $$

$$ w^*(x^*, z^*, t^*) = \left[ \frac{\partial \xi^*}{\partial x^*} u_0^*(x^*, t^*) - (z^* - \xi^*) \frac{\partial u_0^*}{\partial x^*} \right] + \sigma^2(z^* - \xi^*) w_2^*(x^*, t^*) + O(\sigma^4) \qquad (15) $$

showing in addition that $u^*$ and $w^*$ also have even series expansions in $\sigma$ for uneven bathymetries. Using (14), clearly $\overline{u^{*2}} - (\overline{u^*})^2 = O(\sigma^4)$ without having to express $u_0^*$ and $u_2^*$ in terms of $\overline{u^*}$ as done by Dingemans [42].

Left-hand side integrals in (12) are evaluated with the above series expansions expressed in terms of $\overline{u^*}$,

$$ u^*(x^*, z^*, t^*) = \overline{u^*}(x^*, t^*) + O(\sigma^2) \qquad (16) $$

$$ w^*(x^*, z^*, t^*) = \frac{\partial \xi^*}{\partial x^*} \overline{u^*} - (z^* - \xi^*) \frac{\partial \overline{u^*}}{\partial x^*} + O(\sigma^2) \qquad (17) $$
Retaining all terms up to $O(\sigma^2)$ leads to the following set of equations expressed in terms of dimensional variables and unknowns:

$$ \eta_t + (h\bar{u})_x = 0 \qquad (18) $$

$$ h\bar{u}_t + h\bar{u}\bar{u}_x + gh\eta_x + [h^2(\frac{1}{3}\mathcal{P} + \frac{1}{2}\mathcal{Q})]_x + \xi_x h(\frac{1}{2}\mathcal{P} + \mathcal{Q}) = 0 \qquad (19) $$

$$ \mathcal{P} = -h[\bar{u}_{xt} + \bar{u}\bar{u}_{xx} - (\bar{u}_x)^2] \qquad (20) $$

$$ \mathcal{Q} = \xi_x(\bar{u}_t + \bar{u}\bar{u}_x) + \xi_{xx}\bar{u}^2 \qquad (21) $$

where $\mathcal{P} + \mathcal{Q}$ and $\mathcal{Q}$ represent, respectively, the vertical fluid acceleration at the free surface and at the bottom (see Appendix A). This fully nonlinear set of equations contains all $\varepsilon$ order of magnitude terms because no hypothesis have been introduced concerning nonlinearities (i.e. $\varepsilon \sim O(1)$). Some of these $\varepsilon$ terms are related to the convective vertical acceleration neglected in low-order Boussinesq equations (e.g. Reference [4]). One of the interesting properties of Serre equations for horizontal bottom lies in the fact that they have a closed-form solitary wave solution, which happens to be the Rayleigh solitary wave [43], expressed as

$$ h(x,t) = h_0 + a \operatorname{sech}^2[K(x - Ct)] \qquad (22) $$

$$ \bar{u} = C \left( 1 - \frac{h_0}{h} \right) \qquad (23) $$

$$ K = \sqrt{\frac{3a}{4h_0^2(h_0 + a)}}, \quad C = C_0 \sqrt{1 + \frac{a}{h_0}} \quad \text{and} \quad C_0 = \sqrt{gh_0} \qquad (24) $$

where $h_0$ is the water depth at rest, $a$ the amplitude, $C$ the phase speed and $K$ the outskirt decay parameter. This solution will be used later to numerically test the proposed finite volume scheme.

Additionally, the range of validity of this set of equations can be easily extended to wave propagation problems in deeper waters using the dispersion correction technique discussed in References [7, 35]. A Padé (2, 2) expansion for Stokes linear dispersion is achieved by adding the following term on the left-hand side of Equation (19):

$$ B(x,t) = -\alpha h \xi^2 [\bar{u}_{xx} + (\bar{u}\bar{u}_x)_{xx} + g(h + \xi)_{xxx}] \qquad (25) $$

where $\alpha$ is an adjustable parameter ($\sim 1/15$) that must be tuned in order to match (at best) the linear Stokes dispersion relation.

Finally, it is possible to show that continuity and momentum equations, including the dispersion correction term and assuming that bottom variations are only function of the spatial coordinate $x$ can be written in a weak quasi-conservative form by introducing a dependent variable $q$ which aggregates all time derivatives of Equations (19) and (25). This convenient mathematical form first given by Seabra-Santos [36] reads

$$ h_t + F_x = 0 \qquad (26) $$
$$q_t + G_x = S \qquad (27)$$

$$q = (1 + r)\bar{u} - \frac{1}{3h}(h^3\bar{u}_x)_x - \alpha\xi^2\bar{u}_{xx} \qquad (28)$$

with flux functions

$$F = h\bar{u} \qquad (29)$$

$$G = q\bar{u} + g(h + \xi) - \frac{1}{2}\bar{u}^2(1 + \xi_x^2) + (\xi_x\bar{u} - \frac{1}{2}h\bar{u}_x)h\bar{u}_x - \alpha\xi^2[\bar{u}_x^2 + \bar{u}\bar{u}_{xx} + g(h + \xi)_{xx}] \qquad (30)$$

and the non-conservative right-hand side term S (analogous to a source term) in (27) has the following expression:

$$S = -2\alpha\xi\bar{u}\bar{u}_x + g(h + \xi)_{xx} \qquad (31)$$

The latter contains a small correction for the dispersive properties of the system. Moreover, $r = (h_x + \xi_x)\xi_x + \frac{1}{2}h\xi_{xx}$, and thus traduces some bathymetric effects. In what follows, overbars will be dropped for convenience.

### 3. NUMERICAL METHOD

Because dispersive terms embedded in Boussinesq-type equations may have a similar mathematical form than lower-order terms associated to truncation errors of discretized equations, high-order schemes, both in time and space, are needed. Furthermore, nonlinear processes and frequency dispersion are responsible for the broadening/narrowing of initial wavelength scales in time and numerical methods used to solve Boussinesq equations must be able to adequately reproduce interactions between different wave modes, without damping physical ones nor introducing spurious short-wave noise. In this context we will develop a 4th-order compact finite volume scheme for spatial integration of equations derived in the previous section over a uniform 1D staggered grid while time integration will be performed using a standard 4th-order Runge-Kutta method.

#### 3.1. Finite volume method for spatial discretization

The computational domain is schematically described in Figure 2 where $\Omega_i$ and $\Omega_j$ are control volumes centred around i-nodes ($i=0, 1, ..., I-1$) and around j-nodes ($j=0, 1, ..., I$),

<figure><img src="image_3.png" /><figcaption>Figure 2. One-dimensional computational domain meshed with staggered control volumes $\Omega_i$ and $\Omega_j$.</figcaption></figure>
respectively, i.e.

$$ \Omega_i = \{x \in [x_j, x_{j+1}]\} \qquad (32) $$

$$ \Omega_j = \{x \in [x_{i-1}, x_i]\} \qquad (33) $$

where $I$ is the number of $\Omega_i$ control volumes used to discretize the physical domain. Nodal coordinates are defined as $x_i = \Delta x/2 + i\Delta x$ and $x_j = j\Delta x$, where $\Delta x$ is the constant spatial grid resolution (length of control volumes). The governing equations will be integrated over this staggered domain in order to optimize the cell-face interpolation process but also to improve accuracy and minimize short-wave noise generation [22, 44].

Equations (26) and (27) can be easily integrated over control volumes $\Omega_i$ as follows:

$$ \frac{\partial}{\partial t} \int_{x_j}^{x_{j+1}} h \, dx + F(x_{j+1}, t) - F(x_j, t) = 0 \qquad (34) $$

$$ \frac{\partial}{\partial t} \int_{x_j}^{x_{j+1}} q \, dx + G(x_{j+1}, t) - G(x_j, t) = \int_{x_j}^{x_{j+1}} S \, dx \qquad (35) $$

where flux functions are evaluated at cell faces of each control volume (j-nodes).

Unlike previous equations, the auxiliary relation (28) which links variable $q$, to the physical variables $h$ and $u$ will be integrated over staggered control volumes $\Omega_j$:

$$ \int_{x_{i-1}}^{x_i} \left[ (1+r)u - \frac{1}{3h}(h^3 u_x)_x - \alpha \xi^2 u_{xx} \right] dx = \int_{x_{i-1}}^{x_i} q \, dx \qquad (36) $$

Numerical approximations will be introduced in the following to evaluate flux functions at cell faces and volume integrals.

In continuity and momentum equations (34) and (35), volume integrals of $h$ and $q$ can be stored as cell-averaged values and no further approximations need to be done. For example, the cell-averaged value of a generic function $\Phi$ over control volume $\Omega_i$ is written as

$$ \hat{\Phi}_i = \frac{1}{\Delta x} \int_{x_j}^{x_{j+1}} \Phi \, dx \qquad (37) $$

Discretized counterparts of Equations (34) and (35) are then obtained by performing the finite volume integration over the whole domain:

$$ \frac{d\hat{h}_i^n}{dt} = -\frac{1}{\Delta x}(F_{j+1}^n - F_j^n) \quad \text{for } i=0,1,\dots,I-1 \quad \text{and} \quad j=i \qquad (38) $$

$$ \frac{d\hat{q}_i^n}{dt} = \hat{S}_i^n - \frac{1}{\Delta x}(G_{j+1}^n - G_j^n) \quad \text{for } i=0,1,\dots,I-1 \quad \text{and} \quad j=i \qquad (39) $$

where we have introduced for simplicity the notation $\Phi_j^n = \Phi(x_j, t_n)$, with the discrete time coordinate $t_n = n\Delta t$ ($n=0, 1, \dots, N$) and $\Delta t$ is the time step. The non-conservative term will be approximated using a centred-in-space approach where

$$ \hat{S}_i = \frac{1}{\Delta x} \int_{x_j}^{x_{j+1}} S \, dx \approx - \frac{2\alpha}{\Delta x} \xi_i(\xi_x)_i (\tilde{S}_{j+1}^n - \tilde{S}_j^n) \qquad (40) $$
and $\tilde{S} = uu_x + g(h + \xi)_x$. This discretized expression will introduce only slight errors in the overall numerical approximation since this term represents a small correction of $O(\sigma^4)$ for dispersive properties.

On the other hand, evaluation of flux functions at cell faces requires a reconstruction of $h_j$ and $q_j$ values from cell-averaged variables because time integration of Equations (38) and (39) will produce $\hat{h}_i$ and $\hat{q}_i$. Here, cell face values will be estimated at each time step using the implicit 4th-order interpolation technique developed in References [38, 39]. This approach results in the following system of $I-1$ equations for internal nodal points:

$$ \frac{1}{4} \Phi_{j-1} + \Phi_j + \frac{1}{4} \Phi_{j+1} = \frac{3}{4} (\hat{\Phi}_{i-1} + \hat{\Phi}_i) \quad \text{for } j=1,2,\dots,I-1 \text{ and } i=j \quad (41) $$

At boundary nodes, $j=0$ and $j=I$, boundary conditions must provide the correct values for $h$ and $q$.

When using Equation (41), 4th-order accuracy is achieved in the evaluation of cell face values by solving a tridiagonal system. On the contrary, using a pointwise 4th-order method where volume integrals are estimated using Simpson's rule and cell face values are interpolated with piecewise cubic polynomials, it is necessary to solve a pentadiagonal system (see, for instance, Reference [30]). Thus, computational cost is significantly reduced when implementing the present compact scheme.

Evaluation of flux functions requires, in addition, the prescription of depth-averaged velocity values at cell faces. The latter are obtained at the end of each time step, once $h_j$ and $q_j$ have been computed at internal nodes, by solving the auxiliary relation (36). Different integrals appearing in this equation are approximated as follows:

$$ \begin{aligned} & \frac{1}{\Delta x} \int_{x_{i-1}}^{x_i} (1+r)u \, dx \approx (1+r_j)\hat{u}_j \\ & -\frac{1}{\Delta x} \int_{x_{i-1}}^{x_i} \frac{1}{3h}(h^3 u_x)_x \, dx \approx -\frac{1}{3\Delta x h_j}[h_i^3(u_x)_i - h_{i-1}^3(u_x)_{i-1}] \\ & -\frac{1}{\Delta x} \int_{x_{i-1}}^{x_i} \alpha \xi^2 u_{xx} \, dx \approx -\frac{\alpha}{\Delta x} \xi_j^2 [(u_x)_i - (u_x)_{i-1}] \\ & \frac{1}{\Delta x} \int_{x_{i-1}}^{x_i} q \, dx = \hat{q}_j \end{aligned} $$

At this point, some care must be taken for the reconstruction of depth-averaged velocity values at $j$-nodes using Equation (36). Since we want to obtain $u_j$ directly and still preserve a tridiagonal structure in the discrete system, it is necessary to link cell-averaged values $\hat{u}_j$ and diffusive fluxes to pointwise values at $j$-nodes carefully. The latter can be achieved using the following strategy.

On the one hand, combining Simpson's rule to evaluate volume integrals over $\Omega_i$ and $\Omega_j$,

$$ \hat{\Phi}_i = \frac{2}{3} \Phi_i + \frac{1}{6} (\Phi_j + \Phi_{j+1}) + O(\Delta x^4) \quad (42) $$

$$ \hat{\Phi}_j = \frac{2}{3} \Phi_j + \frac{1}{6} (\Phi_{i-1} + \Phi_i) + O(\Delta x^4) \quad (43) $$
and Equation (41), provides the following relation:

$$ \hat{\Phi}_j = \frac{1}{24}(\Phi_{j-1} + 22\Phi_j + \Phi_{j+1}) + O(\Delta x^4) \quad (44) $$

which still has a tridiagonal matrix form. On the other hand, the velocity gradient $(u_x)_i$ must be written in terms of pointwise u-values at j-nodes. Therefore, it is preferable to use an explicit formula to approximate its spatial derivative. However, an explicit 4th-order formula will break down the tridiagonal nature of the system because 5-point stencils will be necessary to compute diffusive fluxes. One way to preserve a tridiagonal structure is to use a deferred-correction approach, which is expressed as [45]

$$ (u_x)_i = \frac{(u_{j+1} - u_j)}{\Delta x} + \left\{ \frac{(u_{j-1} - 27u_j + 27u_{j+1} - u_{j+2})}{24\Delta x} - \frac{(u_{j+1} - u_j)}{\Delta x} \right\}^{\text{old}} \quad (45) $$

where superscript 'old' is used for depth-averaged velocities evaluated at the previous level of the iterative process. Thus, the discretized equation associated to (36) must be solved iteratively because a 2nd-order approximation is used for the implicit part, and a 4th-order formula appears in the explicit contribution in brackets. When some convergence criteria will be satisfied between 'old' and 'new' $u_j$-values, the overall approximation for $u_x$ will be 4th-order accurate. In practice only few iterations will be needed to obtain sufficient accuracy. The convergence criteria for fixed-point iterations is based on a $L_1$ norm and reads

$$ \frac{\|\mathbf{u}_j - \{\mathbf{u}_j\}^{\text{old}}\|_1}{\|\mathbf{u}_j\|_1} \le \tau \quad (46) $$

where $\mathbf{u}_j = [u_{j=2} \ u_{j=3} \ \dots \ u_{j=I-2} \ u_{j=I-1}]^T$ is the vector containing cell face values of the most recently computed depth-averaged velocity and $\tau$ is the tolerance for the relative error between successive iterations. Hence, internal nodal values of $u_j$ are obtained by solving the following implicit tridiagonal system:

$$ \left\{ \begin{aligned} & \left\{ \frac{(1+r_j)}{24} - \frac{1}{\Delta x^2} \left( \frac{h_i^3}{3h_j} + \alpha \xi_j^2 \right) \right\} u_{j-1} + \left\{ \frac{22(1+r_j)}{24} + \frac{1}{\Delta x^2} \left[ \frac{h_{i-1}^3 + h_i^3}{3h_j} + 2\alpha \xi_j^2 \right] \right\} u_j \\ & + \left\{ \frac{(1+r_j)}{24} - \frac{1}{\Delta x^2} \left( \frac{h_i^3}{3h_j} + \alpha \xi_j^2 \right) \right\} u_{j+1} \\ & = \frac{1}{24}(q_{j-1} + 22q_j + q_{j+1}) + \frac{1}{\Delta x^2} \left\{ \frac{h_i^3}{3h_j} + \alpha \xi_j^2 \right\} U_i^{\text{old}} \\ & - \frac{1}{\Delta x^2} \left\{ \frac{h_{i-1}^3}{3h_j} + \alpha \xi_j^2 \right\} U_{i-1}^{\text{old}} \quad \text{for } j=1,2,\dots,I-1 \text{ and } i=j \end{aligned} \right. \quad (47) $$

where

$$ U_i^{\text{old}} = \frac{1}{24} \{u_{j-1} - 3u_j + 3u_{j+1} - u_{j+2}\}^{\text{old}} \quad (48) $$

Near boundaries, evaluation of the last relation will require to compute depth-averaged velocity values outside the physical domain. In the case of non-periodic boundaries, ghost
nodes values will be extrapolated from interior nodes using suitable piecewise polynomial functions. Similarly, boundary conditions must be treated with caution because there is a coupling between the cell face reconstruction of $q_j$ values (Equation (41)) and the right-hand side of Equation (47). However, in the present work only periodic boundary conditions will be considered. The special treatment of boundary and ghost nodes in the framework of wave propagation over a beach will be addressed in a companion paper.

Finally, derivatives of water depth, $h$, and depth-averaged velocity, $u$, must be computed at $j$-nodes in order to complete the evaluation of the flux function $G_j$, the source term $\hat{S}_i$ and the dependent variable $r_j$ at each time step. Similarly, water depth values at $i$-nodes must be interpolated to compute matrix coefficients in Equation (47) while bottom function $\xi$ and its derivatives must be prescribed at $i$ and $j$ nodes. The latter is not difficult because bottom bathymetry constitutes a known input data, which is not time dependent. On the other hand, $i$-values of water depth will require an interpolation which can be performed easily from computed values of $\hat{h}_i$ and $h_j$ using relation (42). It is relevant to note that the strategy of overlapping the control volumes, where Equation (36) is integrated, has allowed us to minimize the number of interpolated nodal values.

Estimates of derivatives are performed using a 4th-order finite difference compact scheme as described in Reference [40]. These compact formulae have better spectral properties for short-length waves than their explicit 4th-order finite difference counterparts and matching conditions at boundary nodes are less tedious. It is important to note that the highest derivative that one needs to estimate here will be of 2nd-order. In fact, the finite volume method applied to the quasi-conservative set of equations derived here has allowed us to lower the highest $x$-derivative order by one when compared to other equivalent set of equations (e.g. References [17, 27]). Derivative estimates at internal $j$-nodes are obtained from cell face values using the following implicit formulae [40]:

$$ \frac{1}{4}(\Phi_x)_{j-1} + (\Phi_x)_j + \frac{1}{4}(\Phi_x)_{j+1} = \frac{3}{4\Delta x}(\Phi_{j+1} - \Phi_{j-1}) \quad \text{for } j=1,2,\dots,I-1 \text{ and } i=j \quad (49) $$

$$ \frac{1}{10}(\Phi_{xx})_{j-1} + (\Phi_{xx})_j + \frac{1}{10}(\Phi_{xx})_{j+1} = \frac{6}{5\Delta x^2}(\Phi_{j+1} - 2\Phi_j + \Phi_{j-1}) \quad \text{for } j=1,2,\dots,I-1 \text{ and } i=j \quad (50) $$

Again, special care must be taken at boundary nodes as discussed in References [40, 46]. It is worth to point out that these compact forms can be verified manipulating relation (44).

Spatial discretization of Equations (34)–(36) using the present finite volume technique allows us to write the original system of partial differential equations (PDE) as the following set of ordinary differential equations (ODE),

$$ \frac{d \hat{\mathbf{x}}_i^n}{dt} = \mathbf{f}_i^n \qquad (51) $$

$$ \mathcal{A}(\mathbf{h}_i^n) \cdot \mathbf{u}_j^n = \mathcal{B} \cdot \mathbf{q}_j^n + \mathcal{C}(\mathbf{h}_j^n) \cdot \{\mathbf{u}_j^n\}^{\text{old}} \quad (52) $$
where

$$ \hat{\mathbf{x}}_i^n = \begin{bmatrix} \hat{\mathbf{h}}_i^n \\ \hat{\mathbf{q}}_i^n \end{bmatrix} \quad \text{and} \quad \mathbf{f}_i^n = \begin{bmatrix} \mathbf{0} \\ \hat{\mathbf{S}}_i^n \end{bmatrix} - \begin{bmatrix} \Delta \mathbf{F}_i^n \\ \Delta \mathbf{G}_i^n \end{bmatrix} $$

Here, bold letters represent vectors containing cell-averaged values of associated variables over the whole discrete domain. For instance, a generic vector $\hat{\mathbf{\Phi}}_i^n$ contains cell-averaged values over control volumes $\Omega_i$ of function $\Phi$ evaluated at time $t = t_n$. Similarly, the notation $\Delta\mathbf{\Phi}_i^n = (1/\Delta x)(\mathbf{\Phi}_i^{n+1} - \mathbf{\Phi}_i^n)$ indicates a cell face flux difference vector, and $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ are matrices obtained from the discretized equation (47). In the next subsection, we will use an explicit 4th-order Runge-Kutta method to numerically integrate the system of ODE just derived.

## 3.2. Numerical time stepping

The implementation of a high-order numerical ODE solver for time integration is justified since spatial derivatives are discretized using a high-order scheme as well. The finite volume discretization technique is theoretically 4th-order accurate in space so we wish to reach a similar discretization error when integrating the system (51)–(52) in time. The 4th-order Runge-Kutta method is an attractive alternative for the solution of wave-type equations and can be quite efficient when compared to other high-order schemes [44]. In the present context, the standard four-stage Runge-Kutta method can be stated as follows:

First estimation at $t = t_{n+1/2}$

$$ \mathbf{p}^{(1)} = \Delta t \mathbf{f}_i^n, \quad \hat{\mathbf{x}}_i^{n+(1/2)} = \hat{\mathbf{x}}_i^n + \frac{1}{2} \mathbf{p}^{(1)} \qquad (53a) $$

$$ \mathbf{u}_j^{n+(1/2)} = \mathcal{A}(\mathbf{h}_j^{n+(1/2)})^{-1} \cdot [\mathcal{B} \cdot \mathbf{q}_j^{n+(1/2)} + \mathcal{C}(\mathbf{h}_j^{n+(1/2)}) \cdot \{\mathbf{u}_j^{n+(1/2)}\}^{\text{old}}] \qquad (53b) $$

Second estimation at $t = t_{n+1/2}$

$$ \mathbf{p}^{(2)} = \Delta t \mathbf{f}_i^{n+(1/2)}, \quad \hat{\mathbf{x}}_i^{n+(1/2)} = \hat{\mathbf{x}}_i^n + \frac{1}{2} \mathbf{p}^{(2)} \qquad (53c) $$

$$ \mathbf{u}_j^{n+(1/2)} = \mathcal{A}(\mathbf{h}_j^{n+(1/2)})^{-1} \cdot [\mathcal{B} \cdot \mathbf{q}_j^{n+(1/2)} + \mathcal{C}(\mathbf{h}_j^{n+(1/2)}) \cdot \{\mathbf{u}_j^{n+(1/2)}\}^{\text{old}}] \qquad (53d) $$

First estimation at $t = t_{n+1}$

$$ \mathbf{p}^{(3)} = \Delta t \mathbf{f}_i^{n+(1/2)}, \quad \hat{\mathbf{x}}_i^{n+1} = \hat{\mathbf{x}}_i^n + \frac{1}{2} \mathbf{p}^{(3)} \qquad (53e) $$

$$ \mathbf{u}_j^{n+1} = \mathcal{A}(\mathbf{h}_j^{n+1})^{-1} \cdot [\mathcal{B} \cdot \mathbf{q}_j^{n+1} + \mathcal{C}(\mathbf{h}_j^{n+1}) \cdot \{\mathbf{u}_j^{n+1}\}^{\text{old}}] \qquad (53f) $$

Second and final estimation at $t = t_{n+1}$

$$ \mathbf{p}^{(4)} = \Delta t \mathbf{f}_i^{n+1}, \quad \hat{\mathbf{x}}_i^{n+1} = \hat{\mathbf{x}}_i^n + \frac{1}{6}(\mathbf{p}^{(1)} + 2\mathbf{p}^{(2)} + 2\mathbf{p}^{(3)} + \mathbf{p}^{(4)}) \qquad (53g) $$

$$ \mathbf{u}_j^{n+1} = \mathcal{A}(\mathbf{h}_j^{n+1})^{-1} \cdot [\mathcal{B} \cdot \mathbf{q}_j^{n+1} + \mathcal{C}(\mathbf{h}_j^{n+1}) \cdot \{\mathbf{u}_j^{n+1}\}^{\text{old}}] \qquad (53h) $$

At each step, j-values of water depth, h, and auxiliary variable q must be reconstructed from cell-averaged quantities using the implicit relation (41) before solving the equation for depth-averaged velocity, u. As explained in the previous subsection, the latter is solved iteratively
and at each stage the last computed value of $u_j$ is used as first guess. A relative tolerance for iteration errors, $\tau$, of $10^{-4}$ has shown to give good accuracy by requiring only 1–3 iterations. On the other hand, the numerical scheme requires at each step the inversion of tridiagonal matrices which can be performed at low computational cost using the Thomas algorithm (see, for instance, Reference [47]). The number of operations needed to solve each matrix equation will then be linearly proportional to the number of unknowns in the system.

## 4. SPECTRAL ANALYSIS AND LINEAR STABILITY

Important information about the general performance of our finite volume resolution can be obtained by analysing its linearized version. In this section we will study in particular the stability, accuracy in spectral space and dispersion relation preserving properties (DRP) of the scheme using spectral analysis [46]. Our main goal here is to find necessary conditions that will ensure that the numerical scheme will behave properly when describing wave propagation phenomena. For simplicity, we will restrict ourselves to the analysis of linear wave propagation over flat bottoms.

When using relation (41) in Equations (38)–(39) and assuming that fixed-point iterations have converged (i.e. {$u_j^n$}old = $u_j^n$) in (47), the linearized version of the discretized system of equations can be recast as

$$ \frac{1}{4} \frac{d h_{j-1}^n}{dt} + \frac{d h_j^n}{dt} + \frac{1}{4} \frac{d h_{j+1}^n}{dt} = - \frac{3}{4 \Delta x} (F_{j+1}^n - F_{j-1}^n) \quad (54) $$

$$ \frac{1}{4} \frac{d q_{j-1}^n}{dt} + \frac{d q_j^n}{dt} + \frac{1}{4} \frac{d q_{j+1}^n}{dt} = - \frac{3}{4 \Delta x} (G_{j+1}^n - G_{j-1}^n) \quad (55) $$

$$ \frac{\beta}{\rho^2} u_{j-2}^n + \left(1 - \frac{28\beta}{\rho^2}\right) u_{j-1}^n + \left(22 + \frac{54\beta}{\rho^2}\right) u_j^n + \left(1 - \frac{28\beta}{\rho^2}\right) u_{j+1}^n + \frac{\beta}{\rho^2} u_{j+2}^n = q_{j-1}^n + 22q_j^n + q_{j+1}^n \quad (56) $$

where

$$ F_j^n = h_0 u_j^n, \quad G_j^n = gh_j^n - \alpha h_0^2 g(h_{xx})_j^n, \quad \beta = (1/3 + \alpha), \quad \rho = \Delta x/h_0 $$

and $h_0 = |\xi|$ is the still water depth of the horizontal bottom flume. This set of equations will be used in what follows in order to analyse the numerical properties of the scheme. It is interesting to note the close resemblance between Equations (54)–(55) and the implicit formula for the first spatial derivative estimate (49).

### 4.1. Semi-discrete system of equations and numerical dispersion

Spectral accuracy of the spatial discretization is studied first by decomposing the numerical solution into finite Fourier series. It is important to point out that spectral accuracy is related
to the ability of the scheme to resolve the different length scales present in the problem [46]. Therefore, it is paramount to assess how well the numerical scheme is able to describe the physical range of wavelengths. In order to perform this analysis we replace the following individual Fourier components into the discretized system of equations:

$$ \begin{bmatrix} h_j^n \\ q_j^n \\ u_j^n \end{bmatrix} = \begin{bmatrix} h_\kappa \\ q_\kappa \\ u_\kappa \end{bmatrix} e^{i(\kappa j \Delta x - \omega_r t)} \quad (57) $$

where $\hat{i} = \sqrt{-1}$ is the imaginary constant, $h_\kappa$, $q_\kappa$ and $u_\kappa$ are Fourier amplitudes, $\kappa$ is the wave number of the individual component and $\omega_r$ is the associated numerical frequency. Numerical dispersion can be studied by substituting (57) in Equations (54), (55) and (56) and solving the resulting system of equations for $\omega_r$. Before doing so, we need to find the spectral expression for the second spatial derivative that appears in the flux function $G_j^n$. The numerical estimation of this derivative has the form $(h_{xx})_j = -\tilde{\kappa}^2 h_\kappa e^{i(\kappa j \Delta x - \omega_r t)}$, where by substitution in the implicit difference formula (50) the associated numerical wave number is written as

$$ (\tilde{\kappa} \Delta x)^2 = 12 \frac{(1 - \cos \kappa \Delta x)}{(5 + \cos \kappa \Delta x)} \quad (58) $$

Similarly, substituting the Fourier components (57) into Equation (56) gives the following relation between $q_\kappa$ and $u_\kappa$:

$$ q_\kappa = \Psi u_\kappa \quad (59) $$

where

$$ \Psi = 1 + \frac{(1/3 + \alpha)}{\rho^2} \frac{(27 + \cos 2\kappa\Delta x - 28 \cos \kappa\Delta x)}{(11 + \cos \kappa\Delta x)} \quad (60) $$

Finally, if we use the Fourier components (57) and relation (59) in continuity and momentum equations (55)–(56), the compact finite volume scheme takes the following form in the Fourier space:

$$ (\mathcal{H} + \mathbf{I} \hat{\omega}_r) \cdot \begin{bmatrix} h_\kappa \\ u_\kappa \end{bmatrix} = 0 \quad (61) $$

where

$$ \mathcal{H} = \frac{3}{\rho} \begin{pmatrix} 0 & K_1 \\ K_2 & 0 \end{pmatrix} $$

is the spatial discretization operator whose components read,

$$ K_1 = -\hat{i} \frac{\sin \kappa \Delta x}{(2 + \cos \kappa \Delta x)} \quad (62) $$
$$ K_2 = \frac{K_1}{\Psi} [\rho^2 + \alpha(\tilde{\kappa}\Delta x)^2] \left(\frac{C_0}{\Delta x}\right)^2 \quad (63) $$

I is the identity matrix and $C_0 = \sqrt{gh_0}$. Non-trivial solution of Equation (61) provides the discrete dispersion relation associated to our compact finite volume method which can be written in a non-dimensional form as

$$ \begin{aligned} \frac{C_{\text{fv}}^2}{C_0^2} &= \frac{1}{C_0^2} \left(\frac{\omega_r}{\kappa}\right)^2 = -\left(\frac{3}{\kappa\Delta x}\right)^2 K_1^2 \frac{[\rho^2 + \alpha(\tilde{\kappa}\Delta x)^2]}{\rho^2 \Psi} \\ &= \left(\frac{3 \sin \kappa\Delta x}{\kappa\Delta x(2 + \cos \kappa\Delta x)}\right)^2 \\ &\quad \times \frac{(11 + \cos \kappa\Delta x)[\rho^2 + \alpha(\tilde{\kappa}\Delta x)^2]}{[\rho^2(11 + \cos \kappa\Delta x) + (1/3 + \alpha)(27 + \cos 2\kappa\Delta x - 28 \cos \kappa\Delta x)]} \end{aligned} \quad (64) $$

with $C_{\text{fv}}$ being the numerical phase speed associated to the finite volume scheme. Therefore, the non-dimensional discrete dispersion relation depends on the wave number, $\kappa\Delta x$, and on the parameter $\rho = \Delta x/h_0$. Although this general expression cannot be reduced to a simple form, some relevant information may be obtained by exploring its asymptotic behaviour as $\rho \to 0$ (which can be viewed as the deep water limit) and as $\rho \to \infty$ (which corresponds to the shallow water limit).

In order to do some comparisons between the discretized system of equations and the original set of PDE we recall that the non-dimensional *physical* dispersion relation that belongs to the linearized Boussinesq model can be written as [41]

$$ \frac{C_{\text{th}}}{C_0} = \left[ \frac{1 + \alpha(\kappa h_0)^2}{1 + (1/3 + \alpha)(\kappa h_0)^2} \right]^{1/2} = \left[ \frac{1 + \alpha \left( \frac{\kappa \Delta x}{\rho} \right)^2}{1 + (1/3 + \alpha) \left( \frac{\kappa \Delta x}{\rho} \right)^2} \right]^{1/2} \quad (65) $$

where $C_{\text{th}}$ stands for the phase speed associated to the theoretical dispersion relation produced by the extended system of Serre equations.

In the shallow water limit it can be shown that the ratio between the *physical* and the discrete dispersion relation is

$$ \lim_{\rho \to \infty} \frac{C_{\text{fv}}}{C_{\text{th}}} = \frac{3 \sin \kappa \Delta x}{\kappa \Delta x (2 + \cos \kappa \Delta x)} = 1 - \frac{1}{180} (\kappa \Delta x)^4 + O[(\kappa \Delta x)^6] \quad (66) $$

proving that the numerical scheme mimics the exact relation with 4th-order accuracy in $\kappa\Delta x$. Similarly, in the deep water limit we have

$$ \begin{aligned} \lim_{\rho \to 0} \frac{C_{\text{fv}}}{C_{\text{th}}} &= \frac{3 \sin \kappa \Delta x}{\kappa \Delta x (2 + \cos \kappa \Delta x)} \left[ \frac{6(11 + \cos \kappa \Delta x)}{(13 - \cos \kappa \Delta x)(5 + \cos \kappa \Delta x)} \right]^{1/2} \\ &= 1 - \frac{11}{2880} (\kappa \Delta x)^4 + O[(\kappa \Delta x)^6] \end{aligned} \quad (67) $$
<figure><img src="image_4.png" /><figcaption>Figure 3. Leading-order truncation error term vs $\sigma = \kappa h_0$ for, (-) $\alpha = 0$, (· · ·) $\alpha = 1/30$, (- -) $\alpha = 1/20$ and (· · -) $\alpha = 1/10$: (a) full 4th-order compact finite volume scheme; and (b) mixed 4th- and 2nd-order discretization strategy.</figcaption></figure>

which shows that the phase speed error is also of 4th-order in $\kappa\Delta x$. In the general case the leading-order truncation error term of $O[(\kappa\Delta x)^4]$ will be a function of the dispersive characteristics of propagating waves, $\sigma = \kappa h_0$. In fact, it can be demonstrated that the maximum value for the leading-order coefficient will always occur at the shallow water limit, where $\sigma \to 0$, independently of the chosen $\alpha$ value as shown in Figure 3(a). It is also interesting to note that the asymptotic behaviour in deep water (i.e. $\sigma \to \infty$) changes when $\alpha = 0$ converging then to $-1/576$.

Even though it is generally acknowledged that numerical dispersion constitutes an important property of any scheme, dispersion relations associated to previously published finite difference or finite element high-order Boussinesq-type models have rarely been reported. Moreover, it is commonly argued that numerical dispersion can be avoided if $O(\sigma^2)$ dispersive terms are discretized with 2nd-order truncation error approximations, while convective or 1st-order terms are approximated through 4th-order centred finite difference formulae (e.g. References [18, 24]). Even if from a practical point of view the latter could be sufficient for well-resolved waves, the numerical scheme will only be 2nd-order accurate in $\kappa\Delta x$ and in the general case more nodal points per wavelength will be needed to reproduce the associated phase speed correctly. In order to illustrate this we will discretize high-order dispersive terms appearing in the flux function $G_j^n$ of Equation (55) and in the linearized version of the auxiliary equation (36) by using explicit 2nd-order finite difference approximations. It is worth emphasizing that
no direct comparison can be made with the finite difference technique used in References [17, 24] or with finite element models based on similar mixed-order approaches [26, 27] since our numerical resolution is based on a compact finite volume integration on a staggered grid. In spite of that, the following analysis can be useful to see what consequences this kind of strategy may have in our particular numerical method. Thus, 2nd-order derivatives are approximated now using the following explicit 2nd-order formula:

$$ (\Phi_{xx})_j = \frac{1}{\Delta x^2} (\Phi_{j+1} - 2\Phi_j + \Phi_{j-1}) \qquad (68) $$

and the corresponding spectral relation is found to be (after substitution of a Fourier component),

$$ (\tilde{\kappa}' \Delta x)^2 = 2(1 - \cos \kappa \Delta x) \qquad (69) $$

Similarly, using 2nd-order approximations for high-order derivatives in the linearized version of auxiliary equation (36) it can be shown that

$$ \Psi' = 1 + \frac{(1/3 + \alpha) 24(1 - \cos \kappa \Delta x)}{\rho^2 (11 + \cos \kappa \Delta x)} \qquad (70) $$

Substituting relations (69) and (70) in (64) provides the following form for the discrete dispersion relation associated to the mixed approach:

$$ \frac{C'_{fv}}{C_{th}} = 1 - \frac{\sigma^2 [(\alpha - 1/3) + \alpha(\alpha + 1/3)\sigma^2]}{48 [1 + \alpha\sigma^2][1 + (\alpha + 1/3)\sigma^2]} (\kappa\Delta x)^2 + O[(\kappa\Delta x)^4] \qquad (71) $$

Therefore, using 2nd-order finite difference formulae to evaluate dispersive terms results in an overall reduction in the spectral resolution of the method. However, the coefficient that multiplies the $O[(\kappa\Delta x)^2]$ term will be relatively small (but positive!) near the shallow water limit (see Figure 3(b)) in such a way that the accuracy of the scheme would not be dramatically affected. The latter constitutes the main argument behind the mixed approach used in several finite difference and finite element Boussinesq-type models. Nevertheless, we will show in the rest of the section that the full 4th-order scheme has better numerical properties and is more robust than its mixed counterpart.

Another important and challenging property to test is the ability of a particular numerical scheme to correctly reproduce group velocities. This property is related to wave energy propagation and usually imposes more severe restrictions on the application of numerical methods than the dispersion relation preserving condition. The wave group velocity associated to the numerical scheme and the linearized set of PDE can be estimated from the non-dimensional dispersion relation as,

$$ \frac{1}{C_0} \frac{\partial \omega}{\partial \kappa} = \frac{C_g}{C_0} = \frac{\partial}{\partial(\kappa \Delta x)} \left[ (\kappa \Delta x) \frac{C}{C_0} \right] \qquad (72) $$

Again, the general discrete relation will depend on $\kappa\Delta x$ and $\rho$ but asymptotic expressions can be found as before. In the shallow water limit the ratio between the numerical and the physical group velocity reads as

$$ \lim_{\rho \to \infty} \frac{C_{gfv}}{C_{gth}} = 1 - \frac{1}{36}(\kappa \Delta x)^4 + O[(\kappa \Delta x)^6] \qquad (73) $$
<figure><img src="image_5.png" /><figcaption>Figure 4. Ratio between computed and physical phase speed vs $\kappa\Delta x$ for $\alpha = 1/15$, when (—) $\rho \to \infty$, (- -) $\rho \to 0$ and (---) $\rho = 0.1$: (a) full 4th-order compact finite volume scheme; and (b) mixed 4th- and 2nd-order discretization strategy.</figcaption></figure>

while in the deep-water case we have

$$ \lim_{\rho \to 0} \frac{C_{gfv}}{C_{gth}} = 1 - \frac{11}{576} (\kappa \Delta x)^4 + O[(\kappa \Delta x)^6] \quad (74) $$

showing that the truncation error is still 4th-order in $\kappa\Delta x$ but with a leading-order term an order of magnitude larger than in the phase speed case. Similarly, it is possible to demonstrate that when the mixed-order strategy is used the overall truncation error for the group velocity estimate is 2nd-order in $\kappa\Delta x$, except in the shallow water limit ($\rho \to \infty$) where it has the same form as (73).

Phase speed estimates using both approaches, the full 4th-order and the mixed one, are plotted using $\alpha = 1/15$ in Figure 4. It can be seen that when $\rho \to \infty$, both strategies provide the same results. This means that for a given water depth, $h_0$, phase speed estimates get better in the mixed model if $\Delta x \to \infty$ which is of course rather problematic. Conversely, for a fixed spatial grid resolution, the modified finite volume method tends to give better phase speed estimates as the depth becomes shallower. However, larger errors could arise in deepest regions of the computational domain exactly where dispersive effects become predominant and the model is intended to perform better. On the contrary, the parameter $\rho$ has only a marginal influence on phase speed errors when the full 4th-order model is used. Thus, the full 4th-order compact finite volume method is more robust, in the sense that phase speed
estimates are less sensitive to changes in the $\rho$ parameter, than the mixed-order version. For comparison, numerical results obtained with $\rho=0.1$, which is a reasonable value for that parameter, are also plotted in Figure 4. It is seen that the numerical estimate converges to the physical solution (within a 0.5% error) when $\kappa\Delta x \approx 0.325\pi$ ($L/\Delta x \approx 6.1$) in the full 4th-order case, and when $\kappa\Delta x \approx 0.202\pi$ ($L/\Delta x \approx 9.9$) in the mixed one. Hence, for a reasonable value of $\rho$, the mixed version requires nearly 62% more nodal points per wavelength than the original one in order to properly describe linear phase speed. As explained before, the phase speed error in the mixed-order model can be reduced if the parameter value $\rho = \Delta x/h_0$ is increased, however, the worst error occurs when $\Delta x/h_0 \to 0$. In that limit the 0.5 % error level is reached for $\kappa\Delta x \approx 0.152\pi$ ($L/\Delta x \approx 13.2$). Compared to the full 4th-order scheme and for a given water depth more than twice nodal points per wavelength will be needed to achieve the same phase speed accuracy in the limit of very good grid resolution. Finally, and in accordance with the leading-order truncation error analysis made before, close examination of Figure 4 shows that when $\rho=0.1$ phase speed of well-resolved waves (i.e. $\kappa\Delta x \ll \pi$) are slightly overestimated by the mixed-order scheme.

Similar qualitative conclusions can be drawn when comparing numerical results for group velocities presented in Figure 5. The main difference here is that one obtains negative values for group velocities in the case of poorly resolved waves ($\kappa\Delta x > 0.67\pi$). On the other hand, in order to describe group velocities within a 0.5% error when $\rho=0.1$, the full 4th-order

<figure><img src="image_6.png" /><figcaption>Figure 5. Ratio between computed and physical group velocity vs $\kappa\Delta x$ for $\alpha = 1/15$, when (—) $\rho \to \infty$, (- -) $\rho \to 0$ and (· · ·) $\rho = 0.1$: (a) full 4th-order compact finite volume scheme; and (b) mixed 4th- and 2nd-order discretization strategy.</figcaption></figure>
finite volume method requires $\kappa\Delta x \approx 0.213\pi$ ($L/\Delta x \approx 9.4$) while the mixed scheme reaches the same accuracy when $\kappa\Delta x \approx 0.117\pi$ ($L/\Delta x \approx 17.2$). It is confirmed then that proper numerical representation of group velocities is a more challenging task because more nodal points per wavelength are needed to converge to the physical solution and poorly resolved waves will transmit their energy in the wrong direction. In this particular test, the full 4th-order numerical scheme shows an even better performance when compared to the mixed version because the latter requires roughly 83% more nodal points per wavelength in order to reach the same level of spectral resolution. In the limit of very good grid resolution for a given water depth ($\Delta x/h_0 \to 0$) the mixed-order model reaches the 0.5% error in group velocities when $\kappa\Delta x \approx 0.088$, thus requiring more than 22 nodal points per wavelength.

The analysis that has been presented not only demonstrates the clear superiority of the full 4th-order version, but also shows that mixed-order approaches must be taken with care since the associated discrete dispersion relation has a misleading dependence on $\rho = \Delta x/h_0$. On the other hand, mixed-order strategies have been introduced in previous numerical models to avoid pentadiagonal or more complicated patterns in matrices that need to be inversed. We believe that compact schemes or deferred-correction strategies constitute a more satisfactory way to overcome this problem in the framework of dispersive wave propagation phenomena.

## 4.2. Analysis of the full discrete system

In order to achieve a comprehensive numerical analysis of the scheme we study in this subsection the full discrete system of equations, which takes into account the Runge-Kutta time stepping. Individual Fourier components that one needs to substitute in linearized discrete Equations (54), (55) and (56) have now the following form:

$$ \begin{bmatrix} h_j^n \\ q_j^n \\ u_j^n \end{bmatrix} = \begin{bmatrix} h_\kappa \\ q_\kappa \\ u_\kappa \end{bmatrix} e^{i(\kappa j \Delta x - (\omega_r + i\omega_i) n \Delta t)} = \begin{bmatrix} h_\kappa \\ q_\kappa \\ u_\kappa \end{bmatrix} |\hat{\lambda}|^n e^{i(\kappa j \Delta x - \omega_r n \Delta t)} \quad (75) $$

where $\omega_r$ and $\omega_i$ correspond to real and imaginary parts of the frequency and $|\hat{\lambda}| = e^{\omega_i \Delta t}$ is the amplification factor of the numerical scheme. It can be seen that the analysis expounded in the previous subsection on the semi-discrete system of equations will be equivalent to the one we are going to present now in the limit of good numerical temporal resolution as $\omega_i \Delta t \to 0$. Indeed, the real part of the frequency will give information about phase speed errors (numerical dispersion), and the associated imaginary part contains information about the stability of the scheme (numerical amplification/dissipation). Numerical properties of the discrete scheme are now driven by three free parameters, $\kappa\Delta x$, $\rho = \Delta x/h_0$ and the Courant number, $C_r = C_0 \Delta t/\Delta x$.

It is worth to point out that analogous linear stability analysis have been performed recently on a finite element resolution of an equivalent Boussinesq set of equations by Woo and Liu [27], and on several finite volume resolutions of a KdV-like system of equations by Bradford and Sanders [29]. In what follows we will compare the numerical performance of our scheme with those previously published. Similarly, some comparisons between the full 4th-order method and the mixed-order version will be given for completion.
As before, we substitute Fourier components in the discrete system of equations (75). After using relation (59) in order to eliminate the auxiliary variable $q$ in the 4th-order Runge-Kutta time stepping process (53), it can be shown that the full discrete scheme takes the following form in the Fourier space,

$$ (\mathcal{G}(\mathcal{H}) - \mathbf{I}|\hat{\lambda}|e^{-\hat{\iota}\omega_r \Delta t}) \cdot \begin{bmatrix} h_{\kappa} \\ u_{\kappa} \end{bmatrix} = 0 \qquad (76) $$

where $\mathcal{G}(s) = \sum_{p=0}^{4} \frac{1}{\rho^p} (\Delta t \cdot s)^p$ is the Runge-Kutta time stepping operator and $\mathcal{H}$ is the spatial discretization operator introduced in the previous subsection. Non-trivial solution of Equation (76) yields

$$ \begin{aligned} & \left[ 1 + \frac{9}{2} \left( \frac{\Delta t}{\rho} \right)^2 K_1 K_2 + \frac{27}{8} \left( \frac{\Delta t}{\rho} \right)^4 (K_1 K_2)^2 - |\hat{\lambda}| e^{-\hat{\iota} \omega_r \Delta t} \right]^2 \\ & = \frac{9}{4} \left( \frac{\Delta t}{\rho} \right)^2 K_1 K_2 \left[ 2 + 3 \left( \frac{\Delta t}{\rho} \right)^2 K_1 K_2 \right]^2 \end{aligned} \qquad (77) $$

Invoking definitions (62) and (63) we write

$$ K = \left(\frac{\Delta t}{\rho}\right)^2 K_1 K_2 = - \frac{\sin^2 \kappa \Delta x}{(2 + \cos \kappa \Delta x)^2} \frac{[\rho^2 + \alpha(\tilde{\kappa} \Delta x)^2]}{\rho^2 \Psi} C_r^2 \qquad (78) $$

Thus, equating real and imaginary parts in Equation (77) results in the following:

$$ |\hat{\lambda}| \cos \omega_r \Delta t = \mathcal{R} = 1 - \frac{9}{2} |K| + \frac{27}{8} |K|^2 \qquad (79a) $$

$$ |\hat{\lambda}| \sin \omega_r \Delta t = \pm \mathcal{I} = \pm \frac{3}{2}(2 - 3|K|)|K|^{1/2} \qquad (79b) $$

Squaring and adding previous relations gives the amplification factor, while dividing (79b) by (79a) provides the discrete dispersion relation of the numerical scheme, i.e.

$$ |\hat{\lambda}| = \sqrt{\mathcal{R}^2 + \mathcal{I}^2} \qquad (80a) $$

$$ \tan \omega_r^{\pm} \Delta t = \pm \frac{\mathcal{I}}{\mathcal{R}} \qquad (80b) $$

Although the dispersion relation (80b) has two roots, in the frequency range of resolvable waves ($\omega_r \Delta t \in [0, \pi]$) only the root with a positive imaginary part is meaningful. Because of the change in sign of the function $\mathcal{I}$ when $|K| = 2/3$, the mode associated with the first root, $\omega_r^+ \Delta t$, remains in the plane of resolvable waves as far as $|K| \le 2/3$. On the other hand, the second mode associated with $\omega_r^- \Delta t$, will be valid when $|K| > 2/3$. It is then concluded that, even though the numerical scheme does not have any spurious computational mode in the frequency range of resolvable waves, it is preferable to remain in the region defined by the inequality $|K| \le 2/3$.
Numerical stability is studied now by looking at the amplification factor defined in (80a). The scheme produces non-amplifying solutions as far as $|\lambda| \le 1$, so it is straightforward to show that if

$$|K| \le \frac{8}{9} \qquad (81)$$

the RK4 compact finite volume scheme is stable. In the general case, the stability condition on the Courant number will be a function of $\alpha$, $\rho$ and $\kappa\Delta x$, but it can be shown that the most stringent situation occurs when $\rho \to \infty$. Replacing relations (58) and (60) in (78) provides the following stability condition in that limit:

$$C_r \le \frac{2\sqrt{2}}{3} \frac{(2 + \cos \kappa \Delta x)}{\sin \kappa \Delta x} \le 2\sqrt{\frac{2}{3}} \qquad (82)$$

The numerical scheme is thus stable if the Courant number is chosen to be less than 1.633 no matter what values of $\alpha$ and $\rho$ are used and for all resolvable wave numbers. Nevertheless, the stability limit can be moved further on by decreasing the parameter value $\rho$, i.e. by improving the spatial grid resolution for a fixed $h_0$. In the limit of very good grid resolution ($\rho \to 0$), the stability condition (81) takes the following form:

$$C_r \le \frac{2\sqrt{3}}{9} \frac{(2 + \cos \kappa \Delta x)}{\sin \kappa \Delta x} \sqrt{\frac{(1/3 + \alpha)(5 + \cos \kappa \Delta x)(13 - \cos \kappa \Delta x)}{\alpha(11 + \cos \kappa \Delta x)}} \le 1.603 \sqrt{\frac{(1/3 + \alpha)}{\alpha}} \qquad (83)$$

Therefore, when taking $\alpha = 1/15$ the method remains stable for all resolvable wave numbers, $\kappa\Delta x$, if the Courant number is less than 3.926 in the limit of very good spatial resolution. It is interesting to note that the high-order dispersion correction term introduces some kind of instability in the system because the stability limit is increased by reducing the weight of this term (i.e. $\alpha \to 0$). On the other hand, for finite values of $\rho$, it can be shown that the stability limit will be bounded by the two asymptotic values given in (82) and (83).

Even though stability is a necessary condition in order to ensure that the numerical scheme will not produce unbounded solutions in time, it does not say anything about the accuracy of the method in terms of dissipative and dispersive errors. For instance, Figure 6 shows numerical dissipation and phase speed errors for different Courant numbers and for a fixed value of $\rho = 0.2$ (keeping $\alpha = 1/15$). This particular $\rho$-value has been chosen in order to perform a comparison with the numerical analysis presented in papers by Woo and Liu [27] and Bradford and Sanders [29]. It is seen in Figure 6(a) that when $C_r = 1.0$ the present compact finite volume scheme is stable and nearly neutral (i.e. $|\lambda| \simeq 1.0$), thus introducing a clear improvement over the Petrov-Galerkin finite element model described in Reference [27] because with the same numerical parameters, amplifying modes were already present in the latter. In fact, those authors were only able to ensure stability when the Courant number was less than 0.5, which is a value more than three times smaller than the most stringent stability limit expressed in Equation (82). Similarly, stability of several linear finite volume schemes (Lax-Wendroff, Fromm and Warming-Beam) and a centred 2nd-order finite difference scheme used to solve a KdV-like system of equations was compared in Reference [29]. For an almost equivalent set of numerical parameters ($C_r = 1.0$ and $\rho = \Delta x/h_0 \simeq 0.18$), Bradford and Sanders showed that only the Fromm and the Warming-Beam methods were stable but they appear
<figure><img src="image_7.png" /><figcaption>Figure 6. Numerical dissipation and phase speed error of the 4th-order compact finite volume scheme for $\alpha=1/15$, $\rho=0.2$ at different Courant numbers: (—) $C_r=1.0$, (- -) $C_r=2.0$, (· · ·) $C_r=3.0$ and (-x-) $C_r=4.0$: (a) amplification factor; and (b) ratio between computed and physical phase speed.</figcaption></figure>

to be dissipative for mid to high-range spatial frequencies. Thus, it can be concluded from this particular comparison that the present numerical scheme has a wider stability region and is less damping than those studied in References [27, 29]. Moreover, using Maclaurin series to expand (80a) it can be shown that for small values of $\kappa\Delta x$,

$$ \lim_{\rho \to \infty} |\lambda| = 1 - \frac{C_r^6}{144} (\kappa \Delta x)^6 + O[(\kappa \Delta x)^8] \quad (84a) $$

$$ \lim_{\rho \to 0} |\lambda| = 1 - \frac{C_r^6}{144} \left[ \frac{\alpha}{(1/3 + \alpha)} \right]^3 (\kappa \Delta x)^6 + O[(\kappa \Delta x)^8] \quad (84b) $$

Therefore, the dissipative error for the present finite volume scheme is $O[(\kappa\Delta x)^6]$, while linear finite volume schemes studied by Bradford and Sanders [29] only shown to be $O[(\kappa\Delta x)^4]$ accurate.

Further information on the damping modes associated to our numerical model can be obtained from Figure 6(a) when looking at amplification factors associated to increasingly high Courant numbers. It is seen that for $C_r=2.0$ and 3.0 some mid-range frequency dissipation is present with a maximum located at $\kappa\Delta x \approx 0.662\pi$. Additional computations show that the numerical scheme is non-amplifying (stable) for $\rho=0.2$ as far as $C_r < 3.71$. However, mid-range frequencies begin to be heavily damped before reaching that limit (see Figure 6(a)).
As discussed in Reference [48], this dissipative behaviour is mostly due to the Runge–Kutta time stepping.

Phase speed errors obtained for $\rho=0.2$ and for different Courant numbers are presented in Figure 6(b). It is seen that phase speed estimates get worse when increasing the Courant number. Again, an unstable rapidly growing mode can be identified and coincides with the one observed for the amplification factor. It is also evident from this figure the passage from one root to the other in the numerical dispersion relation (80b) that takes place between Courant numbers $C_r=3.0$ and 4.0. As before and in the limit of small $\kappa\Delta x$, it can be shown that a Maclaurin expansion of the dispersion relation yields

$$ \lim_{\rho \to \infty} \frac{C_{fv}}{C_{th}} = 1 - \frac{1}{180} \left(1 + \frac{3}{2} C_r^4\right) (\kappa \Delta x)^4 + O[(\kappa \Delta x)^6] \quad (85a) $$

$$ \lim_{\rho \to 0} \frac{C_{fv}}{C_{th}} = 1 - \frac{1}{2880} \left( 11 + 24 \left[ \frac{\alpha}{(1/3 + \alpha)} \right]^2 C_r^4 \right) (\kappa \Delta x)^4 + O[(\kappa \Delta x)^6] \quad (85b) $$

Therefore, phase errors dominate over dissipative errors in the present compact finite volume scheme. Linear finite volume schemes studied by Bradford and Sanders [29] showed similar characteristics but with larger truncation errors. Figure 6 also indicates that in the present example, stable neutral solutions are obtained if $C_r \le 1.5$. It is important to note that phase speed will not be overpredicted in this range of Courant numbers.

In the general case, it is possible to construct stability and phase error contour plots in the $(C_r, \kappa\Delta x)$ plane for different values of the $\rho$-parameter. Before doing so we recall that, to the authors' knowledge, detailed analysis of the spectral resolution of previously published numerical methods used to solve equivalent sets of Boussinesq-type equations have not been reported. However, this property is paramount when dealing with physical problems where different wavelength scales are naturally present or can emerge as a consequence of nonlinear interaction. For example, recent extended versions of Boussinesq-type equations can theoretically match the Stokes linear dispersion relation up to $\sigma \approx 6$ when formulated in surface and horizontal velocity variables (e.g. References [21, 35]), and up to even higher values ($\sigma \approx 25$) if the vertical velocity is in addition retained as unknown (see Reference [6]). Nevertheless, if numerical schemes used to solve the associated systems of PDE are not carefully conceived to correctly reproduce phase speed and group velocity for all physical wavelengths, theoretical improvements concerning the extension of the $\sigma$ range may be unpractical. For instance, if the mixed-order strategy is used to discretize convective and dispersive terms, the asymptotic error analysis given in the previous section showed that the 0.5% phase speed error limit remains between $0.15\pi < (\kappa\Delta x)_{0.5\%} < 0.3\pi$ when $\omega_i\Delta t \to 0$. Therefore, the lower limit of this inequality, which occurs asymptotically in deep waters, is an indication of the numerical resolution of the scheme when dealing with extremely dispersive waves. This means that the discretized model will correctly describe linear phase speed for wave numbers $\kappa\Delta x$ being less than $0.15\pi$ and this value fixes the $\rho$-parameter that will ensure that all physical wavelengths will be well represented (because $\kappa\Delta x = \sigma\rho$). If the physical $\sigma$ range of the discretized Boussinesq-type equations is 6, then a value of $\rho = \Delta x/h_0 \le 0.08$ must be used in the mixed-order finite volume scheme. Similarly, if the physical range is moved further on to reach $\sigma \sim 25$ as it is the case in Reference [6], then the spectral resolution of the numerical scheme imposes the discretization parameter $\rho$ to be less than 0.02. It is evident from results presented in the previous subsection that, for group velocities, the ratio $\Delta x/h_0$ must
<figure><img src="image_8.png" /><figcaption>Figure 7. Contours of 0.5 and 5% dissipative errors for $\alpha = 1/15$ and $\rho = 0.1$ in the $(C_r, \kappa\Delta x)$ plane. $|\lambda| \le 0.95$ in region I (damping), and $|\lambda| \ge 1.05$ in region II (amplifying): (a) full 4th-order compact finite volume scheme; and (b) mixed 4th- and 2nd-order discretization strategy.</figcaption></figure>

be further reduced to achieve an accurate representation of this property. Moreover, for finite Courant number values spectral resolution will drop even more. This simple but important discussion allows us to illustrate the relevancy of the spectral resolution in terms of phase speed and group velocity of numerical schemes. In the framework of numerical resolution of Boussinesq-type equations, this important DRP property must be studied. More specifically, if extended versions of those equations, which may allow for the description of highly dispersive waves are to be discretized, the use of low-resolution explicit finite difference formulae can lead to overwhelming computational efforts if an accurate representation of group velocities is intended (even if high-order schemes are used). Indeed, it is widely acknowledged that such approximations have significantly worse spectral properties than equivalent compact formulae [38–40].

Contour plots for the dissipative errors in the $(C_r, \kappa\Delta x)$ plane are presented in Figure 7 for the full 4th-order scheme and the *mixed*-order approach. The $\rho$-parameter is fixed at 0.1, a value that allows both numerical approaches to describe all physical wavelengths because the linear dispersion relation associated to the PDE system match the theoretical Stokes relation up to $\sigma \sim 3$ when $\alpha = 1/15$ [41]. This figure shows that there is a wide region where both schemes remain stable and nearly neutral (white region). In region I, the numerical scheme starts to damp mid-range frequencies but it remains stable. Amplifying unstable modes begin to disturb the numerical solution when the thin region
<figure><img src="image_9.png" /><figcaption>Figure 8. Contours of 0.5 and 5% phase speed errors for $\alpha = 1/15$ and $\rho = 0.1$ in the $(C_r, \kappa\Delta x)$ plane: (a) full 4th-order compact finite volume scheme; and (b) mixed 4th- and 2nd-order discretization strategy.</figcaption></figure>

located between zones I and II is reached. Therefore, the upper edge of region I roughly represents the stability limit of the scheme. It is seen that the full 4th-order version will be stable for all frequencies if $C_r$ is taken to be less than 3.8, and that the *mixed*-order one will remain non-amplifying if $C_r$ is less than 4.1. Hence, the use of low-order explicit finite difference formulae to discretize dispersive terms results in a slight improvement in stability and similar conclusions were given in Reference [49]. However, for practical computations, spectral resolution in terms of phase speed, group velocity and numerical dissipation will impose more stringent conditions than stability and the Courant number will be in general chosen to be less than 1.5 in order to ensure a proper representation of all physical wavelengths.

As discussed before, spectral resolution and DRP properties of numerical schemes used to solve wave-type equations must be studied as emphasized in Reference [46]. A wide range of stability region is for sure a desired quality for any numerical solver because small unphysical perturbations will not grow unboundedly. However, the use of too high Courant numbers (even if they are below the stability limit) to integrate the equations in time can produce poor results in terms of phase and group velocities and this was already illustrated in Figure 6. Now we generalize this result plotting contours of phase speed and group velocity errors in the $(C_r, \kappa\Delta x)$ plane in Figures 8 and 9. It is seen that numerical errors in phase speeds and group velocities are rather insensitive to changes in the Courant number as long as the latter remains below a value of 1.5. For Courant numbers laying above that limit, the
<figure><img src="image_10.png" /><figcaption>Figure 9. Contours of 0.5 and 5% group velocity errors for $\alpha = 1/15$ and $\rho = 0.1$ in the $(C_r, \kappa\Delta x)$ plane: (a) full 4th-order compact finite volume scheme; and (b) mixed 4th- and 2nd-order discretization strategy.</figcaption></figure>

extent of the DRP region in terms of wave numbers $\kappa\Delta x$ is reduced as the Courant number increased and this is more evident for the full 4th-order scheme. On the other hand, the DRP region for a given Courant number is bigger when using the full 4th-order approach than when discretizing dispersive terms with 2nd-order explicit finite difference formulae. Indeed, if the Courant number is kept below a value of 1.5, the extent of the region where numerical errors in group velocities are less than 0.5% is almost two times bigger in the full 4th-order version than in the *mixed* one. A similar trend is observed in Figure 8 for the comparative analysis of numerical representation of phase speeds. This clearly shows the inconvenience of using *mixed*-order strategies to simplify the numerical treatment of high-order terms because it results in an important overall reduction of the spectral resolution in this particular finite volume scheme. Therefore, significantly more nodal points per wavelength will be needed to properly describe wave kinematics; moreover, if the numerical scheme must also be able to deal with extremely dispersive/short waves, this critical limitation may also obscure theoretical improvements in the $\sigma$-range of application of Boussinesq-type equations.

Finally, it can be concluded from the linear analysis that keeping the Courant number below 1.5 and using $\rho = 0.1$ for numerical computations will ensure that the compact finite volume scheme used to discretize the set of PDE derived here is stable, non-damping and with a maximized DRP region.
## 5. NUMERICAL EXPERIMENTS AND DISCRETIZATION ERRORS

In this last section we study the stability of the proposed numerical scheme using an heuristical approach that allows us to explore how this property is affected by nonlinearity. Besides, the discretization error is estimated numerically by comparing analytical and computed results using the solitary wave solution presented in Section 2. Discretization error is related to the order of approximation of the scheme and is different (although somehow related) from the spectral properties investigated in the previous section. Here, we will instead determine how truncation errors, which belong to numerical approximation of partial derivatives, are reduced when the spatial and temporal grid is refined.

### 5.1. Discretization errors

Without dispersion correction terms (i.e. $\alpha=0$) and for a flat bottom, Equations (18)–(19) have the solitary wave solution expressed by relations (22)–(24). This interesting property is used here to estimate discretization errors by comparing the numerical and the exact solution for this particular PDE system. In order to do that, a solitary wave of relative amplitude $a/h_0=0.2$ is propagated during two characteristic wave periods over a still water depth of $h_0=1$ m. The characteristic period is estimated from the characteristic length of the solitary wave, $L$, following the definition used in Reference [27]. Hence, $L$ is the distance over which $\eta/h_0 \ge 0.001$, where $\eta$ is the location of the free surface. The equivalent solitary wave period is then defined as $T=L/C$ and can be computed invoking the definition for the phase speed of the soliton given in Equation (24). The average numerical error is estimated using the $L_1$ norm and reads

$$L_{\epsilon} = \frac{\sum_{j=0}^{j=L} \|h_j^n - h_j^e\|_1}{\sum_{j=0}^{j=L} \|h_j^e\|_1} \qquad (86)$$

where $h^e$ and $h^n$ are, respectively, the exact and the numerically computed solution, $I$ is the total number of finite volumes used to discretize the spatial domain, and water depth values are evaluated at cell faces of each control volume (j-nodal points). Equations are integrated using the solitary wave solution centred at $x=X_0$ as initial condition. Computations are carried out during two wave periods and results compared with the analytical soliton centred at position $x=X_0+2L$. All numerical tests are performed using a fixed-point-iteration tolerance of $\tau=10^{-4}$ to solve (47). Spatial discretization error is studied using a small time step of $\Delta t=0.002$s in order to ensure that most of the errors belong to the spatial grid resolution used to perform computations. Estimated convergence errors obtained by increasing the number of finite volumes per wavelength (i.e. decreasing $\rho$) are presented in Figure 10(a) where it can be seen that, except for the last two computed errors using the smallest $\rho$, the ideal slope for a 4th-order spatial discretization is fairly well recovered. It is believed that for the two last computed values, round-off errors or some influence of the time step (because it is not small enough for the prescribed $\rho$ value) are responsible for the slight divergence from the ideal slope. Nevertheless, it is widely confirmed that the finite volume discretization of the system of PDE is $O[(\Delta x)^4]$ accurate.

An equivalent numerical test is performed on discretization errors due to time integration for a fixed $\rho=0.1$ and monotonically decreasing the Courant number. Associated numerical errors
<figure><img src="image_11.png" /><figcaption>Figure 10. Discretization errors for a propagating solitary wave. (—) Ideal 4th-order slope; (o) computed errors when $a/h_0 = 0.2$; (Δ) computed errors when $a/h_0 = 0.6$: (a) convergence error in $\Delta x$ using a fixed $\Delta t = 0.002$ s; and (b) convergence error in $\Delta t$ using a fixed $\rho = 0.1$.</figcaption></figure>

are presented in Figure 10(b) and compared against the theoretical 4th-order slope. Again, it is fairly well confirmed that the convergence error is $O[(\Delta t)^4]$, except for the last point where, probably for the same reasons stated before, the computed slope slightly diverges from the theoretical one. Therefore, the proposed RK4 compact finite volume scheme is $O[(\Delta x)^4, (\Delta t)^4]$ accurate.

In addition, convergence errors obtained for a highly nonlinear propagating solitary wave ($a/h_0 = 0.6$) are also plotted in Figure 10. It is seen that the rate of convergence for the discretization error in $\Delta x$ is slightly worse than the one observed in the moderately nonlinear case investigated before. The latter is not surprising since strong nonlinearities may affect the model performance. Nevertheless, the convergence rate is still well above $O[(\Delta x)^3]$ for this challenging test and the computed slope only starts to diverge from the theoretical one when $\Delta x/h_0 \le 0.4$. This trend, which has been also noticed in the previous example, is probably associated to the fixed time step value used to estimate errors since in this highly nonlinear case the equivalent wavelength (and the associated wave period) is much smaller than the previous one. On the contrary, the convergence rate for temporal discretization appears to be rather insensitive to changes in nonlinearity since the ideal slope is well recovered. Thus, the four-stage Runge-Kutta scheme used to integrate the system in time has an excellent performance even for highly nonlinear cases.
<figure><img src="image_12.png" /><figcaption>Figure 11. Strongly nonlinear solitary wave propagating over 50 equivalent wavelengths where $a/h_0 = 0.6$, $L = 14.7$ m, $T = 3.7$ s. Exact solution (−) and numerical computations performed with $C_r = 1.0$, $\rho = 0.1$ (−−), $\rho = 0.2$ (·), and $\rho = 0.5$ (−·−): (a) overview of the whole domain and initial and final locations of the solitary wave; and (b) close inspection of numerical results.</figcaption></figure>

This numerical test also validates the chosen deferred-correction approach used to preserve the tridiagonal structure of the matrix that links the auxiliary variable $q$ with the depth-averaged velocity $u$ in matrix equation (47). Moreover, the convergence rate was rather insensitive to changes in the fixed-point-iteration tolerance as long as $\tau \le 10^{-4}$. Indeed, $\tau = 10^{-4}$ proved to be a sufficient choice in all computations.

Another important property that any numerical scheme used to solve Boussinesq-type equations must fulfil concerns its ability to preserve wave form and phase speed in time. These qualities have been already investigated from an analytical point of view through the linear spectral analysis given in Section 4. However, it is important to test these properties in the case of highly nonlinear waves propagating over long distances. This can be achieved using again the closed-form solitary wave solution of the original Serre equations. Similar tests have been reported in References [24, 26, 27] for moderately nonlinear solitary waves. For instance, Figure 11 shows computed results for the propagation of a solitary wave of relative amplitude $a/h_0 = 0.6$ over 50 equivalent wavelengths using $C_r = 1$ and $\rho = 0.1$. A perfect agreement between the analytical and numerical solution is noticed. It is remarkable how this result confirms the analysis given in previous section where it was shown that these set of numerical parameters may produce stable non-damping solutions for all physical wavelengths. In addition, when decreasing the number of nodal points per wavelength some phase speed and
amplitude errors can be noticed. Nevertheless, even for $ρ=0.2$ the agreement with the analytical solution is almost perfect. Although an objective comparison with previously published models used to solve equivalent set of equations cannot be performed, we recall that a similar grid resolution (~75 elements per wavelength but with $C_r \sim 0.5$) was used in Reference [27] on the basis of a Taylor-Galerkin finite element model and some phase speed errors were already noticed when propagating a moderately nonlinear ($a/h_0 = 0.3$) solitary wave over a shorter distance (~25 wavelengths). A similar trend was also reported in References [24, 26].

Finally, when increasing the $ρ$-value further to reach 0.5 (~29 nodal points per wavelength), bigger discrepancies show up. Nevertheless, the relative error in phase speed over this long distance appears to be less than 2%.

## 5.2. Nonlinear stability

The way by which nonlinearity affects the stability of the proposed scheme is investigated now performing numerical experiments where we estimate the maximum Courant number for which a propagating solitary wave remains stable. Computations are carried out over six equivalent wave periods for different relative amplitudes $a/h_0$ and for monotonically increasing Courant numbers until numerical instability shows up. The influence of the correction dispersion term is also investigated because it was analytically demonstrated in Section 4 that this term may produce some destabilizing effects in the system (see inequality (83)). Therefore, three different $α$-values are investigated, namely, $α = 1/15$, $α = 1/30$ and $α = 0$. It is important to recall that it is expected that when $α = 0$ important numerical (and mathematical) properties of the system may change as it was illustrated in Figure 3 where a totally different behaviour was noticed for this $α$-value. Numerical estimation of the stability limit in terms of the Courant number is given in Figure 12. Theoretical values obtained from linear analysis is also included in the figure (when $a/h_0 \to 0$) for $α = 1/15$ and $α = 1/30$. The limiting Courant number predicted by linear theory when $α = 0$ and $ρ = 0.1$ is roughly $C_r = 16$. Such a high value was not recovered numerically probably because the stability analysis performed in Section 4 may only be valid for finite values of $α$ as discussed before. Nonetheless, for the two finite values of $α$ investigated here, numerical experiments confirm in a satisfactory way the linear stability analysis as depicted in Figure 12. Moreover, when the dispersion correction term is switched off the stability region is considerably increased as it was expected from the numerical analysis of the scheme. It is also noticed that for finite $α$-values, nonlinearity does affect the stability of the scheme because the maximum stable Courant number is reduced. Nevertheless, numerical experiments show that the RK4 compact finite volume scheme remains stable, even for extremely nonlinear waves, as long as the Courant number is kept below a value of 1.8 when $α \le 1/15$. Undoubtedly, this result constitutes a major improvement in stability when compared to previously published numerical solvers for Boussinesq-type equations.

It was noticed that for the finite $α$-values investigated, numerical instability was mainly triggered because the resulting matrix of Equation (47) appeared to be ill conditioned in the presence of high nonlinearity. On the contrary, when dispersion correction terms were set to zero, numerical instability appeared to be mainly due to intrinsic properties associated to the 4th-order Runge-Kutta time stepping, and this may partially explain why the region of stability was only slightly affected by changes in $a/h_0$ when $α = 0$. We believe that this misleading dependence on the chosen $α$ value must be further investigated on a mathematical
<figure><img src="image_13.png" /><figcaption>Figure 12. Stability limits estimated numerically using $ρ=0.1$ for increasing non-linear solitary waves and different values for the dispersion correction parameter: $(-) \alpha=1/15$, $(-−) \alpha=1/30$, and $(-·−) \alpha=0$.</figcaption></figure>

basis because it was demonstrated that the dispersion correction strategy used here may lead to an important reduction of the stability region. Moreover, it was rather surprising to confirm numerically that the stability limit may reach such high Courant numbers when the dispersion correction term was switched off.

## 6. CONCLUSIONS

A novel approach to numerically handle a set of fully nonlinear and weakly dispersive Boussinesq-type equation was presented and deeply investigated using numerical analysis. The chosen set of extended Serre equations could be written in a weak quasi-conservative form, which makes the use of finite volume methods very attractive. Indeed, recently developed compact strategies for cell-face reconstruction in the framework of finite volume methods were borrowed from CFD and adapted to this particular set of PDE. High-order accuracy was achieved in spectral space and numerically estimated discretization errors confirmed that the newly developed scheme is fully 4th-order in space and time while computational efforts were kept at a very reasonable level. The latter constitutes an important feature of the present scheme because extensions in two horizontal dimensions (2D-H), where relatively large scale problems (~ km) are to be discretized, may be achieved at affordable computational costs.
In fact, the use of the present compact approach allowed us to reduce CPU time by more
than three times when compared to our previous explicit 4th-order finite volume model [30].
Moreover, we believe that the particular choice of this compact finite volume method will
allow us to extend the model in 2D-H straightforwardly taking advantage of developments
already available in several space dimensions and structured grids (e.g. References [39, 50]).

On the other hand, the use of linear spectral analysis provides important information con-
cerning the performance of the newly developed scheme. For instance, we have good reasons
to believe that the widely used strategy of discretizing convective and dispersive terms using
*mixed*-order finite difference formulae may lead to phase and group velocity errors which
are difficult to control. This drawback can be specially critical for computations dealing with
extremely dispersive waves since a good description of the whole range of characteristic wave-
lengths will certainly impose the use of a very fine grid resolution. Therefore, some recent
theoretical improvements for Boussinesq-type equations may be unpractical for real world
applications. The latter has a strong analogy with some difficulties encountered in the field
of CFD where, for direct numerical simulations (DNS) of Navier-Stokes equations, the range
of spatial lengths that a particular numerical scheme should be able to correctly represent is
broadened towards smaller scales for increasingly high Reynolds number (thus requiring fur-
ther refining of the grid size). This question has not really been tackled yet in the framework
of coastal engineering. However, if numerical modelling of near-shore hydrodynamics is to
be handled with help of Boussinesq-type equations, efficient, accurate and stable larger scale
models are to be developed.

Additionally, qualitative comparisons with previously published numerical schemes used
to solve equivalent set of equations were performed in Sections 4 and 5. Linear spectral
analysis and some numerical experiments showed that the present RK4 compact finite volume
model possesses a wide stability region, good spectral resolution and an excellent numerical
behaviour even when dealing with high nonlinearities. For instance, stability is only ensured
in some of the previously published numerical strategies used to solve equivalent set of
equations for Courant numbers which should remain below 0.5 (e.g. References [22, 27, 51]).
On the contrary, the present finite volume scheme has shown to be stable (i.e. for all wave
numbers and spatial grid resolutions) if Courant number is fixed to a value below 1.6, thus
introducing fundamental improvements over previous ones. Moreover, this limit could be
moved further on if grid resolution was improved and it was found in Section 5 that, even
in the case of the propagation of an extremely nonlinear solitary wave, solution remained
stable as far as $C_r \le 1.8$ for a grid resolution of $\rho = 0.1$. Indeed, it was found analytically in
Section 4 that the dispersion correction term as given in References [7, 35] introduces some
kind of numerical instability in the system. This undesirable situation was also confirmed by
numerical experiments.

Ongoing developments of the present model concern its extension into the surf zone by
including extra breaking terms. For instance, wave breaking was recently incorporated using
the parametrization given in Reference [52]. Similarly, the numerical treatment of appropri-
ate boundary conditions and the validation of the model for wave propagation over uneven
bottoms will be addressed in a companion paper.

Finally, it is also worth to point out that in more complex wave propagation problems (i.e.
over uneven bottoms, periodic waves and added breaking effects), high-frequency waves can
arise as a consequence of nonlinear interaction. In this context, the use of high-order filters
to damp out spurious wave components is of common usage. The spectral analysis performed
here gives insight of how this can be achieved without affecting the range of resolved physical wavelengths.

APPENDIX A

Series expansion of $u^*$ and $w^*$ are found using Equations (2), (3) and boundary condition (6). As a first step assume the following expansion for both unknowns:

$$u^*(x^*, z^*, t^*) = \sum_{n=0}^{\infty} \frac{(z^* - \xi^*(x^*))^n}{n!} f_n(x^*, t^*) \quad (A1)$$

$$w^*(x^*, z^*, t^*) = \sum_{n=0}^{\infty} \frac{(z^* - \xi^*(x^*))^n}{n!} g_n(x^*, t^*) \quad (A2)$$

where $f_0$ and $g_0$ are, respectively, the bottom horizontal and vertical velocities. Boundary condition (6) imposes that,

$$g_0 = \xi_x^* f_0 \quad (A3)$$

Informations on $f_n$ and $g_n$ are obtained by invoking continuity and irrotationality at each order. At order $n$, Equations (2) and (3) give

$$\begin{aligned} \xi_x^* f_{n+1} - g_{n+1} &= (f_n)_x \\ f_{n+1} + \sigma^2 \xi_x^* g_{n+1} &= \sigma^2 (g_n)_x \end{aligned} \quad (A4)$$

which at the lowest order writes,

$$\begin{aligned} \xi_x^* f_1 - g_1 &= f_{0x} \\ f_1 + \sigma^2 \xi_x^* g_1 &= \sigma^2 g_{0x} \end{aligned} \quad (A5)$$

This linear set of equations has the solution

$$f_1 = \frac{\sigma^2 \xi_x^* f_{0x} + \sigma^2 (f_0 \xi_x^*)_x}{1 + \sigma^2 (\xi_x^*)^2} \quad (A6)$$

$$g_1 = \frac{\sigma^2 \xi_x^* (f_0 \xi_x^*)_x - f_{0x}}{1 + \sigma^2 (\xi_x^*)^2} \quad (A7)$$

Thus, if $\sigma^2(\xi_x^*)^2 \ll 1$, the following expansions in terms of $\sigma^2$ hold

$$f_1 = \sigma^2[\xi_x^* f_{0x} + (f_0 \xi_x^*)_x] + O(\sigma^4) \quad (A8)$$

$$g_1 = -f_{0x} + \sigma^2[\xi_x^* (f_0 \xi_x^*)_x + f_{0x} (\xi_x^*)^2] + O(\sigma^4) \quad (A9)$$

and a justification to (14) and (15) is provided.
Similarly, this functional form for horizontal and vertical velocities allows us to express the vertical flow acceleration function given by Equation (9) as

$$ \Gamma^*(x^*, z^*, t^*) = \left[ \frac{(z^* - \xi^*)}{(\varepsilon\eta^* - \xi^*)} \mathcal{P}^*(x^*, t^*) + \mathcal{Q}^*(x^*, t^*) \right] + O(\sigma^2) \quad (A10) $$

with,

$$ \varepsilon\sigma^2\mathcal{P}^*(x^*, t^*) = \frac{\mathcal{P}}{g}, \quad \varepsilon\sigma^2\mathcal{Q}^*(x^*, t^*) = \frac{\mathcal{Q}}{g} $$

using dimensional functions, $\mathcal{P}$ and $\mathcal{Q}$, which were introduced in relations (20) and (21).

## ACKNOWLEDGEMENTS

The first author would like to thank the French Foreign Office, CROUS Grenoble, the Chilean Research Council (CONICYT) and Pontificia Universidad Católica de Chile for their financial support during his PhD research. Cristián Escauriaza and anonymous reviewers are also acknowledged for their useful suggestions and comments.

## REFERENCES

1. Boussinesq J. 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. *Journal de Mathematiques Pures et Appliquées* 1872; 2:55-108.

2. Korteweg DJ, de Vries D. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves. *Philosophical Magazine* 1895; 39:422-443.

3. Serre F. Contribution à l'étude des écoulements permanents et variables dans les canaux. *Houille Blanche* 1953; 8:374-388.

4. Peregrine DH. Long waves on beaches. *Journal of Fluid Mechanics* 1967; 27(4):815-827.

5. Lynett P, Liu L-F. A two-layer approach to wave modelling. *Proceedings of the Royal Society of London A* 2004; 460:2637-2669.

6. Madsen PA, Bingham HB, Schäffer HA. Boussinesq formulation for fully non-linear and extremely dispersive water waves. *Proceedings of the Royal Society of London A* 2003; 459:1075-1104.

7. Madsen PA, Murray R, Sørensen OR. A new form of the Boussinesq equations with improved linear dispersion characteristics. *Coastal Engineering* 1991; 15:371-388.

8. Meftah K, Sergent P, Gomi P. Linear analysis of a new type of extended Boussinesq model. *Coastal Engineering* 2004; 51(2):185-206.

9. Nwogu O. Alternative form of Boussinesq equations for nearshore wave propagation. *Journal of Waterway, Port, Coastal, and Ocean Engineering* 1993; 119(6):618-638.

10. Witting JM. A unified model for the evolution of nonlinear water waves. *Journal of Computational Physics* 1984; 56:203-236.

11. Kennedy AB, Chen Q, Kirby JT, Dalrymple RA. Boussinesq modelling of wave transformation, breaking and runup. I: 1D. *Journal of Waterway, Port, Coastal, and Ocean Engineering* 2000; 126(1):39-48.

12. Schäffer HA, Madsen PA, Deigaard RA. A Boussinesq model for waves breaking in shallow water. *Coastal Engineering* 1993; 20:185-202.

13. Veeramony J, Svendsen IA. The flow in surf-zone waves. *Coastal Engineering* 2000; 39:93-122.

14. Zelt JA. The run-up of nonbreaking and breaking solitary waves. *Coastal Engineering* 1991; 15:205-246.

15. Seabra-Santos FJ, Renouard DP, Temperville AM. Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle. *Journal of Fluid Mechanics* 1987; 176:117-134.

16. Seabra-Santos FJ, Renouard DP, Temperville AM. Etude théorique et expérimentale des domaines de validité des théories d'évolution des ondes en eau peu profonde. *Annales de Geophysique* 1988; 6:671-680.

17. Wei G, Kirby JT, Grilli ST, Subramanya R. A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves. *Journal of Fluid Mechanics* 1995; 294:71-92.

18. Abbott MB, McCowan AD, Warren IR. Accuracy of short-wave numerical model. *Journal of Hydraulic Research* 1984; 110(10):1287-1301.
19. Antunes do Carmo JS, Seabra-Santos FJ, Almeida AB. Numerical solution of the generalized Serre equations with the McCormack finite-difference scheme. *International Journal for Numerical Methods in Fluids* 1993; **16**:725-738.

20. Fuhram DR, Bingham HB. Numerical solution of fully non-linear and highly dispersive Boussinesq equations in two horizontal dimensions. *International Journal for Numerical Methods in Fluids* 2004; **44**:231-255.

21. Gobbi MF, Kirby JT, Wei G. A fully non-linear Boussinesq model for surface waves. Part 2. Extension to $O(kh)^4$. *Journal of Fluid Mechanics* 2000; **405**:181-210.

22. Shi F, Dalrymple RA, Kirby JT, Kennedy A. A fully nonlinear Boussinesq model in generalized curvilinear coodinates. *Coastal Engineering* 2001; **42**:337-358.

23. Su CH, Mirie RM. On head-on collisions between two solitary waves. *Journal of Fluid Mechanics* 1980; **98**: 509-525.

24. Wei G, Kirby JT. Time-dependent numerical code for extended Boussinesq equations. *Journal of Waterway, Port, Coastal, and Ocean Engineering* 1995; **121**(5):251-261.

25. Antunes do Carmo JS, Seabra-Santos FJ, Barthélemy E. Surface waves propagation in shallow water: a finite element model. *International Journal for Numerical Methods in Fluids* 1993; **16**:447-459.

26. Walkley M, Berzins M. A finite element method for the one-dimensional extended Boussinesq equations. *International Journal for Numerical Methods in Fluids* 1999; **29**:143-157.

27. Woo S-B, Liu PL-F. A Petrov-Galerkin finite element model for one-dimensional fully non-linear and weakly dispersive wave propagation. *International Journal for Numerical Methods in Fluids* 2001; **37**:541-575.

28. Woo S-B, Liu PL-F. Finite-element model for modified Boussinesq equations. II: applications to nonlinear harbor oscillations. *Journal of Waterway, Port, Coastal, and Ocean Engineering* 2004; **130**(1):17-28.

29. Bradford SF, Sanders BF. Finite-volume models for unidirectional, nonlinear, dispersive waves. *Journal of Waterway, Port, Coastal, and Ocean Engineering* 2002; **128**(4):173-182.

30. Cienfuegos R, Barthélemy E, Bonneton P. Résolution numérique en volumes finis d'un système d'équations de Serre étendu. *Révue Française de Génie Civil* 2005; **9**(7-8):889-902.

31. Erduran KS, Ilic S, Kutija V. Hybrid finite-volume finite-difference scheme for the solution of Boussinesq equations. *International Journal for Numerical Methods in Fluids* 2005; **49**(11):1213-1232.

32. Stansby PK. Solitary wave run up and overtopping by a semi-implicit finite-volume shallow-water Boussinesq model. *Journal of Hydraulic Research* 2003; **41**(6):639-647.

33. Toro EF. *Shock-Capturing Methods for Free-Surface Flows* (1st edn). Wiley: New York, 2001.

34. Bernetti R, Toro EF, Brocchini M. An operator-splitting method for long waves. *Long Waves Symposium*, vol. 1, Tessaloniki, Greece, 2003; 49-56.

35. Madsen PA, Schäffer HA. Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis. *Philosophical Transactions of the Royal Society of London A* 1998; **356**:3123-3184.

36. Seabra-Santos FJ. Contribution à l'étude des ondes de gravité bidimensionnelles en eau peu profonde. *PhD Thesis*, Institut National Polytechnique de Grenoble, France, 1985.

37. Whitham GB. *Linear and Nonlinear Waves* (1st edn). Wiley Inter-Science: New York, 1974.

38. Kobayashi MH. On a class of Padé finite volume methods. *Journal of Computational Physics* 1999; **156**: 137-180.

39. Lacor C, Smirnov S, Baelmans M. A finite volume formulation of compact central schemes on arbitrary structured grids. *Journal of Computational Physics* 2004; **198**:535-566.

40. Lele SK. Compact finite difference schemes with spectral-like resolution. *Journal of Computational Physics* 1992; **103**:16-42.

41. Barthélemy E. Nonlinear shallow water theories for coastal waves. *Surveys in Geophysics* 2004; **25**(3-4): 315-337.

42. Dingemans MW. *Water Wave Propagation Over Uneven Bottoms* (1st edn). World Scientific Publication: Singapore, 1997.

43. Guizien K, Barthélemy E. Accuracy of solitary wave generation by a piston wave maker. *Journal of Hydraulic Research* 2002; **40**(3):321-331.

44. Durran DR. *Numerical Methods for Wave Equations in Geophysical Fluid Dynamics* (1st edn). Springer: Berlin, 1999.

45. Ferziger JH, Perić M. *Computational Methods for Fluid Dynamics* (3rd edn). Springer: Berlin, 2002.

46. Sengupta TK. *Fundamentals of Computational Fluid Dynamics* (1st edn). Universities Press (India) Private Limited, 2004.

47. Press WH, Flannery BP, Teukolsky SA, Vetterling WT. *Numerical Recipies in Fortran* (2nd edn). Cambridge University Press: Cambridge, U.K., 1992.

48. Hu FQ, Hussaini MY, Manthey JL. Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics. *Journal of Computational Physics* 1996; **124**:177-191.

49. Fuhram DR, Bingham HB, Madsen PA, Thomsen PG. Linear and non-linear stability analysis for finite difference discretizations of high-order Boussinesq equations. *International Journal for Numerical Methods in Fluids* 2004; **45**:751-773.
50. Shukla RK, Zhong X. Derivation of high-order compact finite difference schemes for non-uniform grid using polynomial interpolation. *Journal of Computational Physics* 2005; **204**:404–429.

51. Gobbi MF, Kirby JT. Wave evolution over submerged sills: tests of a high-order Boussinesq model. *Coastal Engineering* 1999; **37**:57–96.

52. Cienfuegos R, Barthélemy E, Bonneton P. A new wave-breaking parametrization for Boussinesq-type equations. *5th International Symposium on Ocean Wave Measurement and Analysis*, Madrid, Spain, 2005 (CD-ROM).