Ralph Lteif

St Éphane Gerbi

A NEW CLASS OF HIGHER-ORDERED/EXTENDED BOUSSINESQ SYSTEM FOR EFFICIENT NUMERICAL SIMULATIONS BY SPLITTING OPERATORS

Keywords: Water waves, Boussinesq system, Higher-order asymptotic model, Splitting scheme, Hybrid finite volume/finite difference scheme

In this work, we numerically study the higher-ordered/extended Boussinesq system describing the propagation of water-waves over flat topography. A reformulation of the same order of precision that avoids the calculation of high order derivatives on the surface deformation is proposed. We show that this formulation enjoys an extended range of applicability while remaining stable. Moreover, a significant improvement in terms of linear dispersive properties in high frequency regime is made due to the suitable adjustment of a dispersion correction parameter. We develop a second order splitting scheme where the hyperbolic part of the system is treated with a high-order finite volume scheme and the dispersive part is treated with a finite difference approach. Numerical simulations are then performed under two main goals: validating the model and the numerical methods and assessing the potential need of such higher-order model. The applicability of the proposed model and numerical method in practical problems is illustrated by a comparison with experimental data.

1. Introduction 1.1. Motivation. In this paper, we numerically study higher-ordered/extended Boussinesq type models. These equations describe the one-dimensional flow of the free surface of a homogeneous, immiscible fluid moving above a flat topography. They are obtained from the free surface waterwaves (Euler) equations [START_REF] Lannes | The water waves problem[END_REF] (see Eq. (1.3) therein) for an irrotational and incompressible fluid.

The derivation classically relies on a re-scaling of the system in order to reveal small dimensionless parameters which allow to perform asymptotic expansions of non-local operators (Dirichlet-Neumann), thus ignoring the terms whose influence is minimal. We start by introducing respectively the commonly known nonlinear and shallowness parameters:

0 ≤ ε = a h 0 = wave amplitude reference depth ≤ 1 , 0 ≤ √ µ = h 0 λ = reference depth wave-length of the wave < 1 .
The order of magnitude of these parameters makes it possible to identify the considered physical regime.

In a specific long wave regime, ε is considered of the same order as µ (ε ∼ µ). In this regime, Boussinesq derived in [START_REF] Boussinesq | Théorie de l'intumescence liquide appelée onde solitaire ou de translation se propageant dans un canal rectangulaire[END_REF][START_REF] Boussinesq | Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond[END_REF] a weakly nonlinear model bearing his name. In what follows we refer to it as the "original" or "standard" Boussinesq system. Using the horizontal depth-mean velocity v and the free surface parametrization ζ, the standard Boussinesq (sB) equations reads:

(1)

   ∂ t ζ + ∂ x (1 + εζ)v = 0 , 1 -ε 1 3 ∂ 2 x ∂ t v + ∂ x ζ + εv∂ x v = O(ε 2 ) .
This model can be derived from the Green-Naghdi (GN) equations (see [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF]) by neglecting all terms of order O(ε 2 , µε, µ 2 ). Equivalent Boussinesq systems enjoying a better mathematical structure or physical properties have been studied and derived extensively in the literature, see for instance [START_REF] Linares | Well-posedness of strongly dispersive two-dimensional surface wave Boussinesq systems[END_REF][START_REF] Ming | Long-time existence of solutions to Boussinesq systems[END_REF][START_REF] Saut | The Cauchy problem on large time for surface-waves-type Boussinesq systems ii[END_REF][START_REF] Saut | The Cauchy problem on large time for surface waves Boussinesq systems[END_REF]. The sB equations are restricted by containing only weak dispersion and non-linearity (only O(µ, ε) terms are retained). This normally limits precise applications to a small zone moderately exterior to the surf zone. Significant improvement have been made in recent years to expand the application range and cover the range fully from deep water into the surf zone. Madsen et al. [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF] reached this goal by rearranging the dispersive terms (or O(µ) terms) in order to improve linear dispersion properties. On the other hand, Nwogu [START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF] achieved the same result by redefining the dependent velocity variable. These models have been extensively examined for their utility in the prediction of near-shore problems (wave breaking, run-up, wave-induced circulation) as detailed in [START_REF] Madsen | Surf zone dynamics simulated by a Boussinesq type model. Part I. Model description and cross-shore motion of regular waves[END_REF][START_REF] Madsen | Surf zone dynamics simulated by a Boussinesq type model. Part II: surf beat and swash oscillations for wave groups and irregular waves[END_REF][START_REF] Sørensen | Surf zone dynamics simulated by a Boussinesq type model. III. Wave-induced horizontal nearshore circulations[END_REF]. However, these Boussinesq-type models which assume the velocity profiles to be second order polynomials in the vertical coordinate induce inaccuracies near wave breaking [START_REF] Gobbi | A fully nonlinear Boussinesq model for surface waves. II. Extension to O(kh) 4[END_REF]. Several mechanisms exists in order to handle wave breaking that occurs as waves approach the shore. For instance, a hybrid method consisting in suppressing the dispersive terms in breaking regions was initially suggested by Tonelli and Petti [START_REF] Tonelli | Hybrid finite volume -finite difference scheme for 2DH improved Boussinesq equations[END_REF]. Another strategy consists of on an eddy viscosity approach based on the solution of a turbulent kinetic energy following early work by Nwogu [START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF]. The interested reader is referred to [START_REF] Kazolea | On wave breaking for Boussinesq-type models[END_REF] for a comparison between the hybrid and eddy viscosity strategies. Other efficient mechanisms that allow solving the wave breaking problem are worth mentioning, see for instance [START_REF] Kirby | Boussinesq models and their application to coastal processes across a wide range of scales[END_REF][START_REF] Escalante | An efficient hyperbolic relaxation system for dispersive nonhydrostatic water waves and its solution with high order discontinuous Galerkin schemes[END_REF] and references therein.

Many attempts have been made to extend Boussinesq-type models in order to offer better dispersive properties. To this end, in order to incorporate high-order dispersive and nonlinear effects in Boussinesq-type models, one should include some high-order terms in the asymptotic expansion of the velocity potential. The first attempt to derive higher-order Boussinesq-like equations (retaining O(µ 2 ) terms) was performed by Dingemans in [START_REF] Dingemans | Water waves over an uneven bottom; a discussion of long-wave equations[END_REF]. Two versions of equations were given, one based on the depth-averaged velocity and one based on the velocity at the still water level. Dingemans did not provide analyses or computations based on these equations. One can see also the review papers written by Kirby [START_REF] Kirby | Nonlinear, dispersive long waves in water of variable depth[END_REF] and Madsen & Schäffer [START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF][START_REF] Madsen | A review of Boussinesq-type equations for surface gravity waves[END_REF] where Boussinesqtype equations of higher order in dispersion as well as in non-linearity are derived, intensified and analyzed with emphasis on linear dispersion, shoaling and nonlinear properties for large wave numbers. A fully non-linear Boussinesq-type model retaining O(µ 2 ) terms was derived for an horizontal bottom in [START_REF] Gobbi | A fully nonlinear Boussinesq model for surface waves. II. Extension to O(kh) 4[END_REF] and examined for its ability to represent weakly nonlinear wave evolution in intermediate depth and its numerical properties of solitary wave solutions in shallow water. Algorithms for the numerical solution of the latter model are described in [START_REF] Gobbi | Wave evolution over submerged sills: tests of a high-order Boussinesq model[END_REF], where the model is applied to the study of wave shoaling and harmonic generation in the problem of waves propagating over an isolated step. Including higher order terms in the model is not the only method to improve the dispersion relation. In fact, there are other strategies to improve the disperive properties, for instance, improving the degree of freedom of the velocity unknown [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF][START_REF] Duchêne | A new class of two-layer Green-Naghdi systems with improved frequency dispersion[END_REF][START_REF] Kazakova | A new model of shoaling and breaking waves: one-dimensional solitary wave on a mild sloping beach[END_REF]. The interested reader is also referred to [START_REF] Fernández-Nieto | A hierarchy of dispersive layeraveraged approximations of Euler equations for free surface flows[END_REF] where a hierarchy of new models is derived with a layer-wise approach incorporating non-hydrostatic effects to approximate the Euler equations. The linear dispersion relation of these models is analyzed therein and proved to converge to the dispersion relation for the Euler equations when the number of layers goes to infinity. One can see also [START_REF] Escalante | An efficient two-layer non-hydrostatic approach for dispersive water waves[END_REF] where a two-layer non-hydrostatic model with improved dispersive properties is derived and [START_REF] Sánchez | Numerical simulations of a dispersive model approximating free-surface Euler equations[END_REF] where a numerical scheme is designed for theses models. Other higher-order asymptotic shallow-water models were derived in the literature. For instance, the extended Green-Naghdi (eGN) equations (accurate up to the order O(µ 3 ) while the full non-linearity is preserved) were firstly derived in their Hamiltonian formulation by Matsuno in [START_REF] Matsuno | Hamiltonian formulation of the extended Green-Naghdi equations[END_REF][START_REF] Matsuno | Hamiltonian structure for two-dimensional extended Green-Naghdi equations[END_REF]. More recently, Khorbatly et al. derived the eGN equations in [START_REF] Khorbatly | Derivation and well-posedness of the extended Green-Naghdi equations for flat bottoms with surface tension[END_REF] by performing an asymptotic analysis of the Dirichlet-Neumann operator that originates from the formulation of the water wave problem [START_REF] Lannes | The water waves problem[END_REF]. In the aforementioned papers of Khorbatly et al., the mathematical analysis is addressed which is mainly devoted to the well-posedness of the equations.

Note, however, that less determined efforts were made to numerically study higher-ordered asymptotic equations for the water-waves problem. Due to their extensive length, the high-order equations incorporating very high order derivatives (see for example fifth-order derivatives in [START_REF] Gobbi | A fully nonlinear Boussinesq model for surface waves. II. Extension to O(kh) 4[END_REF]) may not seem viable as a basis for a numerical model. At this point, it is worth mentioning several recent advances presenting novel hyperbolic reformulation of Serre-Green-Naghdi [START_REF] Bassi | A hyperbolic reformulation of the Serre-Green-Naghdi model for general bottom topographies[END_REF] and Boussinesq-type [START_REF] Escalante | A general non-hydrostatic hyperbolic formulation for Boussinesq dispersive shallow flows and its numerical approximation[END_REF] models, see also [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF] for the first derivation of hyperbolic reformulation of a dispersive system from variational principles in the flat bottom case. Those new first-order reformulations are based on a relaxed augmented system in which the divergence constraints of the velocity flow variables are coupled with the other conservation laws via an evolution equation for the depth-averaged non-hydrostatic pressures. They avoid the use of high order derivatives which are not easy to treat numerically due to the large stencil usually needed. Moreover allow to overcome the numerical difficulties and the severe time step restrictions arising from higher order terms [START_REF] Busto | On high order ADER discontinuous Galerkin schemes for first order hyperbolic reformulations of nonlinear dispersive systems[END_REF].

Bearing these facts in mind, we develop in this paper a numerical model solving a class of higher-ordered/extended Boussinesq type models where we discuss their usefulness for practical applications. We propose a new reformulation of the model with improved linear dispersive properties and an extended range of applicability. This new reformulation is suitable to the implementation of a hybrid scheme splitting the hyperbolic and dispersive parts of the equations. This strategy has been initially introduced for Boussinesq-like and Green-Naghdi equations in order to handle correctly wave breaking, see [START_REF] Tonelli | Hybrid finite volume -finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF]. The hyperbolic part of the system is treated with a high-order finite volume scheme whereas the dispersive part is treated with a finite difference method at the same order. This splitting strategy allowed us to overcome the severe time step restriction induced due to the presence of high order derivatives by calculating the time step in the first finite-volume sub-step. The numerical investigations show that the time step restriction from the CFL condition (according to which the time step must be chosen proportional to the mesh spacing) of the finite-volume step is enough to ensure stability for the whole numerical method. As a preliminary step, we treat in this paper the case of flat bottom topography, investigating the variable topography case can be based on the results of this paper and will be faced in a forthcoming work. We would like to emphasize that this paper provides sufficient grounds for treating the more complex variable topography case. The goal of this paper is to give a numerical assessment of the potential need of such a new formulation of higher-order model that avoids the calculation of higher order derivatives and evaluate the effect of adding factorized high order terms to standard models. 1.2. Higher-ordered/extended Boussinesq equations. Let us start first by introducing the higherordered extended Boussinesq (eB) equations that naturally show up through the asymptotic approximation of the Dirichlet-Neumann operator. These equations are a straightforward extension of the sB equations [START_REF] Bassi | A hyperbolic reformulation of the Serre-Green-Naghdi model for general bottom topographies[END_REF]. Neglecting the terms of order O(ε 3 ) while keeping the O(ε 2 ) terms in the equations one gets the eB equations. Alternatively, one can easily recover the eB equations from the eGN equations derived in [START_REF] Khorbatly | Derivation and well-posedness of the extended Green-Naghdi equations for flat bottoms with surface tension[END_REF], by considering weak non-linearity (ε ∼ µ) and dropping all terms of order O(µε 2 , µ 2 ε, ε 3 ) therein. Thus, one can write the weakly nonlinear Boussinesq system including higher order dispersive effects as follows:

(2)

∂ t ζ + ∂ x (hv) = 0 , (1 + εT [εζ] + ε 2 T)∂ t v + ∂ x ζ + εv∂ x v + ε 2 Qv = O(ε 3 ) ,
where h = 1 + εζ is the non-dimensionalised height of the fluid and denote by

T [εζ]w = - 1 3h ∂ x h 3 ∂ x w , Tw = - 1 45 ∂ 4 x w, Qv = - 1 3 ∂ x vv xx -v 2 x .
Actually, neglecting terms of order O(ε 2 ) in ( 2), the second order differential operator becomes

T [εζ]w = -1 3 ∂ 2
x w and one can easily recover [START_REF] Bassi | A hyperbolic reformulation of the Serre-Green-Naghdi model for general bottom topographies[END_REF]. The eB model ( 2) can be found in [START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF] (see Eq. (3.11) therein after neglecting all terms of order O(εµ 4 , ε 2 µ 2 )). Unfortunately, the eB model ( 2) seems to suffer from instabilities that turn out to be fatal for any practical use. This is due to the positive sign in front of the elliptic fourth-order linear operator T which also prevent the invertibility of the factorized operator, see [START_REF] Khorbatly | Mathematical modeling and numerical analysis for the higher order Boussinesq system[END_REF]Section 3.1]. In fact, linearizing the eB model ( 2) around some rest state solution, one gets the following unstable dispersion relation:

(3)

w 2 = k 2 1 + 1 3 εk 2 - 1 45 ε 2 k 4 ,
where k is the spatial wave number and w represents the time frequency. Note that one can recover the dispersion relation associated to (1) by neglecting ε 2 terms in (3). As expected, Madsen & Schäffer [START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF] noticed an improved accuracy for small wave number values k when comparing with the lower-order equations (see Figure 1 therein). However, a quick functional study shows that the denominator of the right hand side of (3) becomes negative whenever k 1, preventing any hope concerning the well-posedness of the initial value problem. In fact, a singularity occurs in (3) for k 2 = 1 2 (15 + 9 √ 5), i.e. k ≈ 4.2. Despite being a large wave number value, this singularity proves to be inoperable for any reasonable application of the eB equations [START_REF] Benjamin | Model equations for long waves in nonlinear dispersive systems[END_REF], see [START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF].

To overcome this problem and in order to gain some confidence into the model validity, an enhanced set of equations of same order of precision, without instabilities, has been derived in [START_REF] Khorbatly | Mathematical modeling and numerical analysis for the higher order Boussinesq system[END_REF] by replacing the left most term of the second equation of ( 2

), (1 + εT [εζ] + ε 2 T)∂ t v, by (1 + εT [εζ] -ε 2 T)(∂ t v) + 2ε 2 T(∂ t v)
and making use of a BBM trick (Benjamin-Bona-Mahony) [START_REF] Benjamin | Model equations for long waves in nonlinear dispersive systems[END_REF] represented in the following approximate equation [START_REF] Khorbatly | Derivation and well-posedness of the extended Green-Naghdi equations for flat bottoms with surface tension[END_REF] for more details. In light of these remarks and after setting ±ε 2 T [εζ](vv x ) in the second equation of (2), one obtains the following model:

∂ t v = -∂ x ζ + O(ε), see
(4)      ∂ t ζ + ∂ x (hv) = 0, J (∂ t v + εv∂ x v) + ∂ x ζ + 2 45 ε 2 ∂ 5 x ζ + 2 3 ε 2 ∂ x ((∂ x v) 2 ) = O(ε 3 ).
where h = 1 + εζ and (5)

J = 1 + εT [εζ] -ε 2 T.
w 2 = k 2 1 + 2 45 ε 2 k 4 1 + 1 3 εk 2 + 1 45 ε 2 k 4 .
This has an important side effect of removing every singularity found in (3), thus making a set of useless high-order equations applicable. An equivalent formulation of model ( 4) 1 was fully justified recently in [START_REF] Khorbatly | Mathematical modeling and numerical analysis for the higher order Boussinesq system[END_REF]. In fact, a unique solution of the model ( 4) exist over the time scale of order 1 √ ε and stay close the solution of the full Euler system. Although every singularity found in ( 3) is now removed in (5), the phase velocity associated to (5) has the same classical velocity limit of long waves as k → 0, but a finite limit equal to √ 2 as k → ∞ instead of expected zero limit. In order to improve the model linear dispersive properties in intermediate regime of wave numbers and extend the range of applicability without changing the basic limit properties of the dispersion relation, we derive in this paper a new reformulation [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] of the same order of precision of (4) that allows the adjustment of a dispersion correction parameter α (see Section 2.1) and prevents the calculation of high order derivatives on ζ (see Section 2.2). The dispersion relation obtained around some rest state solution associated to the new formulation ( 14) is the following (see detailed calculation in Appendix A):

(6) w 2 = k 2 1 + ε(α -1)k 2 3 + ε 2 (α -1)k 4 45 + (7 -5α)ε 2 k 4 45(1 + εα 3 k 2 ) 1 + εα 3 k 2 + ε 2 α 45 k 4
.

The corresponding phase velocity has in the limit when k → 0 a classical velocity of long waves, and in the limit of short waves when k → ∞ the phase velocity has a finite limit equal to α -1 α . Consequently, setting α = 1, the phase velocity vanishes in the limit of short waves when k → ∞ as in the case of the exact linear dispersion relation [START_REF] Duchêne | A new class of two-layer Green-Naghdi systems with improved frequency dispersion[END_REF] for the full Euler equation. This paper is devoted to the numerical study of the eB system [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF]. We will show that thanks to the proper choice of the dispersion correction parameter α (see Section 2.3) the eB model with factorized high order derivatives [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] have better dispersive properties in intermediate regime of wave numbers than the eB model without factorization [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF]. Moreover, we will show that this newly derived formulation is stable with respect to high frequency perturbations (see Section 2.4).

The paper is organized as follows, we firstly propose a reformulation of the same order of precision of the extended Boussinesq model (4) up to the third order. This reformulation makes the model more appropriate for the numerical implementation and significantly improved in terms of linear dispersive properties due to the suitable adjustment of a dispersion correction parameter. The reformulation is performed then via the factorization of high order derivatives on the surface deformation ζ. We will show that the improvement is significant in the dispersive properties of the model with factorization of high order derivatives on the surface deformation together with an appropriate choice of an optimal value of the dispersion correction parameter α. We then study the stability of two models, with and without factorization, and we will show that factorizing only the fifth order derivative presented in the second model equation induces a destabilizing effect : we need to factor every high order derivative on ζ.

Secondly, we propose a suitable Strang splitting of operators to solve the improved model : a hyperbolic part representing the Nonlinear Shallow Water system and a dispersive part representing the high order derivatives. The hyperbolic part of the system is treated with a high-order finite volume scheme whereas the dispersive part is treated with a finite difference method at the same order. This splitting strategy allowed us to calculate the time step in the first finite-volume sub-step allowing to overcome the time step restrictions induced by high order derivatives existing 1 The equivalent formulation is obtained by multiplying both sides of the second equation of system (4) by the water height function, h.

in the dispersive part of the model and the numerical investigations show that the time step restriction from the CFL condition (according to which the time step must be chosen proportional to the mesh spacing) of the finite-volume step is enough to ensure stability for the whole numerical method. Moreover, a reconstruction of nodal unknowns and centered unknowns is presented.

Finally, numerical validations are presented under two main goals: showing the interest of the proposed formulation of the extended Boussinesq model as well as the good behavior of the numerical scheme and assessing the potential need of such higher-order models.

Reformulation of the extended Boussinesq system

The system (4) is much easier to solve numerically than the standard formulation [START_REF] Benjamin | Model equations for long waves in nonlinear dispersive systems[END_REF]. In fact, the operator J has an appropriate structure allowing its inversion. Using straightforward asymptotic expansions, the left-most term of the second equation of ( 4) can be written under the form:

J(∂ t v + εv∂ x v) = (1 + εT [0] -ε 2 T)(∂ t v + εv∂ x v) - 2 3 ε 2 ζ∂ 2 x (∂ t v) -ε 2 ∂ x ζ∂ x (∂ t v) + O(ε 3 ). Now using the fact that ∂ t v = -∂ x ζ + O(ε)
, one can deduce that the above equation can be recast under the following form:

J(∂ t v + εv∂ x v) = (1 + εT [0] -ε 2 T)(∂ t v + εv∂ x v) + 2 3 ε 2 ζ∂ 3 x ζ + ε 2 ∂ x ζ∂ 2 x ζ + O(ε 3 ).
where

T [0]w = - 1 3 ∂ 2
x w and Tw = -

1 45 ∂ 4 x w. Hence, system (4) becomes:
(7)

           ∂ t ζ + ∂ x hv = 0, 1 + εT [0] -ε 2 T ∂ t v + εv∂ x v + ∂ x ζ + 2 45 ε 2 ∂ 5 x ζ + 2 3 ε 2 ∂ x ((∂ x v) 2 ) + 2 3 ε 2 ζ∂ 3 x ζ + ε 2 ∂ x ζ∂ 2 x ζ = O(ε 3 ),
where h = 1 + εζ. The left-most factorized operator of the second equation of (7) also enjoys a structure allowing its inversion. For the proof of the invertibility, one has to apply a Lax-Milgram theorem where the coercivity condition of the bilinear form is satisfied (see [START_REF] Khorbatly | Derivation and well-posedness of the extended Green-Naghdi equations for flat bottoms with surface tension[END_REF]Lemma 1]). Moreover, this operator can be inverted once for all numerical time steps because of its timeindependent form. The simple one-dimensional structure of the model with a time-independent operator reduce slightly the computational time. In fact, the strategy of removing time-dependency from the left-most factorized operator was originally initiated for numerical simulations of the fully nonlinear and weakly dispersive GN models in the two-dimensional case [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF], in order to reduce significantly the computational time.

2.1.

A one-parameter family of extended Boussinesq equations. The eB equations are significantly improved in terms of linear dispersive properties due to the higher-order terms existing in these equations, see [START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF]. Additional improvement providing a finer characterization in high frequency regimes can be brought by adjusting a dispersion correction parameter α. Following the lines in [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. II. Boundary conditions and validation[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF] and without affecting the accuracy of the model, we improve the frequency dispersion of problem [START_REF] Boussinesq | Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond[END_REF]. This is possible, if one adds to the second equation of (7) some terms of the same order as the equation precision and adjusts the parameter α in an appropriate way. See section 2.3 for the discussion on the choice of the parameter α. From the second equation of [START_REF] Boussinesq | Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond[END_REF], one deduce the following approximation:

(8) ∂ t v + εv∂ x v + ∂ x ζ + εT [0](∂ t v) = O(ε 2 ),
where

T [0]w = - 1 3 ∂ 2 x w. Using again the fact that ∂ t v = -∂ x ζ + O(ε)
, thus approximation (8) can be written as:

∂ t v + εv∂ x v + ∂ x ζ + ε 3 ∂ 3 x ζ = O(ε 2 ),
and hence, for any α ∈ R * + :

(9)

∂ t v = α∂ t v + (α -1)[∂ x ζ + εv∂ x v + ε 3 ∂ 3 x ζ] + (1 -α)O(ε 2 ).
The second equation of ( 7) can be recast after substituting ∂ t v by its approximation given in (9):

1 + εT [0] -ε 2 T α∂ t v + αεv∂ x v + (α -1)∂ x ζ + ε 3 (α -1)∂ 3 x ζ + (1 -α)O(ε 2 ) + α -1 α + 1 α ∂ x ζ + 2 45 ε 2 ∂ 5 x ζ + 2 3 ε 2 ∂ x ((∂ x v) 2 ) + 2 3 ε 2 ζ∂ 3 x ζ + ε 2 ∂ x ζ∂ 2 x ζ = O(ε 3 ).
After neglecting the terms (εT [0] -

ε 2 T)(1 -α)O( 2 ) and -ε 2 T ε 3 (α -1)∂ 3 x ζ of order O(ε 3
), one has:

α∂ t v + αεv∂ x v + (α -1)(∂ x ζ + ε 3 ∂ 3 x ζ) + (1 -α)O(ε 2 ) + εT [0] -ε 2 T α∂ t v + αεv∂ x v + (α -1)∂ x ζ + ε 2 3 (α -1)T [0]∂ 3 x ζ + α -1 α + 1 α ∂ x ζ + 2 45 ε 2 ∂ 5 x ζ + 2 3 ε 2 ∂ x ((∂ x v) 2 ) + 2 3 ε 2 ζ∂ 3 x ζ + ε 2 ∂ x ζ∂ 2 x ζ = O(ε 3 ).
Using the fact that

α∂ t v + αεv∂ x v = ∂ t v + εv∂ x v + (α -1)(∂ t v + εv∂ x v) one has: ∂ t v + εv∂ x v + (α -1)(∂ t v + εv∂ x v + ∂ x ζ + ε 3 ∂ 3 x ζ) + (1 -α)O(ε 2 ) + εαT [0] -ε 2 αT ∂ t v + εv∂ x v + (α -1) α ∂ x ζ + ε 2 3 (α -1)T [0]∂ 3 x ζ + α -1 α + 1 α ∂ x ζ + 2 45 ε 2 ∂ 5 x ζ + 2 3 ε 2 ∂ x ((∂ x v) 2 ) + 2 3 ε 2 ζ∂ 3 x ζ + ε 2 ∂ x ζ∂ 2 x ζ = O(ε 3 ).
Following straightforward computations and using the fact that

ε 2 3 (α-1)T [0]∂ 3 x ζ = - ε 2 9 (α-1)∂ 5 x ζ, one has: 1+εαT [0]-ε 2 αT ∂ t v +εv∂ x v + α -1 α ∂ x ζ +(α-1)(∂ t v + v∂ x v +∂ x ζ + ε 3 ∂ 3 x ζ)+(1-α)O(ε 2 ) + 1 α ∂ x ζ + (7 -5α) 45 ε 2 ∂ 5 x ζ + 2 3 ε 2 ∂ x ((∂ x v) 2 ) + 2 3 ε 2 ζ∂ 3 x ζ + ε 2 ∂ x ζ∂ 2 x ζ = O(ε 3 ).
Finally using [START_REF] Busto | On high order ADER discontinuous Galerkin schemes for first order hyperbolic reformulations of nonlinear dispersive systems[END_REF], one has (α-1)

(∂ t v + v∂ x v +∂ x ζ + ε 3 ∂ 3 x ζ)+(1-α)O(ε 2
) = 0 and thus system (7) with improved frequency dispersion can be written as:

(10)              ∂ t ζ + ∂ x hv = 0, 1 + εαT [0] -ε 2 αT ∂ t v + εv∂ x v + α -1 α ∂ x ζ + 1 α ∂ x ζ + (7 -5α) 45 ε 2 ∂ 5 x ζ + 2 3 ε 2 ∂ x ((∂ x v) 2 ) + 2 3 ε 2 ζ∂ 3 x ζ + ε 2 ∂ x ζ∂ 2 x ζ = O(ε 3 ).
Similarly, a significant improvement of the dispersive properties has been attained in the derivation of a three-parameter family of GN equations, see [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF]. In here, we will limit ourselves to the oneparameter family of eB equations [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] for the sake of simplicity.

2.2.

Reformulation of the extended Boussinesq equations [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF]. In what follows, we derive an equivalent model to [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] (in the sense of precision), i.e. O(ε 3 ), that prevents the calculation of high order derivatives on ζ. To this effect, we call such a model eB with factorized high order derivatives. Certainly, the model enclose high order derivatives on ζ, but we make it possible not to compute them by factoring them out by (1 + εαT [0]). The price to pay is an increase in computational cost, since one needs to solve an extra linear system but the gain is significant in extending the range of applicability.

Using the fact that (1

+ εαT [0])(∂ x ζ) = ∂ x ζ + O(ε), one has ∂ x ζ = (1 + εαT [0]) -1 (∂ x ζ) + O(ε), and thus the terms ∂ 2 x ζ, ∂ 3 x ζ and ∂ 5 x ζ become respectively: (11) ∂ 2 x ζ = ∂ x (1 + εαT [0]) -1 (∂ x ζ) + O(ε), (12)
∂ 3 x ζ = ∂ 2 x (1 + εαT [0]) -1 (∂ x ζ) + O(ε), (13)
∂ 5 x ζ = ∂ 4 x (1 + εαT [0]) -1 (∂ x ζ) + O(ε). Replacing ∂ 2 x ζ, ∂ 3 x ζ and ∂ 5
x ζ by their expression obtained in ( 11), ( 12) and ( 13) respectively in the second equation of [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF], one can write the eB equations with improved dispersion and factorized high order derivatives on ζ as:

(14)                    ∂ t ζ + ∂ x hv = 0, 1 + εαT [0] -ε 2 αT ∂ t v + εv∂ x v + α -1 α ∂ x ζ + 1 α ∂ x ζ + (7 -5α) 45 ε 2 ∂ 4 x (1 + εαT [0]) -1 (∂ x ζ) + 2 3 ε 2 ∂ x ((∂ x v) 2 ) + 2 3 ε 2 ζ∂ 2 x (1 + εαT [0]) -1 (∂ x ζ) + ε 2 ∂ x ζ∂ x (1 + εαT [0]) -1 (∂ x ζ) = O(ε 3
), where h = 1 + εζ. Note, that this formulation avoids the calculation of high order derivatives, in particular on the surface deformation ζ, but there are still fourth and second order derivatives (this time not on ζ) in the second equation of ( 14) and a large stencil is still needed though. Significant interest is behind the derivation of the eB formulation [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF]. In fact, we believe that factorizing high order derivatives, namely on the surface deformation ζ, will extend the range of applicability to high frequency regimes (see discussion in section 2.3), while remaining stable (see section 2.4.2). In the following section, we will highlight the advantages of working with the factorized eB model ( 14) rather than [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF].

Choice of the parameter α. The main comparison between any asymptotic model and the full

Euler equations is performed at the stage of linear periodic plane wave solutions. At this point, a part of the model's properties [START_REF] Stoker | Water Waves: The mathematical theory with applications[END_REF] are summed up in the dispersion relation, relating the spatial wave number k and the time frequency w. It comes from the earlier linearisation of the system around some rest state. Improving the dispersive characteristics of our model require a suitable choice of the parameter α so that the dispersion characteristics of the full Euler system corresponds with those of the improved eB systems at the dispersion relation level. Following [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], we adjust this parameter so that both phase and group velocities are minimized over a range of values of k ∈ [0, K]. This can be done by minimizing a weighted averaged error (see for instance [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. II. Boundary conditions and validation[END_REF]) introduced for this reason.

In what follows, we will show that the eB model with factorized high order derivatives have better dispersive properties when compared to models including the eB model without factorization and other lower-order models. The dispersion relation corresponding to [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] can be derived by investigating the linear behavior of small perturbation to a constant state solution (ζ, v) and then looking for the corresponding plane wave solutions of the form (ζ 0 , v 0 )e i(kx-wt) :

(15) (w -εkv) 2 hk 2 = 1 + ε(α -1)k 2 3 + (6 -4α)ε 2 k 4 45 - 2ε 2 k 2 ζ 3 1 + εα 3 k 2 + ε 2 α 45 k 4 .
The choice of α is classically made to obtain a good matching with the dispersion relation of the full Euler equations around the rest state (ζ, v) = (0, 0). The exact dispersion relation for the full Euler system is recalled below:

(16) w 2 S = |k| √ ε tanh( √ ε|k|).
For small wave numbers, the Taylor expansions of (15) (with (ζ, v) = (0, 0)) and ( 16) are equivalent and the choice of α does not play any role in the leading terms (see numerical test 4.5). Indeed, one has:

w 2 α,eB = w 2 S ≈ k 2 - εk 4 3 + 2 15 ε 2 k 6 + O(ε 3 k 8 ) .
Classically, finding an optimal value of α for a range of values of k requires the minimization of the squared relative weighted error defined below:

(17) Err = K 0 1 k C p eB -C p S C p S + C g eB -C g S C g S 2 dk, over some range k ∈ [0, K],
where C p S (k) and C g S (k) are respectively the reference phase and group velocities associated with the Stokes linear theory. The division by k will emphasize the importance of keeping errors to a minimum in shallow water. The linear phase and group velocities associated to [START_REF] Dingemans | Water waves over an uneven bottom; a discussion of long-wave equations[END_REF] are defined as:

C p eB (k) = w α,eB (k) |k| and C g eB (k) = dw α,eB (k) dk .
The weighted averaged error ( 17) has an absolute minimum of (60%) for w α,eB defined in (15) (with (ζ, v) = (0, 0)) in the dispersive range 0 ≤ k ≤ 10. The optimal value for α in this case does not play any role and α = 1 is set. This very big error shows that the eB model [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] without factorizing high order derivatives on ζ has a limited range of applicability and thus poor dispersion properties in large wave numbers regime. C g S (group velocity, blue solid line) with α = 1 when compared with the reference from Stokes theory (black solid line) when k > 2. On the contrary, a very good correspondence is observed in small wave numbers regime (i.e when k ≤ 1). This shows that the eB model [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] without factorizing high order derivatives on ζ has good dispersion properties in small wave numbers regime (k ≤ 1). However, this is not the case in larger wave numbers regime (k > 2). Now, we discuss the dispersive properties of model [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF]. The dispersion relation associated to [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] is the following (see Appendix A for detailed calculation): [START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF] ( wα,eB -εkv)

2 hk 2 = 1 + ε(α -1)k 2 3 + ε 2 (α -1)k 4 45 + (7 -5α)ε 2 k 4 45(1 + εα 3 k 2 ) - 2ε 2 k 2 ζ 3(1 + εα 3 k 2 ) 1 + εα 3 k 2 + ε 2 α 45 k 4 .
At this stage, one has to minimize the error function [START_REF] Dutykh | Efficient computation of steady solitary gravity waves[END_REF] for wα,eB defined in ( 18) with (ζ, v) = (0, 0). In Figure 2 we plot the associated error in terms of α. One can clearly see that the weighted average error has an absolute minimum (≈ 1%) in the dispersive range 0 ≤ k ≤ 10. The optimal value for α is 1.0610. Meanwhile, the absolute minimum of the weighted averaged error associated with the lower-order Green-Naghdi model in the Camassa-Holm regime (GN-CH) (precise up to O(µ 2 , µε 2 )) [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] is much larger (≈ 10%) with an optimal value α = 1.0800. In Figure 3, errors on linear phase (red) and group (blue) velocities are plotted. The ratio C p eB,αopt C p S with an optimal choice α = 1.0610 (red solid line) is very close to the reference from Stokes theory (black solid line) which shows a very good correspondence between the dispersion relation obtained using the eB model [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] with factorized high order derivatives on ζ with α = 1.0610 and the theoretical one over 0 ≤ k ≤ 10. Larger difference exists between the group velocity (blue solid line) and the reference from Stokes theory (black solid line). This difference starts to proliferate when K > 10, showing an overestimation of this property. When the modeled dispersion relation is obtained using the lower-order GN-CH model [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] with α = 1.0800, a clear discrepancy exists between both linear group and phase velocity errors (blue and red dot lines) when compared with the reference from Stokes theory (black solid line). We also highlight that comparisons of errors on both linear phase and group velocities of the eB model ( 14) against other lower ordered models (for instance, sB or GN type models) are successfully reproduced and a clear superiority of our higher order model is proved when it comes to approximating properly phase and group speeds. We do not include these comparisons in the present study for the sake of shortness. In fact similar results are obtained since the improved formulation (with dispersion correction parameter) of all second order models (sB, GN-CH and GN) is the same when linearized around some rest state and thus one expects a similar linear dispersion relation.

In conclusion, the eB model ( 14) contains factorized higher-order dispersive terms that are neglected in lower-order models which enlarge the application scope remarkably to cover the area from deep water (long wavelength regime) into the area of breaking waves (short wavelength regime i.e. high frequency), see numerical test in section 4.4. This is the reason why, in the numerical experiments, we choose to work with the model [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF], where high order derivatives are factorized, which as seen above, has an extended range of applicability and good dispersive properties in large wave numbers regime.

Remark 1. In the case of variable topography, the study of dispersive properties must be supplemented with some hints about linear shoaling. In fact, the minimization of phase and group velocities errors is quite problematic in the variable topography case as a better phase and group velocities involves a deterioration of shoaling properties, which are a paramount for near-shore oceanography. For this reason, and in order to consider stiff configurations that include high harmonics while keeping the improved dispersive properties of the model, a three-parameter family can be derived by adding two additional parameters θ and γ using a change of variables for the velocity [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF]. The range of validity of the three parameter family of equations is extended to deeper water and considerable improvements in challenging configurations are obtained [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF]. This issue is out of the scope of this article and the variable topography case will be considered in a forthcoming paper.

2.4. Stability of the extended Boussinesq models. The high frequency instabilities of an improved GN-CH model are studied in [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF]. These instabilities are due to the third order derivative existing in the equation. One of the advantages of the eB is its stability in high frequency regime due to the presence of higher order derivatives, namely derivatives of order five in ζ. These terms seems to have a stabilizing effect. In what follows, we discuss qualitatively the stability of both models [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] and [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF]. Since the parameter ε do not play a direct role in the stability results and for notational convenience we set ε = 1 throughout this section.

2.4.1. Stability of the extended Boussinesq model with high order derivatives. Before discussing the stability of the eB models in high frequency regime, we would like to mention that the choice of α in the model [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] in high frequency regime does not play an important role. In fact the latter model has poor dispersive properties in intermediate and large wave numbers regime, see discussion in Section 2.3. Therefore, one has to choose α = 1.

The dispersion relation [START_REF] Dingemans | Water waves over an uneven bottom; a discussion of long-wave equations[END_REF] associated to [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] with α = 1 reads:

(19) (w -kv) 2 hk 2 = 1 + 2k 4 45 - 2k 2 ζ 3 1 + 1 3 k 2 + 1 45 k 4
.

The perturbations of the rest state (ζ, v) = (0, 0) are always stable as per the below dispersion relation:

(20)

w 2 = k 2 1 + 2 45 k 4 1 + 1 3 k 2 + 1 45 k 4 .
However, a quick functional study shows that the numerator of the right-hand side of ( 4). However, as we have mentioned in the introduction, from the dispersion relation [START_REF] Escalante | An efficient two-layer non-hydrostatic approach for dispersive water waves[END_REF] one can remark that the phase velocity associated to the eB model with high order derivatives (10) has the same classical velocity limit of long waves as k → 0, but a finite limit equal to √ 2 instead of expected zero limit.

Remark 2. Replacing ∂ 5 x ζ by ∂ 4 x (1 + T [0]) -1 (∂ x ζ)
in the second equation of (10) with α = 1, one gets the following model:

(21)            ∂ t ζ + ∂ x hv = 0, 1 + T [0] -T ∂ t v + v∂ x v + ∂ x ζ + 2 45 ∂ 4 x (1 + T [0]) -1 (∂ x ζ) + 2 3 ∂ x ((∂ x v) 2 ) + 2 3 ζ∂ 3 x ζ + ∂ x ζ∂ 2 x ζ = 0.
This replacement modifies the dispersion relation [START_REF] Escalante | An efficient hyperbolic relaxation system for dispersive nonhydrostatic water waves and its solution with high order discontinuous Galerkin schemes[END_REF] into:

(22) (w -kv) 2 hk 2 = 1 - 2 3 k 2 ζ + k 4 45 2 1 + 1 3 k 2 1 + 1 3 k 2 + 1 45 k 4
.

A similar functional study to the previous one shows that the r.h.s numerator of ( 22) is negative whenever ζ > 2k 4 + 15k 2 + 45 10k 4 + 30k 2 . In high frequency regime, the upper bound of the overall surface deformation ζ is approximately close to 0.2, hence reducing the range of values of ζ for which the condition is satisfied, namely -1 < ζ < 0.2 (keeping mind that h = 1 + ζ should remain always positive). Therefore, if this condition is not satisfied the complex square root of w will generate an instability in the model. Actually, a stability in high frequency regime is ensured if the condition -1 < ζ < 0.2 is satisfied which we believe is a limitation for the applications that we have in mind. Thus, one can deduce that factorizing only the fifth order derivative present in the second equation of [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] does not stabilize the model, at least for a big range of values of the overall surface deformation ζ (see Figure 4).

At this stage, one may wonder if a factorization of only the third order derivative on ζ may be enough to control the incriminated sign in the dispersion relation [START_REF] Escalante | An efficient hyperbolic relaxation system for dispersive nonhydrostatic water waves and its solution with high order discontinuous Galerkin schemes[END_REF]. Indeed, this is true and w remains real for large wave numbers, under the condition that ζ < 2k 6 + 6k 4 + 45k 2 + 135 90k 2 (for k ≈ 10 the upper bound of the latter inequality becomes very large). However, factoring only the third order derivative leads to a model with the same poor linear dispersive properties as [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF]. This was expected because the third order derivative term is nonlinear. This is why we suggested in section 2.2 to factorize second, third and fifth order derivatives present in the second equation of [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF]. The stability of the eB model with factorized high order derivatives [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] is discussed in the next subsection.

Stability of the extended Boussinesq model with factorized high order derivatives.

In what follows, we explore the stability of the eB model with factorized high order derivatives [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF]. We recall here the dispersion relation [START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF] associated to the eB model ( 14):

(23) (w -kv) 2 hk 2 = 1 + (α -1)k 2 3 - 2k 2 ζ 3(1 + α 3 k 2 ) + (α -1)k 4 45 + (7 -5α)k 4 45(1 + α 3 k 2 ) 1 + α 3 k 2 + α 45 k 4 .
The r.h.s numerator in ( 23) is positive if and only if

ζ < k 2 α (α -1) k 2 + 15α -12 + (18 -15α) + 90α -45 90 + 3 2k 2 .
We recall that α > 1. Indeed, improving the dispersive properties of the model ( 14) in large frequency regime requires an appropriate choice α = 1.0610 (see discussion in section 2.3). With this choice of α and in high frequency regime, the r.h.s of the above inequality becomes very large, namely same order as α(α -1)k 4 90 . Thus, relaxing the stability condition on ζ. This is another reason why, for the rest of the paper, we choose to work with the model [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF], where high order derivatives are factorized, which as seen above is stable for the majority of applications we have in mind. Of course, one could use model ( 10) which seems to be stable, but at the price of losing the improved dispersive properties that model ( 14) enjoys, see section 2.3 and test case of section 4.4.

As we have already mentioned in the introduction, one can remark from the dispersion relation [START_REF] Favre | Ondes de translation dans les canaux découverts[END_REF] around some rest state solution (ζ, v) = (0, 0) that the phase velocity associated to ( 14) has the same classical velocity limit of long waves as k → 0, namely 1, and a finite limit equal to α -1 α as k → ∞ which is equal 0 (as in the case of the exact linear dispersion relation for the full Euler equation ( 16)) when the dispersion correction parameter is set to α = 1.

Figure 4 shows the left to right propagation of a solitary wave initially centered at x 0 = 15 of amplitude a = 0.6. The computational domain length is L = 30 and discretized with 480 cells. The water surface profiles of our numerical solutions provided by the models (10) (blue line), ( 14) (red line) and (21) (green line) are compared at t = 0.5 and t = 0.7 using the fifth order discretization "WENO5-DF4-RK4" (see sections 3.2 and 3.3). One can clearly see the stability of the eB models ( 10) and ( 14) while the model ( 21) seems to be unstable in high frequency regime. In fact, when implementing in model [START_REF] Escalante | A general non-hydrostatic hyperbolic formulation for Boussinesq dispersive shallow flows and its numerical approximation[END_REF] an initial solution that does not satisfy the limiting stability condition discussed in remark 2, more precisely when choosing ζ = 0.6 > 0.2, one observes a high frequency instability. We would like to mention that we tried the same test but starting with an initial solution where the overall surface deformation is ζ = 0.1. All models seems to be stable but we do not include this test here for the sake of simplicity.

Numerical methods

In what follows, we will just introduce the numerical scheme devoted to solve the eB model with factorized high order derivatives [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] in order to ease the reading. A similar numerical scheme is adopted when the eB model with high order derivatives (10) is concerned.

The remarkable structure of the eB models make them suitable for the implementation of a hybrid scheme splitting the hyperbolic and dispersive parts of the equations. This strategy has been initially introduced for Boussinesq-like and Green-Naghdi equations in order to handle correctly wave breaking that occurs as waves approach the shore, see [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF]. A computation of a half-time step of the hyperbolic part is used as a sensor to evaluate the energy loss occurring during wave breaking (accurate detection of wave fronts), see [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF]. Near the breaking points, the dynamics of the waves are described correctly using the hyperbolic part but the dispersive components of the equation become very singular. In order to handle wave breaking, switching from the dispersive part to the hyperbolic part is indispensable. In this paper, we do not investigate breaking waves. In fact, our work is limited to the flat topography case and we leave for future research works the treatment of breaking waves in the variable bottom configuration. However, we stick here to the splitting strategy since it is computationally efficient, stable and cheap. Moreover, this splitting strategy allowed us to overcome the severe time step restriction induced due to the presence of high order derivatives by calculating the time step in the first finite-volume sub-step. The numerical investigations show that the time step restriction from the CFL condition (according to which the time step must be chosen proportional to the mesh spacing) of the finite-volume step is enough to ensure stability for the whole numerical method.The splitting scheme following the lines in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF][START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] is presented in the section below.

3.1. The splitting scheme. We recall first the eB system ( 14) under consideration:

(24)                    ∂ t ζ + ∂ x hv = 0, 1 + εαT [0] -ε 2 αT ∂ t v + εv∂ x v + α -1 α ∂ x ζ + 1 α ∂ x ζ + (7 -5α) 45 ε 2 ∂ 4 x (1 + εαT [0]) -1 (∂ x ζ) + 2 3 ε 2 ∂ x ((∂ x v) 2 ) + 2 3 ε 2 ζ∂ 2 x (1 + εαT [0]) -1 (∂ x ζ) + ε 2 ∂ x ζ∂ x (1 + εαT [0]) -1 (∂ x ζ) = O(ε 3 ), where h = 1 + εζ, T [0]w = - 1 3 ∂ 2
x w and Tw = -

1 45 ∂ 4
x w. The solution operator S(.) related to [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF] is decomposed at each time step ∆t following a hybrid Strang splitting scheme:

S(∆t) = S 1 (∆t/2)S 2 (∆t)S 1 (∆t/2).
• S 1 (t) is the solution operator related to the hyperbolic nonlinear shallow water equations, NSWE:

(25)      ∂ t ζ + ∂ x hv = 0, ∂ t v + εv∂ x v + α -1 α ∂ x ζ + 1 α ∂ x ζ = 0.
The NSWE system (25) can be written in the following conservative form:

(26)      ∂ t ζ + ∂ x hv = 0, ∂ t v + ∂ x ε 2 v 2 + ζ = 0,
where h = 1 + εζ.

• S 2 (t) is the solution operator related to the remaining (dispersive) part of the equations.

(

)                    ∂ t ζ = 0, 1 + εαT [0] -ε 2 αT ∂ t v - 1 α ∂ x ζ + 1 α ∂ x ζ + (7 -5α) 45 ε 2 ∂ 4 x (1 + εαT [0]) -1 (∂ x ζ) + 2 3 ε 2 ∂ x ((∂ x v) 2 ) + 2 3 ε 2 ζ∂ 2 x (1 + εαT [0]) -1 (∂ x ζ) + ε 2 ∂ x ζ∂ x (1 + εαT [0]) -1 (∂ x ζ) = 0 27
The hyperbolic conservative structure of system (26) allows a computation of S 1 following a finitevolume method. Whereas, a classical finite difference method is used to compute S 2 .

Remark 3. Treating the two-dimensional case is not the objective of this paper, however, the extension of the splitting approach to two-dimensional surface waves does not raise theoretical difficulty. Following the steps of [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF], the eB model [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF] in the two-dimensional case can be reformulated in a way that is suitable to the implementation of the splitting strategy with the benefit of removing numerical obstructions. Indeed, the computation of the dispersive part in the above splitting scheme requires the inversion of a fourth order differential operator. This operator is a matricial operator that can be replaced by a new one having a diagonal structure whose inversion is not computationally demanding; numerically this is equivalent to the resolution of two sparse linear systems. Moreover, time dependency can be removed from this operator while keeping its diagonal structure so that it has not to be modified at each time step.

3.2.

Finite volume scheme. The hyperbolic system ( 26) is conveniently rewritten with conservative variables and a flux function:

(28)

∂ t U + ∂ x (F (U )) = 0,
where,

U = ζ v , F (U ) = hv ε 2 v 2 + gζ , (29)
with h = 1 + εζ. The Jacobian matrix is given by: (30) A(U ) = d(F (U )) = εv h g εv .

The homogeneous system (28) is strictly hyperbolic if inf x∈R h > 0 that is to say the domain of the fluid must remain strictly connected.

The Cauchy problem associated to [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF] is the following:

(31)    ∂ t U + ∂ x (F (U )) = 0, t ≥ 0, x ∈ R. U (0, x) = U 0 (x), x ∈ R.
The finite volume method used to the approximation of (31) imposes conservation laws in a onedimensional control volume

[x i-1/2 , x i+1/2 ] × [t n , t n+1 ] of dimensions ∆x = x i+1/2 -x i+1/2 and ∆t = t n+1 -t n .
Figure 5. The space discretization.

The approximate cell average of U on the cell m i = [x i-1/2 , x i+1/2 ] at time t is denoted by U i and given by:

U i = 1 ∆x mi U (t, x) dx .
The approximate cell average of U on the cell

m i = [x i-1/2 , x i+1/2 ] at time t n is denoted by U n i
and given by:

U n i = 1 ∆x mi U (t n , x) dx .
Integrating (25) over the computational cell m i , the semi-discrete form can be represented as:

(32) dU i (t) dt + 1 ∆x F i+1/2 -F i-1/2 = 0
where F i±1/2 are the numerical fluxes defined at each cell interface as:

(33) F i+1/2 = F (U i , U i+1 ) ≈ 1 ∆x mi F (U (t, x i+1/2 ))dx.
VFRoe method. In what follows, we consider the numerical approximation of the hyperbolic system of conservation laws in the form of [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF]. To this end, we adopt the VFRoe method (see [START_REF] Buffard | A sequel to a rough Godunov scheme: application to real gases[END_REF][START_REF] Gallouët | Some recent finite volume schemes to compute Euler equations using real gas EOS[END_REF][START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF]) which is an approximate Godunov scheme. It relies on the exact resolution of the following linearized Riemann problem:

           ∂ t U + A(U n i , U n i+1 )∂ x U = 0, U (0, x) =    U n i if x < x i+1/2 , U n i+1 if x > x i+1/2 , where A(U n i , U n i+1 ) = A U n i + U n i+1 2 . (34)
By solving the linearized Riemann problem we obtain U * i+1/2 = U (x = x i+1/2 , t = t n ), the interface value between two neighboring cells. CFL condition. It is always necessary to impose what is called a CFL condition (for Courant, Friedrichs, Levy) on the time step to prevent the blow up of the numerical values. It comes usually under the form [START_REF] Kazolea | On wave breaking for Boussinesq-type models[END_REF] a i+1/2 ∆t ≤ ∆x, i = 1, . . . , N,

where a i+1/2 = max i∈[1,N ] (j = 1, 2, |λ j ( U i )|) and λ j ( U i ) are the eigenvalues of A U i = U n i + U n i+1
2 .

The restriction [START_REF] Kazolea | On wave breaking for Boussinesq-type models[END_REF] enables in practice to compute the time step at each time level t n , in order to determine the new time level t n+1 = t n + ∆t (within this view, ∆t is not constant, it is computed in an adaptive fashion). Consistency. The numerical flux F (U l , U r ) is called consistent with ( 28) if (36) F (U, U ) = F (U ) for all U.

3.2.1. First order finite-volume scheme. The semi-discrete equation ( 32) is discretized by an explicit Euler (in time) method to obtain:

(37) U n+1 i = U n i - ∆t h i (F n i+1/2 -F n i-1/2 ),
where the numerical flux is defined directly as the value of the exact flux at the interface value, namely:

F n i+1/2 = F (U n i , U n i+1 ) = F (U * i+1/2 ) ( 38
)
F n i-1/2 = F (U n i-1 , U n i ) = F (U * i-1/2 ).
Let us remark that by construction the numerical flux given by [START_REF] Kirby | Nonlinear, dispersive long waves in water of variable depth[END_REF] ensures the consistency property.

In the sequel, we will suppose that the space discretization is uniform. Algorithm. In the following, we state the algorithm for computing the discrete values U n+1 i at t n+1 . Given the initial data and boundary conditions and the number CF L ≤ 1, we start with the known discrete averaged values (U n i ) for i = 0, ..., N +1 at t n . As long as (t < T ) one has to do: 1) Computation of A i for i = 1, ..., N + 1 where

A i = A U n i-1 + U n i 2 .
2) Computation of r 1 i ,r 2 i and λ 1 i , λ 2 i set respectively as the eigenvectors and eigenvalues of A i .

3) Computation of ∆t, such that ∆t ∆x ≤ CF L a i+1/2 .

4) Computation of U *

i-1/2 for i = 1, ..., N + 1 by solving the linearized Riemann problem.

In fact we have 3 cases:

• if λ 1 i ,λ 2 i <0 then U * i-1/2 = U n i . • if λ 1 i ,λ 2 i >0 then U * i-1/2 = U n i-1 . • if λ 1 i < 0, λ 2 i > 0 then for: (39)
             x < λ 1 i t one has U * i-1/2 = U n i-1 , x > λ 1 i t or x < λ 2 i t one has U * i-1/2 = U n i -(R -1 [U ]) 2 r 2 i = U n i-1 + (R -1 [U ]) 1 r 1 i , x > λ 2 i t one has U * i-1/2 = U n i , with R = (r 1 i |r 2 i ) and [U ] = U n i -U n i-1 . 5) Computation of F (U * i-1/2 ). 6) Computation of U n+1 i = U n i - ∆t ∆x (F n i+1/2 -F n i-1/2 ) for i = 1, ..., N .
We repeat this algorithm for the new level of time (t n+1 + ∆t), until we reach the required final time T . In what follows, the computation of high-order accurate numerical fluxes is reached by reconstructing left and right constant averaged values using a fifth-order WENO scheme, before applying the numerical flux. The only change is in the computation of the interface values U * i+1/2 which will depend on the high order reconstructed right and left values when solving the linearized Riemann problem.

3.2.2.

High order finite-volume scheme: WENO5-RK4. Considering numerical approximations of the hyperbolic system (26), we seek a numerical scheme that reach high order accuracy in smooth regions while avoiding the spurious oscillations around discontinuity. In fact, we aim at dealing with discontinuous initial data (dam-break problem) generating dispersive shock waves [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF][START_REF] Mitsotakis | On the Galerkin/finite-element method for the Serre equations[END_REF] which needs a special treatment at the numerical scheme level. This can be achieved by using a fifth-order accuracy WENO reconstruction for hyperbolic conservation laws, following [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF][START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF]. Second-order schemes tend to alter the dispersive properties of the model due to dispersive truncation errors. To prevent this in the study of dispersive waves, high order schemes are imperative [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF][START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF]. For the sake of simplicity, the reader is referred to [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] for more details concerning the high order discretization of the hyperbolic system.

Regarding time discretization, fourth-order explicit Runge-Kutta "RK4" method is used and one gets the "WENO5-RK4" scheme.

3.3. Finite difference scheme for the dispersive part. The splitting scheme is a mix between a finite volume discretization and a finite difference method. This mix induces a switching between cell-averaged values defined by the finite volume discretization and nodal values used by the finite difference discretization for each unknown and at each time step. Using fifth-order accuracy WENO reconstruction, one can approximate the nodal values (i.e finite difference unknowns) (U n i ) i=1,N +1 in terms of the cell-averaged values (i.e finite volume unknowns) (U n i ) i=1,N by the following relation:

(40) U n i = 1 30 U n i-2 - 13 60
U n i-1 + 47 60 U n i + 9 20 U n i+1 - 1 20 U n i+2 + O(∆x 5 ), 1 ≤ i ≤ N + 1,
One can easily recover the relation that allows to determine the cell-averaged values (U n i ) i=1,N in terms of the nodal values (U n i ) i=1,N +1 by inverting the relation [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]. The global order of the scheme is preserved. In fact, the formula is precise up to order O(∆x 5 ) terms. Before proceeding by the computation of S 2 , we recall first the remaining (dispersive part) of the equations, given in section 3.1. ( 41)

                   ∂ t ζ = 0, 1 + εαT [0] -ε 2 αT ∂ t v - 1 α ∂ x ζ + 1 α ∂ x ζ + (7 -5α) 45 ε 2 ∂ 4 x (1 + εαT [0]) -1 (∂ x ζ) + 2 3 ε 2 ∂ x ((∂ x v) 2 ) + 2 3 ε 2 ζ∂ 2 x (1 + εαT [0]) -1 (∂ x ζ) + ε 2 ∂ x ζ∂ x (1 + εαT [0]) -1 (∂ x ζ) = 0.
Inverting the operator (1 + εαT [0]) requires the discretization of (1 + εαT [0])X = B (i.e resolution of linear systems). The linear system is solved using the matrix division operator of Matlab, X = (1 + εαT [0]) \ B. This produces the solution using Gaussian elimination, without forming the inverse. Using an explicit Euler in time scheme, the finite discretization of the system (41) using classical finite difference methods leads to the following discrete problem:

(42)                              ζ n+1 -ζ n ∆t = 0, v n+1 -v n ∆t - 1 α D 1 (ζ n ) + 1 - εα 3 D 2 + ε 2 α 45 D 4 -1 1 α D 1 (ζ n ) + (7 -5α) 45 ε 2 D 4 1 - εα 3 D 2 -1 (D 1 (ζ n )) + 2 3 ε 2 D 1 ((D 1 v n ) 2 ) + 2 3 ε 2 ζ n D 2 1 - εα 3 D 2 -1 (D 1 (ζ n )) +ε 2 D 1 (ζ n )D 1 1 - εα 3 D 2 -1 (D 1 (ζ n )) = 0.
The matrices D 1 , D 2 and D 4 are respectively the classical centered discretizations of the derivatives ∂ x , ∂ 2 x , and ∂ 4

x given below. The spatial derivatives are discretized using fourth-order formulas, "DF4":

(∂ x U ) i = 1 12∆x (-U i+2 + 8U i+1 -8U i-1 + U i-2 ), (∂ 2 x U ) i = 1 12∆x 2 (-U i+2 + 16U i+1 -30U i + 16U i-1 -U i-2 ), (∂ 4 x U ) i = 1 6∆x 4 (-U i+3 + 12U i+2 -39U i+1 + 56U i -39U i-1 + 12U i-2 -U i-3 ).
A standard extension to fourth-order classical Runge-Kutta "RK4" scheme is used, and thus one obtains the "DF4-RK4" scheme.

Remark 4. At this stage, it is worth mentioning that the numerical scheme of the eB system with high order derivatives [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] is similar to the one developed for model [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF]. In fact, the hyperbolic part of the system is the same as in [START_REF] Fernández-Nieto | A hierarchy of dispersive layeraveraged approximations of Euler equations for free surface flows[END_REF], but the high order derivatives involved in the second equation of the remaining (dispersive) part should be treated accordingly. More precisely, third and fifth order derivatives are discretized using the following fourth-order formulas:

(∂ 3 x U ) i = 1 8∆x 3 (-U i+3 + 8U i+2 -13U i+1 + 13U i-1 -8U i-2 + U i-3 ), (∂ 5 x U ) i = 1 6∆x 5 (-U i+4 + 9U i+3 -26U i+2 + 29U i+1 -29U i-1 + 26U i-2 -9U i-3 + U i-4 ).
Boundary conditions are imposed using the method presented in Section 3.4.

Remark 5. An investigation of the dispersive properties of the splitting numerical scheme adopted here is done in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] for a Green-Naghdi type model. The splitting in time of the hyperbolic and dispersive parts is the main originality of this approach. Considering the semi-discretized (in time) version of the splitting scheme adopted in this paper, we believe that the corresponding dispersion relation can be determined following the same classical steps as in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF]. An extension to the fullydiscretized scheme is possible but very technical and would not lead to any important insight on the dispersive properties of the hyperbolic/dispersive splitting. To avoid repetition, we do not include in this paper the analysis of the discrete dispersion relation, the interested reader is referred to [4, Section 3.4.2].

3.4.

Boundary conditions. To close the differential problems, boundary conditions need to be imposed. Boundary conditions for both the hyperbolic and dispersive parts of the splitting scheme are treated by imposing suitable relations on both cell-averaged and nodal quantities. In this paper, we only consider periodic boundary conditions as it was already done in [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] for the study of internal waves.

For the hyperbolic part, "ghosts" cells are introduced respectively at the western and eastern boundaries of the domain. The imposed relations on the cell-averaged quantities are the following: • U -k+1 = U N -k+1 , and U N +k = U k , k ≥ 1, for periodic conditions on western and eastern boundaries.

For the dispersive part, we simply impose the boundary conditions on the nodal values located outside of the domain. In this way, we maintain centered formula at the boundaries, while keeping a regular structure in the discretized model:

• U -k+1 = U N -k+1 , and U N +k = U k , k ≥ 1, for periodic conditions on western and eastern boundaries. Remark 6. As we have already said, in the proposed model we do not try to totally avoid the computations of high order derivatives as in [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF][START_REF] Bassi | A hyperbolic reformulation of the Serre-Green-Naghdi model for general bottom topographies[END_REF][START_REF] Escalante | A general non-hydrostatic hyperbolic formulation for Boussinesq dispersive shallow flows and its numerical approximation[END_REF], instead we avoid the direct computation of high order derivatives on the flow variable ζ. There are still fourth order derivatives (but this time not directly on ζ) and a large stencil is still needed though. Approximating high order derivatives is one drawback of the proposed numerical scheme due to the large stencil needed and the development of high-order schemes is not an easy task. In the dispersive part of our model, the spatial derivatives are discretized using fourth-order formulas and one disadvantage is that the treatment of boundary conditions is more complex in that case. Moreover, in the case of variable topography, several challenging situations arise involving the design of robust, well-balanced and positive preserving numerical schemes for high-order partial differential equations (PDEs).

Numerical validations

This part is devoted to the numerical validations of the model and the numerical scheme. We begin by examining the numerical solution of the eB model for the case of solitary waves and show that the latter enjoys better approximate solution with respect to lower-order models.

Secondly, we study the propagation of a solitary wave solution with correctors of order O(ε 3 ) established in [START_REF] Khorbatly | Mathematical modeling and numerical analysis for the higher order Boussinesq system[END_REF]. We compare our numerical solution with an analytic one (up to an O(ε 3 ) remainder) at several times and show that our numerical scheme is very efficient and accurate.

As usual when dealing with a model in oceanographical science, one has to test the ability of the model and the numerical method for the test of a head on collision of counter propagating solitary waves. A very good agreement is observed.

An important fact to reveal is whether or not the improved model is pertinent. The third and fourth numerical tests reveal that in presence of large wave number, the choice of the parameter α and the high order derivatives factorization are crucial. In the fifth test, we study the dam-break problem supplemented by a comparison between the standard and extended Boussinesq models. This test is build to test the ability of the eB model and the numerical method to deal with irregular solution. A very good behavior is observed. In the last numerical test, we consider the well-known "Favre waves" resulting after the impact of a wave on a vertical wall. The numerical results obtained by the eB model are compared with experimental results of Favre [START_REF] Favre | Ondes de translation dans les canaux découverts[END_REF] and Treske [START_REF] Treske | Undular bores (Favre-waves) in open channels -experimental studies[END_REF] and a good agreement is observed.

In all the test cases, we use the "WENO5-DF4-RK4" discretization and a CFL number of 1. Our numerical investigations highlight that the time step restriction from the CFL condition of the finite-volume step is enough to ensure stability for the whole numerical method. 4.1. Numerical solitary wave solution of the extended Boussinesq system. Solitary waves consist of steadily translating disturbances where the nonlinear and dispersive effects counterbalance each other to create a permanent-form with a single crest solution. In this section, we numerically study solitary waves solutions of the eB model and compare it to other models including the accurate full solution of the water-waves model [START_REF] Dutykh | Efficient computation of steady solitary gravity waves[END_REF][START_REF] Tanaka | The stability of solitary waves[END_REF] which will be considered as an "exact" solution in our comparisons. We will show that the eB model enjoys a better approximate solitary wave solution (when compared to the "exact" solution) with respect to other lower-order models.

In the last decades, there have been several works on nonlinear PDEs modeling solitary waves. The famous Korteweg-de Vries (KdV) scalar equation or the coupled Boussinesq and Green-Naghdi evolution equations describe the shallow water waves and admit explicit families of solitary wave solutions [START_REF] Boussinesq | Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond[END_REF][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][START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF][START_REF] Chen | Exact solutions of various Boussinesq systems[END_REF]. The calculation of the exact analytic expressions for solitary wave solutions can be done in many ways, one of which, is the direct integration method. This method seeks traveling wave solutions. The PDE is replaced by an ordinary differential equation (ODE) by working in a traveling frame of reference. Hence, one looks for closed-form solutions in terms of special functions.

In [START_REF] Khorbatly | Mathematical modeling and numerical analysis for the higher order Boussinesq system[END_REF], a careful examination of the direct integration method revealed that the third order non linear ODE associated to the eB model does not admit an explicit analytic solution. The explicit solution of the eB system remains an open problem. To this effect, we employ the Matlab solver ode45 to compute numerically the solution of the ODE. The obtained numerical solutions are compared with the accurate full solution of the water-waves model computed using the Matlab script of Clamond and Dutykh [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] where a fast and precise approach for computing solitary waves solution is introduced. The fast and accurate Matlab script in [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] is limited to relatively small velocities. To this end, three values of velocity are used, namely c = 1.025, c = 1.01 and c = 1.002. We compare the obtained solutions with other models including the accurate full solution of the water-waves model, the original Green-Naghdi system (ζ GN ) and the Boussinesq system (ζ B ). The explicit solution of the original Green-Naghi model [START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF][START_REF] Su | Korteweg-de Vries equation and generalizations. III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF] reads: [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF] εζ

GN (x) = (c 2 -1) sech 2 3(c 2 -1) 4c 2 ε x = εc 2 ζ B (x) .
The waves are re-scaled so that the Boussinesq solution do not depend on c. Consistently, we set ε = 1. After re-scaling, Figure 6 shows clearly that the solitary waves tend towards the Boussinesq solution (ζ B ) as c -1 → 0. A zoom-in on the crest of half of the solitary waves shows that, of all models, the eB model has the best match with the full Euler system (water-waves) solution. The normalized l 2 -norm of the error as a function of c-1 is plotted in a log-log scale in Figure 7. The convergence rate is indeed quadratic for the extended Boussinesq model whereas it is only linear for the original Green-Naghdi model. Results are gathered in Table 1. This highlights the fact that the higher-order extended Boussinesq model have a better approximate solution when compared to lower-order models. 4.2. Propagation of a solitary wave solution with correctors. A careful examination revealed that the extended Boussinesq system (2) does not admit an exact solitary wave solution, see [START_REF] Khorbatly | Mathematical modeling and numerical analysis for the higher order Boussinesq system[END_REF]Section 4]. In order to validate our numerical scheme we use the explicit solution with correctors of order O(ε 3 ) found in [37, Section 5] that we disclose below. Such solitary waves are analytical solutions of the extended Boussinesq system (2) up to O(ε 3 ) remainders. Therefore, this family of solutions can be used as a validation tool for our present numerical scheme and its given by (ζ, v) with

(44) ζ = ζ 1 + ε 2 2 (ζ 0 2 + v 0 2 )(x -t) + (ζ 0 2 -v 0 2 )(x + t) + t 0 f (s, x -t + s)ds - t 0 f (s, x + t -s)ds , and (45)
v = v 1 + ε 2 2 (ζ 0 2 + v 0 2 )(x -t) -(ζ 0 2 -v 0 2 )(x + t) + t 0 f (s, x -t + s)ds + t 0 f (s, x + t -s)ds ,
where (ζ 1 , v 1 ) is the well known explicit solution of solitary traveling wave of the sB system (1) given by:

     ζ 1 (t, x) = a sech 2 k (x -ct) , v 1 (t, x) = cζ 1 (t, x) 1 + εζ 1 (t, x) ,
where k = 3a 4 and c = 1 1 -aε and a is an arbitrary chosen constant. The initial conditions

ζ 0 2 and v 0 2 are both given in C ∞ (R) and set ζ 0 2 = v 0 2 = exp - 3πx 10
2
. The function f (t, x) is defined by:

f (ζ 1 , v 1 ) = ∂ x ζ 1 ∂ x ∂ t v 1 + 2 3 ζ 1 ∂ 2 x ∂ t v 1 + 1 45 ∂ 4 x ∂ t v 1 + 1 3 ∂ x v 1 (v 1 ) xx -(v 1 ) 2 x .
In this test, we investigate the left to right propagation of a solitary wave initially centered at x 0 = 20, of amplitude a = 0.2. The computational domain length is L = 100 and discretized with 1600 cells. The solitary wave is initially far from boundaries, thus the periodic boundary conditions do not affect the computation. The water surface profile of our numerical solution provided by the model ( 14) with α = 1, is compared with the analytical one given by ( 44)-( 45) at several times using the fifth order discretization "WENO5-DF4-RK4". An excellent agreement between numerical and analytical solutions is observed in Figure 8. The amplitude and shape of the computed solitary wave are accurately preserved during the propagation, indicating an accurate discretization of the governing equations in both space and time. To complete the picture and assess the convergence of our numerical scheme, we compute the numerical solution for this particular test case for an increasing number of cells N , over a duration T = 1. We start with N = 400 number of cells (δ x = L N = 0.25) and successively multiply the number of cells by two. The relative errors E L 2 (ζ) and E L 2 (v) on the water surface deformation and the averaged velocity are computed at t = 1, using the discrete L 2 norm . 2 :

(46) E L 2 (ζ) = ζ num -ζ sol 2 ζ sol 2 ; E L 2 (v) = v num -v sol 2 v sol 2 ,
where (ζ num , v num ) are the numerical solutions and (ζ sol , v sol ) are the analytical ones coming from ( 44)- [START_REF] Madsen | Surf zone dynamics simulated by a Boussinesq type model. Part I. Model description and cross-shore motion of regular waves[END_REF]. Results are presented in Table 2 and Figure 9 where E L 2 (ζ) and E L 2 (v) are plotted against δ x in log scales, for the considered relative amplitude a = 0.2. Very accurate results are obtained, indicating that the employed numerical method is capable of computing in a stable way the propagation of a solitary wave. Moreover, computing a linear regression on all points yields a slope equal to 2.33 for ζ and 2.34 for the averaged velocity v. This result sounds rational because the global (time and space) order of our scheme may be limited by the order of the splitting method used here, which is of order two as already discussed by Bonneton et al.

in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF]. 4.3. Head on collision of counter propagating solitary waves. A standard nonlinear test case for numerical methods is the interaction of solitary wave. In this numerical test, we study an important phenomenon in the study of nonlinear dispesive waves, the head on collision of two counter propagating waves with different amplitudes. We used solitary wave solutions with correctors of order O(ε 3 ) for the eB system established in [START_REF] Khorbatly | Mathematical modeling and numerical analysis for the higher order Boussinesq system[END_REF] and defined in ( 44)- [START_REF] Madsen | Surf zone dynamics simulated by a Boussinesq type model. Part I. Model description and cross-shore motion of regular waves[END_REF]. We consider two solitary waves centered at x = -50 and x = 50 at t = 0 on a spatial domain L = 200 with a constant depth h 0 = 2, see Figure 10. The solitary wave centered at x = -50 travels to the right with a speed c s,1 = 1.0206 and an initial amplitude a 1 = 0.4 while the one centered at x = 50 travels to the left with a speed c s,2 = 1.0102 and an initial amplitude a 2 = 0.2. The domain is discretized using 1200 cells and periodic boundary conditions are imposed. The numerical solutions are computed using model [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] with α = 1. The collision of the two waves starts at about t = 43, see Figure 10.

After the interaction, each wave continue moving in its own direction and turn up to be unaffected by the collision, see Figure 11. A proper description of the distinctive nature of nonlinear interactions is illustrated when zooming at the oscillating dispersive tails of very small amplitude appearing at the center of the domain at t = 70 in Figure 12. One can also observe two dispersive tails with smaller amplitudes located to the left and right boundaries. The generation of such dispersive tails is due to the O(ε 3 ) remainder terms as mentioned in the beginning of Section 4.2.

The high precision of our numerical scheme is verified after accurately capturing this phenomenon and inducing similar observations to earlier works [START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF][START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF] where the head-on collision is carried out. Breaking of a regular heap of water with a large wave number. In this numerical test, we highlight the importance of factorizing high order derivatives present in the improved eB model ( 14) together with the appropriate choice of the parameter α in improving the frequency dispersion in high frequency regimes. To this end, we consider a sufficiently regular heap of water with a large wave number represented by the initial data:

ζ(0, x) = 0.7e -80x 2 , v(0, x) = 0, (dashed lines) with a domain of computation x ∈ (-2, 2) discretized with 512 cells and under periodic boundary conditions. The non-linearity parameter is set as follow: ε = 0.1 (non-dimensional setting). Our numerical solutions are computed using models (10) (without factorization) and ( 14) (with factorization). We compare our numerical solutions with the numerical solutions computed using the Matlab script of Duchêne, Israwi and Talhouk [START_REF] Duchêne | A new class of two-layer Green-Naghdi systems with improved frequency dispersion[END_REF] and with the lower-order GN-CH model obtained in [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF]. In [START_REF] Duchêne | A new class of two-layer Green-Naghdi systems with improved frequency dispersion[END_REF], the original Green-Naghdi (GN) model describing a two-layer flow 2 is improved in terms of frequency dispersion by introducing a new class of tailored GN models with a slight modification of the dispersion components using a class of Fourier multipliers. In particular, an "improved" GN model is derived sharing the same dispersion relation as the full Euler (FE) system. These type of models are commonly called full dispersion models [START_REF] Bona | Asymptotic models for internal waves[END_REF]. The "improved" GN model is a full dispersion model and its numerical solution will be used as a reference solution. Eventually, the goal of this test case is to show that the optimized eB model ( 14) may provide some results which are in the same league of full dispersion models.

2 One can easily recover the one-layer configuration by setting γ = 0 and δ = 1. 14) and the eB model [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] with the "improved" GN model and the GN-CH model [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] at t = 3.

Figure 13(a) shows when α is chosen appropriately as discussed in section 2.3, namely α opt = 1.0610, our numerical solution computed over a sufficient duration t = 3 using model [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] behave similarly to the one computed with the "improved" GN model (sharing the same dispersion relation as the FE system). In contrary, when choosing α opt = 1 or when using the model [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] to compute the solution (Figure 13 (b)), the behavior is different than the "improved" GN model. Note that the GN-CH model has an improved frequency dispersion due to the careful choice of the parameter α. Nevertheless, the numerical solution computed using the GN-CH model is far from the numerical solution of the "improved" GN model. The observed agreement in high frequency regime between the numerical solutions of the "improved" GN model and the eB model ( 14) rather than the GN-CH model is due to the factorized high order dispersion terms existing in [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF]. In fact, the eB model ( 14) is precise up to O(ε 3 ) order and thus contains factorized high-order dispersive terms that do not exist in the GN-CH model. As already stated in Section 2.3, one can see that the choice of an optimal value of α when using the model [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] has no beneficial effect due to the high frequency regime setting. This numerical test confirms that the model [START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] has a range of applicability limited to k ≤ 1 and thus has poor dispersion properties in intermediate and large wave numbers regime. 4.5. Breaking of a regular heap of water with a small wave number. In this numerical test, we consider the breaking of a sufficiently regular heap of water with a small wave number whose initial data is:

ζ(0, x) = 0.7e -0.4x 2 , v(0, x) = 0, (dashed lines) within a domain of computation x ∈ (-2, 2) discretized with 512 cells and under periodic boundary conditions. The non-linearity parameter is set as follow: ε = 0.5 (non-dimensional setting). Our numerical solutions are computed using models ( 10) and ( 14) and compared with the numerical solution of the "improved" Green-Naghdi model and the numerical solution of the GN-CH model [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] over a sufficient duration t = 3. The parameter α is fixed as 1 since varying α does not yield significant improvements. In fact, both eB models ( 10) and ( 14) have an equivalent dispersion relation to the one of the full Euler system for small wave numbers and the choice of α does not play any role in the leading terms. Indeed, Figure 14 shows a fairly good agreement between the solutions of the eB models (10) (yellow line) and ( 14) (blue line), and the solution of the GN-CH model [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] (red line) and the one of the "improved" Green-Naghdi model (green line). This confirms the fact that, in small wave numbers regime, every aforementioned model behave similarly and enjoy similar dispersive properties as the one of the full Euler system. 4.6. Dam-break problem. Dealing with non-regular solutions needs a special treatment at the numerical scheme level. Earlier works [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF] have shown that the use of high-order schemes in dispersive waves study is necessary to prevent the corruption of the dispersive properties of the model. In general, dispersive shock waves are generated due to the dispersive effects [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF][START_REF] Mitsotakis | On the Galerkin/finite-element method for the Serre equations[END_REF] when considering discontinuous initial data. In this numerical test, we implement a dam break problem in order to investigate the performance of our numerical scheme in handling non-regular solutions. We study the dam-break problem in the extended (eB) and standard (sB) Boussinesq models. We consider the following initial data: [START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] ζ(0, x) = a[1 + tanh(250 -|x|)], v(0, x) = 0, with a = 0.2091 m defined on the computational domain x ∈ (-700, 700) discretized using 2800 cells and imposed under periodic boundary conditions. The solutions of the eB model are computed using the dimensional version of the eB model with factorized high order derivatives knowing that similar results were obtained when using the model without factorization. To this end, we set ε = 1 and add the gravity term to the equations as needed. We make this choice in order to test the scheme in the challenging conditions of the strongly dispersive regime. In this test, the choice of α is not important, thus we choose α = 1. The solutions of the sB model are computed accordingly using the same splitting scheme. The dam break wave shapes of the sB and eB models are shown at different times in Figure 15. The initial data break and generate dispersive shock waves. A closeup on the profiles at t = 30s and t = 65s shows two dispersive shock waves counter-propagating on both sides of the "dam", and two rarefaction waves moving in the direction of the center. The dispersive tail generated by the eB model is larger and have higher amplitude oscillations when compared to the respective tail generated by the sB model. This result is consistent with the fact that the eB model contains high order nonlinear dispersive terms not present in the sB model. This numerical test shows that our high-order numerical scheme was able to capture accurately the dispersive shock waves phenomenon. Dispersive shock waves in a large class of dispersive shallow water models were carried out in several works [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF][START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF][START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF][START_REF] Mitsotakis | On the Galerkin/finite-element method for the Serre equations[END_REF] and show a good agreement with our numerical simulation. Remark 7. The splitting strategy may limit the whole method to second order accuracy, however, for the study of dispersive waves, it is necessary to use high-order schemes to prevent the corruption of the dispersive properties of the model by some dispersive truncation errors linked to second-order schemes. The reader is referred to [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF], where the study of dam-break problem in the GN-CH model is supplemented by a comparison between second (MUSCL-RK2) and fifth (WENO5-RK4) order accuracy methods. Although a splitting scheme of order 2 was adopted therein, the dispersive properties were clearly corrupted by the second-order method while the higher-order method was able to deal with discontinuous initial data and capture the rapid oscillations in dispersive shock waves.

4.7. Favre waves. In this numerical test, we consider the well-known "Favre waves" resulting after the impact of a wave on a vertical wall. Experimental investigations of those waves were addressed for the first time by Favre [START_REF] Favre | Ondes de translation dans les canaux découverts[END_REF] in a rectangular channel. Similar experiments have been performed by Treske in [START_REF] Treske | Undular bores (Favre-waves) in open channels -experimental studies[END_REF] and Soares Frazao and Zech in [START_REF] Frazao | Undular bores and secondary waves -experiments and hybrid finite-volume modelling[END_REF] in open channels. This problem was numerically studied more recently in [START_REF] Soares-Frazão | A second-order semi-implicit hybrid scheme for one-dimensional Boussinesq-type waves in rectangular channels[END_REF][START_REF] Tissier | Nearshore dynamics of Tsunami-like undular bores using a fully nonlinear Boussinesq model[END_REF]. Due to dispersion, the uniform free surface flow impacting a wall reflects and free surface undulations (called "Favre waves", see Figure 16) appear. The leading wave has a maximum amplitude a max and a minimum amplitude a min and is followed by waves of decreasing heights. The jump height is denoted by a m and D denotes the constant velocity of the wave front. In this kind of experiments, the Froude number F is defined as the ratio between the wave speed v 0 -D and the celerity √ gh 0 . Consequently, according to conservation of mass and momentum, a relation between the Froude number and the upstream and downstream water depths h 0 and h 0 + a m can be obtained (see the Bélanger formula Eq. (4.10) in [START_REF] Gavrilyuk | Spilling breakers in shallow water: applications to Favre waves and to the shoaling and breaking of solitary waves[END_REF]): In what follows, we will compare the numerical results obtained by the extended Boussinesq model for the "Favre waves" problem with the experimental results of Favre [START_REF] Favre | Ondes de translation dans les canaux découverts[END_REF] and Treske [START_REF] Treske | Undular bores (Favre-waves) in open channels -experimental studies[END_REF]. On a computational domain x ∈ (0, 300), we consider a uniform initial profile defined by ζ(0, x) = 0 and an impact velocity v(0, x) = v 0 related to the relative Froude number F by the formula [START_REF] Gavrilyuk | Spilling breakers in shallow water: applications to Favre waves and to the shoaling and breaking of solitary waves[END_REF]:

(48) h 0 + a m h 0 = √ 1 + 8F 2 -1 2 .
v 0 = gh 0 F - 1 + √ 1 + 8F 2 4F ,
where g = 10 m/s and h 0 = 1 m. Reflecting boundary conditions are used on the right wall at x = 300 m. The solutions of the eB model are computed using the dimensional version of the eB model ( 14) with factorized high order derivatives. To this end, we set ε = 1 and add the gravity term to the equations as needed. In this test, the choice of α is not important, thus we choose α = 1. The numerical scheme adopted in this paper for the extended Boussinesq equations can be used to solve this problem until some critical impact velocity determined in terms of the relative Froude number F . For higher Froude numbers (F >≈ 1.3), the transition from the undular bore to the bore consisting of a steep front (wave breaking) occurs. At this stage, a special treatment of the numerical scheme based on the adopted splitting strategy is required in order to handle wave breaking (see introduction of Section 3). This is not investigated in this paper and is left to a future work (see [START_REF] Tissier | Nearshore dynamics of Tsunami-like undular bores using a fully nonlinear Boussinesq model[END_REF] for modeling of breaking waves).

In Figure 17, we compare the numerical results obtained by the eB model for different mesh sizes (N = 1000, 2000, and 3000) with the numerical results of the Serre Green-Naghdi (SGN) equations obtained by the method in [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF]. The comparison is done at time t = 54 s for a Froude number F = 1.16 that corresponds to the following upstream velocity v 0 = 0.6490 m/s. The solid red line corresponds to the numerical solution of the SGN equations obtained by the method [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF] on a 2000 cell mesh. A good agreement is observed. Our results for 2000 and 3000 cells are almost superposed, thus the convergence is guaranteed. The first wave amplitude is well estimated with a finer mesh, however a small discrepancy is observed in the prediction of the jump height a m between the numerical solutions of the eB model and the SGN model. For a Froude number F = 1.16, the corresponding jump height obtained using equation ( 48) is a m = 0.215 m which corresponds exactly to the one predicted by the eB model, see Figure 17. Thus, an accurate prediction of the jump height a m is provided by the eB model rather than the SGN equations. In addition we compare the amplitudes of undular bores obtained by the eB model with experimental data of Favre [START_REF] Favre | Ondes de translation dans les canaux découverts[END_REF] and Treske [START_REF] Treske | Undular bores (Favre-waves) in open channels -experimental studies[END_REF]. The computations are performed for different Froude numbers from the interval F ∈ [1.02, 1.36]. The maximum a max , the minimum a min amplitudes of the leading wave and the jump height a m are taken at t = 54 s with N = 2000. Figure 18 shows that the results obtained by the eB model are in good agreement with experimental data until the wave breaking occurs corresponding to the Froude number about 1.25. After this critical value, the transition from the undular bore to the bore consisting of a steep front (breaking bore) occurs and our numerical scheme is no more valid since it does not handle wave breaking.

Conclusion

In this work, a numerical model is developed for a class of higher-ordered/extended Boussinesq (eB) equations describing the propagation of water-waves over a flat bottom. A reformulation with the same order of precision that avoids the calculation of high order derivatives on the surface deformation is proposed including a dispersion correction parameter α to be precisely chosen which improve the linear dispersion properties. Insights about linear stability and dispersion optimization are provided showing that the improvement is significant when factorizing every high order derivative. In fact, factorizing only the fifth order derivative admits a destabilizing effect. The eB model with factorized high order derivatives provide linear dispersion characteristics which are accurate for wave numbers k up to 10 which are superior to lower-order models.

A second order time splitting scheme is then proposed relying on a combination of finite volume and finite difference methods. The hyperbolic part of the equations is discretized using a highorder finite volume WENO scheme, while the dispersive part is treated using classical high-order finite differences.

Finally several numerical computations are exhibited validating the model and the numerical methods. We began by examining the numerical solution of the eB model for the case of solitary waves and compared it to other lower-order models. The eB model is found to have a better approximate solution. The propagation of a solitary wave solution with correctors is then considered allowing to study the accuracy and convergence properties of the proposed numerical scheme. In the following case, the interaction of two counter propagating solitary waves is considered as a standard nonlinear test case showing the high precision of our scheme. The optimized higherorder eB formulation with factorization of high order derivatives provide some results which are on par with full dispersion models. In particular, the breaking of a regular heap of water with both large and small wave numbers are studied highlighting the importance of factorizing high order derivatives in improving the frequency dispersion in high frequency regimes. The dam-break problem supplemented by a comparison between the standard and extended Boussinesq models is studied. The eB model containing high order nonlinear dispersive terms not present in the sB model generates larger dispersive tails with higher amplitude oscillations. This test shows that the dispersive properties of the model are well captured thanks to the high order accuracy of the numerical scheme. In the last numerical test case, a comparison with experimental data was presented for the study of "Favre waves". The proposed scheme reproduces the measurements with a good agreement.

Following the steps of this work, next steps may cover the occurrence of variable topography to seriously discriminate high order models and two dimensional extension to study more real-life cases.

Appendix A. Dispersion relation of the eB model [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF] This appendix is devoted to the computation of the dispersion relation associated with the new eB model [START_REF] Clamond | Fast accurate computation of the fully nonlinear solitary surface gravity waves[END_REF]. We start by investigating the linear behavior of small perturbation ( ζ, ṽ) to a constant state solution (ζ, v). The linear equations governing these perturbations are (49) The second equation of [START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF] .

             ∂ t ζ + vε∂ x ζ + h∂ x ṽ = 0, 1 + εαT [0] -ε 2 αT ∂ t ṽ + εv∂ x ṽ + α -1 α ∂ x ζ + 1 α ∂ x ζ + (7 -5α)ε 2 45 ∂ 4 x (1 + εαT [0]) -1 (∂ x ζ) +

Figure 1 .

 1 Figure 1. Errors on linear phase velocity (red) and group velocity (blue) for the eB model (10).

Figure 2 .

 2 Figure 2. Phase and group velocities weighted averaged error as a function of α for 0 ≤ k ≤ 10. The eB model (14) is in solid line, the Green-Naghdi model in the Camassa-Holm regime (GN-CH)[START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] is in dashes.

Figure 3 .

 3 Figure 3. Errors on linear phase velocity (red) and group velocity (blue). The reference from Stokes theory (black solid line), the eB model (14) (α = 1.0610) in solid lines, the eB model (14) (α = 1) in dashes, the GN-CH model[START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] (α = 1.0800) in dots.

Figure 4 .

 4 Figure 4. Comparison at different times between the solutions of the models (10) (blue line), (14) (red line) and (21) (green line) in high frequency regime.

Figure 6 .

 6 Figure 6. Comparison of the solitary waves solutions.

Figure 7 .

 7 Figure 7. Errors as a function of c -1 (log-log plot).

Figure 8 .

 8 Figure 8. Propagation of a solitary wave: water surface profiles at t = 0, 10, 30, 50 and 70.

Figure 9 .

 9 Figure 9. Propagation of a solitary wave: L 2 -error on the water surface deformation and the averaged velocity for a = 0.2.

Figure 10 .

 10 Figure 10. Head on collisions: surface wave shape at t = 0, 43, 46 and t = 49.

Figure 11 .

 11 Figure 11. Head on collisions: surface wave shape at t = 53, 55, 58 and 60.

Figure 12 .

 12 Figure 12. Head on collisions: surface wave shape at t = 70.

( a )

 a Comparison of the numerical solutions of the eB model (14) (blue) with the "improved" GN model (green) and the GN-CH model (red). (b) Comparison of the numerical solutions of the eB model (10) (blue) with the "improved" GN model (green) and the GN-CH model (red).

Figure 13 .

 13 Figure 13. Comparison of the numerical solutions of the eB model (14) and the eB model[START_REF] Castro | Fully nonlinear long-wave models in the presence of vorticity[END_REF] with the "improved" GN model and the GN-CH model[START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] at t = 3.

Figure 14 .

 14 Figure 14. Comparison of the numerical solutions of model (14) (blue line) and model (10) (yellow line) and the GN-CH model [5] (red line) with the "improved" Green-Naghdi model (green line) at t = 3.

Figure 15 .

 15 Figure 15. Dam break: wave shape at different times, comparison between the numerical solution of the eB model (solid red line) and sB model (dashed blue lines).

Figure 16 .

 16 Figure 16. A sketch of Favre waves.

Figure 17 .

 17 Figure 17. Comparison of Favre waves at time t = 54 s for the Froude number F = 1.16. Red solid line corresponds to the numerical solution of the SGN equations obtained by the method [24]. The results obtained with the eB model are shown for different mesh sizes: 1000 (blue solid line), 2000 (black dashed line), 3000 (green dashed line).

Figure 18 .

 18 Figure 18. Amplitude of undular bores for different Froude numbers.

2ε 2 3 ζ∂ 2 x ( 1 +

 21 εαT [0]) -1 (∂ x ζ) = O(ε 3 ), where h = 1 + εζ, T [0]w = -Looking for the corresponding plane wave solutions of the form (ζ 0 , v 0 )e i(kx-wt) with k the spatial wave number and w the time pulsation, one finds the dispersion relation. From the first equation of (49) one obtains: (50) -iw ζ + v(ikε ζ) + h(ikṽ) = 0 ⇒ ṽ = ζ(w -εkv) kh .

  , hence extending the range of values of the overall surface deformation ζ for which the condition is satisfied. This ensures a numerical stability in most of the situations considered for applications (see Figure

			19) becomes
	negative whenever ζ > numbers, under the condition that ζ < 2k 4 + 45 30k 2 . Thus, as provided by (19), w remains real for large wave 2k 4 + 45 30k 2 . In the majority of the applications we have in mind, the overall surface deformation ζ does not go beyond 2k 4 + 45 30k 2 , and this condition is satisfied. Actually, as k gets large (for instance k ≈ 10) in high frequency regime, the upper
	bound of ζ gets also large (≈	k 2 15	= 6.67)

Table 1 .

 1 The normalized l 2 -norm of the error for the Boussinesq, original GN and eB models.

		Boussinesq		original GN		ex-Bouss
	c -1	Error	Conv. rate	Error	Conv. rate	Error	Conv. rate
	0.025	0.0459	-	0.0058	-	0.0019	-
	0.01	0.0188	0.9743	0.0023	1.0361	3.07 × 10 -4	1.9984
	0.002	0.0038	0.9864	4.45 × 10 -4	1.0186	1.23 × 10 -5	1.9990

Table 2 .

 2 Propagation of a solitary wave: relative L 2 -error table for the conservative variables.

	N	E L 2 (ζ)	Conv. rate E L 2 (v)	Conv. rate
	400	3.50 × 10 -3 -	3.33 × 10 -3 -
	800	9.32 × 10 -4 1.9070	8.29 × 10 -4 2.0064
	1600 2.05 × 10 -4 2.0466	1.70 × 10 -4 2.1433
	3200 3.23 × 10 -5 2.2449	2.48 × 10 -5 2.3479
	6400 4.79 × 10 -6 2.3873	3.50 × 10 -6 2.4845
	12800 1.44 × 10 -6 2.3344	1.49 × 10 -6 2.3436