R Cienfuegos

E Barth Elemy
email: eric.barthelemy@hmg.inpg.fr

P Bonneton
email: p.bonneton@epoc.u-bordeax1.fr

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

Keywords: Boussinesq-type equations, Serre equations, nite volume method, compact schemes

A high-order nite 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.

INTRODUCTION

In recent years, important research eorts 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 [START_REF] Boussinesq | Th eorie 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[END_REF], Korteweg and de Vries [START_REF] Korteweg | On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves[END_REF], Serre [START_REF] Serre | Contribution al ' etude des ecoulements permanents et variables dans les canaux[END_REF] and Peregrine [START_REF] Peregrine | Long waves on beaches[END_REF], 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 [START_REF] Lynett | A two-layer approach to wave modelling[END_REF][START_REF] Madsen | Boussinesq formulation for fully non-linear and extremely dispersive water waves[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF][START_REF] Meftah | Linear analysis of a new type of extended Boussinesq model[END_REF][START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF][START_REF] Witting | A unied model for the evolution of nonlinear water waves[END_REF], and to describe surf-zone dynamics by adding extra wave-breaking terms [START_REF] Kennedy | Boussinesq modelling of wave transformation, breaking and runup. I: 1D[END_REF][START_REF] Sch Aer | A Boussinesq model for waves breaking in shallow water[END_REF][START_REF] Veeramony | The ow in surf-zone waves[END_REF][START_REF] Zelt | The run-up of nonbreaking and breaking solitary waves[END_REF] and removing the originally embedded weakly nonlinear hypothesis [START_REF] Seabra-Santos | Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle[END_REF][START_REF] Seabra-Santos | Etude th eorique et exp erimentale des domaines de validit e des th eories d' evolution des ondes en eau peu profonde[END_REF][START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]. Even though those important improvements have provided a powerful set of mathematical tools to describe the hydrodynamics of coastal zones, the question of constructing ecient, stable and accurate numerical schemes to solve them have received considerably less attention.

Numerical solutions of Boussinesq-type equations have been mostly tackled with nite difference methods [START_REF] Abbott | Accuracy of short-wave numerical model[END_REF][START_REF] Do Carmo | Numerical solution of the generalized Serre equations with the McCormack nite-dierence scheme[END_REF][START_REF] Fuhram | Numerical solution of fully non-linear and highly dispersive Boussinesq equations in two horizontal dimensions[END_REF][START_REF] Gobbi | A fully non-linear Boussinesq model for surface waves. Part 2. Extension to O(kh) 4[END_REF][START_REF] Shi | A fully nonlinear Boussinesq model in generalized curvilinear coodinates[END_REF][START_REF] Su | On head-on collisions between two solitary waves[END_REF][START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF], and to some extent, using nite element techniques [START_REF] Do Carmo | Surface waves propagation in shallow water: a nite element model[END_REF][START_REF] Walkley | A nite element method for the one-dimensional extended Boussinesq equations[END_REF][START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF][START_REF] Woo | Finite-element model for modied Boussinesq equations. II: applications to nonlinear harbor oscillations[END_REF]. On the contrary, the application of the nite volume method to discretize them has been only reported very recently [START_REF] Bradford | Finite-volume models for unidirectional, nonlinear, dispersive waves[END_REF][START_REF] Cienfuegos | R esolution num erique en volumes nis d'un syst eme d' equations de Serre etendu[END_REF][START_REF] Erduran | Hybrid nite-volume nite-dierence scheme for the solution of Boussinesq equations[END_REF][START_REF] Stansby | Solitary wave run up and overtopping by a semi-implicit nite-volume shallow-water Boussinesq model[END_REF]. Even though nite volume techniques have been successfully applied to the strictly hyperbolic dispersionless set of nonlinear shallow water equations (NSWE) using shock-capturing solvers (see Reference [START_REF] Toro | Shock-Capturing Methods for Free-Surface Flows (1st edn)[END_REF] 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 dicult. However, shock-capturing methods may provide an ecient and elegant way to numerically integrate Boussinesq equations but novel strategies are to be developed in order to adequately handle the mixed hyperbolicparabolic mathematical nature of the system. Several authors are working in that direction, developing in particular hybrid or operator splitting techniques where convective hyperbolictype terms are treated using shock-capturing approaches while high-order dispersive terms are handled with help of standard nite dierence approximations (e.g. References [START_REF] Erduran | Hybrid nite-volume nite-dierence scheme for the solution of Boussinesq equations[END_REF][START_REF] Bernetti | An operator-splitting method for long waves[END_REF]). 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 ux reconstruction strategy would be the most appropriate in the framework of dispersive wave propagation problems (see, for instance, Reference [START_REF] Erduran | Hybrid nite-volume nite-dierence scheme for the solution of Boussinesq equations[END_REF]).

In the present work, a nite volume technique is applied to the fully nonlinear and weakly dispersive set of Serre equations rst derived by Seabra-Santos et al. [START_REF] Seabra-Santos | Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle[END_REF]. The inclusion of a dispersion correction term following the lines given by Madsen and Sch aer [START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF] allows for an extension of the range of application of this set into moderately deep waters. The choice of integrating equations with a nite 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 [START_REF] Fj | Contribution al ' etude des ondes de gravit e bidimensionnelles en eau peu profonde[END_REF]. 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 ux. We recall that recent applications of nite 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 [START_REF] Bradford | Finite-volume models for unidirectional, nonlinear, dispersive waves[END_REF][START_REF] Bernetti | An operator-splitting method for long waves[END_REF] and a low-order Boussinesq model in References [START_REF] Erduran | Hybrid nite-volume nite-dierence scheme for the solution of Boussinesq equations[END_REF][START_REF] Stansby | Solitary wave run up and overtopping by a semi-implicit nite-volume shallow-water Boussinesq model[END_REF]). Similarly, most of discretization strategies implemented in previous nite volume approaches appears to be 2nd-order accurate in space. On the contrary, our aim here is to develop an ecient high-order nite 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 [START_REF] Whitham | Linear and Nonlinear Waves (1st edn)[END_REF]. Nevertheless, our numerical approach borrows some recently developed ideas from modern computational uid dynamics (CFD) which make high-order accuracy and high spectral resolution possible at an aordable computational cost. Such techniques are based on a compact reconstruction of cell face values at a high-order by solving implicit algebraical equations [START_REF] Kobayashi | On a class of Pad e nite volume methods[END_REF][START_REF] Lacor | A nite volume formulation of compact central schemes on arbitrary structured grids[END_REF]. Indeed, 4th-order accuracy is achieved by solving an implicit tridiagonal system. Similarly, spatial derivatives are estimated with a 4thorder truncation error using compact dierencing [START_REF] Lele | Compact nite dierence schemes with spectral-like resolution[END_REF]. 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 nite volume scheme. A comprehensive linear numerical analysis of the method is given in Section 4 and additional information on the performance of the nite volume scheme is obtained with help of numerical experiments in Section 5. In particular, the inuence that strong nonlinearities may have on stability and accuracy is investigated.

GOVERNING EQUATIONS

Serre equations describe the 2D, irrotational and shallow water ow of an incompressible and inviscid uid 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 elemy in Reference [START_REF] Barth | Nonlinear shallow water theories for coastal waves[END_REF].

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

w u ∼ (1)
where u and w are, respectively, the horizontal and vertical uid velocities (∼ stands for order of magnitude). The bottom boundary condition suggests that bottom slopes x should be of order , since w=u = x at the bottom. Obviously, bottom horizontal scales are not disconnected from the scales of the waves. In other words, the ow 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 [START_REF] Seabra-Santos | Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle[END_REF]. This settled, the scaling detailed in Reference [START_REF] Barth | Nonlinear shallow water theories for coastal waves[END_REF] leads to the following dimensionless forms of continuity, irrotationality, momentum equations and boundary conditions:

@u * @x * + @w * @z * = 0 (2) @u * @z * -2 @w * @x * = 0 (3) @u * @t * + 2 @ @x * (u * 2 )+ 2 @ @z * (u * w * )=- @p * @x * (4)
2 @w * @t * + 2 2 u * @w * @x * + 2 2 w * @w * @z * = -@p * @z * -1 ( 5)

w * = u * @ * @x * at the bottom z * = * (6)
w * = @ * @t * + u * @ * @x * at the free surface z * = * (7)

p * = 0 at the free surface z * = * (8)
where * indicates dimensionless variables. Here, p * represents the pressure, * is the free surface elevation and = 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: * (x * ;z * ;t * )= @w * @t * + u * @w * @x * + w * @w * @z * (9)

which is related to the vertical acceleration of uid particles.

Recall, we seek a set of equations for the depth-averaged horizontal velocity u * and the water depth h * = * - * . Depth averaging a function f * results from applying the following operator:

f * = 1 h * * * f * dz * (10)
Depth-averaging continuity equation yields

@ * @t * + @ @x * [h * u * ] = 0 (11)
Depth averaging the x-momentum equation ( 4) and following the derivation along similar lines as in Reference [START_REF] Barth | Nonlinear shallow water theories for coastal waves[END_REF] gives

h * @u * @t * + h * u * @u * @x * + h * @ * @x * + 2 @ @x * * * z * * dz * + 2 @ * @x * = - @ @x * * * (u * 2 -(u * ) 2 )dz * (12)
So far no approximation has been invoked, but the integrals in ( 12) cannot be reduced to partial dierential expressions without approximating the integrands. We aim at obtaining a set of equations comprising terms at most of O( 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 * *

(u * 2 -(u * ) 2 )dz * = h * (u * 2 -(u * ) 2 ) ( 13
)
This term has given rise in the past to discussions as to whether it is of O( 4 ) or not for uneven bathymetries. Dingemans [START_REF] Dingemans | Water Wave Propagation Over Uneven Bottoms (1st edn)[END_REF]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 4 . Using series expansions for u * and w * given in Appendix A we have,

u * (x * ;z * ;t * )=u * 0 (x * ;t * )+ 2 (z * - * )u * 2 (x * ;t * )+O( 4 ) (14)
w * (x * ;z * ;t * )= @ * @x * u * 0 (x * ;t * ) -(z * - * ) @u * 0 @x * + 2 (z * - * ) w * 2 (x * ;t * )+O( 4 ) (15)
showing in addition that u * and w * also have even series expansions in for uneven bathymetries. Using [START_REF] Zelt | The run-up of nonbreaking and breaking solitary waves[END_REF], clearly u * 2 -(u * ) 2 = O( 4 ) without having to express u * 0 and u * 2 in terms of u * as done by Dingemans [START_REF] Dingemans | Water Wave Propagation Over Uneven Bottoms (1st edn)[END_REF].

Left-hand side integrals in [START_REF] Sch Aer | A Boussinesq model for waves breaking in shallow water[END_REF] are evaluated with the above series expansions expressed in terms of u * ,

u * (x * ;z * ;t * )=u * (x * ;t * )+O( 2 ) ( 16
)
w * (x * ;z * ;t * )= @ * @x * u * -(z * - * ) @u * @x * + O( 2 ) (17)
Retaining all terms up to O( 2 ) leads to the following set of equations expressed in terms of dimensional variables and unknowns:

t +(h u) x = 0 ( 18
)
h u t + h u u x + gh x +[h 2 ( 1 3 P + 1 2 Q)] x + x h( 1 2 P + Q) = 0 ( 19
) P = -h[ u xt + u u xx -( u x ) 2 ] ( 20
) Q = x ( u t + u u x )+ xx u 2 (21)
where P + Q and Q represent, respectively, the vertical uid acceleration at the free surface and at the bottom (see Appendix A). This fully nonlinear set of equations contains all order of magnitude terms because no hypothesis have been introduced concerning nonlinearities (i.e. ∼O( 1)). Some of these terms are related to the convective vertical acceleration neglected in low-order Boussinesq equations (e.g. Reference [START_REF] Peregrine | Long waves on beaches[END_REF]). 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 [START_REF] Guizien | Accuracy of solitary wave generation by a piston wave maker[END_REF], expressed as

h(x; t)=h 0 + a sech 2 [K(x -Ct)] (22) u = C 1 - h 0 h (23)
K = 3a 4h 2 0 (h 0 + a) ;C = C 0 1+ a h 0 and C 0 = gh 0 ( 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 nite 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 [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF][START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF]. A Pad e (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)= -h 2 [ u xxt +( u u x ) xx + g(h + ) xxx ] ( 25
)
where is an adjustable parameter (∼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 [START_REF] Do Carmo | Surface waves propagation in shallow water: a nite element model[END_REF]. This convenient mathematical form rst given by Seabra-Santos [START_REF] Fj | Contribution al ' etude des ondes de gravit e bidimensionnelles en eau peu profonde[END_REF] reads

h t + F x = 0 ( 26
)
q t + G x = S (27) q =(1+r) u - 1 3h (h 3 u x ) x -2 u xx (28)
with ux functions

F = h u (29)
G = q u + g(h + ) -1 2 u 2 (1 + 2 x )+( x u -1 2 h u x )h u x -2 [ u 2 x + u u xx + g(h + ) xx ] (30)
and the non-conservative right-hand side term S (analogous to a source term) in [START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF] has the following expression:

S = -2 x [ u u x + g(h + ) x ] x (31)
The latter contains a small correction for the dispersive properties of the system. Moreover, r =(h x + x ) x + 1 2 h xx , and thus traduces some bathymetric eects. In what follows, overbars will be dropped for convenience.

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 dierent wave modes, without damping physical ones nor introducing spurious short-wave noise. In this context we will develop a 4th-order compact nite 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.

Finite volume method for spatial discretization

The computational domain is schematically described in Figure 2 where i and j are control volumes centred around i-nodes (i =0; 1;:::;I -1) and around j-nodes (j =0; 1;:::;I),

Left Boundary

Right Boundary respectively, i.e.

i = {x ∈ [x j ;x j+1 ]} (32) j = {x ∈ [x i-1 ;x i ]} ( 33
)
where I is the number of i control volumes used to discretize the physical domain. Nodal coordinates are dened as x i =x=2+ix and x j = jx, where 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 [START_REF] Shi | A fully nonlinear Boussinesq model in generalized curvilinear coodinates[END_REF][START_REF] Durran | Numerical Methods for Wave Equations in Geophysical Fluid Dynamics[END_REF]. Equations ( 26) and ( 27) can be easily integrated over control volumes i as follows:

@ @t xj+1 xj h dx + F(x j+1 ;t) -F(x j ;t) = 0 (34) @ @t xj+1 xj q dx + G(x j+1 ;t) -G(x j ;t)= xj+1 xj S dx (35)
where ux functions are evaluated at cell faces of each control volume (j-nodes). Unlike previous equations, the auxiliary relation [START_REF] Woo | Finite-element model for modied Boussinesq equations. II: applications to nonlinear harbor oscillations[END_REF] which links variable q, to the physical variables h and u will be integrated over staggered control volumes j :

xi xi-1 (1 + r)u - 1 3h (h 3 u x ) x -2 u xx dx = xi xi-1 q dx (36)
Numerical approximations will be introduced in the following to evaluate ux functions at cell faces and volume integrals.

In continuity and momentum equations [START_REF] Bernetti | An operator-splitting method for long waves[END_REF] and [START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF], 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 over control volume i is written as

i = 1 x xj+1 xj dx (37)
Discretized counterparts of Equations ( 34) and [START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF] are then obtained by performing the nite volume integration over the whole domain:

d ĥn i dt = - 1 x (F n j+1 -F n j ) for i =0; 1;:::;I -1 and j = i ( 38
)
dq n i dt = Ŝn i - 1 x (G n j+1 -G n j ) for i =0; 1;:::;I -1 and j = i (39)
where we have introduced for simplicity the notation n j =(x j ;t n ), with the discrete time coordinate t n = n t (n =0; 1;:::;N) and t is the time step. The non-conservative term will be approximated using a centred-in-space approach where Ŝi = 1 x

xj+1 xj S dx ≃- 2 x i ( x ) i ( Sn j+1 -Sn j ) (40)
and S = uu x + g(h + ) x . This discretized expression will introduce only slight errors in the overall numerical approximation since this term represents a small correction of O( 4 ) for dispersive properties. On the other hand, evaluation of ux 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 ĥi and qi . Here, cell face values will be estimated at each time step using the implicit 4th-order interpolation technique developed in References [START_REF] Kobayashi | On a class of Pad e nite volume methods[END_REF][START_REF] Lacor | A nite volume formulation of compact central schemes on arbitrary structured grids[END_REF]. This approach results in the following system of I -1 equations for internal nodal points: 1 4 j-1 + j + 1 4 j+1 = 3 4 ( i-1 + i ) for j =1; 2;:::;I -1 and i = j

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 [START_REF] Cienfuegos | R esolution num erique en volumes nis d'un syst eme d' equations de Serre etendu[END_REF]). Thus, computational cost is signicantly reduced when implementing the present compact scheme.

Evaluation of ux 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 [START_REF] Fj | Contribution al ' etude des ondes de gravit e bidimensionnelles en eau peu profonde[END_REF]. Dierent integrals appearing in this equation are approximated as follows:

1 x xi xi-1 (1 + r)u dx ≃ (1 + r j )û j - 1 x xi xi-1 1 3h (h 3 u x ) x dx ≃- 1 3xh j [h 3 i (u x ) i -h 3 i-1 (u x ) i-1 ] - 1 x xi xi-1 2 u xx dx ≃- x 2 j [(u x ) i -(u x ) i-1 ] 1 x xi xi-1 q dx =q j
At this point, some care must be taken for the reconstruction of depth-averaged velocity values at j-nodes using Equation [START_REF] Fj | Contribution al ' etude des ondes de gravit e bidimensionnelles en eau peu profonde[END_REF]. 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 ûj and diusive uxes 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 i and j , i = 2 3 i + 1 6 ( j + j+1 )+O(x 4 ) (42)

j = 2 3 j + 1 6 ( i-1 + i )+O(x 4 ) (43)
and Equation ( 41), provides the following relation: j = 1 24 ( j-1 + 22 j + j+1 )+O(x 4 ) [START_REF] Durran | Numerical Methods for Wave Equations in Geophysical Fluid Dynamics[END_REF] 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 diusive uxes. One way to preserve a tridiagonal structure is to use a deferredcorrection approach, which is expressed as [START_REF] Ferziger | Computational Methods for Fluid Dynamics (3rd edn)[END_REF] (u

x ) i = (u j+1 -u j ) x + (u j-1 -27u j +27u j+1 -u j+2 ) 24x - (u j+1 -u j ) x old ( 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 [START_REF] Fj | Contribution al ' etude des ondes de gravit e bidimensionnelles en eau peu profonde[END_REF] 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 satised between 'old' and 'new' u j -values, the overall approximation for u x will be 4thorder accurate. In practice only few iterations will be needed to obtain sucient accuracy.

The convergence criteria for xed-point iterations is based on a L 1 norm and reads

u j -{u j } old 1 u j 1 6 (46)
where

u j =[u j =2 u j =3 ••• 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 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:

(1 + r j ) 24

- 1 x 2 h 3 i-1 3h j + 2 j u j-1 + 22(1 + r j ) 24 + 1 x 2 (h 3 i-1 + h 3 i ) 3h j +2 2 j u j + (1 + r j ) 24 - 1 x 2 h 3 i 3h j + 2 j u j+1 = 1 24 (q j-1 +22q j + q j+1 )+ 1 x 2 h 3 i 3h j + 2 j U old i - 1 x 2 h 3 i-1 3h j + 2 j U old i-1
for j =1; 2;:::;I -1 and i = j

where

U old i = 1 24 {u j-1 -3u j +3u j+1 -u j+2 } old (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 [START_REF] Press | Numerical Recipies in Fortran (2nd edn)[END_REF]. 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 ux function G j , the source term Ŝi and the dependent variable r j at each time step. Similarly, water depth values at i-nodes must be interpolated to compute matrix coecients in Equation ( 47) while bottom function and its derivatives must be prescribed at i and j nodes. The latter is not dicult 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 ĥi and h j using relation [START_REF] Dingemans | Water Wave Propagation Over Uneven Bottoms (1st edn)[END_REF]. 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 nite dierence compact scheme as described in Reference [START_REF] Lele | Compact nite dierence schemes with spectral-like resolution[END_REF]. These compact formulae have better spectral properties for shortlength waves than their explicit 4th-order nite dierence 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 nite volume method applied to the quasi-conservative set of equations derived here has allowed us to lower the highest xderivative order by one when compared to other equivalent set of equations (e.g. References [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF][START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF]). Derivative estimates at internal j-nodes are obtained from cell face values using the following implicit formulae [START_REF] Lele | Compact nite dierence schemes with spectral-like resolution[END_REF]:

1 4 ( x ) j-1 +( x ) j + 1 4 ( x ) j+1 = 3 4x
( j+1 -j-1 ) for j =1; 2;:::;I-1 and i = j (49)

1 10 ( xx ) j-1 +( xx ) j + 1 10 ( xx ) j+1 = 6
5x 2 ( j+1 -2 j + j-1 ) for j =1; 2;:::;I -1 and i = j

Again, special care must be taken at boundary nodes as discussed in References [START_REF] Lele | Compact nite dierence schemes with spectral-like resolution[END_REF][START_REF] Sengupta | Fundamentals of Computational Fluid Dynamics (1st edn)[END_REF]. It is worth to point out that these compact forms can be veried manipulating relation [START_REF] Durran | Numerical Methods for Wave Equations in Geophysical Fluid Dynamics[END_REF]. Spatial discretization of Equations ( 34)-( 36) using the present nite volume technique allows us to write the original system of partial dierential equations (PDE) as the following set of ordinary dierential equations (ODE),

d xn i dt = f n i (51) A(h n j ) • u n j = B • q n j + C(h n j ) •{u n j } old (52)
where

xn i = ⎡ ⎣ ĥn i qn i ⎤ ⎦ and f n i = ⎡ ⎣ 0 Ŝn i ⎤ ⎦ - F n i G n i
Here, bold letters represent vectors containing cell-averaged values of associated variables over the whole discrete domain. For instance, a generic vector n i contains cell-averaged values over control volumes i of function evaluated at time t = t n . Similarly, the notation

n i =(1=x)( n j+1 -n j )
indicates a cell face ux dierence vector, and A, B and C are matrices obtained from the discretized equation [START_REF] Press | Numerical Recipies in Fortran (2nd edn)[END_REF]. In the next subsection, we will use an explicit 4th-order Runge-Kutta method to numerically integrate the system of ODE just derived.

Numerical time stepping

The implementation of a high-order numerical ODE solver for time integration is justied since spatial derivatives are discretized using a high-order scheme as well. The nite 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 ecient when compared to other high-order schemes [START_REF] Durran | Numerical Methods for Wave Equations in Geophysical Fluid Dynamics[END_REF]. In the present context, the standard four-stage Runge-Kutta method can be stated as follows:

First estimation at t = t n+1=2 p (1) =t

f n i ; xn+(1=2) i = xn i + 1 2 p (1) (53a) u n+(1=2) j = A(h n+(1=2) j ) -1 • [B • q n+(1=2) j + C(h n+(1=2) j ) •{u n+(1=2) j } old ] (53b) Second estimation at t = t n+1=2 p (2) =t f n+(1=2) i ; xn+(1=2) i = xn i + 1 2 p (2) (53c) u n+(1=2) j = A(h n+(1=2) j ) -1 • [B • q n+(1=2) j + C(h n+(1=2) j ) •{u n+(1=2) j } old ] (53d) First estimation at t = t n+1 p (3) =tf n+(1=2) i ; xn+1 i = xn i + 1 2 p (3) (53e) u n+1 j = A(h n+1 j ) -1 • [B • q n+1 j + C(h n+1 j ) •{u n+1 j } old ] (53f)
Second and nal estimation at t = t n+1

p (4) =t f n+1 i ; xn+1 i = xn i + 1 6 (p (1) +2p (2) +2p (3) + p (4) ) (53g) u n+1 j = A(h n+1 j ) -1 • [B • q n+1 j + C(h n+1 j ) •{u n+1 j } old ] (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 [START_REF] Barth | Nonlinear shallow water theories for coastal waves[END_REF] before solving the equation for depthaveraged 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 rst guess. A relative tolerance for iteration errors, ,o f1 0 -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 [START_REF] Press | Numerical Recipies in Fortran (2nd edn)[END_REF]). The number of operations needed to solve each matrix equation will then be linearly proportional to the number of unknowns in the system.

SPECTRAL ANALYSIS AND LINEAR STABILITY

Important information about the general performance of our nite 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 [START_REF] Sengupta | Fundamentals of Computational Fluid Dynamics (1st edn)[END_REF]. Our main goal here is to nd 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 at bottoms.

When using relation [START_REF] Barth | Nonlinear shallow water theories for coastal waves[END_REF] in Equations ( 38)-( 39) and assuming that xed-point iterations have converged (i.e. {u n j } old = u n j ) in [START_REF] Press | Numerical Recipies in Fortran (2nd edn)[END_REF], the linearized version of the discretized system of equations can be recasted as

1 4 dh n j-1 dt + d h n j dt + 1 4 d h n j+1 dt = - 3 4x (F n j+1 -F n j-1 ) (54) 1 4 d q n j-1 dt + d q n j dt + 1 4 d q n j+1 dt = - 3 4x (G n j+1 -G n j-1 ) (55) 2 u n j-2 + 1 - 28 2 u n j-1 + 22 + 54 2 u n j + 1 - 28 2 u n j+1 + 2 u n j+2 = q n j-1 +22q n j + q n j+1 (56)
where

F n j = h 0 u n j ;G n j = gh n j -h 2 0 g(h xx
) n j ; =(1=3+); =x=h 0 and h 0 = || is the still water depth of the horizontal bottom ume. 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 rst spatial derivative estimate (49).

Semi-discrete system of equations and numerical dispersion

Spectral accuracy of the spatial discretization is studied rst by decomposing the numerical solution into nite Fourier series. It is important to point out that spectral accuracy is related to the ability of the scheme to resolve the dierent length scales present in the problem [START_REF] Sengupta | Fundamentals of Computational Fluid Dynamics (1st edn)[END_REF]. 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:

⎡ ⎢ ⎢ ⎣ h n j q n j u n j ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ h q u ⎤ ⎥ ⎥ ⎦ e (jx-!rt) (57)
where = (-1) is the imaginary constant, h , q and u are Fourier amplitudes, is the wave number of the individual component and ! 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 ! r . Before doing so, we need to nd the spectral expression for the second spatial derivative that appears in the ux function G n j . The numerical estimation of this derivative has the form (h xx ) j = -2 h e (jx-!r t) , where by substitution in the implicit dierence formula (50) the associated numerical wave number is written as

(x) 2 =12 (1 -cos x) (5 + cos x) (58)
Similarly, substituting the Fourier components (57) into Equation (56) gives the following relation between q and u :

q = u (59)
where

=1+ (1=3+) 2 (27 + cos 2x -28 cos x) (11 + cos x) (60)
Finally, if we use the Fourier components (57) and relation (59) in continuity and momentum equations ( 55)-(56), the compact nite volume scheme takes the following form in the Fourier space:

(K + I ! r ) • h u = 0 (61)
where

K = 3 0 K 1 K 2 0
is the spatial discretization operator whose components read,

K 1 = - sin x (2 + cos x) (62)
K 2 = K 1 [ 2 + (x) 2 ] C 0 x 2 (63)
I is the identity matrix and C 0 = gh 0 . Non-trivial solution of Equation (61) provides the discrete dispersion relation associated to our compact nite volume method which can be written in an non-dimensional form as

C 2 fv C 2 0 = 1 C 2 0 ! r 2 = - 3 x 2 K 2 1 [ 2 + (x) 2 ] 2 = 3 sin x x(2 + cos x) 2 × (11 + cos x)[ 2 + (x) 2 ] [ 2 (11 + cos x)+(1=3+)(27 + cos 2x -28 cos x)] (64)
with C fv being the numerical phase speed associated to the nite volume scheme. Therefore, the non-dimensional discrete dispersion relation depends on the wave number, x, and on the parameter =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 → 0 (which can be viewed as the deep water limit) and as →∞ (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 [START_REF] Barth | Nonlinear shallow water theories for coastal waves[END_REF]

C th C 0 = 1+(h 0 ) 2 1+(1=3+)(h 0 ) 2 1=2 = ⎡ ⎢ ⎢ ⎢ ⎣ 1+ x 2 1+(1=3+) x 2 ⎤ ⎥ ⎥ ⎥ ⎦ 1=2 ( 65
)
where C 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 →∞ C fv C th = 3 sin x x(2 + cos x) =1- 1 180 (x) 4 + O[(x) 6 ] (66)
proving that the numerical scheme mimics the exact relation with 4th-order accuracy in x.

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

lim →0 C fv C th = 3 sin x x (2 + cos x) 6(11 + cos x) (13 -cos x)(5 + cos x) 1=2 =1- 11 2880 (x) 4 + O[(x) 6 ] (67)
Even though it is generally acknowledged that numerical dispersion constitutes an important property of any scheme, dispersion relations associated to previously published nite dierence or nite element high-order Boussinesq-type models have rarely been reported. Moreover, it is commonly argued that numerical dispersion can be avoided if O( 2 ) dispersive terms are discretized with 2nd-order truncation error approximations, while convective or 1st-order terms are approximated through 4th-order centred nite dierence formulae (e.g. References [START_REF] Abbott | Accuracy of short-wave numerical model[END_REF][START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF]). Even if from a practical point of view the latter could be sucient for well-resolved waves, the numerical scheme will only be 2nd-order accurate in 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 ux function G n j of Equation ( 55) and in the linearized version of the auxiliary equation ( 36) by using explicit 2nd-order nite dierence approximations. It is worth emphasizing that no direct comparison can be made with the nite dierence technique used in References [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF][START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF] or with nite element models based on similar mixed-order approaches [START_REF] Walkley | A nite element method for the one-dimensional extended Boussinesq equations[END_REF][START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF] since our numerical resolution is based on a compact nite 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:

( xx ) j = 1 x 2 ( j+1 -2 j + j-1 ) (68)
and the corresponding spectral relation is found to be (after substitution of a Fourier component),

( ′ x) 2 = 2(1 -cos x) (69)
Similarly, using 2nd-order approximations for high-order derivatives in the linearized version of auxiliary equation [START_REF] Fj | Contribution al ' etude des ondes de gravit e bidimensionnelles en eau peu profonde[END_REF] it can be shown that

′ =1+ (1=3+) 2 24(1 -cos x) (11 + cos x) (70)
Substituting relations (69) and ( 70) in (64) provides the following form for the discrete dispersion relation associated to the mixed approach:

C ′ fv C th =1- 2 48 [( -1=3) + ( +1=3) 2 ] [1 + 2 ][1+( +1=3) 2 ] (x) 2 + O[(x) 4 ] (71)
Therefore, using 2nd-order nite dierence formulae to evaluate dispersive terms results in an overall reduction in the spectral resolution of the method. However, the coecient that multiplies the O (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 aected. The latter constitutes the main argument behind the mixed approach used in several nite dierence and nite 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,

1 C 0 @! @ = C g C 0 = @ @(x) (x) C C 0 (72)
Again, the general discrete relation will depend on x and 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 while in the deep-water case we have

lim →∞ C gfv C gth =1- 1 36 (x) 4 + O[(x) 6 ] (73)
lim →0 C gfv C gth =1- 11 576 (x) 4 + O[(x) 6 ] (74)
showing that the truncation error is still 4th-order in 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 x, except in the shallow water limit ( →∞) 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 =1=15 in Figure 4. It can be seen that when →∞, 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 x →∞ which is of course rather problematic. Conversely, for a xed spatial grid resolution, the modied nite 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 eects become predominant and the model is intended to perform better. On the contrary, the parameter has only a marginal inuence on phase speed errors when the full 4th-order model is used. Thus, the full 4th-order compact nite volume method is more robust, in the sense that phase speed estimates are less sensitive to changes in the parameter, than the mixed-order version. For comparison, numerical results obtained with =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 x≃0:325 (L=x≃6:1) in the full 4th-order case, and when x≃0:202 (L=x≃9:9) in the mixed one. Hence, for a reasonable value of , 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 =x=h 0 is increased, however, the worst error occurs when x=h 0 → 0. In that limit the 0.5 % error level is reached for x≃0:152 (L=x≃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 =0:1 phase speed of well-resolved waves (i.e. x ≪ ) 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 dierence here is that one obtains negative values for group velocities in the case of poorly resolved waves (x¿0:67). On the other hand, in order to describe group velocities within a 0.5% error when =0:1, the full 4th-order nite volume method requires x≃0:213 (L=x≃9:4) while the mixed scheme reaches the same accuracy when x≃0:117 (L=x≃17:2). It is conrmed 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 (x=h 0 → 0) the mixed-order model reaches the 0.5% error in group velocities when x≃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 =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.

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:

⎡ ⎢ ⎢ ⎣ h n j q n j u n j ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ h q u ⎤ ⎥ ⎥ ⎦ e (jx-(!r+! i ) nt) = ⎡ ⎢ ⎢ ⎣ h q u ⎤ ⎥ ⎥ ⎦ || n e (jx-!r nt) (75)
where ! r and ! i correspond to real and imaginary parts of the frequency and || =e ! i t is the amplication 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 ! i t → 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 amplication=dissipation). Numerical properties of the discrete scheme are now driven by three free parameters, x, =x=h 0 and the Courant number, C r = C 0 t=x.

It is worth to point out that analogous linear stability analysis have been performed recently on a nite element resolution of an equivalent Boussinesq set of equations by Woo and Liu [START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF], and on several nite volume resolutions of a KdV-like system of equations by Bradford and Sanders [START_REF] Bradford | Finite-volume models for unidirectional, nonlinear, dispersive waves[END_REF]. 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,

(G(K) -I||e -!r t ) • h u = 0 ( 76
)
where G(s)= 4 p =0 1 p! (t • s) p is the Runge-Kutta time stepping operator and K is the spatial discretization operator introduced in the previous subsection. Non-trivial solution of Equation (76) yields

1+ 9 2 t 2 K 1 K 2 + 27 8 t 4 (K 1 K 2 ) 2 -||e -!r t 2 = 9 4 t 2 K 1 K 2 2+3 t 2 K 1 K 2 2 ( 77
)
Invoking denitions ( 62) and (63) we write

K = t 2 K 1 K 2 = - sin 2 x (2 + cos x) 2 [ 2 + (x) 2 ] 2 C 2 r (78)
Thus, equating real and imaginary parts in Equation (77) results in the following:

|| cos ! r t = R =1-9 2 |K| + 27 8 |K| 2 (79a) || sin ! r t = ± I = ± 3 2 (2 -3|K|)|K| 1=2 (79b)
Squaring and adding previous relations gives the amplication factor, while dividing (79b) by (79a) provides the discrete dispersion relation of the numerical scheme, i.e.

|| = R 2 + I 2 (80a) tan ! ± r t = ± I R (80b)
Although the dispersion relation (80b) has two roots, in the frequency range of resolvable waves (! r t ∈ [0;]) only the root with a positive imaginary part is meaningful. Because of the change in sign of the function I when |K| =2=3, the mode associated with the rst root, ! + r t, remains in the plane of resolvable waves as far as |K|62=3. On the other hand, the second mode associated with ! - r 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 dened by the inequality |K|62=3. Numerical stability is studied now by looking at the amplication factor dened in (80a). The scheme produces non-amplifying solutions as far as ||61, so it is straightforward to show that if

|K|6 8 9 (81)
the RK4 compact nite volume scheme is stable. In the general case, the stability condition on the Courant number will be a function of , and x, but it can be shown that the most stringent situation occurs when →∞. Replacing relations (58) and ( 60) in (78) provides the following stability condition in that limit: C r 6 2

√ 2 3 (2 + cos x) sin x 62 2 3 ( 82
)
The numerical scheme is thus stable if the Courant number is chosen to be less than 1.633 no matter what values of and are used and for all resolvable wave numbers. Nevertheless, the stability limit can be moved further on by decreasing the parameter value , i.e. by improving the spatial grid resolution for a xed h 0 . In the limit of very good grid resolution ( → 0), the stability condition (81) takes the following form:

C r 6 2 √ 3 9 (2 + cos x) sin x (1=3+) (5 + cos x)(13 -cos x) (11 + cos x) 61:603 (1=3+) (83)
Therefore, when taking =1=15 the method remains stable for all resolvable wave numbers, 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. → 0). On the other hand, for nite values of , 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 dierent Courant numbers and for a xed value of =0:2 (keeping =1=15). This particular -value has been chosen in order to perform a comparison with the numerical analysis presented in papers by Woo and Liu [START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF] and Bradford and Sanders [START_REF] Bradford | Finite-volume models for unidirectional, nonlinear, dispersive waves[END_REF]. It is seen in Figure 6(a) that when C r =1:0 the present compact nite volume scheme is stable and nearly neutral (i.e. ||≃1:0), thus introducing a clear improvement over the Petrov-Galerkin nite element model described in Reference [START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF] 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 nite volume schemes (Lax-Wendro, Fromm and Warming-Beam) and a centred 2nd-order nite dierence scheme used to solve a KdV-like system of equations was compared in Reference [START_REF] Bradford | Finite-volume models for unidirectional, nonlinear, dispersive waves[END_REF]. For an almost equivalent set of numerical parameters (C r =1:0 and =x=h 0 ≃0:18), Bradford and Sanders showed that only the Fromm and the Warming-Beam methods were stable but they appear 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 [START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF][START_REF] Bradford | Finite-volume models for unidirectional, nonlinear, dispersive waves[END_REF]. Moreover, using Maclaurin series to expand (80a) it can be shown that for small values of x,

lim →∞ || =1- C 6 r 144 (x) 6 + O[(x) 8 ] (84a) lim →0 || =1- C 6 r 144 (1=3+) 3 (x) 6 + O[(x) 8 ] (84b)
Therefore, the dissipative error for the present nite volume scheme is O[(x) 6 ], while linear nite volume schemes studied by Bradford and Sanders [START_REF] Bradford | Finite-volume models for unidirectional, nonlinear, dispersive waves[END_REF] only shown to be O[(x) 4 ] accurate.

Further information on the damping modes associated to our numerical model can be obtained from Figure 6(a) when looking at amplication 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 x≃0:662. Additional computations show that the numerical scheme is non-amplifying (stable) for =0:2 as far as C r ¡3:71. However, midrange frequencies begin to be heavily damped before reaching that limit (see Figure 6(a)).

As discussed in Reference [START_REF] Hu | Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics[END_REF], this dissipative behaviour is mostly due to the Runge-Kutta time stepping.

Phase speed errors obtained for =0:2 and for dierent 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 identied and coincides with the one observed for the amplication factor. It is also evident from this gure 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 x, it can be shown that a Maclaurin expansion of the dispersion relation yields

lim →∞ C fv C th =1- 1 180 1+ 3 2 C 4 r (x) 4 + O[(x) 6 ] (85a) lim →0 C fv C th =1- 1 2880 11+24 (1=3+) 2 C 4 r (x) 4 + O[(x) 6 ] (85b)
Therefore, phase errors dominate over dissipative errors in the present compact nite volume scheme. Linear nite volume schemes studied by Bradford and Sanders [START_REF] Bradford | Finite-volume models for unidirectional, nonlinear, dispersive waves[END_REF] 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 61: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 ;x) plane for dierent values of the -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 dierent 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 ≃6 when formulated in surface and horizontal velocity variables (e.g. References [START_REF] Gobbi | A fully non-linear Boussinesq model for surface waves. Part 2. Extension to O(kh) 4[END_REF][START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF]), and up to even higher values (≃25) if the vertical velocity is in addition retained as unknown (see Reference [START_REF] Madsen | Boussinesq formulation for fully non-linear and extremely dispersive water waves[END_REF]). 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 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¡(x) 0:5% ¡0:3 when ! i t → 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 x being less than 0:15 and this value xes the -parameter that will ensure that all physical wavelengths will be well represented (because x = ). If the physical range of the discretized Boussinesq-type equations is 6, then a value of =x=h 0 60:08 must be used in the mixed-order nite volume scheme. Similarly, if the physical range is moved further on to reach ∼ 25 as it is the case in Reference [START_REF] Madsen | Boussinesq formulation for fully non-linear and extremely dispersive water waves[END_REF], then the spectral resolution of the numerical scheme imposes the discretization parameter to be less than 0:02. It is evident from results presented in the previous subsection that, for group velocities, the ratio x=h 0 must be further reduced to achieve an accurate representation of this property. Moreover, for nite 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 specically, 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 nite dierence formulae can lead to overwhelming computational eorts 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 signicantly worse spectral properties than equivalent compact formulae [START_REF] Kobayashi | On a class of Pad e nite volume methods[END_REF][START_REF] Lacor | A nite volume formulation of compact central schemes on arbitrary structured grids[END_REF][START_REF] Lele | Compact nite dierence schemes with spectral-like resolution[END_REF].

Contour plots for the dissipative errors in the (C r ;x) plane are presented in Figure 7 for the full 4th-order scheme and the mixed-order approach. The -parameter is xed 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 ∼ 3 when =1=15 [START_REF] Barth | Nonlinear shallow water theories for coastal waves[END_REF]. This gure 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 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 nite dierence formulae to discretize dispersive terms results in a slight improvement in stability and similar conclusions were given in Reference [START_REF] Fuhram | Linear and non-linear stability analysis for nite dierence discretizations of high-order Boussinesq equations[END_REF]. 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 [START_REF] Sengupta | Fundamentals of Computational Fluid Dynamics (1st edn)[END_REF]. 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 ;x) plane in Figures 8 and9. 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 extent of the DRP region in terms of wave numbers 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 nite dierence 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 nite volume scheme. Therefore, signicantly 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 -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 =0:1 for numerical computations will ensure that the compact nite volume scheme used to discretize the set of PDE derived here is stable, non-damping and with a maximized DRP region.

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 aected 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 dierent (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 rened.

Discretization errors

Without dispersion correction terms (i.e. = 0) and for a at bottom, Equations ( 18)-( 19) have the solitary wave solution expressed by relations ( 22)- [START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF]. 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 denition used in Reference [START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF]. Hence, L is the distance over which =h 0 ¿0:001, where is the location of the free surface. The equivalent solitary wave period is then dened as T = L=C and can be computed invoking the denition for the phase speed of the soliton given in Equation [START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF]. The average numerical error is estimated using the L 1 norm and reads

L = j = I j =0 h n j -h e j 1 j = I j =0 h e j 1 ( 86
)
where h e and h n are, respectively, the exact and the numerically computed solution, I is the total number of nite 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 xed-point-iteration tolerance of =10 -4 to solve [START_REF] Press | Numerical Recipies in Fortran (2nd edn)[END_REF]. Spatial discretization error is studied using a small time step of t =0:002s 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 nite volumes per wavelength (i.e. decreasing ) are presented in Figure 10(a) where it can be seen that, except for the last two computed errors using the smallest , the ideal slope for a 4th-order spatial discretization is fairly well recovered. It is believed that for the two last computed values, round-o errors or some inuence of the time step (because it is not small enough for the prescribed value) are responsible for the slight divergence from the ideal slope. Nevertheless, it is widely conrmed that the nite volume discretization of the system of PDE is O[(x) 4 ] accurate.

An equivalent numerical test is performed on discretization errors due to time integration for a xed =0:1 and monotonically decreasing the Courant number. Associated numerical errors are presented in Figure 10(b) and compared against the theoretical 4th-order slope. Again, it is fairly well conrmed that the convergence error is O[(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 nite volume scheme is O[(x) 4 ; (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 x is slightly worse than the one observed in the moderately nonlinear case investigated before. The latter is not surprising since strong nonlinearities may aect the model performance. Nevertheless, the convergence rate is still well above O[(x) 3 ] for this challenging test and the computed slope only starts to diverge from the theoretical one when x=h 0 60:4. This trend, which has been also noticed in the previous example, is probably associated to the xed 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. 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 depthaveraged velocity u in matrix equation [START_REF] Press | Numerical Recipies in Fortran (2nd edn)[END_REF]. Moreover, the convergence rate was rather insensitive to changes in the xed-point-iteration tolerance as long as 610 -4 . Indeed, =10 -4 proved be a sucient choice in all computations.

Another important property that any numerical scheme used to solve Boussinesq-type equations must full 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 [START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF][START_REF] Walkley | A nite element method for the one-dimensional extended Boussinesq equations[END_REF][START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF] 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 =0:1. A perfect agreement between the analytical and numerical solution is noticed. It is remarkable how this result conrms 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 ∼ 0:5) was used in Reference [START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF] on the basis of a Taylor-Galerkin nite 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 [START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF][START_REF] Walkley | A nite element method for the one-dimensional extended Boussinesq equations[END_REF].

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

Nonlinear stability

The way by which nonlinearity aects 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 dierent relative amplitudes a=h 0 and for monotonically increasing Courant numbers until numerical instability shows up. The inuence of the correction dispersion term is also investigated because it was analytically demonstrated in Section 4 that this term may produce some destabilizing eects in the system (see inequality (83)). Therefore, three dierent -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 dierent 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 gure (when a=h 0 → 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 nite values of as discussed before. Nonetheless, for the two nite values of investigated here, numerical experiments conrm in a satisfactory way the linear stability analysis as depicted in Figure 12. Moreover, when the dispersion correction term is switched o the stability region is considerably increased as it was expected from the numerical analysis of the scheme. It is also noticed that for nite -values, nonlinearity does aect the stability of the scheme because the maximum stable Courant number is reduced. Nevertheless, numerical experiments show that the RK4 compact nite volume scheme remains stable, even for extremely nonlinear waves, as long as the Courant number is kept below a value of 1.8 when 61=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 nite -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 aected 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 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 conrm numerically that the stability limit may reach such high Courant numbers when the dispersion correction term was switched o.

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 nite volume methods very attractive. Indeed, recently developed compact strategies for cell-face reconstruction in the framework of nite 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 conrmed that the newly developed scheme is fully 4th-order in space and time while computational eorts 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 aordable 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 nite volume model [START_REF] Cienfuegos | R esolution num erique en volumes nis d'un syst eme d' equations de Serre etendu[END_REF]. Moreover, we believe that the particular choice of this compact nite 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 [START_REF] Lacor | A nite volume formulation of compact central schemes on arbitrary structured grids[END_REF][START_REF] Shukla | Derivation of high-order compact nite dierence schemes for non-uniform grid using polynomial interpolation[END_REF]).

On the other hand, the use of linear spectral analysis provides important information concerning 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 nite dierence formulae may lead to phase and group velocity errors which are dicult 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 wavelengths will certainly impose the use of a very ne 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 diculties encountered in the eld 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 further rening 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, ecient, 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 nite 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 [START_REF] Shi | A fully nonlinear Boussinesq model in generalized curvilinear coodinates[END_REF][START_REF] Woo | A Petrov-Galerkin nite element model for one-dimensional fully non-linear and weakly dispersive wave propagation[END_REF][START_REF] Gobbi | Wave evolution over submerged sills: tests of a high-order Boussinesq model[END_REF]). On the contrary, the present nite volume scheme has shown to be stable (i.e. for all wave numbers and spatial grid resolutions) if Courant number is xed 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 61:8 for a grid resolution of =0:1. Indeed, it was found analytically in Section 4 that the dispersion correction term as given in References [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF][START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF] introduces some kind of numerical instability in the system. This undesirable situation was also conrmed 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 [START_REF] Cienfuegos | A new wave-breaking parametrization for Boussinesq-type equations[END_REF]. Similarly, the numerical treatment of appropriate 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 eects), high-frequency waves can arise as a consequence of nonlinear interaction. In this context, the use of high-order lters to damp out spurious wave components is of common usage. The spectral analysis performed here gives insight of how this can be achieved without aecting the range of resolved physical wavelengths.

APPENDIX A

Series expansion of u * and w * are found using Equations ( 2), (3) and boundary condition [START_REF] Madsen | Boussinesq formulation for fully non-linear and extremely dispersive water waves[END_REF]. As a rst step assume the following expansion for both unknowns: where f 0 and g 0 are, respectively, the bottom horizontal and vertical velocities. Boundary condition [START_REF] Madsen | Boussinesq formulation for fully non-linear and extremely dispersive water waves[END_REF] imposes that,

g 0 = * x f 0 (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 *

x f n+1 -g n+1 =(f n ) x f n+1 + 2 *

x g n+1 = 2 (g n ) x

(A4) which at the lowest order writes, * x f 1 -g 1 =f 0x f 1 + 2 *

x g 1 = 2 g 0x (A5)

This linear set of equations has the solution

f 1 = 2 * x f 0x + 2 (f 0 * x ) x 1+ 2 ( * x ) 2
(A6)

g 1 = 2 * x (f 0 * x ) x -f 0x 1+ 2 ( * x ) 2 (A7)
Thus, if 2 ( * x ) 2 ≪ 1, the following expansions in terms of 2 hold f 1 = 2 [ * x f 0x +(f 0 * x ) x ]+O( 4 ) (A8)

g 1 = -f 0x + 2 [ * x (f 0 * x ) x + f 0x ( * x ) 2 ]+O( 4 ) (A9)
and a justication to ( 14) and ( 15) is provided.

Figure 1 .

 1 Figure 1. Denition sketch for wave propagation over uneven bathymetries. z = is the free surface elevation, z = is the bottom location, and h = -is the local water depth (here has a negative value).

Figure 2 .

 2 Figure 2. One-dimensional computational domain meshed with staggered control volumes i and j .

Figure 3 .

 3 Figure 3. Leading-order truncation error term vs = h 0 for, (-) =0, (•••) =1=30, (--) =1=20 and (-• -) =1=10: (a) full 4th-order compact nite volume scheme; and (b) mixed 4th-and 2nd-order discretization strategy.

Figure 4 .

 4 Figure 4. Ratio between computed and physical phase speed vs x for =1=15, when (-) →∞, (--) → 0 and (-• -) =0:1: (a) full 4th-order compact nite volume scheme; and (b) mixed 4th-and 2nd-order discretization strategy.

Figure 5 .

 5 Figure 5. Ratio between computed and physical group velocity vs x for =1=15, when (-) →∞, (--) → 0 and (-•-) =0:1: (a) full 4th-order compact nite volume scheme; and (b) mixed 4th-and 2nd-order discretization strategy.

Figure 6 .

 6 Figure 6. Numerical dissipation and phase speed error of the 4th-order compact nite volume scheme for =1=15, =0:2 at dierent Courant numbers: (-) Cr =1:0, (--) Cr =2:0, (-• -) Cr =3:0 and (-×-) Cr =4:0: (a) amplication factor; and (b) ratio between computed and physical phase speed.

Figure 7 .

 7 Figure 7. Contours of 0.5 and 5% dissipative errors for =1=15 and =0:1 in the (Cr;x) plane. ||60:95 in region I (damping), and ||¿1:05 in region II (amplifying): (a) full 4th-order compact nite volume scheme; and (b) mixed 4th-and 2nd-order discretization strategy.

Figure 8 .

 8 Figure 8. Contours of 0.5 and 5% phase speed errors for =1=15 and =0:1 in the (Cr;x) plane: (a) full 4th-order compact nite volume scheme; and (b) mixed 4th-and 2nd-order discretization strategy.

Figure 9 .

 9 Figure 9. Contours of 0.5 and 5% group velocity errors for =1=15 and =0:1i nt h e (Cr;x) plane: (a) full 4th-order compact nite volume scheme; and (b) mixed 4th-and 2nd-order discretization strategy.

Figure 10 .

 10 Figure 10. Discretization errors for a propagating solitary wave. (--) Ideal 4th-order slope; (•) computed errors when a=h 0 =0:2; (M) computed errors when a=h 0 =0:6: (a) convergence error in x using a xed t =0:002 s; and (b) convergence error in t using a xed =0:1.

Figure 11 .

 11 Figure 11. Strongly nonlinear solitary wave propagating over 50 equivalent wavelengths where a=h 0 =0:6, L =14:7m , T =3:7 s. Exact solution (-) and numerical computations performed with Cr =1:0, =0:1( --), =0:2( ••), and =0:5( -•-): (a) overview of the whole domain and initial and nal locations of the solitary wave; and (b) close inspection of numerical results.

Figure 12 .

 12 Figure 12. Stability limits estimated numerically using =0:1 for increasing nonlinear solitary waves and dierent values for the dispersion correction parameter:(-) =1=15, (--) =1=30, and (-•-) =0.

u

  * (x * ;z * ;t * )= ∞ n =0 (z * - * (x * )) n n! f n (x * ;t * ) (A1) w * (x * ;z * ;t * )= ∞ n =0 (z * - * (x * )) n n! g n (x * ;t * ) (A2)

ACKNOWLEDGEMENTS

The rst author would like to thank the French Foreign Oce, CROUS Grenoble, the Chilean Research Council (CONICYT) and Ponticia Universidad Cat olica de Chile for their nancial support during his PhD research. Cristi an Escauriaza and anonymous reviewers are also acknowledged for their useful suggestions and comments.

Contract=grant sponsor: French Foreign Oce Contract=grant sponsor: Chilean Research Council (CONICYT)

Similarly, this functional form for horizontal and vertical velocities allows us to express the vertical ow acceleration function given by Equation ( 9

using dimensional functions, P and Q, which were introduced in relations ( 20) and ( 21).