# A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: boundary conditions and validation

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 II: boundary conditions and validation. International Journal for Numerical Methods in Fluids, 2007, 53 (9), pp.1423-1455. 10.1002/fld.1359. hal-00268892

HAL Id: hal-00268892

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

Submitted on 23 Sep 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 II: Boundary conditions and validation

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 de Chile, Chile

This paper supplements the validation of the fourth-order compact finite volume Boussinesq-type model presented by Cienfuegos et al. (Int. J. Numer. Meth. Fluids 2006, in press). We discuss several issues related to the application of the model for realistic wave propagation problems where boundary conditions and uneven bathymetries must be considered. We implement a moving shoreline boundary condition following the lines given by Lynett et al. (Coastal Eng. 2002; 46:89–107), while an absorbing-generating seaward boundary and an impermeable vertical wall boundary are approximated using a characteristic decomposition of the Serre equations. Using several benchmark tests, both numerical and experimental, we show that the new finite volume model is able to correctly describe nonlinear wave processes from shallow waters and up to wavelengths which correspond to the theoretical deep water limit. The results compare favourably with those reported using former fully nonlinear and weakly dispersive Boussinesq-type solvers even when time integration is conducted with Courant numbers greater than 1.0. Furthermore, excellent nonlinear performance is observed when numerical computations are compared with several experimental tests on solitary waves shoaling over planar beaches up to breaking. A preliminary test including the wave-breaking parameterization described by Cienfuegos (Fifth International Symposium on Ocean Wave Measurement Analysis, Madrid, Spain, 2005) shows that the Boussinesq model can be extended to deal with surf zone waves. Finally, practical aspects related to the application of a high-order implicit filter as given by Gaitonde et al. (Int. J. Numer. Methods Engng 1999; 45:1849–1869) to damp out unphysical wavelengths, and the numerical robustness of the finite volume scheme are also discussed.

<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 de Chile, Chile.</footnote>

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

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

<footnote>Contract/grant sponsor: CROUS Grenoble</footnote>

<footnote>Contract/grant sponsor: Chilean Research Council (CONICYT)</footnote>

<footnote>Contract/grant sponsor: Pontificia Universidad Católica de Chile</footnote>
# 1. INTRODUCTION

In the companion paper [1], analytical properties of a new fourth-order compact finite volume scheme used to integrate Boussinesq-type equations were investigated. In particular, linear stability, spectral resolution and nonlinear performance of the scheme were analysed in the framework of wave propagation over flat bottoms. In the present work, practical considerations intended to apply model equations to realistic problems of wave propagation over uneven bathymetries are considered. Experimental test cases will provide additional insight concerning model performance when used to describe complex wave processes.

The set of mass and momentum conservation equations has already been presented in Reference [1]. Here we only recall the weak quasi-conservative form that has been discretized with the compact finite volume method. This particular fully nonlinear and weakly dispersive set of Boussinesq-type equations reads,

$$h_t + F_x = 0 \qquad (1)$$

$$q_t + G_x = S \qquad (2)$$

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

where the different flux functions, source term and dependent variable are defined, respectively, as

$$F = hu \qquad (4)$$

$$G = qu + g(h + \xi) - \frac{1}{2}u^2(1 + \xi_x^2) + (\xi_x u - \frac{1}{2}hu_x)hu_x - \alpha\xi^2[u_x^2 + uu_{xx} + g(h + \xi)_{xx}] \qquad (5)$$

$$S = -2\alpha\xi\xi_x[uu_x + g(h + \xi)_x]_x \qquad (6)$$

$$r = (h_x + \xi_x)\xi_x + \frac{1}{2}h\xi_{xx} \qquad (7)$$

Here h is the total water depth, u is the depth-averaged horizontal velocity of the fluid, z = ξ is the vertical coordinate of the bottom, g is gravitational acceleration, and α is an adjustable parameter associated with the dispersion correction term as given in References [2, 3]. The latter provides an extension of this shallow water set to deeper water wave propagation problems. The cartesian system of coordinates and a sketch of the problem are presented in Figure 1.

In the companion paper [1], details on the application of a high-order compact finite volume scheme to this set of equations have been given and thoroughly analysed for wave propagation over horizontal bottoms. For the reader's convenience, the main features of the numerical scheme are summarized here: (i) cell-face values are reconstructed from calculated cell-averaged ones using a compact strategy developed in References [4, 5]; (ii) high-order dispersive terms are discretized using compact finite difference formulae [6]; (iii) a standard four-stage Runge-Kutta method is used to integrate the discrete system in time; (iv) fourth-order accuracy in both space and time
<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 (here $\xi$ has a negative value), and $h = \eta - \xi$ is the total water depth.</figcaption></figure>

is achieved; and (v) stable non-damping solutions are obtained if the Courant number is taken to be less than 2.0 when using a grid size $\Delta x = 0.1 h_0$ ($h_0$ being a representative value for the still water depth).

In the following, we will pursue the validation of the proposed model using both numerical tests and experimental measurements for wave propagation over uneven bathymetries. The important question of conceiving and discretizing suitable boundary conditions for this system of equations will also be addressed. We are specifically interested in developing boundary conditions able to accurately describe several nonlinear wave processes taking place in a nearshore context.

The present paper is organized as follows: in Section 2, dispersive properties of the model are studied, while Section 3 deals with the numerical implementation of different types of boundary conditions. In Section 4 a high-order filtering technique is discussed and applied to damp out non-physical wavelengths not resolved by the model, and in Section 5 a complete validation of the proposed scheme is presented using both analytical and experimental examples.

## 2. OPTIMUM VALUE FOR THE DISPERSION CORRECTION PARAMETER

Improvements in dispersive characteristics embedded in the Boussinesq-type model can be achieved by adjusting the $\alpha$-parameter associated with the high-order dispersion correction term. From an analytical point of view it was shown by Barthélemy [7] that a Padé (2, 2) approximation for the exact linear Stokes dispersion relation is obtained by choosing $\alpha = \frac{1}{15}$. Therefore, this extended system of Serre equations can adequately match the Stokes relation up to $\sigma = \kappa h_0 \approx \pi$, where $\kappa$ is a characteristic wave number and $h_0$ is a typical water depth. Comparisons between dispersion characteristics for the partial differential equation (PDE) set and theoretical Stokes phase speed and group velocity are presented in Figure 2 for different values of $\alpha$.

Similarly, a Stokes expansion analysis was performed in the same Reference [7] to investigate the nonlinear properties of this particular system and this will not be repeated here. Nevertheless, bearing in mind that linear wave propagation properties of the set of PDE will theoretically allow for the description of wavelengths up to, say $\sigma_m = \kappa h_0 = 3.0$, one can try to find an optimum value for $\alpha$ by minimizing the joint root mean square (RMS) error of the phase and group velocity estimates. Therefore, instead of choosing the $\alpha$-value, which produces the exact Padé (2, 2) Stokes expansion, we will tune this parameter in order to minimize phase and group velocity errors over
<figure><img src="image_3.png" /><figcaption>Figure 2. Comparisons between linear dispersion characteristics of the model and exact Stokes relations for different values of $\alpha$: (-) $\alpha=0$; (---) $\alpha=1/10$; (·−) $\alpha=1/15$; (··) $\alpha=1/20$; (−△−) $\alpha=1/25$: (a) phase speeds; and (b) group velocities.</figcaption></figure>

the dispersive range $0 \le \sigma \le 3$. This can be done by introducing the following quadratic form for the joint error:

$$Err = \sqrt{\frac{\int_{0}^{\sigma_m} (C_{\text{th}} - C_{\text{Stokes}})^2 d\sigma}{\int_{0}^{\sigma_m} C_{\text{Stokes}}^2 d\sigma} + \frac{\int_{0}^{\sigma_m} (C_{g \text{th}} - C_{g \text{Stokes}})^2 d\sigma}{\int_{0}^{\sigma_m} C_{g \text{Stokes}}^2 d\sigma}} \quad (8)$$

where $C_{\text{th}}$ and $C_{\text{Stokes}}$ represent, respectively, linear phase speeds associated with the PDE system and the Stokes theory, while $C_{g \text{th}}$ and $C_{g \text{Stokes}}$ correspond to group velocities. They read,

$$\frac{C_{\text{th}}}{C_0} = \sqrt{\frac{1 + \alpha\sigma^2}{1 + (1/3 + \alpha)\sigma^2}} \quad (9)$$

$$\frac{C_{\text{Stokes}}}{C_0} = \sqrt{\frac{\tanh \sigma}{\sigma}} \quad (10)$$

$$\frac{C_{g \text{th}}}{C_0} = \frac{\alpha(1/3 + \alpha)\sigma^4 + 2\alpha\sigma^2 + 1}{\sqrt{[(1/3 + \alpha)\sigma^2 + 1]^3(1 + \alpha\sigma^2)}} \quad (11)$$

$$\frac{C_{g \text{Stokes}}}{C_0} = \frac{\tanh \sigma + \sigma(1 - \tanh^2 \sigma)}{2\sqrt{\sigma \tanh \sigma}} \quad (12)$$

with $C_0 = \sqrt{gh_0}$.
<figure><img src="image_4.png" /><figcaption>Figure 3. Optimum adjustment of dispersion correction parameter $\alpha$ over the range $0 \le \sigma \le 3.0$: (a) joint RMS error for phase and group velocities as a function of $\alpha$; and (b) comparison between model phase and group velocities and theoretical Stokes relations for the optimum value $\alpha = 0.053$: $(-) C_{\text{th}}/C_{\text{Stokes}}$, and $(- -) C_{g\text{th}}/C_{g\text{Stokes}}$.</figcaption></figure>

Hence, the optimization process consists in minimizing the error function (8) over $0 \le \sigma \le 3$. The associated error as a function of $\alpha$ is plotted in Figure 3(a). It is clearly seen that the joint RMSE has an absolute minimum in this dispersive range. The optimum value for $\alpha$ is 0.053 ($\approx \frac{1}{19}$), which is smaller than the analytical result given in Reference [7]. In Figure 3(b) the phase and group velocities obtained with this $\alpha$-value are plotted and compared with the corresponding Stokes relations. The agreement between the modelled dispersion relation and the theoretical one is fairly good over this $\sigma$-range. Larger discrepancies exist for the group velocity but the maximum error remains below 6% for all the physical wavelengths that the model is able to describe. Nevertheless, the error in group velocities starts to grow rapidly when $\sigma \ge 2.3$, showing that this property will be overestimated in this range.

Unless explicitly stated, all numerical examples that we will present in the following sections will have the dispersion correction parameter fixed at the optimum value determined here, i.e. $\alpha = 0.053$.

### 3. BOUNDARY CONDITIONS

It is worth emphasizing that boundary conditions (BCs) affect the numerical properties of the scheme being implemented. Stability and accuracy of high-order methods can be dramatically reduced if boundary treatment is not performed carefully [8, 9]. In the case of hyperbolic systems
it is generally acknowledged that the formal spatial accuracy of an inner Nth-order scheme is preserved if the system is closed by discretizing BCs with at least an (N - 1)th-order numerical approximation (see for instance Reference [10]). In the present work no analytical investigation of the numerical treatment of BCs is performed because of the complexity of the high-order method implemented. Instead, this important property is studied using a heuristic approach based on several benchmark tests. Our objective is to ensure that the discretized counterpart of BCs does not introduce any numerical instability in the system and at the same time preserves the spatial accuracy of the numerical scheme used in the inner region.

The numerical strategies intended to reproduce the physical properties of the system of equations at the boundaries must closely follow its intrinsic mathematical nature. Hence, analysis of the BCs will begin by studying the set of PDEs from a mathematical standpoint. Writing Serre equations in characteristic directions provides useful information for developing suitable BC's for this set. We focus on BCs able to reproduce important physical processes in a nearshore context, namely, an open sea BC, a moving shoreline strategy for the swash zone and a wall-type BC.

### 3.1. Decomposition of Serre equations in characteristic directions

Serre equations, which are recovered from (1)–(3) when $\alpha = 0$ and bottom variations are neglected, can be recast into a first-order system of equations by introducing auxiliary variables following the lines given in Reference [11], i.e.

$$h_t + uh_x + u_x h = 0 \qquad (13)$$

$$u_t + uu_x + gh_x + \frac{2}{3}h_x\Gamma + \frac{1}{3}h\Gamma_x = 0 \qquad (14)$$

$$\Gamma = w_t + uw_x \qquad (15)$$

$$w = -hu_x \qquad (16)$$

where $w$ is the vertical velocity evaluated at the free surface, and $\Gamma$ is related to the vertical acceleration of the fluid evaluated at the same location. Thus, all dispersive terms associated with the free surface curvature are contained in $\Gamma$. The mathematical properties of this first-order set of PDEs can be studied by writing it in characteristic form; the development is not reproduced here but further details can be found in References [11, 12]. This analysis shows that Serre equations possess two double characteristic directions expressed by

$$\frac{dx}{dt} = u, u \quad \text{and} \quad \frac{dx}{dt} = \infty, \infty \qquad (17)$$

and consequently the system is not of strictly hyperbolic type. From a mathematical point of view, this means that it is not possible to find Riemann invariants for this PDE set. Nevertheless, in order to write BCs which would reflect this hyperbolic–parabolic mathematical nature, it is useful to recast the Serre equations in the following quasi-hyperbolic form,

$$\frac{dR^+}{dt} = -\frac{1}{3h} \frac{\partial}{\partial x} (h^2 \Gamma) \quad \text{along} \quad \frac{dx}{dt} = u + \sqrt{gh} \qquad (18)$$

$$\frac{dR^-}{dt} = -\frac{1}{3h} \frac{\partial}{\partial x} (h^2 \Gamma) \quad \text{along} \quad \frac{dx}{dt} = u - \sqrt{gh} \qquad (19)$$
introducing positive and negative Riemann functions, $R^+ = u + 2\sqrt{gh}$ and $R^- = u - 2\sqrt{gh}$. Hence, vertical accelerations, which are neglected in nonlinear shallow water equations (NSWE), are responsible for the loss of hyperbolicity in Boussinesq-type equations by introducing a horizontal dependence in the characteristic plane $(x, t)$. It is important to recall that in this form, Equations (18)–(19) can only be applied at topographically uniform boundaries.

It is worth noting that, from a physical point of view, time scales associated to dispersive effects are likely to be larger than those associated with nonlinearities from intermediate to shallow waters. It may therefore be assumed that, over short distances/times, Riemann variables might be locally conserved along the characteristics. We will use this physical argument in what follows to develop suitable BCs for the finite volume resolution of the extended system of Serre equations.

## 3.2. Open sea boundary condition

The cross-shore Boussinesq model will be used to describe waves propagating over non-uniform beach profiles. Hence, an absorbing/generating BC is required at the seaward boundary. We need to prescribe an incident wave field there, allowing residual reflected waves to leave the computational domain without causing spurious oscillations. The numerical strategy consists in separating incident and outgoing waves on the basis of the quasi-hyperbolic form derived in the previous subsection. This kind of approach has been applied to Boussinesq-type equations in References [13, 14] following the lines given by Van Dongeren and Svendsen [15] for fairly long waves. We will use a similar approach but without assuming that incident and outgoing waves are of small amplitude (weakly nonlinear) as is the case in the aforementioned works. Instead, we locally neglect dispersive effects to estimate depth-averaged velocity and water depth at the seaward boundary using characteristic equations. We hypothesize that, over short times/lengths, dispersive effects do not influence boundary values. This heuristic treatment will be validated later by numerical tests. Therefore, assuming that $\sqrt{gh}>u$, the left and right-going characteristic curves are expressed by Equations (18)–(19), where vertical acceleration is neglected (see Figure 4 for a sketch of the positive and negative characteristic curves).

It should be recalled that a four stage Runge–Kutta time stepping is applied to integrate equations at interior nodes. This implies that the numerical solution is evaluated twice at $t + \Delta t/2$ and twice at $t + \Delta t$ [1]. On the other hand, seaward depth-averaged velocity and water depth values are obtained by solving systems (18)–(19). At each half time step, the right-going Riemann variable is estimated assuming that the incident wave field propagates from $-\infty$ (a fictitious domain) to the seaward boundary with a constant wave form. Hence, if the free surface elevation time series of incident waves at the seaward boundary is prescribed as the function $\eta_A(t)$, the corresponding Riemann variable should read,

$$R_L^+ = u_L + 2\sqrt{gh} = \sqrt{g[\eta_A(t + \Delta t/2) - \xi_A]} \left[ \frac{\eta_A(t + \Delta t/2)}{\eta_A(t + \Delta t/2) - \xi_A} + 2 \right] \quad (20)$$

where a long wave approximation is used to relate the flow velocity to the local water depth, i.e. $hu = (h+\xi)\sqrt{gh}$. Estimation of the negative Riemann variable requires more effort since the initial location at time $t$ of the outgoing characteristic, $x_R$, is not known beforehand. A first estimate for this location is obtained by integrating the left characteristic trajectory (19) and thus solving the following nonlinear equation for $x_R$:

$$x_R - x_A + \frac{\Delta t}{2}[u(x_R, t) - \sqrt{gh(x_R, t)}] = O(\Delta t^2) \quad (21)$$
<figure><img src="image_5.png" /><figcaption>Figure 4. Definition sketch of characteristic decomposition of incident and outgoing waves. C+ and C- represent respectively positive and negative characteristic curves in the (x, t) plane.</figcaption></figure>

A linear interpolation can be used to compute variables at location $x_R$ from known nodal depth-averaged velocities and water depths at time $t$. This implicit equation is easily solved by implementing, for example, a standard Newton-Raphson algorithm. Once $x_R$ is computed, the outgoing Riemann variable is estimated as

$$R_R^- = u(x_R, t) - 2\sqrt{gh(x_R, t)} \quad (22)$$

Then, solving systems (18)–(19) yields the first estimate for boundary values of $h$ and $u$ at $t + \Delta t/2$,

$$u_{A'} = \frac{1}{2}(R_L^+ + R_R^-) + O(\sigma^2 \Delta t, \Delta t) \quad (23)$$

$$h_{A'} = \frac{1}{16g}(R_L^+ - R_R^-)^2 + O(\sigma^2 \Delta t, \Delta t) \quad (24)$$

Cell-face values of water depth are readily computed for interior nodal points at $t + \Delta t/2$ from the finite volume integration of continuity equation (1). However, determining depth-averaged velocity from the momentum equation also means prescribing the dependent variable $q$ at the seaward boundary. This variable reflects the parabolic nature of the system and its boundary value must be fixed taking into account the horizontal dependence in the (x, t) plane given by Equation (3). Indeed, the mathematical form of this equation imposes an implicit treatment for this BC because the $q_A$ and internal j-values of the depth-averaged velocity are coupled. To cope with this additional constraint, Equation (3) is used to estimate $q_A$ from interior nodal values and
previously computed seaward variables $h_A$ and $u_A$. For this purpose, Equation (3) is discretized over the cell $\Omega_A = \{x \in [x_i=-1, x_i=0]\}$ (using the definitions given in the companion paper [1]) where nodal points located outside the physical domain are extrapolated with cubic piecewise functions, which provide fourth-order accuracy [16]. This additional equation must be solved simultaneously with those obtained at the interior nodes by implementing an iterative procedure. This is a straightforward matter because a deferred-correction approach is used to estimate $u_j$ values from the discretized counterpart of Equation (3).

The strategy used to advance $u$ and $q$ in time taking into account the BC for $u$ and $h$ is summarized as follows:

1. Reconstruct cell-face values $q_j$ from cell-averaged ones using the additional relation obtained from discretizing Equation (3) at the seaward boundary and the last computed $u_j$ values for interior nodes.

2. Using $q_j$ values estimated at the previous step, compute depth-averaged velocities at the interior nodes solving Equation (3).

3. If a convergence criterion for the relative error between the new and old values for $q_A$ and $u_j$ is not satisfied, then come back to 1. Otherwise, the iterative process stops.

Errors are estimated using the same expression based on a $L_1$ norm given in Reference [1]. The relative tolerance is fixed at $\tau = 10^{-4}$ for all computations.

At the end of this process all cell-face values for variables $h$, $u$ and $q$ are known at $t + \Delta t/2$. Nevertheless, the Runge-Kutta scheme will produce a second estimate of the variables at the same time level. In this second stage, characteristic equation (19) can be integrated with greater accuracy by taking advantage of the previous guess. The $x_R$ location is now computed using the following implicit equation:

$$x_R - x_A + \frac{\Delta t}{4}[u(x_R, t) + u_{A'} - \sqrt{gh(x_R, t)} - \sqrt{gh_{A'}}] = O(\Delta t^3) \quad (25)$$

and the second estimation for the boundary values reads,

$$u_{A'} = \frac{1}{2}(R_L^+ + R_R^-) + O(\sigma^2 \Delta t, \Delta t^2) \quad (26)$$

$$h_{A'} = \frac{1}{16g}(R_L^+ - R_R^-)^2 + O(\sigma^2 \Delta t, \Delta t^2) \quad (27)$$

where Riemann variables are defined as before. Thus, at each intermediate time step of the Runge-Kutta scheme, boundary values for water depth and depth-averaged velocity are estimated with an error of $O(\sigma^2 \Delta t, \Delta t^2)$. In the general case, the time step $\Delta t$ will be small enough and the error associated with dispersion would be small. Similarly, the boundary value for $q$ and velocities at the interior nodes are obtained following the same iterative strategy expounded before. The solution is then advanced to the next half step $t + \Delta t$ using an analogous procedure starting from the computed values at $t + \Delta t/2$.

It is worth emphasizing that even though a restrictive assumption regarding dispersion has been adopted in deriving seaward BCs, our two-step approach takes into account the parabolic nature of Boussinesq-type equations. Indeed, by using Equation (3) to compute the boundary value of the dependent variable $q$, dispersive effects are considered. Otherwise, $q$ would have been taken to be exactly the depth-averaged horizontal velocity $u$. However, this treatment will certainly not be
sufficient for highly dispersive waves (i.e. $kh_0 \gg 1$) where an internal source function in combination with a sponge layer will probably be the method of choice (see for instance Reference [17]).

### 3.3. Vertical wall

At a vertical and impermeable wall, horizontal velocity must vanish. Again, the characteristic decomposition of Serre equations provides a way to fix the value for the variable $h$. Therefore, the numerical strategy for estimating boundary values of $u$, $h$ and $q$ is equivalent to the one used for the open sea condition. The only difference is that the fictitious right-going Riemann variable, $R_L^+$, must be chosen in order to ensure that the horizontal velocity is zero at the corresponding boundary.

This particular BC is written for a situation where the vertical wall is located on a sloping bottom. Integration of the characteristic equations now provides the following relations:

$$R_{A'}^{+} - R_{L}^{+} = -g(\xi_{x})_{A} \frac{\Delta t}{2} + O(\sigma^{2}\Delta t) \quad (28)$$

$$R_{A'}^{-} - R_{R}^{-} = -g(\xi_{x})_{A} \frac{\Delta t}{2} + O(\sigma^{2}\Delta t) \quad (29)$$

where the associated bottom slope term is included as a source in the right-hand side. In order to impose a vanishing horizontal velocity at the wall we write,

$$R_{L}^{+} = 2g(\xi_{x})_{A} \frac{\Delta t}{2} - R_{R}^{-} \quad (30)$$

$$u_{A'} = -g(\xi_x)_A \frac{\Delta t}{2} + \frac{1}{2}(R_L^+ + R_R^-) + O(\sigma^2 \Delta t, \Delta t) \quad (31)$$

and the water depth at the boundary $A'$ is estimated with the same expression as before (Equation (24)). Once boundary values of $u$ and $h$ are computed in this way, the remaining variable $q$ is estimated seaward using the procedure outlined before. Therefore, it is possible to obtain greater accuracy in the second stage of the Runge-Kutta time stepping.

### 3.4. Moving shoreline

An accurate representation of swash motions at the shoreline is very important since infragravity waves are generated in nature by nonlinear interactions, by breaking and by partial reflection produced in particular by the run-up and run-down process that takes place in the swash zone (e.g. References [18, 19]). Several authors have developed numerical strategies intended to reproduce the dynamic movement of the shoreline in the framework of Boussinesq-type models. These schemes have been mainly applied using fixed grids, where the moving boundary problem is only solved in an approximate way. Although it is generally acknowledged that this kind of strategy can sometimes produce loss of mass and numerical instability (see for instance Reference [20]), it has been widely used in nearshore problems because the treatment of the moving shoreline is fundamentally simplified. Some recent examples of such approaches range from ad-hoc or essentially numerical approaches (e.g. References [21–23]) to more sophisticated ones where a Riemann wet-dry problem is explicitly solved [24].

In the present work, the simple extrapolation technique given by Lynett et al. [21] has been adapted in our finite volume resolution. Lynett et al.'s [21] moving shoreline technique is based
on a linear extrapolation of free surface and velocity through the wet-dry boundary. In the present finite volume scheme, wet cells are included in the computation of Equations (1)-(2) if,

$$ \frac{1}{\Delta x} \int_{\Omega_j} h \, dx \geq h_{\text{tol}} \qquad (32) $$

where $\Omega_j$ denotes a finite volume, $\Delta x$ is the grid size and $h_{\text{tol}}$ is a threshold value. Cell-averaged values for water depth $1/\Delta x \int h \, dx$ and cell-averaged variable $1/\Delta x \int q \, dx$ are linearly extrapolated through dry cells. Following Lynett et al. [21], a four-point filter may be passed through the extrapolated region in order to avoid slope discontinuities. Equations (1)-(2) are integrated in time and the location of the wet-dry interface is determined once every half time step ($\Delta t/2$) in the Runge-Kutta time stepping. Spatial derivatives are computed in the physical and extrapolated domain, where five extra dry nodal points are considered. Depth-averaged velocity is finally computed from Equation (3) in the extended domain (wet + dry).

The numerical implementation of this moving shoreline will be validated in the next section, where particular care must be paid to the issue of mass conservation. The linear extrapolation at the shoreline locally reduces the accuracy of the scheme to second order as pointed out in Reference [21] but allows larger time steps to be used in the present context. This is an important feature because our finite volume method may be integrated in time without necessarily respecting the CFL condition (i.e. Courant numbers larger than 1.0).

### 3.5. Compact differencing at boundary nodes

Boundary estimates for first and second-order spatial derivatives appearing in flux functions (4)-(5) are performed using compact approximations, which preserve fourth-order accuracy for first derivatives and third-order for second derivatives. The resulting matrix for the whole domain will be tri-diagonal even when non-periodic boundaries are considered. Matching conditions at the left boundary read [6],

$$ (\Phi_x)_A + 3(\Phi_x)_{j=1} = \frac{1}{\Delta x} \left( -\frac{17}{6}\Phi_A + \frac{3}{2}\Phi_{j=1} + \frac{3}{2}\Phi_{j=2} - \frac{1}{6}\Phi_{j=3} \right) \qquad (33) $$

$$ (\Phi_{xx})_A + 11(\Phi_{xx})_{j=1} = \frac{1}{\Delta x^2} (13\Phi_A - 27\Phi_{j=1} + 15\Phi_{j=2} - \Phi_{j=3}) \qquad (34) $$

where $\Phi_j$ is a generic function evaluated at node j and A represents the left boundary (see Reference [1] for nodal definitions). Similarly, matching conditions at the right boundary can be easily derived.

## 4. SPATIAL FILTERING

In the case of complex wave propagation problems over uneven bathymetries, dispersion and nonlinearity may interact to produce higher harmonics. These shorter waves may eventually lie out of the range of physical validity of the model equations and a high-order filtering technique is often considered to damp them out. This is also justified because numerical schemes used to discretize PDEs are not able to describe properly all wavelengths as was shown for our finite volume method in Reference [1].
Regarding the latter, the spectral analysis of the numerical scheme performed in Reference [1] provides an important theoretical background for improving the performance of the method. Indeed, since the Runge-Kutta finite volume scheme used here is based on centred-in-space approximations, it was shown that there is almost no damping of different wave modes when appropriate discretization parameters are used. However, not all wavelengths are properly described by the numerical method and examples of dispersion relation preserving (DRP) regions were given in the companion paper [1]. From a practical point of view this means that high frequency waves may introduce destabilizing effects as numerical integration advances in time. These short waves can arise as a consequence of nonlinear interaction or even be introduced by approximated boundary conditions, which do not match PDEs exactly. One effective way to eliminate non-physical high frequencies is to use high-order filtering techniques [25]. This approach is of widespread usage in the numerical treatment of wave-type equations (e.g. References [6, 26, 27]) and has also been used in the context of Boussinesq-type modelling (e.g. References [14, 22, 28]).

Consequently, our main concern here is to damp out frequencies that cannot be resolved by the numerical scheme in order to eliminate non-physical wavelengths. Since spectral analysis provides reliable information on the numerical properties of the finite volume scheme, an appropriate high-order filter can be chosen. The eighth-order implicit filter given by Gaitonde et al. [26] is used in the model. It can be applied to cell-averaged variables,

$$ \alpha_f \tilde{\Phi}_{i-1} + \tilde{\Phi}_i + \alpha_f \tilde{\Phi}_{i+1} = \sum_{n=0}^{4} \frac{a_n}{2} (\tilde{\Phi}_{i+n} + \tilde{\Phi}_{i-n}) \quad (35) $$

where $\tilde{\alpha}$ denotes filtered values, $\alpha_f$ is a free parameter ($-0.5 < \alpha_f < 0.5$) and the right-hand side coefficients read,

$$ \begin{aligned} a_0 &= \frac{93 + 70\alpha_f}{128}, & a_1 &= \frac{7 + 18\alpha_f}{16}, & a_2 &= \frac{-7 + 14\alpha_f}{32} \\ a_3 &= \frac{1 - 2\alpha_f}{16}, & a_4 &= \frac{-1 + 2\alpha_f}{128} \end{aligned} $$

When $\alpha_f = 0.5$ there are no damped frequencies and as $\alpha_f$ is reduced a wider range of frequencies is partially filtered as shown in Reference [26]. The filter can be applied once every $N_f$ Runge-Kutta time steps in order to eliminate high frequencies when computations are carried out for long periods. When high-order filtering is required, a fixed value of $\alpha_f = 0.4$ can be used, thus ensuring that only unresolved frequencies are suppressed.

Even if from a theoretical point of view the chosen high-order filter should not act on well resolved wavelengths, it is important to assess whether the numerical representation of wave kinematics and energy are affected. Indeed, the use of a filter to damp out high frequencies is delicate since it removes energy from the system. It is paramount then to evaluate the associated energy dissipation rate. This issue will be investigated first with the help of the following numerical example.

Consider a periodic domain with horizontal bottom where a motionless Gaussian hump located at its centre is prescribed as initial free surface condition (see Figure 5). The computational domain has a length of 20 m, and the still water depth is fixed at $h_0 = 1$ m. The initial condition reads,

$$ h(x, t=0) = H_0 \exp \left[ -\frac{1}{2s^2} (x - x_0)^2 \right] - \zeta(x) \quad (36) $$

$$ u(x, t=0) = 0 \quad (37) $$
<figure><img src="image_6.png" /><figcaption>Figure 5. Initial free surface configuration (−) for the numerical test, free surface profile after 3 s (−−) and free surface at $t = 500$ s (−·−).</figcaption></figure>

where $H_0$ is the height of the perturbed free surface, $x_0$ is its location and $s$ is a shape coefficient which sets its initial spatial length. The model is run over a long period (say ~10 000 × $\Delta$t) under circumstances where nonlinear interactions may generate higher harmonics. For instance, using $H_0/h_0 = 1.2$ and $s = 0.3$ m ensures that nonlinearities are relatively important and that dispersion will also have an influence. From this initial situation we compute the wave evolution in the periodic domain choosing $\Delta x = 0.05$ m and $\Delta t = 0.024$ s, which corresponds to a Courant number $C_r = 1.5$, and for 500 s (≈20 000 × $\Delta$t). A spatial snapshot of the initial perturbation, the free surface position few seconds after and at the end of the computation are depicted in Figure 5. It can be seen that at early stages there are several different wavelengths interacting with each other and that wave amplitudes are not small ($a/h_0 \sim 0.3$). During the computation, nonlinearities, dispersion and wave-wave interaction should favour the redistribution of energy at different wavenumbers and at $t = 500$ s different wavelengths are visible. Thus, this numerical example may be used to illustrate how filtering may affect computations.

Power density spectra estimated from free surface time series evaluated at $x = 10$ m, i.e. at the centre of the periodic domain are compared, and the results presented in Figure 6. It can be seen that there are no noticeable differences in the spectral distribution of energy when no filter is applied and when it is used once every 100 time steps or even once every time step. It is important to point out that, for this particular example, fairly short waves were generated since the spectral signature shows that some energy has been transferred to frequencies ranging between 0.6 and 1.0 Hz. Using the associated linear phase speed to compute wavelengths, these frequencies roughly correspond to $h_0/L \sim 0.2–0.3$. Such ratios are close to the theoretical deep water limit given by $h_0/L \sim 0.5$, thus proving that the present numerical example may be useful for testing the filtering technique when used in practical applications. An additional test concerning the filtering technique will be presented in Section 5.4 for wave propagation over a submerged bar.

Another important feature that should be investigated is how the total mass and energy contained in the periodic domain is affected by the application of the high-order filter. In theory, these quantities must be conserved for the entire computation and differences that may arise should be attributed to the numerical method used to integrate PDEs or to the filtering technique. In order to compute the total energy evolution contained in the computational domain, we use theoretical
<figure><img src="image_7.png" /><figcaption>Figure 6. Power density spectra of free surface time series evaluated at $x_0 = 10$ m: (a) applying the filter once per time step ($\alpha_f = 0.4$); (b) applying the filter once every 100 time steps ($\alpha_f = 0.4$); and (c) without using the eighth-order filter.</figcaption></figure>

information on horizontal and vertical velocity profiles for Serre equations (i.e. $\alpha = 0$) already presented in the companion paper [1]. It should be recalled that an *ad hoc* high-order dispersion correction term has been added to improve linear dispersive properties of the model. Therefore, when $\alpha \neq 0$ there is no analytical expression for the velocity profile. On the contrary, when fixing $\alpha = 0$ an exact expression (up to $O(\sigma^2)$) for the total kinetic and potential energy exists for the Serre equations. It follows that,

$$E_k^T(t) = \int_{x_l}^{x_r} \frac{1}{2} \rho h(u^2 + h^2 u_x^2) \, dx + O(\sigma^4) \approx \frac{1}{2} \rho \Delta x \sum_{j \in J} h_j (u^2 + h^2 u_x^2)_j \quad (38)$$

$$E_p^T(t) = \int_{x_l}^{x_r} \frac{1}{2} \rho g \eta^2 \, dx \approx \frac{1}{2} \rho g \Delta x \sum_{j \in J} \eta_j^2 \quad (39)$$

where $x_l$ and $x_r$ correspond, respectively, to the spatial coordinates of the left and right domain boundaries, and summation is carried out over the whole discrete domain, where $j$ is a nodal index and $J$ represents the set of discrete nodes. The total mass contained in the domain can be computed in a similar way. Time evolution of mass and total energy ($E_k^T + E_p^T$) in the case where the eighth-order filter is applied once per time step is plotted in Figure 7. It is confirmed that the filter does not alter mass nor energy since the relative errors on both quantities at the end of the computation are negligible, being less than 0.2% at $t = 500$ s.

This numerical example demonstrates that the use of the implicit eighth-order filter with $\alpha_f = 0.4$, as given in Reference [26], does not significantly affect energy distribution even when it is applied once per time step. Hence the filter is only acting on unresolved wavelengths and may be useful for improving numerical stability if computations are carried out for longer periods.
<figure><img src="image_8.png" /><figcaption>Figure 7. Time evolution of mass and energy in the computational domain when the eighth-order filter is applied once per time step ($\alpha_f = 0.4$): (a) relative error on mass; and (b) relative error on total energy (initial mass and energy are taken as references).</figcaption></figure>

# 5. MODEL VALIDATION

In the following subsections the model's capabilities will be evaluated by comparing its predictions with several experimental and numerical benchmark tests. Different numerical examples have been chosen in order to test in particular the physical validity of the extended system of Serre equations and the adequacy of the numerical implementation of BCs.

## 5.1. Swash motion and shoreline movement

Numerical implementation of the moving shoreline BC is validated using Carrier and Greenspan's [29] analytical solution for long wave propagation over a planar beach. Dispersive terms are switched off in order to use the analytical solution obtained for the dispersionless set of non-linear shallow water equations. The model is tested using the same parameters as in References [21–23], i.e. an incident sinusoidal wave with height $H = 0.006$ m and period $T = 10$ s propagating over a beach of constant slope 1:25, where the depth in the horizontal part of the channel is $h_0 = 0.5$ m.

The RK4 finite volume scheme theoretically allows larger time steps compared to previously published numerical Boussinesq-type models. The computation is therefore performed using a Courant number of 1.5 (i.e. $\Delta t/\Delta x\sqrt{gh_0} = 1.5$) and a grid size of $\Delta x = 0.05$ m. No spatial filtering is used in this particular case and comparisons between computed and analytical solutions are presented in Figure 8, where the threshold value for wet cells is fixed at $h_{tol} = 10^{-4}$ m. It can be seen that the agreement between the numerical and theoretical solutions is very good even when using a time step which is more than three times larger than the one used in the results reported for validating the original moving shoreline approach in Reference [21]. Only slight discrepancies arise in vertical excursion at the shoreline but the maximum and minimum wave run-up and run-down are very well predicted. Additionally, mass loss was almost negligible in this numerical example, being less than 0.005% over 100 s of computation.
<figure><img src="image_9.png" /><figcaption>Figure 8. Run-up and run-down at the shoreline. Comparison between Carrier and Greenspan’s [29] analytical solution (−) and numerical (−−) prediction using $Δx = 0.05$ m and $Δt = 0.034$ s: (a) vertical excursion at the shoreline; and (b) spatial free surface profiles at different times.</figcaption></figure>

It may be concluded that the moving shoreline BC proposed by Lynett *et al.* [21] can be successfully adapted to our finite volume scheme. Indeed, the numerical results compare favourably with the ones reported by these authors, even if larger time steps are used. However it is important to go further in the validation and test this approach in situations where dispersive terms are also included.

The run-up and run-down of a solitary wave propagating over a planar beach of 1:19.85 slope are used for this purpose. The model results are compared with laboratory measurements carried out by Synolakis [30] for an incident solitary wave of relative amplitude $a_0/h_0 = 0.28$, which broke strongly on run-up. Numerical computation is performed using a grid size $Δx = 0.025$ m, $C_r = 1.0$ and $h_{tol} = 10^{-4}$ m; the still water depth at the horizontal part of the channel is fixed at $h_0 = 0.25$ m and a four-point filter is passed through the extrapolated region following Lynett *et al.* [21]. The eddy viscosity-type breaking parameterization proposed by Cienfuegos *et al.* [31] is implemented in the Boussinesq model (further details can be found in Reference [32]).

Video measurements are available for this experiment and spatial snapshots for measured and computed free surface profiles are presented in Figure 9. As pointed out in Reference [33], incipient breaking takes place slightly before the dimensionless time $t^* = t\sqrt{g/h_0} = 20$ is reached. Even if some slight discrepancies may be noticed, the overall agreement is quite good. Some differences may be attributed to the lack of a friction term in the model. This phenomenon can be important when water becomes very shallow as in the run-down stage, where the computed water tongue is clearly thinner than the measured one. However, this example confirms that the moving shoreline condition can provide accurate predictions even in the presence of dispersion and breaking.
<figure><img src="image_10.png" /><figcaption>Figure 9. Spatial snapshots of free surface elevation for Synolakis [30] breaking solitary wave run-up and run-down experiment. $a_0/h_0 = 0.28$, $h_0 = 0.25$ m and $t^* = t\sqrt{g/h_0}$. (x): experimental data; (-) numerical computation using $\Delta x = 0.025$ m and $\Delta t = 0.0160$ s.</figcaption></figure>

## 5.2. Total reflection of a solitary wave at a wall

Validation of the wall BC can be performed on the basis of experimental measurement reported by Walkley and Berzins [34]. In this experimental set-up, two different tests where solitary waves propagate over a gentle slope of 1:50 collapsing into a vertical wall, were studied. The bathymetry for this problem is presented in Figure 10, where, following the numerical investigation performed
in Reference [34], at time $t=0$ s computations are initialized with the solitary wave centred at $x=50$ m. The vertical wall is located in our case at $x=0.0$ m and measured time series of free surface elevation are available at $x=2.25$ m. The still water depth in the horizontal part of the channel is $h_0=0.7$ m and incident solitary waves have a relative amplitude of $a_0/h_0=0.10$ and 0.174, respectively.

Numerical integration of the equations is conducted using the wall BC at $x=0.0$ m, a grid size $\Delta x = 0.07$ m and $C_r = 1.0$ in both runs and no filter is applied. Measured and computed free surface elevation time series are compared at the section $x = 2.25$ m in Figure 11. The numerical results are similar to those reported in Reference [34], where the first peak represents the incident wave and the second one corresponds to the reflected wave. While the first peak is very well reproduced, the second one is slightly overpredicted by the numerical model for both solitary waves. Nevertheless, this overprediction is significantly less than the results obtained with the weakly nonlinear and weakly dispersive finite element Boussinesq code developed in Reference [34]. Indeed, for the

<figure><img src="image_11.png" /><figcaption>Figure 10. Bathymetric configuration for the solitary wave test (total reflection on a vertical wall) and initial location of the soliton.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Comparison between measured (–) and computed (—) free surface elevation time series at $x = 2.25$ m; $\Delta x = 0.07$ m and $\Delta t = 0.0267$ s: (a) incident wave amplitude: $a_0/h_0 = 0.10$; and (b) incident wave amplitude: $a_0/h_0 = 0.174$.</figcaption></figure>
incident wave of relative amplitude $a_0/h_0 = 0.10$ the maximum value of free surface elevation is roughly overpredicted by 3% in the present computation while Walkley and Berzins [34] reported an overprediction of 8% for this case. More importantly, for the second solitary wave of amplitude $a_0/h_0 = 0.174$, our model overpredicts the second peak by roughly 5% while the weakly nonlinear code produced an error of nearly 38%. Hence, error discrepancies must be attributed to the weakly nonlinear hypothesis embedded in Boussinesq-type equations solved with the finite element model. In our case, the slight discrepancies between measured and computed secondary peaks may be reasonably attributed to the inviscid hypothesis used to derive Boussinesq-type equations since there are no losses at the reflecting wall as pointed out by Walkley and Berzins [34].

The present numerical test thus validates the wall BC developed here and confirms the importance of including the so-called full nonlinearity in Boussinesq-type equations.

### 5.3. Model robustness and open sea boundary condition

In this subsection the ability of the open sea BC to evacuate different wavelengths is tested. Additionally, important properties of the model are also evaluated using a challenging numerical test. Our aim is to investigate the stability of the model when used to compute wave propagation phenomena over sharp bathymetries producing highly discontinuous bottom gradients.

The first example aims at exploring the influence of the source term by running the model at rest over a non-trivial bathymetry. It may be recalled from the governing equations that the source-like term appearing in the right-hand side of (2) contains bottom derivatives up to second order. This non-conservative term must be *well-balanced* in order to ensure the stability and accuracy of the numerical scheme (see for instance References [35–38]). The bathymetry chosen for this first example was proposed in the *working group on dam break modelling* [39] to test shallow water equation solvers. The bed configuration is presented in Figure 12(a) and represents a very challenging geometry that is useful for testing conservative properties of discretized source terms in numerical models. The channel length is 1500 m and the still water depth in the horizontal part is taken to be $h_0 = 15$ m. In Figure 12(b) and (c) first and second bottom spatial derivatives are plotted as estimated from fourth-order compact finite difference formulae [6]. It is clearly seen that this non-trivial bathymetry produces sharp gradients and consequently the numerical capability of the discretized source term may be studied from this example.

We perform a computation using periodic BCs and without applying the filter. Numerical integration is carried out with $\Delta x = 3.0$ m and $\Delta t = 0.371$ s for 7500 s running the model at rest. It is worth noting that this time step corresponds to a Courant number $C_r = 1.5$ in the present example and that computation is thus carried out over more than 20 000 time steps. These numerical values for space and time grid resolution were chosen because the RK4 finite volume model is expected to be non-damping in this case (see Reference [1]). The free surface profile obtained at $t = 7500$ s is presented in Figure 13(a), while depth-averaged velocity and flow rate at the same time coordinate are depicted in Figure 13(b) and (c). From these figures it can be seen that, even if some small disturbances ($\sim 10^{-5}$ m) are visible at the free surface, the flow velocity is negligibly small over the whole domain at the end of the computation. This result proves that the model is able to deal with steady-state conditions and that small free surface disturbances will not be amplified when running the model at rest. Therefore, it appears that the centred in space strategy used to discretize the source term does not introduce any numerical unbalance in fluxes and no spurious behaviour has been noticed even when computations are carried out for long periods over this challenging benchmark test.
<figure><img src="image_13.png" /><figcaption>Figure 12. Geometrical configuration with *j*-nodes denoted by dots: (a) bottom bathymetry and free surface at the end of the computation; (b) numerically estimated first spatial derivative of the bottom; and (c) numerically estimated second spatial derivative of the bottom.</figcaption></figure>

In order to investigate the adequacy of the absorbing BC, an additional numerical test with a motionless Gaussian hump with $H_0 = 0.10$ m and $s = 0.5$ m as initial condition is considered. The geometrical configuration for this example is sketched in Figure 14 where a step with height $\Delta h = 0.4$, width $\Delta L = 5$ m and side slopes of 1:10 is centred in the computational domain whose still water depth is fixed at $h_0 = 0.5$ m in the horizontal part of the channel. At left and right boundaries the open sea condition is applied in order to evacuate the extra volume of water associated with the Gaussian hump.

Numerical integration of the equations is carried using $\Delta x = 0.025$ m and a Courant number of 1.5. The eighth-order implicit filter is passed over cell-averaged variables once per time step using $\alpha_f = 0.4$. The numerical computation is shown in Figure 14, where it can be seen that short waves arise as the resulting perturbations travel in both directions and enter into the deep part of the channel. It can be seen that different wavelengths are adequately evacuated at left and right boundaries by the open sea condition, but some small residual waves are radiated back into the computational domain. It is worth noting that these short waves are beyond the physical validity of the model equations. Indeed, partially reflected waves have wavelengths $L \sim 0.5-1$ m, thus $h_0/L \sim 0.5-1$ in this example ($h_0/L = 0.5$ corresponds to the deep water limit). This test shows that the absorbing BC partially reflects short wavelengths. This situation is consistent with the approximation made when invoking the characteristic equations and can be improved by
<figure><img src="image_14.png" /><figcaption>Figure 13. Computed results at $t = 7500$ s using $\Delta x = 0.3$ m and $\Delta t = 0.371$ s: (a) free surface, $z = \eta$; (b) depth-averaged velocity, $u$; and (c) flow rate ($=hu$).</figcaption></figure>

reducing the time step. Nevertheless, for practical purposes involving wave propagation over gently sloping beaches, reflected waves generated by the bottom bathymetry and the moving shoreline are expected to be of smaller amplitude and longer period than incoming waves [18, 19] because the breaking process must dissipate most of the wave energy. If highly dispersive waves are to be evacuated, sponge layers and internal source functions for prescribing incident wave fields should be considered (see for instance Reference [17]). Hence, the present example is also useful for illustrating the limitations of the open sea BC implemented.

### 5.4. Regular waves propagating over a submerged bar

A very important test for wave propagation models was presented by Dingemans [40] in the framework of the MAST II G8-M project. In this work the performance of several Boussinesq-type models was analysed by comparing numerical predictions with experimental measurements of regular waves propagating over a submerged bar [41, 42]. The experimental set-up was conceived to investigate the frequency dispersion characteristics and nonlinear interaction of complex wave propagation phenomena and has become an unavoidable and challenging benchmark test for any phase resolving numerical code. Major difficulties arise in this experimental set-up because higher bounded harmonics, which are generated by nonlinear interaction as waves shoal in the
<figure><img src="image_15.png" /><figcaption>Figure 14. Spatial snapshots of free surface at different times using a motionless Gaussian hump of water as initial condition. $H_0 = 0.10 \text{ m}$, $s = 0.5 \text{ m}$, $\Delta x = 0.025 \text{ m}$, $\Delta t = 0.0169 \text{ s}$ and passing a eighth-order implicit filter at interior domain with $\alpha_f = 0.4$.</figcaption></figure>

upstream slope of the bar, are freely released as they encounter a sudden increase in water depth downslope. Thus, numerical models, which are usually derived using long wave approximations, are not able to describe the wave kinematics of shorter waves correctly [40]. Both nonlinear and dispersive properties embedded in the model equations contribute to the successful respresentation of measured time series of free surface elevations as demonstrated by Gobbi and Kirby [28], thus making this test a very demanding one.

Even though research efforts in recent years have successfully led to major improvements in the linear dispersive characteristics and nonlinear properties of Boussinesq-type equations (see for instance Reference [43] for a review), this has been achieved mainly by incorporating higher-order terms, which introduce numerical complexity. An alternative way of improving model capabilities has been recently investigated by several authors who have developed multi-layer approaches for deriving extended Boussinesq equations, which provide a better description of short wavelengths
without necessarily introducing very high-order derivatives but at higher computational cost [44, 45]. Therefore, any extension of Boussinesq models to deal with highly nonlinear and extremely dispersive waves increases the CPU time. Nevertheless the aim here is to develop an efficient numerical scheme that is able to describe wave propagation phenomena over natural beaches, and consequently nonlinearity has been considered a more important property to enhance than dispersion. However, wave propagation over a submerged bar constitutes a demanding test, which can be used in practice to evaluate the theoretical limitations of the chosen set of Boussinesq-type equations.

Having the latter in mind, we compare the numerical predictions with the non-breaking experimental data available in References [41, 42] and summarized by Dingemans [40]. Two sets of regular waves propagating over a bar and generated in a small wave-flume are considered (cases A and C from Reference [40]). Free surface elevation was measured at several wave gauges located before, above and after the bar. A sketch of the bottom configuration and the location of the wave gauges is depicted in Figure 15. The reader is referred to Dingemans [40] and Gobbi and Kirby [28] for further information on the experimental set-up.

The first test (case A) was run for periodic waves with incident amplitude $a_0 = 0.01$ m and period $T = 2.02$ s. The absorbing-generating condition is considered at the left boundary in order to prescribe measured time series of free surface elevation at the location $x = 5.7$ m, just before the toe of the bar. The remaining wave gauges located in the shoaling region, above and behind the bar, are used to compare the model predictions with laboratory measurements. The equations are integrated using a grid size $\Delta x = 0.04$ m and a Courant number $C_r = 1.5$, which corresponds to a $\Delta t = 0.303$ s. The computational domain is made long enough in order to avoid any interference with reflected waves that might be generated at the right boundary during the computation. In Figure 16, the numerical results are compared with measured time series at the same time windows taking care to synchronize the input data because there was a mismatch error for the wave gauge located at $x = 5.7$ m, as pointed out by Gobbi and Kirby [28].

The model performs very well up to the gauge located at $x = 15.7$ m downslope of the bar while, as expected from the theoretical limitations of the present Boussinesq set, some discrepancies show up as higher harmonics are released behind the bar. Nevertheless, the agreement is still quite

<figure><img src="image_16.png" /><figcaption>Figure 15. Bottom bathymetry for tests reported in Reference [40] and spatial location of wave gauges (all distances in meters).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16. Comparison of free surface elevation time series for test case A from Dingemans [40] at several locations. $a_0 = 0.01$ m, $T = 2.02$, $\Delta x = 0.04$ m and $\Delta t = 0.303$ s. (-): measured, (---) computed.</figcaption></figure>

reasonable and is found to be comparable to the results reported by Gobbi and Kirby [28] using the fully nonlinear and weakly dispersive model developed by Wei *et al.* [46]. However, it is worth emphasizing that in our case, time integration was performed using a Courant number of 1.5 while the numerical results reported in Reference [28] were conducted with a Courant number that remained below 0.3 in order to ensure stability. Therefore, it is confirmed that the present compact finite volume scheme is more efficient than the finite difference scheme used to discretize the PDE system in Reference [28].
Experiment C reported by Dingemans [40] is considered next. Here, incident waves have an amplitude $a_0 = 0.0205$ m and period $T = 1.01$ s. This is a far more difficult test for Boussinesq models since the wave amplitude is twice the value prescribed before, and the incident wavelength is nearly 2/5 shorter. As in the previous example, the incident wave field is prescribed at the left boundary using free surface time series measured by the wave gauge located at $x = 5.7$ m. The numerical model is run with $\Delta x = 0.04$ m and $\Delta t = 0.303$ s, i.e. a Courant number of 1.5 and the computational domain is again made long enough to avoid residual reflection. In this example, the dispersion parameter for incident waves is $kh_0 = 1.69$, thus more than half the theoretical dispersive limit for the present Boussinesq-type equations (see Section 2). This means that the higher harmonic decomposition that takes place after the bar will rapidly generate wavelengths that are beyond the range of validity of the model. Consequently, phase speed estimates for freely released short waves may not be accurately reproduced.

The numerical results are reported and compared with experimental data at several locations in Figure 17. On the shoal and above the bar crest up to gauge location $x = 13.5$ m, the computed free surface elevation shows good agreement with the measurements. This result proves that the Boussinesq model is able to adequately reproduce the nonlinear energy transfer, from the leading wave component to higher harmonics, that takes place in this region. After waves pass over the bar, higher harmonics are released as free waves, thus propagating at different phase speeds. Discrepancies then arise between the computed and measured free surface because dispersive properties of the Boussinesq model are not accurate enough for waves with $kh_0 > \pi$. In this experiment, the second harmonic is well above this theoretical limit and the numerical results are similar to those reported by other authors using equivalent sets of fully nonlinear and weakly dispersive Boussinesq-type equations (e.g. the one layer model in Reference [44] and Wei et al.'s [46] model in Reference [28]). However, the temporal and spatial grid sizes used to integrate the equations in the aforementioned works were smaller than the ones used in the present example.

In order to analyse further the effects of the high-order spatial filtering technique on the spectral properties of the computed results, a wave propagation test over Dingemans [40] bathymetry is performed. This example offers the opportunity of investigating the filter's performance under forced conditions, where energy is continuously added in the system. The model is run using the second-order Stokes theory to prescribe the incident wave field, fixing $a_0 = 0.01$ m and $T = 1.01$ s. As previously noted, after waves pass the bar, higher harmonics are freely released and, as in Dingemans [40] test case C, for $T = 1.01$ s, the second harmonic is still beyond the range of validity of the equations in terms of linear dispersion. To investigate the filter's performance, power density spectra evaluated from the numerical model at location $x = 21$ m are compared, by applying the filter once per time step, once every 100 time steps and without using it. The numerical computations are carried out with $\Delta x/h_0 = 0.1$ and $C_r = 1.5$, which should produce non-damping solutions [1]. The results are presented in Figure 18 taking $\alpha_f = 0.4$ for the eighth-order filter described in Section 4. In this plot the fundamental frequency is located at $f = 1$ Hz, and higher harmonics are visible at $f = 2, 3, 4, ...$ Hz. It can be seen that up to the fourth harmonic there are no noticeable differences in the power density spectra for computations performed with or without applying the filter. It is only at the fifth harmonic that its effect can be seen. Indeed, the small amount of energy present at $f = 5$ Hz is suppressed. In addition, it appears that applying the filter once per time step or once every 100 time steps has no great influence on the overall result. It is thus confirmed that the eighth-order filter does not act on resolved wavelengths. In this particular example, it extracts energy from the system only at the fifth harmonic while the
<figure><img src="image_18.png" /><figcaption>Figure 17. Comparison of free surface elevation time series for test case C from Dingemans [40] at several locations. $a_0 = 0.0205$ m, $T = 1.01$, $\Delta x = 0.04$ m and $\Delta t = 0.303$ s. (-): measured, (---) computed.</figcaption></figure>

physical wavelength that the model can accurately represent roughly corresponds to the second harmonic.

5.5. *Nonlinear shoaling of solitary waves propagating over a beach*

In this subsection, the nonlinear shoaling characteristics of the model are investigated by comparing
the numerical results with laboratory measurements of solitary waves propagating over planar
beaches. This is an important check for the nonlinear properties of the chosen set of Boussinesq-
type equations in the framework of nearshore wave propagation because on natural beaches, waves
<figure><img src="image_19.png" /><figcaption>Figure 18. Power density spectra for free surface time series computed at $x_0 = 21$ m, test case C from Dingemans [40]: (a) applying the filter once per time step ($\alpha_f = 0.4$); (b) applying the filter once every 100 time steps ($\alpha_f = 0.4$); and (c) without using the eighth-order filter.</figcaption></figure>

running up gentle slopes behave in a quasi independent manner, thus showing close resemblances to solitary waves [47]. Moreover, accurate prediction of the limiting amplitude of nearly breaking waves is paramount in the context of morphodynamics since the location at which wave-breaking occurs is often related to sandbar formation.

Experiments were conducted in the LEGI's wave-flume (36 m) at Grenoble with two different beach slope configurations of 1:30 and 1:60. In the first case, the still water depth was fixed at $h_0 = 0.25$ m in the horizontal part of the channel, while in the second configuration, the still water depth was $h_0 = 0.18$ m [48]. Four solitary waves were generated for each bathymetry using a piston wavemaker and following the experimental procedure given by Guizien and Barthélemy [49]. Free surface displacements were measured at several wave gauges located just before breaking using resistive probes. Measurement locations and relative wave amplitudes in the experiments are summarized in Table I for the beach slope of 1:30 and in Table II for the planar beach with a slope of 1:60.

Four solitary wave amplitudes were studied in order to investigate the performance of the Boussinesq model when used to describe the shoaling of waves with different lengths. Indeed, the smaller the solitary wave amplitude, the longer the equivalent wavelength. This is an important test, which can illustrate the physical relevance of the extended system of Serre equations in real world applications.

The model is initialized using Serre's solitary wave solution [50] with the same measured incident amplitude. For both beach configurations, the computations are carried out with a grid size $\Delta x/h_0 = 0.1$ and a Courant number $C_r = 1.0$. The adequacy of the numerically generated solitary wave is checked using measurements from an additional resistive probe located in the
<table><caption>Table I. Location of wave gauges for solitary waves shoaling on a planar beach with slope 1:30 and measured and computed wave amplitudes.</caption><tbody><tr><td colspan="6"><strong>Incident wave amplitude: a<sub>0</sub>/h<sub>0</sub> = 0.096</strong></td></tr><tr><td>Probe location relative to the shoreline (m)</td><td>2.430</td><td>2.215</td><td>1.960</td><td>1.740</td><td>1.502</td></tr><tr><td>Measured wave amplitude (a/|ξ|)</td><td>0.368</td><td>0.415</td><td>0.495</td><td>0.582</td><td>0.725</td></tr><tr><td>Computed wave amplitude (a/|ξ|)</td><td>0.365</td><td>0.408</td><td>0.472</td><td>0.547</td><td>0.659</td></tr><tr><td>Relative error (%)</td><td>-0.7</td><td>-1.5</td><td>-4.6</td><td>-6.0</td><td>-9.2</td></tr><tr><td colspan="6"><strong>Incident wave amplitude: a<sub>0</sub>/h<sub>0</sub> = 0.298</strong></td></tr><tr><td>Probe location relative to the shoreline (m)</td><td>3.980</td><td>3.765</td><td>3.510</td><td>3.290</td><td>3.052</td></tr><tr><td>Measured wave amplitude (a/|ξ|)</td><td>0.631</td><td>0.699</td><td>0.766</td><td>0.843</td><td>0.940</td></tr><tr><td>Computed wave amplitude (a/|ξ|)</td><td>0.649</td><td>0.702</td><td>0.773</td><td>0.854</td><td>0.943</td></tr><tr><td>Relative error (%)</td><td>2.7</td><td>0.3</td><td>0.9</td><td>1.3</td><td>0.4</td></tr><tr><td colspan="6"><strong>Incident wave amplitude: a<sub>0</sub>/h<sub>0</sub> = 0.456</strong></td></tr><tr><td>Probe location relative to the shoreline (m)</td><td>4.910</td><td>4.695</td><td>4.440</td><td>4.220</td><td>3.982</td></tr><tr><td>Measured wave amplitude (a/|ξ|)</td><td>0.724</td><td>0.800</td><td>0.852</td><td>0.921</td><td>0.983</td></tr><tr><td>Computed wave amplitude (a/|ξ|)</td><td>0.746</td><td>0.789</td><td>0.852</td><td>0.907</td><td>0.962</td></tr><tr><td>Relative error (%)</td><td>3.0</td><td>-1.4</td><td>0.0</td><td>-1.5</td><td>-2.1</td></tr><tr><td colspan="6"><strong>Incident wave amplitude: a<sub>0</sub>/h<sub>0</sub> = 0.534</strong></td></tr><tr><td>Probe location relative to the shoreline (m)</td><td>5.180</td><td>4.965</td><td>4.710</td><td>4.490</td><td>4.252</td></tr><tr><td>Measured wave amplitude (a/|ξ|)</td><td>0.818</td><td>0.867</td><td>0.944</td><td>1.012</td><td>1.067</td></tr><tr><td>Computed wave amplitude (a/|ξ|)</td><td>0.807</td><td>0.861</td><td>0.920</td><td>0.986</td><td>1.063</td></tr><tr><td>Relative error (%)</td><td>-1.3</td><td>-0.7</td><td>-2.5</td><td>-2.6</td><td>-0.4</td></tr></tbody></table>

<table><caption>Table II. Location of wave gauges for solitary waves shoaling on a planar beach with slope 1:60 and measured and computed wave amplitudes.</caption><tbody><tr><td colspan="6"><strong>Incident wave amplitude: a<sub>0</sub>/h<sub>0</sub> = 0.091</strong></td></tr><tr><td>Probe location relative to the shoreline (m)</td><td>4.500</td><td>4.285</td><td>4.030</td><td>3.810</td><td>3.572</td></tr><tr><td>Measured wave amplitude (a/|ξ|)</td><td>0.272</td><td>0.297</td><td>0.321</td><td>0.356</td><td>0.400</td></tr><tr><td>Computed wave amplitude (a/|ξ|)</td><td>0.275</td><td>0.295</td><td>0.324</td><td>0.352</td><td>0.392</td></tr><tr><td>Relative error (%)</td><td>0.8</td><td>-0.6</td><td>1.0</td><td>-1.0</td><td>-2.2</td></tr><tr><td colspan="6"><strong>Incident wave amplitude: a<sub>0</sub>/h<sub>0</sub> = 0.286</strong></td></tr><tr><td>Probe location relative to the shoreline (m)</td><td>5.450</td><td>5.235</td><td>4.980</td><td>4.760</td><td>4.522</td></tr><tr><td>Measured wave amplitude (a/|ξ|)</td><td>0.715</td><td>0.783</td><td>0.832</td><td>0.917</td><td>1.007</td></tr><tr><td>Computed wave amplitude (a/|ξ|)</td><td>0.724</td><td>0.773</td><td>0.847</td><td>0.909</td><td>0.993</td></tr><tr><td>Relative error (%)</td><td>1.3</td><td>-1.4</td><td>1.9</td><td>-0.9</td><td>-1.4</td></tr><tr><td colspan="6"><strong>Incident wave amplitude: a<sub>0</sub>/h<sub>0</sub> = 0.479</strong></td></tr><tr><td>Probe location relative to the shoreline (m)</td><td>7.950</td><td>7.735</td><td>7.480</td><td>7.260</td><td>7.022</td></tr><tr><td>Measured wave amplitude (a/|ξ|)</td><td>0.722</td><td>0.770</td><td>0.809</td><td>0.858</td><td>0.914</td></tr><tr><td>Computed wave amplitude (a/|ξ|)</td><td>0.722</td><td>0.748</td><td>0.789</td><td>0.828</td><td>0.871</td></tr><tr><td>Relative error (%)</td><td>0.0</td><td>-2.8</td><td>-2.5</td><td>-3.5</td><td>-4.7</td></tr><tr><td colspan="6"><strong>Incident wave amplitude: a<sub>0</sub>/h<sub>0</sub> = 0.558</strong></td></tr><tr><td>Probe location relative to the shoreline (m)</td><td>9.150</td><td>8.935</td><td>8.680</td><td>8.460</td><td>8.222</td></tr><tr><td>Measured wave amplitude (a/|ξ|)</td><td>0.681</td><td>0.721</td><td>0.758</td><td>0.798</td><td>0.851</td></tr><tr><td>Computed wave amplitude (a/|ξ|)</td><td>0.678</td><td>0.701</td><td>0.728</td><td>0.758</td><td>0.793</td></tr><tr><td>Relative error (%)</td><td>-0.4</td><td>-2.8</td><td>-3.9</td><td>-5.0</td><td>-6.8</td></tr></tbody></table>
<figure><img src="image_20.png" /><figcaption>Figure 19. Time series of free surface elevation for solitary waves shoaling on a planar beach of slope 1:30 (where $t^* = t\sqrt{g/h_0}$). (–) Measured values, and (---) computed using $\Delta x = 0.025$ m and $\Delta t = 0.0160$ s.</figcaption></figure>

horizontal part of the flume. Similarly, this time series is used to synchronize the remaining probes located in the shoaling region.

Comparisons of free surface elevations between measured and computed waves are presented in Figures 19 and 20 for the 1:30 and 1:60 slopes, respectively. Overall agreement between the computed and measured wave amplitudes and shapes is found to be adequate for the different cases investigated. Larger amplitude errors occur for the incident solitary wave $a_0/h_0 = 0.096$ with the beach slope of 1:30, and for $a_0/h_0 = 0.558$ with the 1:60 slope, thus showing the important influence of the bottom bathymetry in terms of model performance. Nevertheless, the error is never larger than 10% as shown in Tables I and II, and the Serre model is performing outstandingly in predicting nonlinear shoaling even for relative amplitudes as high as $a/|\xi| \sim 1.0$ and for both beach slopes. This is an extremely nonlinear situation where waves are nearly breaking. It should be recalled that the breaking criterion for horizontal bottoms based on solitary wave theory by Munk [51], predicts a limiting amplitude of $a/|\xi| \simeq 0.78$. This confirms then that the present experimental test is useful for investigating the performance of the model in strongly nonlinear situations.

## 6. CONCLUSIONS

This paper deals with the validation of the 1D finite volume Boussinesq-type model described in Reference [1]. BCs intended to reproduce an absorbing/generating seaward boundary and an
<figure><img src="image_21.png" /><figcaption>Figure 20. Time series of free surface elevation for solitary waves shoaling on a planar beach of slope 1:60 (where $t^* = t\sqrt{g/h_0}$). (–) Measured values, and (---) computed using $\Delta x = 0.018$ m and $\Delta t = 0.0135$ s.</figcaption></figure>

impermeable vertical wall located on a sloping bottom is developed with the help of a quasi-hyperbolic decomposition of Serre equations in characteristic directions. Similarly, the numerical strategy for shoreline movement which consists of a linear extrapolation for fictitious dry nodes through the wet-dry interface proposed by Lynett *et al.* [21] is successfully implemented in the model. Even though the numerical approximations used to close the system of equations at the boundary nodes reduce the formal accuracy of the scheme, several benchmark tests demonstrate that accurate predictions of modelled water wave processes can be obtained without affecting the stability of the high-order finite volume method. Indeed, numerical examples chosen to validate the different BCs are conducted with large time steps using Courant numbers ranging between 1.0 and 1.5. The numerical computations performed with these high Courant numbers compare favourably with the results reported in former investigations using much smaller time steps (e.g. References [21, 28, 34]).

The experimental tests for regular waves propagating over a bar reported by Dingemans [40] are used to check the dispersive properties of the Boussinesq-type model. The numerical results confirm several theoretical limitations of the set of PDEs, in particular the physical range of wavelengths that the model is able to describe. Freely released shorter wavelengths propagating after the bar are not accurately captured but the computed results appear to be in agreement with those reported in References [28, 34, 44] using equivalent sets of Boussinesq-type equations.
Comparisons with numerical results previously reported by Walkley and Berzins [34] for total reflection and collapsing of a highly nonlinear soliton on a vertical wall clearly confirm the importance of including the so-called full nonlinearity in Boussinesq-type equations. Moreover, several experimental measurements of solitary waves propagating over two different beach slopes up to nearly breaking conditions indicate that Serre equations can accurately predict nonlinear shoaling evolution and the associated limiting wave amplitudes. The relative error between predicted and computed amplitudes are never larger than 10% in all test cases. This important result suggests that this particular set of Boussinesq-type equations may be used to describe the propagation of nearly breaking or breaking waves. A numerical computation including both, dispersion and breaking, is performed for the run-up, breaking and run-down of an experimental solitary wave [30]. This test confirms that the moving shoreline BC performs well even for a realistic situation; moreover, it is a good opportunity to show that the inclusion of the wave-breaking parameterization described in References [31, 32] allows for the extension of the model into the surf zone.

Finally, a practical and efficient way of suppressing spurious short-wave noise that can arise as a consequence of nonlinear interaction or mismatches between discretized BCs and PDEs is presented. The high-order implicit filtering technique developed in Reference [26] is implemented and numerical examples demonstrate that the latter may be used to improve model stability in practical applications without affecting the spectral resolution of physical wavelengths.

## 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. Professor Maurizio Brocchini from Universita' di Genova is also acknowledged for his useful suggestions and comments.

## REFERENCES

1. Cienfuegos R, Barthélemy E, Bonneton P. 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):1217-1253.

2. 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.

3. 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, Series A* 1998; **356**(1749):3123-3181.

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

5. 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.

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

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

8. Carpenter MH, Gottlieb D, Abarbanel S. The stability of numerical boundary treatments for compact high-order finite difference schemes. *Journal of Computational Physics* 1993; **108**:272-295.

9. Sengupta TK. *Fundamentals of Computational Fluid Dynamics* (1st edn). Universities Press: India, 2004.

10. Gustafsson B. The convergence rate for difference approximations to mixed initial boundary value problems. *Mathematics of Computation* 1975; **29**(130):396-406.

11. Staub R. Contribution à l'étude des ondes de gravité en eau peu profonde. *Ph.D. Thesis*, Université de Grenoble, France, 1972.

12. Dingemans MW. *Water Wave Propagation Over Uneven Bottoms* (1st edn). World Scientific Publishers: Singapore, 1997.
13. Musumeci RE. The flow in the surf zone: a fully nonlinear Boussinesq model for breaking waves. *Ph.D. Thesis*, University of Catania, Italy, 2002.

14. Musumeci RE, Svendsen IA, Veeramony J. The flow in the surf zone: a fully nonlinear Boussinesq-type approach. *Coastal Engineering* 2005; **52**:565–598.

15. Van Dongeren AR, Svendsen IA. Absorbing-generating boundary conditions for shallow water models. *Journal of Waterway, Port, Coastal and Ocean Engineering (ASCE)* 1997; **123**(6):303–313.

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

17. Wei G, Kirby JT, Sinha A. Generation of waves in boussinesq models using a source function method. *Coastal Engineering* 1999; **36**:271–299.

18. Baldock TE, Holmes P, Horn DP. Low frequency swash motion induced by wave grouping. *Coastal Engineering* 1997; **32**:197–222.

19. Mase H. Frequency down-shift on swash oscillation compared to incident waves. *Journal of Hydraulic Research* 1995; **33**(3):397–411.

20. Balzano A. Evaluation of methods for numerical simulation of wetting and drying in shallow water flow models. *Coastal Engineering* 1998; **33**:83–107.

21. Lynett PJ, Wu TR, Liu PLF. Modeling wave runup with depth-integrated equations. *Coastal Engineering* 2002; **46**:89–107.

22. 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 (ASCE)* 2000; **126**(1):39–48.

23. Madsen PA, Sørensen OR, Schäffer HA. Surf zone dynamics simulated by a Boussinesq-type model. Part I. Model description and cross-shore motion of regular waves. *Coastal Engineering* 1997; **32**:255–288.

24. Bellotti G, Brocchini M. On the shoreline boundary conditions for Boussinesq-type models. *International Journal for Numerical Methods in Fluids* 2001; **37**:479–500.

25. Vichnevetsky R, Bowles JB. *Fourier Analysis of Numerical Approximations of Hyperbolic Equations* (2nd edn). SIAM: Philadelphia, PA, 1982.

26. Gaitonde DV, Shang JS, Young JL. Practical aspects of higher-order numerical schemes for wave propagation phenomena. *International Journal for Numerical Methods in Engineering* 1999; **45**:1849–1869.

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

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

29. Carrier GF, Greenspan HP. Water waves of finite amplitude on a sloping beach. *Journal of Fluid Mechanics* 1958; **4**:97–109.

30. Synolakis CE. The runup of solitary waves. *Journal of Fluid Mechanics* 1987; **185**:523–545.

31. Cienfuegos R, Barthélemy E, Bonneton Ph. A new wave-breaking parametrization for Boussinesq-type equations. *Fifth International Symposium on Ocean Wave Measurements and Analysis (Waves'05)*, Madrid, Spain, 2005, CD-ROM.

32. Cienfuegos R. Numerical modelling of two dimensional water wave propagation processes and topographically induced breaking. *Ph.D. Thesis*, Institut National Polytechnique de Grenoble, France, 2005.

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

34. 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.

35. García-Navarro P, Vazquez-Cendón ME. On the numerical treatment of source terms in the shallow water equations. *Computers and Fluids* 2000; **29**:951–979.

36. Greenberg J-M, Leroux A-Y. A well-balanced scheme for the numerical processing of source terms in hyperbolic equations. *SIAM Journal on Numerical Analysis* 1996; **33**:1–16.

37. LeVeque RJ. Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm. *Journal of Computational Physics* 1998; **146**(1):346–365.

38. Vázquez-Cendón ME. Improved treatment of source terms in upwind schemes for the shallow water equations in channel with irregular geometry. *Journal of Computational Physics* 1999; **148**:497–526.

39. Goutal N, Maurel F. *Proceedings of the 2nd Workshop on Dam-Break Wave Simulation; Technical Report HE-43/97/016/a*, Électricité de France, LNHE, 1997.

40. Dingemans MW. Comparison of computations with Boussinesq-like models and laboratory measurements. *Report H-1684.12*, Delft Hydraulics, 1994; 32.

41. Beji S, Battjes JA. Experimental investigation of wave propagation over a bar. *Coastal Engineering* 1993; **19**:151–162.
42. Luth HR, Klopman G, Kitou N. Kinematics of waves breaking partially on an offshore bar; LDV measurements of waves with and without a net onshore current. *Report H-1573*, Delft Hydraulics, 1994; 40.

43. 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, Series A* 2003; **459**(2033):1075–1104.

44. Lynett P, Liu L-F. A two-layer approach to wave modelling. *Proceedings of the Royal Society of London, Series A* 2004; **360**(2049):2637–2669.

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

46. 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.

47. Stiassnie M, Peregrine DH. Shoaling of finite amplitude surface waves on water of slowly varying depth. *Journal of Fluid Mechanics* 1980; **97**:783–805.

48. Guibourg S. Modélisation numérique et expérimentale des houles bidimensionnelles en zone côtière. Ph.D. Thesis, Université Joseph Fourier—Grenoble I, France, 1994.

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

50. 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.

51. Munk WH. The solitary wave theory and its application to surf problems. *Annals of the New York Academy of Sciences* 1949; **51**:376–424.