A new class of higher-ordered/extended Boussinesq system for efficient numerical simulations by splitting operators
AA NEW CLASS OF HIGHER-ORDERED/EXTENDED BOUSSINESQ SYSTEM FOREFFICIENT NUMERICAL SIMULATIONS BY SPLITTING OPERATORS
ST´EPHANE GERBI AND RALPH LTEIF
Abstract.
In this work, we numerically study the higher-ordered/extended Boussinesq system de-scribing the propagation of water-waves over flat topography. An equivalent suitable reformulationis proposed, making the model more appropriate for the numerical implementation and significantlyimproved in terms of linear dispersive properties in high frequency regimes due to the suitable adjust-ment of a dispersion correction parameter. Moreover, we show that a significant interest is behindthe derivation of a new formulation of the higher-ordered/extended Boussinesq system that avoids thecalculation of high order derivatives existing in the model. We show that this formulation enjoys anextended range of applicability while remaining stable. We develop a second order splitting schemewhere the hyperbolic part of the system is treated with a high-order finite volume scheme and thedispersive part is treated with a finite difference approach. Numerical simulations are then performedto validate the model and the numerical methods.
Contents
1. Introduction 11.1. The water-wave equations 11.2. Boussinesq systems 31.3. Outline of the paper 41.4. A fully justified extended Boussinesq system. 42. Reformulation of the extended Boussinesq system 52.1. A one-parameter family of extended Boussinesq equations 52.2. Reformulation of the extended Boussinesq equations (2.6) 62.3. Choice of the parameter α Introduction
The water-wave equations.
In this paper, we investigate the one-dimensional flow of the free surfaceof a homogeneous, immiscible fluid moving above a flat topography z = − h . The horizontal and verticalvariables are denoted respectively by x ∈ R and z ∈ R and t ≥ ζ ( t, x ) denoting the variation with respect to its Mathematics Subject Classification.
Key words and phrases.
Water waves, Boussinesq system, higher-order asymptotic model, splitting scheme, hybridfinite volume/finite difference scheme. a r X i v : . [ m a t h . A P ] F e b est state z = 0 (see Figure 1). The fluid occupies the strictly connected ( ζ ( t, x ) + h >
0) domain Ω t at time t ≥ t = { ( x, z ) ∈ R ; − h ≤ z ≤ ζ ( t, x ) } . Figure 1.
One-dimensional flat bottom fluid domain.The fluid is considered to be perfect, that is with no viscosity and only affected by the force of gravity.We also assume the fluid to be incompressible and the flow to be irrotational so that the velocity field isdivergence and curl free. We denote by ( ρ, V ) the constant density and velocity field of the fluid. Thefirst boundary condition at the free surface expresses a balance of forces. Kinematic boundary conditionsare considered assuming that both the surface and bottom are impenetrable, that is no particle of fluidcan cross. The set of equations describing the flow is now complete and is commonly known as the fullEuler equations:(1.1) ∂ t V + V · ∇ x,z V = − g −→ e z − ∇ x,z P in ( x, z ) ∈ Ω t , t ≥ , ∇ x,z · V = 0 in ( x, z ) ∈ Ω t , t ≥ , ∇ x,z × V = 0 in ( x, z ) ∈ Ω t , t ≥ ,P | z = ζ ( t,x ) = 0 for t ≥ , x ∈ R ,∂ t ζ − (cid:112) | ∂ x ζ | n ζ · V | z = ζ ( t,x ) = 0 for t ≥ , x ∈ R , − V · −→ e z = 0 at z = − h , t ≥ , where n ζ = 1 (cid:112) | ∂ x ζ | ( − ∂ x ζ, T denotes the upward normal vector to the free surface. The theoret-ical study of the above system of equations is extremely difficult due to its large number of unknownsand its time-dependent moving domain Ω t . In fact, we have a free boundary problem, in other wordsthe domain is itself one of the unknowns.Using the assumption of irrotational velocity field, one can express the latter as the gradient ofa potential function ϕ . This potential satisfies the Laplace equation inside the fluid, ∆ x,z ϕ = 0 in( x, z ) ∈ Ω t . Consequently, the evolution of the velocity potential is written now using Bernoulli’sequation. Although the system now is simpler, a free boundary problem still exists. To get over thisobstacle, Craig and Sulem [11, 12] had an interesting idea following Zakharov work [45], consisting of areformulation of the system of equations (1.1) using the introduction of a Dirichlet-Neumann operator,thus reducing the dimension of the considered space and the unknowns number. Denoting by ψ the traceof the velocity potential at the free surface, ψ ( t, x ) = ϕ ( t, x, ζ ( t, x )) = ϕ | z = ζ , the Dirichlet-Neumannoperator is introduced G [ ζ ] ψ = − (cid:0) ∂ x ζ (cid:1) · (cid:0) ∂ x ϕ (cid:1) | z = ζ + (cid:0) ∂ z ϕ (cid:1) | z = ζ = (cid:113) (cid:12)(cid:12) ∂ x ζ (cid:12)(cid:12) (cid:0) ∂ n ϕ (cid:1) | z = ζ here ϕ is defined uniquely from ( ζ, ψ ) as a solution of the following Laplace problem: ∂ x ϕ + ∂ z ϕ = 0 in − h < z < ζ ( t, x ) ,∂ z ϕ | z = − h = 0 ,ϕ | z = ζ = ψ ( t, x ) . with ∂ n = n. ∇ x,z the normal derivative in the direction of the concerned vector n . Thus, the evolution ofonly the two variables ( ζ, ψ ) located at the free surface characterize the flow. Although Zakharov’s refor-mulation resulted in a reduced system of equations, the description of these solutions from a qualitativeand quantitative point of view remains very complex. A remedy for this situation requires the construc-tion of simplified asymptotic models whose solutions are approximate solutions of the full system. Theseapproximate models allow to describe in a fairly precise way the behaviour of the complete system in aspecific physical regime. This requires a rescaling of the system in order to reveal small dimensionlessparameters which allow to perform asymptotic expansions of non-local operators (Dirichlet-Neumann),thus ignoring the terms whose influence is minimal. The order of magnitude of these parameters makes itpossible to identify the considered physical regime. We start by introducing respectively the commonlyknown nonlinear and shallowness parameters:0 ≤ ε = ah = amplitude of the wavereference depth ≤ , ≤ √ µ = h λ = reference depthwave-length of the wave < , The introduction of these dimensionless parameters allows to write the dimensionless form of the fullEuler system:(1.2) ∂ t ζ − µ G µ [ εζ ] ψ = 0 ,∂ t ψ + ζ + ε | ∂ x ψ | − εµ ( µ G µ [ εζ ] ψ + ∂ x ( εζ ) · ∇ ψ ) ε µ | ∂ x ζ | ) = 0 . where the dimensionless Dirichlet-Neumann operator G µ [ εζ ] is now defined by: G µ [ εζ ] ψ = − µ (cid:0) ∂ x ζ (cid:1) · (cid:0) ∂ x ϕ (cid:1) | z = εζ + (cid:0) ∂ z ϕ (cid:1) | z = εζ = (cid:113) µε (cid:12)(cid:12) ∂ x ζ (cid:12)(cid:12) (cid:0) ∂ n ϕ (cid:1) | z = εζ with ϕ solving the dimensionless Laplace problem µ∂ x ϕ + ∂ z ϕ = 0 in − < z < εζ ( t, x ) ,∂ n ϕ | z = − = 0 ,ϕ | z = εζ = ψ ( t, x ) . The essential ingredient in constructing approximate models consists in the asymptotic expansion of thenon-local operators. Replacing these operators in the full Euler system (1.2) by their approximations,one obtains the different asymptotic models. Many asymptotic models can be derived in the shallowwater regime ( µ (cid:28)
1) based on the smallness assumptions made on the nonlinear parameter ε .1.2. Boussinesq systems.
In what follows, we are interested in a specific long wave regime where ε isconsidered of the same order as µ ( ε ∼ µ ). In this regime, Boussinesq derived in [5, 6] a weakly nonlinearmodel bearing his name. In what follows we refer to it as the “original” or “standard” Boussinesq system.Using the horizontal depth-mean velocity(1.3) v ( t, x ) := 11 + εζ ( t, x ) (cid:90) εζ ( t,x ) − ∇ ϕ ( t, x, z ) dz , the standard Boussinesq equations reads:(1.4) ∂ t ζ + ∂ x (cid:0) (1 + εζ ) v (cid:1) = 0 , (cid:0) − ε ∂ x (cid:1) ∂ t v + ∂ x ζ + εv∂ x v = O ( ε ) . This model can be derived from the
Green-Naghdi (GN) equations (see [20]) by neglecting all terms oforder O ( ε , µε, µ ). Equivalent Boussinesq systems enjoying a better mathematical structure or physiscalpropreties have been studied and derived extensively in the litterature, see for instance [1, 7, 8, 26, 29,35, 38, 39].To improve dispersion characteristics, shoaling and nonlinear properties, some efforts have been doneto extend the Boussinesq equations to a higher order of dispersion terms, namely to include dispersiveterms of order µ (higher-order dispersive effects), see for instance [19, 25, 32]. Neglecting the terms f order O ( µ ) while keeping the O ( µ ) terms in the equations prompt a more accurate system. Thestandard Boussinesq equations are restricted by containing only weak dispersion and nonlinearity. Thisnormally limits precise applications to a small zone moderately exterior to the surf zone. Significantimprovement can be made in the sense of expanding the application range and covering the rangefully from deep water into the surf zone by adding higher order terms. Extended Green-Naghdi (eGN)equations were derived and studied in [22, 24, 33, 34]. These equations are accurate up to the order O ( µ ) while the full nonlinearity is preserved. In fact, no assumption is made on the nonlinearityparameter. The extensive length of the eGN equations due to the incorporation of very high orderderivatives introduce numerical complexity and have a high computational cost. In order to bring somesimplifications with respect to the eGN equations, we decide to make a smallness assumption on thenonlinearity parameter, namely we consider weak nonlinearity ε ∼ µ . Following the work in [22, 24,33, 34], the extended Boussinesq (eB) equations, can be easily obtained by dropping all terms of order O ( µε , µ ε, ε ). Thus, one can write the weakly nonlinear Boussinesq system including higher orderdispersive effects as follows(1.5) (cid:40) ∂ t ζ + ∂ x ( hv ) = 0 , (1 + ε T [ h ] + ε T ) ∂ t v + ∂ x ζ + εv∂ x v + ε Q v = O ( ε ) , where h = 1 + εζ is the non-dimensionalised height of the fluid and denote by T [ h ] w = − h ∂ x (cid:0) h ∂ x w (cid:1) , T w = − ∂ x w, Q v = − ∂ x (cid:0) vv xx − v x (cid:1) . This paper is devoted to the numerical study of the eB system (1.5). To this end, an equivalent (in termsof precision) suitable reformulation is proposed, making the model more appropriate for the numericalimplementation.1.3.
Outline of the paper.
In this work, we firstly propose an equivalent reformulation of the extendedBoussinesq model (1.5) up to the third order, that makes the model more appropriate for the numericalimplementation and significantly improved in terms of linear dispersive properties in high frequencyregimes due to the suitable adjustment of a dispersion correction parameter. The reformulation isperformed via two methods : • with the factorization of high order derivatives, • without factorization of high order derivatives.We will show that the improvement is significant in the dispersive properties of the model with factoriza-tion of high order derivatives with an appropriate choice of optimal alpha in the high frequency regime.We then study the stability of the two models with and without factorization and we will show thatfactorizing only the fifth order derivative presented in the second model equation admits a destabilizingeffect : we need to factor all high order derivatives.Secondly, we propose a suitable Strang splitting of operators to solve the improved model : a hyper-bolic part representing the Nonlinear Shallow Water system and a dispersive part representing the highorder derivatives. The hyperbolic part of the system is treated with a high-order finite volume schemewhereas the dispersive part is treated with a finite difference method at the same order. To this end, areconstruction of nodal unknowns and centred unknowns is presented.Finally, numerical validations are presented, showing the interest of the extended Boussinesq modelas well as the good behaviour of the numerical scheme.1.4. A fully justified extended Boussinesq system.
A major drawback arise in the left-most term of thesecond equation of (1.5) due to the positive sign in front of the elliptic fourth-order linear operator T preventing the invertibility of the factorized operator, (see [23, Section 3.1]). A remedy to this situationwould be replacing (1 + ε T [ h ] + ε T ) ∂ t v by (1 + ε T [ h ] − ε T )( ∂ t v ) + 2 ε T ( ∂ t v ) and using a BBM trickrepresented in the following approximate equation ∂ t v = − ∂ x ζ + O ( ε ), see [22, 24] for more details. Inlight of these remarks and after setting ± ε T [ h ]( vv x ) in the second equation of (1.5), one obtains the One may realize that some components in the second order operator T are of size O ( ε ). Actually they have beenkept to maintain the good properties of the operator J (1.7), otherwise these properties would have been disrupt, see [23,Remark 1]. ollowing model:(1.6) ∂ t ζ + ∂ x ( hv ) = 0 , J ( ∂ t v + εv∂ x v ) + ∂ x ζ + 245 ε ∂ x ζ + + ε ∂ x (( ∂ x v ) ) = O ( ε ) . where h = 1 + εζ and(1.7) J = 1 + ε T [ h ] − ε T . The benefit of the new formulation (1.6) is in replacing (1 + ε T [ h ] + ε T ) by a new operator J . Thisreplacement induce a fifth order derivative term on ζ , namely 245 ε ∂ x ζ , but the invertibility of theoperator J is now earned, see remark 1 in [24] for more details. An equivalent formulation of model (1.6) was fully justified recently in [23]. In fact, a unique solution of the model (1.6) exist over the time scaleof order 1 √ ε and stay close the solution of the full Euler system. Remark 1.8.
Provided that the effect of surface tension is taken into consideration, the well-posednessresult for system (1.6) can be directly deduced from [22, 24]. The existence time scale is up to order 1 /ε .The control of higher order derivatives resulting from the performance of the BBM trick is now possibledue to the presence of surface tension which seems to be very essential (see remarks in [22]). If surfacetension is neglected, which is the case in this paper, the aforementioned strategy has to be adapted.In fact, the capillary terms are replaced by a vanishing term ± ε ζ xxx where the negative sign term isused for a suitable definition of the energy space in such a way that the other term can be controlled.Consequently, the existence time scale reached is up to order 1 / √ ε , see [23], which is smaller comparedto the one obtained in the presence of surface tension.2. Reformulation of the extended Boussinesq system
The system (1.6) is much easier to solve numerically than the standard formulation (1.5). In fact, theoperator J has an appropriate structure allowing its inversion. However, the second order dispersive op-erator T [ h ] present in J depends on ζ ( t, x ) and thus on time. In fact, during the numerical computationsthe operator J has to be inverted at each time step so one can solve equation (1.6). Following [15, 27],the time dependency has to be amended in order to reduce the computational time. Using straightfor-ward asymptotic expansions, the left-most term of the second equation of (1.6) can be written underthe form: J ( ∂ t v + εv∂ x v ) = (1 + ε T [0] − ε T )( ∂ t v + εv∂ x v ) − ε ζ∂ x ( ∂ t v ) − ε ∂ x ζ∂ x ( ∂ t v ) + O ( ε ) , where T [0] w = − ∂ x w and T w = − ∂ x w . Hence the system (1.6) is equivalent to the following model:(2.1) ∂ t ζ + ∂ x (cid:0) hv (cid:1) = 0 , (cid:16) ε T [0] − ε T (cid:17)(cid:16) ∂ t v + εv∂ x v (cid:17) + ∂ x ζ + 245 ε ∂ x ζ + ε ∂ x (( ∂ x v ) ) − ε ζ∂ x ( ∂ t v ) − ε ∂ x ζ∂ x ( ∂ t v ) = O ( ε ) , where h = 1 + εζ . The left-most factorized operator of the second equation of (2.1) is now time-independent and enjoys a structure allowing its inversion, thus can be inverted once for all numericaltime steps.2.1. A one-parameter family of extended Boussinesq equations.
The eB equations are significantly im-proved in terms of linear dispersive properties due to the higher-order terms existing in these equations,see [30]. Additional improvement providing a finer characterization in high frequency regimes can bebrought by adjusting a dispersion correction parameter α . Following the lines in [10, 31, 44] and with-out affecting the accuracy of the model, we improve the frequency dispersion of problem (2.1). Thisis possible, if one adds to the second equation of (2.1) some terms of the same order as the equation The equivalent formulation is obtained by multiplying both sides of the second equation of system (1.6) by the waterheight function, h . recision and adjusts the parameter α in an appropriate way. See section 2.3 for the discussion on thechoice of the parameter α .From the second equation of (2.1), one deduce the following approximation:(2.2) ∂ t v + εv∂ x v + ∂ x ζ = O ( ε ) , and hence, for any α ∈ R ∗ + :(2.3) ∂ t v = α∂ t v − (1 − α )[ ∂ x ζ + εv∂ x v − O ( ε )] . The second equation of (2.1) can be recast after substituting ∂ t v by its approximation given in (2.3)and neglecting all terms of order O ( ε ) (cid:16) ε T [0] − ε T (cid:17)(cid:16) α∂ t v − (1 − α )[ ∂ x ζ + εv∂ x v − O ( ε )] + εv∂ x v (cid:17) + ∂ x ζ + 245 ε ∂ x ζ + ε ∂ x (( ∂ x v ) ) − ε εζ∂ x ( ∂ t v ) − ε ∂ x ζ∂ x ( ∂ t v ) = O ( ε ) . After straightforward computations,(2.4) (cid:16) εα T [0] − ε α T (cid:17)(cid:16) ∂ t v + εv∂ x v + α − α ∂ x ζ (cid:17) + 1 α ∂ x ζ + 245 ε ∂ x ζ + ε ∂ x (( ∂ x v ) ) − ε εζ∂ x ( ∂ t v ) − ε ∂ x ζ∂ x ( ∂ t v ) = O ( ε ) . One can also deduce from (2.2) that,(2.5) ∂ t v = − ∂ x ζ + O ( ε ) . Using (2.5) in the last two terms of the equation (2.4) and dropping O ( ε ) terms one gets:(2.6) ∂ t ζ + ∂ x (cid:0) hv (cid:1) = 0 , (cid:16) εα T [0] − ε α T (cid:17)(cid:16) ∂ t v + εv∂ x v + α − α ∂ x ζ (cid:17) + 1 α ∂ x ζ + 245 ε ∂ x ζ + ε ∂ x (( ∂ x v ) )+ ε ζ∂ x ζ + ε ∂ x ζ∂ x ζ = O ( ε ) . Similarly, a significant improvement of the dispersive properties has been attained in the derivation ofa three-parameter family of GN equations, see [9]. In here, we will limit ourselves to the one-parameterfamily of eB equations (2.6) for the sake of simplicity.2.2.
Reformulation of the extended Boussinesq equations (2.6) . In what follows, we derive an equivalentmodel to (2.6) at the same order of precision, i.e. O ( ε ), that prevents the calculation of high orderderivatives on ζ . In here, we call such a model eB with factorized high order derivatives . Certainly, themodel enclose high order derivatives, but we make it possible not to compute them by factoring themout by (1 + εα T [0]). The price to pay is an increase in computational cost, since one needs to solve anextra linear system . In fact, dropping all terms of order O ( ε ) from the second equation of (2.6), onecan easily check the following: (1 + εα T [0])( ∂ t v ) = − ∂ x ζ + O ( ε ) , thus we have ∂ t v = − (1 + εα T [0]) − ( ∂ x ζ ) + O ( ε ) , which may be written as ∂ t v = − ∂ x ζ + O ( ε ) . Using the approximation ∂ t v = − ∂ x ζ + O ( ε ), the terms ∂ x ζ , ∂ x ζ and ∂ x ζ become respectively:(2.7) ∂ x ζ = ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + O ( ε ) , (2.8) ∂ x ζ = ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + O ( ε ) , (2.9) ∂ x ζ = ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + O ( ε ) . Inverting the operator (1 + εα T [0]) requires the discretization of (1 + εα T [0]) A = B (i.e resolution of linear systems). eplacing ∂ x ζ , ∂ x ζ and ∂ x ζ by their expression obtained in (2.7), (2.8) and (2.9) respectively in thesecond equation of (2.6), one can write the eB equations with improved dispersion and factorized highorder derivatives as:(2.10) ∂ t ζ + ∂ x (cid:0) hv (cid:1) = 0 , (cid:16) εα T [0] − ε α T (cid:17)(cid:16) ∂ t v + εv∂ x v + α − α ∂ x ζ (cid:17) + 1 α ∂ x ζ + 245 ε ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + ε ∂ x (( ∂ x v ) ) + ε ζ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + ε ∂ x ζ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) = O ( ε ) , where h = 1 + εζ . Significant interest is behind the derivation of the eB formulation (2.10). In fact,this formulation avoid the calculation of high order derivatives which we believe will extend the rangeof applicability to high frequency regimes (see discussion in section 2.3), while remaining stable (seesection 2.4.2).2.3. Choice of the parameter α . The main comparison between any asymptotic model and the fullEuler equations is performed at the stage of linear periodic plane wave solutions. At this point, allcharacteristics of the model [41] are summed up in the dispersion relation, relating the spatial wavenumber k and the time frequency w . It comes from the earlier linearisation of the system around somerest 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 improvedeB system (2.6) at the dispersion relation level. Hence, instead of choosing the α -value that reduce theerror on the phase velocity for a specific value of k (limiting the improvement to monochromatic waves),we adjust this parameter so that both phase and group velocity are minimized over a range of valuesof k ∈ [0 , K ]. This can be done by minimizing a weighted averaged error introduced for this reason.Firstly, we discuss the dispersive properties of the dimensionalized version of model (2.6) and secondlywe discuss the dispersive properties of the dimensionalized version of model (2.10).First, lets go back to variables with dimensions. Setting ε = 1 and adding the gravity term g asneeded. The system of equations (2.6) reads(2.11) ∂ t ζ + ∂ x (cid:0) hv (cid:1) = 0 , (cid:16) α T [0] − α T (cid:17)(cid:16) ∂ t v + v∂ x v + α − α g∂ x ζ (cid:17) + 1 α g∂ x ζ + 245 g∂ x ζ + 23 ∂ x (( ∂ x v ) )+ 23 gζ∂ x ζ + g∂ x ζ∂ x ζ = 0 , where h = h + ζ (dimensionalized height). The dispersion relation corresponding to (2.11) can bederived by linearising first the equations about the rest state (˜ ζ, ˜ v ) = (0 ,
0) and then looking for thecorresponding plane wave solutions of the form ( ζ, v ) = ( ζ, v ) e i ( kx − wt ) . From the first equation of (2.11)one obtains:(2.12) − iwζ + ikh v = 0 ⇒ v = wkh ζ. The second equation of (2.11) becomes:(2.13) (cid:16) α k + α k (cid:17)(cid:16) − iwv + α − α igkζ (cid:17) + 1 α igkζ + 245 igk ζ = 0 . Substituting v in (2.13) by its expression given in (2.12) and multiplying (2.13) by − kiζ yields thedispersion relation for (2.11):(2.14) w α,eB = gh k (cid:16) α − k + ( α + 1)45 k (cid:17)(cid:16) α k + α k (cid:17) . The exact dispersion relation for the dimenisonalized full Euler system is recalled below:(2.15) w F.E = gh | k | tanh( | k | ) . or small wavenumbers, the Taylor expansions of (2.14) and (2.15) are respectively: w α,eB ≈ gh (cid:16) k − k k (cid:0)
145 + α (cid:1) + O ( k ) (cid:17) ,w F E ≈ gh (cid:16) k − k k + O ( k ) (cid:17) . The two previous Taylor expansions are equivalent provided that α = 1. (See numerical test 4.4). Thisequivalence doest not hold any more for larger wavenumbers which makes the choice of α in (2.14)essential.We describe hereafter a classical approach to find an optimal value of α when studying monochromaticwaves of wavenumber k . Let us define firstly the linear phase and group velocities associated to (2.14)as C peB ( k ) = w α,eB ( k ) | k | and C geB ( k ) = dw α,eB ( k ) dk . After determining the value of k , we choose α such that the phase velocity C peB ( k ) and the referencephase velocity C pS ( k ), coming from Stokes linear theory are uniform.This strategy is only appropriate when considering a particular wavelength. However, finding anoptimal value of α for a range of values of k requires the minimization the squared relative weightederror defined below:(2.16) (cid:90) K k (cid:16) C peB − C pS C pS + C geB − C gS C gS (cid:17) dk, over some range k ∈ [0 , K ], where C gS ( k ) is the reference group velocity associated with the Stokes lineartheory.It is clearly seen that the weighted averaged error has an absolute minimum of (0 . ≤ k ≤
1. The optimum value for α is 0.8351. Nevertheless, the error starts to grow rapidlywhen K ≥
2, showing that the model (2.11) has a limited range of applicability and thus poor dispersionproperties in intermediate and large wave numbers regime (see Figures 2 and 3). a) 0 ≤ k ≤ ≤ k ≤ Figure 2.
Phase and group velocities weighted averaged error as a function of α . a) α = 0 . α = 1 in dashes.(b) α = 0 . α = 1 in dashes. Figure 3.
Errors on linear phase velocity (red) and group velocity (blue) for the eBmodel (2.11).Now, we discuss the dispersive properties of the dimensionalized version of model (2.10). The dimen-sionalized version of equations (2.10) read:(2.17) ∂ t ζ + ∂ x (cid:0) hv (cid:1) = 0 , (cid:16) α T [0] − α T (cid:17)(cid:16) ∂ t v + v∂ x v + α − α g∂ x ζ (cid:17) + 1 α ∂ x ζ + 245 ∂ x (cid:16) (1 + α T [0]) − ( g∂ x ζ ) (cid:17) + 23 ∂ x (( ∂ x v ) ) + 23 ζ∂ x (cid:16) (1 + α T [0]) − ( g∂ x ζ ) (cid:17) + ∂ x ζ∂ x (cid:16) (1 + α T [0]) − ( g∂ x ζ ) (cid:17) = 0 , where h = h + ζ . Following the same techniques as for the derivation of (2.14), one obtains thedispersion relation for (2.17):(2.18) ˜ w α,eB = gh k (cid:16) α − k + k
45 ( α − αk ) (cid:17)(cid:16) α k + α k (cid:17) . igure 4. Phase and group velocities weighted averaged error as a function of α for0 ≤ k ≤
10. The model (2.17) is in solid line, the Green-Naghdi model in the Camssa-Holm regime (GN-CH) [3] is in dots.
Figure 5.
Errors on linear phase velocity (red) and group velocity (blue). The refer-ence from Stokes theory (black solid line), the eB model (2.17) ( α = 1 . α = 1) in dashes, the GN-CH model [3] ( α = 1 . w α,eB defined in (2.18). In Figure 4 weplot the related error in terms of α . One can see clearly that that the smallest value of the weightedaveraged error ( ≈ ≤ k ≤
10 is reached at α = 1 . O ( µ , µε )) [3] is much larger ( ≈ C peB,α opt C pS with an optimal choice α = 1 . α = 1 . ≤ k ≤
10. Larger difference exist between the group velocity (blue solidline) and the reference from Stokes theory (black solid line). This difference starts to proliferate when >
10, showing an overestimation of this property. When the modelled dispersion relation is obtainedusing the GN-CH model [3] with α = 1 . α opt in the large frequency regime.In conclusion, the Green-Naghdi equations in the Camassa-Holm regime (precise up to O ( µ , µε )) arerestrictive due to the sole inclusion of weak dispersion and nonlinearity. This usually restrain preciseapplications inside the surf area. The eB model (2.17) contains higher-order dispersive terms and anadjustable parameter α which enlarge the application scope remarkably to cover the area from deepwater (long wavelength regime) into the area of breaking waves (short wavelength regime i.e. highfrequency), see numerical test in section 4.3.2.4. Stability of the extended Boussinesq models.
The standard Boussinesq system is not stable in highfrequency regime due to the third order derivatives in ζ that are involved in the model. One can seefor instance [3], where the high frequency instabilities of an improved GN-CH model are studied. Theseinstabilities are due to the third order derivative existing in the equation. One of the advantages ofthe eB is its stability in high frequency regime due to the presence of higher order derivatives, namelyderivatives of order five in ζ . These terms seems to have a stabilizing effect. In what follows, we discussqualitatively the stability of both models (2.11) and (2.17).2.4.1. Stability of the extended Boussinesq model with high order derivatives.
Before discussing the sta-bility of the eB models in high frequency regime, we would like to mention that the choice of α in themodel (2.11) in high frequency regime does not play an important role. In fact the latter model haspoor dispersive properties in intermediate and large wave numbers regime, see discussion in section 2.3.Therefore, one has to choose α = 1.The system of equations (2.11) with α = 1 reads:(2.19) ∂ t ζ + ∂ x (cid:0) hv (cid:1) = 0 , (cid:16) T [0] − T (cid:17)(cid:16) ∂ t v + v∂ x v (cid:17) + g∂ x ζ + 245 g∂ x ζ + 23 ∂ x (( ∂ x v ) )+ 23 gζ∂ x ζ + g∂ x ζ∂ x ζ = 0 , where h = h + ζ . We carry out the linear behaviour examination of slight perturbation (˜ ζ, ˜ v ) of a steadystate solution ( ζ, v ). These perturbations are governed by the following linear system of equations:(2.20) ∂ t ˜ ζ + ( h + ζ ) ∂ x ˜ v + v∂ x ˜ ζ = 0 , (cid:16) T [0] − T (cid:17)(cid:16) ∂ t ˜ v + v∂ x ˜ v (cid:17) + g∂ x ˜ ζ + 245 g∂ x ˜ ζ + 23 gζ∂ x ˜ ζ = 0 . Going in search for plane wave solution of the above system, one can derive the corresponding dispersionrelation. These solutions are of the form (˜ ζ ; ˜ v ) = e i ( kx − wt ) ( ζ, v ) and the dispersion relation is thefollowing:(2.21) ( w − kv ) g ( h + ζ ) k = (cid:16) − k ζ + 245 k (cid:17)(cid:16) k + 145 k (cid:17) . The perturbations of the rest state ( ζ, v ) = (0 ,
0) are always stable as per the below dispersion relation:(2.22) w = gh k (cid:16) k (cid:17)(cid:16) k + 145 k (cid:17) . However, a quick functional study shows that the numerator of the right-hand side of (2.21) becomesnegative whenever ζ > k + 4530 k . Thus, as provided by (2.21), w remains real for large wavenumbers,under the condition that ζ < k + 4530 k . In the majority of the applications we have in mind, the overallsurface deformation ζ does not go beyond 2 k + 4530 k , and this condition is satisfied. Actually, as k gets arge (i.e in high frequency regime), the upper bound of ζ gets also large ( ≈ k
15 ), hence extending therange of values of the overall surface deformation ζ for which the condition is satisfied. This ensure anumerical stability in most of the situations considered for applications (see Figure 6). Remark 2.23.
Replacing ∂ x ζ by ∂ x (cid:16) (1 + T [0]) − ( ∂ x ζ ) (cid:17) in the second equation of (2.19), one gets thefollowing model:(2.24) ∂ t ζ + ∂ x (cid:0) hv (cid:1) = 0 , (cid:16) T [0] − T (cid:17)(cid:16) ∂ t v + v∂ x v (cid:17) + g∂ x ζ + 245 g∂ x (cid:16) (1 + T [0]) − ( ∂ x ζ ) (cid:17) + 23 ∂ x (( ∂ x v ) )+ 23 gζ∂ x ζ + g∂ x ζ∂ x ζ = 0 . This replacement modifies the dispersion relation (2.21) into:(2.25) ( w − kv ) g ( h + ζ ) k = (cid:16) − k ζ + k (cid:0)
21 + k (cid:1)(cid:17)(cid:16) k + 145 k (cid:17) . A similar functional study to the previous one shows that the r.h.s numerator of (2.25) is negativewhenever ζ > k + 15 k + 4510 k + 30 k . In high frequency regime, the upper bound of the overall surfacedeformation ζ is approximately close to 0 .
2, hence reducing the range of values of ζ for which thecondition is satisfied, namely − h < ζ < . h = h + ζ should remain alwayspositive). Therefore, if this condition is not satisfied the complex square root of w will generate aninstability in the model. Actually, a stability in high frequency regime is ensured if the condition − h < ζ < . ζ (see Figure 6). This is why we suggested in section 2.2 to factorize second, third and fifth orderderivatives present in the second equation of (2.19). The stability of the eB model with factorized highorder derivatives (2.17) is discussed in the next subsection.2.4.2. Stability of the extended Boussinesq model with factorized high order derivatives.
In what fol-lows, we explore the stability of the eB model with factorized high order derivatives (2.17). Whenlinearizing (2.17) about a non steady state, one gets:(2.26) ( w − kv ) g ( h + ζ ) k = (cid:16) α − k − k ζ α k ) + ( α − k
45 + 2 k α k ) (cid:17)(cid:16) α k + α k (cid:17) . The r.h.s numerator in (2.26) is positive if and only if ζ < k (cid:0) α (cid:0) ( α − k + 15 α − (cid:1) + 3 (cid:1) + 90 α − k . We recall that α >
1. Indeed, improving the dispersive properties of the model (2.17) in large frequencyregime requires an appropriate choice α = 1 . α and in high frequency regime, the r.h.s of the above inequality becomes very large, namely same orderas α ( α − k
90 . Thus, relaxing the stability condition on ζ . For the numerical tests, we choose to workwith the model (2.17) which appears to be stable. igure 6. Comparison at different times between the solutions of the models (2.19)(blue line), (2.17) (red line) and (2.24) (green line) in high frequency regime.Figure 6 shows clearly the stability of the eB models (2.19) and (2.17) while the model (2.24) seemsto be unstable in high frequency regime. In fact, when implementing in model (2.24) an initial solution hat does not satisfy the limiting stability condition discussed in remark 2.23, more precisely whenchoosing ζ = 0 . > .
2, one observes a high frequency instability. We would like to mention that wetried the same test but starting with an initial solution where the overall surface deformation is ζ = 0 . ζ . This isthe reason why, for the rest of the paper, we choose to work with the model (2.17), where high orderderivatives are factorized which, as seen above, is always stable. Of course, one could use model (2.19)which seems to be stable, but at the price of loosing the improved dispersive properties that model (2.17)enjoys, see section 2.3 and test case of section 4.3.3. Numerical methods
In what follows, we will just introduce the numerical scheme devoted to solve the eB model withfactorized high order derivatives (2.10) in order to ease the reading. A similar numerical scheme isadopted when the eB model with high order derivatives (2.6) is concerned.This class of models has two main components, hyperbolic and dispersive. This remarkable structuremakes them suitable for the implementation of a hybrid scheme splitting the two main components of theequations. This strategy has been initially introduced for Boussinesq-like and Green-Naghdi equationsin order to handle correctly wave breaking that occurs as waves approach the shore, see [2, 9, 17, 18, 27].A computation of a half-time step of the hyperbolic part is used as a sensor to evaluate the energyloss occurring during wave breaking (accurate detection of wave fronts), see [42]. Near the breakingpoints, the dynamics of the waves are described correctly using the hyperbolic part but the dispersivecomponents of the equation become very singular. In order to handle wave breaking, switching fromthe dispersive part to the hyperbolic part is indispensable. In this paper, we do not investigate breakingwaves. In fact, our work is limited to the flat topography case and we leave for future research thetreatment of breaking waves in the variable bottom configuration. However, we stick here to the splittingstrategy since it is computationally efficient, stable and cheap. The splitting scheme following the linesin [2, 27, 3] is presented in the section below.3.1.
The splitting scheme.
We recall first the eB system (2.10) under consideration:(3.1) ∂ t ζ + ∂ x (cid:0) hv (cid:1) = 0 , (cid:16) εα T [0] − ε α T (cid:17)(cid:16) ∂ t v + εv∂ x v + α − α ∂ x ζ (cid:17) + 1 α ∂ x ζ + 245 ε ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + ε ∂ x (( ∂ x v ) ) + ε ζ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + ε ∂ x ζ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) = O ( ε ) , where h = 1 + εζ . The solution operator S ( . ) related to (3.1) is decomposed at each time step ∆ t following a hybrid Strang splitting scheme: S (∆ t ) = S (∆ t/ S (∆ t ) S (∆ t/ . • S ( t ) is the solution operator related to the hyperbolic nonlinear shallow water equations, NSWE:(3.2) ∂ t ζ + ∂ x (cid:0) hv (cid:1) = 0 ,∂ t v + εv∂ x v + α − α ∂ x ζ + 1 α ∂ x ζ = 0 . The NSWE system (3.2) can be written in the following conservative form:(3.3) ∂ t ζ + ∂ x (cid:0) hv (cid:1) = 0 ,∂ t v + ∂ x (cid:16) ε v + ζ (cid:17) = 0 , where h = 1 + εζ . S ( t ) is the solution operator related to the remaining (dispersive) part of the equations.(3.4) ∂ t ζ = 0 , (cid:16) εα T [0] − ε α T (cid:17)(cid:16) ∂ t v − α ∂ x ζ (cid:17) + 1 α ∂ x ζ + 245 ε ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + ε ∂ x (( ∂ x v ) ) + ε ζ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + ε ∂ x ζ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) = 0 . The hyperbolic conservative structure of system (3.3) allows a computation of S following a finite-volume method. Whereas, a classical finite difference method is used to compute S .3.2. Finite volume scheme.
The hyperbolic system (3.3) is conveniently rewritten with conservativevariables and a flux function:(3.5) ∂ t U + ∂ x ( F ( U )) = 0 , where,(3.6) U = (cid:18) ζv (cid:19) , F ( U ) = (cid:32) hvε v + ζ (cid:33) , with h = 1 + εζ . The Jacobian matrix is given by:(3.7) A ( U ) = d ( F ( U )) = (cid:18) εv h εv (cid:19) . The homogeneous system (3.5) is strictly hyperbolic if inf x ∈ R h > ∂ t U + ∂ x ( F ( U )) = 0 , t ≥ , x ∈ R .U (0 , x ) = U ( x ) , x ∈ R . The finite volume method used to the approximation of (3.8) imposes conservation laws in a one-dimensional control volume [ x i − / , x i +1 / ] × [ t n , t n +1 ] of dimensions ∆ x = x i +1 / − x i +1 / and ∆ t = t n +1 − t n . Figure 7.
The space discretization.The approximate cell average of U on the cell m i = [ x i − / , x i +1 / ] at time t is denoted by U i and givenby: U i = 1∆ x (cid:90) m i U ( t, x ) dx . The approximate cell average of U on the cell m i = [ x i − / , x i +1 / ] at time t n is denoted by U ni andgiven by: U ni = 1∆ x (cid:90) m i U ( t n , x ) dx . Integrating (3.2) over the computational cell m i , the semi-discrete form can be represented as:(3.9) dU i ( t ) dt + 1∆ x (cid:16) F i +1 / − F i − / (cid:17) = 0 here F i ± / are the numerical fluxes defined at each cell interface as:(3.10) F i +1 / = (cid:101) F ( U i , U i +1 ) ≈ x (cid:90) m i F ( U ( t, x i +1 / )) dx. The computation of high-order accurate numerical fluxes is reached by reconstructing left and rightconstant averaged values using a fifth-order WENO scheme, before applying the numerical flux.3.2.1.
High order finite-volume scheme: WENO5-RK4.
Fixed stencil interpolation of high order accu-racy (second or higher) may not be adequate and thus oscillatory near a discontinuity. In nonlinearproblems containing discontinuities, such oscillations usually leads to numerical instabilities. Reachinghigh order accuracy in smooth regions while avoiding the spurious ossicllations around discontinuitycan be achieved by using a fifth-order accuracy WENO reconstruction for hyperbolic conservation laws,following [21, 40]. Problems with piecewise smooth solutions containing discontinuities can be solvedusing the appropriately-designed high order accurate WENO schemes. The main idea lies at the ap-proximation stage, where a nonlinear adaptive stencil is used in the reconstruction procedure following aconvex combination of all candidate stencils instead of just one to achieve the essentially non-oscillatoryproperty. A proper weight is assigned to each of the candidate stencils determining its contributionto the final approximation. In fact, a nearly zero weight is assigned to the stencils which contain thediscontinuities. Second-order schemes tend to alter the dispersive properties of the model due to dis-persive truncation errors. To prevent this in the study of dispersive waves, high order schemes areimperative [2, 9, 27]. Using the same reconstruction method proposed in [2, 3], we consider a cell m i ,and the corresponding constant averaged value U ni = ( ζ ni , v ni ). Consequently, high order reconstructedleft and right values U n, − i and U n, + i , built following the five points stencil, are introduced as follows:(3.11) U n, + i = U ni + 12 δU n, + i and U n, − i = U ni − δU n, − i , where δU n, + i and δU n, − i are defined as follows: δU n, + i = 23 ( U ni +1 − U ni ) + 13 ( U ni − U ni − ) −
110 ( − U ni − + 3 U ni − U ni +1 + U ni +2 ) −
115 ( − U ni − + 3 U ni − − U ni + U ni +1 )(3.12) δU n, − i = 23 ( U ni − U ni − ) + 13 ( U ni +1 − U ni ) −
110 ( − U ni − + 3 U ni − − U ni + U ni +1 ) −
115 ( − U ni − + 3 U ni − U ni +1 + U ni +2 ) . (3.13)The coefficients 23 , 13 , −
110 and −
115 are set for better dissipation and dispersion properties in the trun-cature error. As a result, the modified scheme becomes the following:(3.14) U n +1 i = U ni − ∆ t ∆ x (cid:16) (cid:101) F ( U n, + i , U n, − i +1 ) − (cid:101) F ( U n, + i − , U n, − i ) (cid:17) . At this point, it is important to apply a limitation procedure to prevent oscillations near discontinuities.To this end, we use the same limitation strategy as in [2, 3]. Consequently, scheme (3.14) becomes(3.15) U n +1 i = U ni − ∆ t ∆ x (cid:16) (cid:101) F ( L U n, + i , L U n, − i +1 ) − (cid:101) F ( L U n, + i − , L U n, − i ) (cid:17) . The limited high-order reconstructed values are defined as follows:(3.16) L U n, + i = U ni + 12 L + i ( U n ) and L U n, − i = U ni − L − i ( U n ) . These reconstructions must be performed on both conservative variables U ni = ( ζ ni , v ni ). L + i ( U n ) and L − i ( U n ) are defined, using the following limiter,(3.17) L ( u, v, w ) = min ( | u | , | v | , | w | ) sgn ( u ) if sgn ( u ) = sgn ( v ) , . epending on (3.17), we define the limiting process as, L + i ( U n ) = L ( δU ni , δU ni +1 , δU n, + i ) and L − i ( U n ) = L ( δU ni +1 , δU ni , δU n, − i ) , with δU ni +1 = U ni +1 − U ni and δU ni = U ni − U ni − are upstream and downstream variations, and δU n, + i and δU n, − i taken from (3.12) and (3.13). Remark 3.18.
At this stage, it is worth mentioning that steady states solution ζ = cst and v = 0 arepreserved at the discrete level.Regarding time discretization, fourth-order explicit Runge–Kutta “RK4” method is used, we describeit in what follows. Given the ODE dydt = f ( t, y ), one has, k = f ( t n , y n ) ,k = f ( t n + ∆ t , y n + ∆ t k ,k = f ( t n + ∆ t , y n + ∆ t k ,k = f ( t n + ∆ t, y n + ∆ tk ) ,y n +1 = y n + ∆ t (cid:0) k + 2 k + 2 k + k (cid:1) . (3.19)Applying (3.19) to (3.15), one gets the “WENO5-RK4” scheme.3.3. Finite difference scheme for the dispersive part.
The splitting scheme is a mix between a finite vol-ume discretization and a finite difference method. This mix induces a switching between cell-averagedvalues defined by the finite volume discretization and nodal values used by the finite difference discretiza-tion for each unknown and at each time step. Using fifth-order accuracy WENO reconstruction, one canapproximate the nodal values (i.e finite difference unknowns) ( U ni ) i =1 ,N +1 in terms of the cell-averagedvalues (i.e finite volume unknowns) ( U ni ) i =1 ,N by the following relation:(3.20) U ni = 130 U ni − − U ni − + 4760 U ni + 920 U ni +1 − U ni +2 + O (∆ x ) , ≤ i ≤ N + 1 , The global order of the scheme is preserved. In fact, the formula is precise up to order O (∆ x ) terms.Before proceeding by the computation of S , we recall first the remaining (dispersive part) of theequations, given in section 3.1.(3.21) ∂ t ζ = 0 , (cid:16) εα T [0] − ε α T (cid:17)(cid:16) ∂ t v − α ∂ x ζ (cid:17) + 1 α ∂ x ζ + 245 ε ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + ε ∂ x (( ∂ x v ) ) + ε ζ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) + ε ∂ x ζ∂ x (cid:16) (1 + εα T [0]) − ( ∂ x ζ ) (cid:17) = 0 . Using an explicit Euler in time scheme, the finite discretization of the system (3.21) using classical finitedifference methods leads to the following discrete problem:(3.22) ζ n +1 − ζ n ∆ t = 0 ,v n +1 − v n ∆ t − α D ( ζ n ) + (cid:16) I − εα D + ε α D (cid:17) − (cid:104) α D ( ζ n ) + 245 ε D (cid:16)(cid:0) − εα D (cid:1) − ( D ζ n )) (cid:17) + ε D (( D v n ) ) + ε ζ n D (cid:16)(cid:0) − εα D (cid:1) − ( D ζ n )) (cid:17) + ε D ( ζ n ) D (cid:16)(cid:0) − εα D (cid:1) − ( D ζ n )) (cid:17)(cid:105) = 0 . The matrices D , D and D are respectively the classical centered discretizations of the derivatives ∂ x , ∂ x , and ∂ x given below. The spatial derivatives are discretized using fourth-order formulae, “DF4”:( ∂ x U ) i = 112∆ x ( − U i +2 + 8 U i +1 − U i − + U i − ) , ( ∂ x U ) i = 112∆ x ( − U i +2 + 16 U i +1 − U i + 16 U i − − U i − ) , ∂ x U ) i = 16∆ x ( − U i +3 + 12 U i +2 − U i +1 + 56 U i − U i − + 12 U i − − U i − ) . A standard extension to fourth-order classical Runge-Kutta “RK4” scheme is used, and thus oneobtains the “DF4-RK4” scheme.
Remark 3.23.
At this stage, it is worth mentioning that the numerical scheme of the eB system withhigh order derivatives (2.6) is similar to the one developed for model (2.10). In fact, the hyperbolic partof the system is the same as in (3.2), but the high order derivatives involved in the second equationof the remaining (dispersive) part should be treated accordingly. More precisely, third and fifth orderderivatives are discretized using the following fourth-order formulae:( ∂ x U ) i = 18∆ x ( − U i +3 + 8 U i +2 − U i +1 + 13 U i − − U i − + U i − ) , ( ∂ x U ) i = 16∆ x ( − U i +4 + 9 U i +3 − U i +2 + 29 U i +1 − U i − + 26 U i − − U i − + U i − ) . 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 byimposing suitable relations on both cell-averaged and nodal quantities. In this paper, we only considerperiodic boundary conditions as it was already done in [4] for the study of internal waves.For the hyperbolic part, “Ghosts” cells are introduced respectively at the western and eastern bound-aries 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 outsideof the domain. In this way, we maintain centered formula at the boundaries, while keeping a regularstructure 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 bound-aries. 4.
Numerical validations
This part is devoted to the numerical validations of the model and the numerical scheme.We begin by studying the propagation of a solitary wave solution with correctors of order O ( ε )established in [23]. We compare our numerical solution with an analytic one (up to an O ( ε ) remainder)at several times and show that our numerical scheme is very efficient and accurate.Secondly, as usual when dealing with a model in oceanographical science, one has to test the abilityof the model and the numerical method for the test of a head on collision of counter propagating solitarywaves. A very good agreement is observed.An important fact to reveal is whether or not the improved model is pertinent. The third numericaltest reveals that in presence of large wave number, the choice of the parameter α and the high orderderivatives factorization are crucial.The fourth numerical test shows a very good agreement of the proposed extended Boussinesq modeland the Green-Naghdi model proposed by Duchˆene et al. [13]The fifth test is build to test the ability of the model and the numerical method to deal with irregularsolution. A very good behaviour is observed.4.1. Propagation of a solitary wave solution with correctors.
A careful examination revealed that theextended Boussinesq system (1.5) does not admit an exact solitary wave solution, see [23, Section 4]. Inorder to validate our numerical scheme we use the explicit solution with correctors of order O ( ε ) foundin [23, Section 5] that we disclose below. Such solitary waves are analytical solutions of the extendedBoussinesq system (1.5) up to O ( ε ) remainders. Therefore, this family of solutions can be used as avalidation tool for our present numerical scheme and its given by ( ζ, v ) with(4.1) ζ = ζ + ε (cid:104) ( ζ + v )( x − t ) + ( ζ − v )( x + t ) + (cid:90) t f ( s, x − t + s ) ds − (cid:90) t f ( s, x + t − s ) ds (cid:105) , and(4.2) v = v + ε (cid:104) ( ζ + v )( x − t ) − ( ζ − v )( x + t ) + (cid:90) t f ( s, x − t + s ) ds + (cid:90) t f ( s, x + t − s ) ds (cid:105) , here ( ζ , v ) is the well known explicit solution of solitary travelling wave of the standard Boussinesqsystem (1.4) given by: ζ ( t, x ) = a sech (cid:16) k ( x − ct ) (cid:17) ,v ( t, x ) = cζ ( t, x )1 + εζ ( t, x ) , where k = (cid:114) a c = (cid:114) − aε and a is an arbitrary chosen constant. The initial conditions ζ and v are both given in C ∞ ( R ) and set ζ = v = exp (cid:16) − (cid:16) πx (cid:17) (cid:17) . The function f ( t, x ) is defined by: f ( ζ , v ) = ∂ x ζ ∂ x ∂ t v + 23 ζ ∂ x ∂ t v + 145 ∂ x ∂ t v + 13 ∂ x (cid:0) v ( v ) xx − ( v ) x (cid:1) . In this test, we investigate the left to right propagation of a solitary wave initially centred at x = 20,of relative amplitude a = 0 .
2. The computational domain length is L = 100 and discretized with 1600cells. The solitary wave is initially far from boundaries, thus the periodic boundary conditions do notaffect the computation. We compare the water surface profile of our numerical solution provided bythe model (2.10) with ε = 0 .
01 and α = 1, with the analytical one given by (4.1)-(4.2) at several timesusing the fifth order discretization “WENO5-DF4-RK4”. An excellent agreement between numerical andanalytical solutions is observed in Figure 8. The amplitude and shape of the computed solitary waveare accurately preserved during the propagation, indicating an accurate discretization of the governingequations in both space and time. Figure 8.
Propagation of a solitary wave: water surface profiles at t = 0 , , , N , over a duration T = 1. Westart with N = 400 number of cells and successively multiply the number of cells by two. The relativeerrors E L ( ζ ) and E L ( v ) on the water surface deformation and the averaged velocity are computed at t = 1, using the discrete L norm (cid:107) . (cid:107) :(4.3) E L ( ζ ) = (cid:107) ζ num − ζ sol (cid:107) (cid:107) ζ sol (cid:107) ; E L ( v ) = (cid:107) v num − v sol (cid:107) (cid:107) v sol (cid:107) , where ( ζ num , v num ) are the numerical solutions and ( ζ sol , v sol ) are the analytical ones coming from (4.1)-(4.2). E L ( ζ ) E L ( v )400 3 . × − . × −
800 9 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − Table 1.
Propagation of a solitary wave: relative L -error table for the conservativevariables. Figure 9.
Propagation of a solitary wave: L -error on the water surface deformationand the averaged velocity for a = 0 . E L ( ζ ) and E L ( v ) are plotted against δ x = LN inlog scales, for the considered relative amplitude a = 0 .
2. Very accurate results are obtained, indicatingthat the employed numerical method is capable of computing in a stable way the propagation of the asolitary 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 . Remark 4.4.
We believe that the main reason for not obtaining the predicted order in space discretizationis due to the fact that the analytic solution given in (4.1)-(4.2) satisfies the model (2.10) up to O ( ε )remainders, that is to say it is an approximate solution. A remedy for this situation could be throughcalculating explicitly the O ( ε ) remainders by injecting the solitary wave (4.1)-(4.2) into system (2.10).These residuals are therefore given by explicit formulae. This technique has been used in [27]. In thispaper, we do not try to give some optimal convergence result and calculating the residuals explicitly isleft to future work. Let us also mention that the global (time and space) order of our scheme may belimited by the order of the splitting method used here, which is of order two as already discussed byBonneton et al. in [2].4.2. Head on collision of counter propagating solitary waves.
A standard nonlinear test case for nu-merical methods is the interaction of solitary wave. In this numerical test, we study an importantphenomenon in the study of nonlinear dispesive waves, the head on collision of two counter propagatingwaves with different amplitudes. We used solitary wave solutions with correctors of order O ( ε ) forthe eB system established in [23] and defined in (4.1)-(4.2). We consider two solitary waves centred at x = −
50 and x = 50 at t = 0 on a spatial domain L = 200, see Figure 10. The solitary wave centred at x = −
50 travels to the right with a speed c s, = 1 . a = 0 . x = 50 travels to the left with a speed c s, = 1 . a = 0 .
2. The omain is discretized using 1200 cells and periodic boundary conditions are imposed. The numericalsolutions are computed using model (2.10) with ε = 0 . α = 1. The collision of the two waves startsat about t = 43, see Figure 10. After the interaction, each wave continue moving in its own directionand turn up to be unaffected by the collision, see Figures 11 and 12. A proper description of the dis-tinctive nature of nonlinear interactions is illustrated when zooming at the oscillating dispersive tailsof very small amplitude at t = 70 in Figure 13. The high precision of our numerical scheme is verifiedafter accurately capturing this phenomenon and inducing similar observations to earlier works [3, 16, 28]where the head-on collision is carried out. Figure 10.
Head on collisions: surface wave shape at t = 0 and t = 43. igure 11. Head on collisions: surface wave shape at t = 46 ,
49 and 53. igure 12. Head on collisions: surface wave shape at t = 55 ,
58 and 60. igure 13. Head on collisions: surface wave shape at t = 70.4.3. Breaking of a regular heap of water with a large wave number.
In this numerical test, we highlightthe importance of factorizing high order derivatives present in the improved eB model (2.10) togetherwith an appropriate choice of the parameter α in improving the frequency dispersion in high frequencyregimes. To this end, we consider a sufficiently regular heap of water with a large wavenumber repre-sented by the initial data: ζ (0 , x ) = 0 . e − x , v (0 , x ) = 0 , (dashed lines) with a domain of computation x ∈ ( − ,
2) discretized with 512 cells and under periodicboundary conditions. The nonlinearity parameter is set as follow: ε = 0 . γ = 0 and δ = 1. In [14], the authorscompute the solution of the obtained Green-Naghdi systems by using spectral methods [43] for spa-tial discretization and the Matlab solver ode45 for temporal discretization. In this numerical test, wecompare th numerical solutions of the proposed model eB, with the solutions corresponding to the “im-proved” Green-Naghi model. This model enjoys the same dispersion relation as the one of the full Eulersystem . In [3], the authors improve the frequency dispersion of the original Green–Naghdi model in he Camassa–Holm regime (GN-CH) describing a two-layer flow thanks to a new reformulation with asuitably chosen parameter α keeping the same order of precision as the original one, namely O ( µ , µε ).One can easily recover the one-layer configuration by setting γ = 0 and δ = 1. Using similar numericalmethods to the one adopted in this paper, in [3] the authors solve numerically the GN-CH model. (a) Comparison of the numerical solutions of model (2.10) (blue) with the “improved” GN model (green) andthe GN-CH model (red).(b) Comparison of the numerical solutions of model (2.6) (blue) with the “improved” GN model (green) andthe GN-CH model (red). Figure 14.
Comparison of the numerical solutions of model (2.10) and model (2.6)with the “improved” GN model and the GN-CH model [3] at t = 3.Figure 14(a) shows when α is chosen appropriately as discussed in section 2.3, namely α opt = 1 . t = 3 using model (2.10) behave similarlyto the one computed with the “improved” Green-Naghdi model. In contrary, when choosing α opt = 1or when using the model (2.6) to compute the solution (Figure 14 (b)), the behaviour is different thanthe “improved” Green-Naghdi model. The GN-CH model has an improved frequency dispersion due tothe careful choice of the parameter α . Nevertheless, the numerical solution computed using the GN-CH model is far from the numerical solution of the “improved” Green-Naghdi model. In fact, the eBmodel (2.10) is precise up to O ( ε ) order and thus contains high-order dispersive terms that are neglectedin the GN-CH model. This explains the fair agreement in high frequency regime between the numerical olutions of the “improved” Green-Naghdi model and the eB model (2.10) rather than the GN-CH model.One can see that the choice of an optimal value of α when using the model (2.6) has no beneficial effectdue to the high frequency regime setting. This numerical test confirms that the model (2.6) has a rangeof applicability limited to k ≤ full Euler system when α is appropriately chosen (see discussionin section 2.3). One can deduce the importance of using the model (2.10) containing factorized highorder derivatives along with an optimal choice of α in the implementation of smooth data with highfrequencies.4.4. Breaking of a regular heap of water with a small wave number.
In this numerical test, we considerthe breaking of a sufficiently regular heap of water with a small wavenumber whose initial data is: ζ (0 , x ) = 0 . e − . x , v (0 , x ) = 0 , (dashed lines) within a domain of computation x ∈ ( − ,
2) discretized with 512 cells and under periodicboundary conditions. The nonlinearity parameter is set as follow: ε = 0 . t = 3. The parameter α is fixed as 1. In fact, both eB models (2.6) and (2.10) havean equivalent dispersion relation to the one of the full Euler system for small wave numbers whenever α = 1. Indeed, Figure 15 shows a fairly good agreement between the solutions of the eB models (2.6)(yellow line) and (2.10) (blue line), and the solution of the GN-CH model [3] (red line) and the one ofthe “improved” Green-Naghdi model (green line). This confirms the fact that all the aforementionedmodels behave similarly in small wavenumbers regime, hence one can deduce that these models enjoyssimilar dispersive properties as the one of the full Euler system in small wavenumbers regime. Figure 15.
Comparison of the numerical solutions of model (2.10) (blue line) andmodel (2.6) (yellow line) and the GN-CH model [3] (red line) with the “improved”Green-Naghdi model (green line) at t = 3.4.5. Dam-break problem.
Dealing with non-regular solutions needs a special treatment at the numericalscheme level. Earlier works [2, 9, 27] have shown that the use of high-order schemes in dispersivewaves 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 [28, 37] when considering discontinuousinitial data. In this numerical test, we implement a dam break problem in order to investigate theperformance of our numerical scheme in handling non-regular solutions. We consider the followinginitial data:(4.5) ζ (0 , x ) = a [1 + tanh(250 − | x | )] , v (0 , x ) = 0 , ith a = 0 . m defined on the computational domain x ∈ ( − , α isnot important, thus we choose α = 1. Figure 16.
Dam break: wave shape at t = 0 s, s and 30 s . igure 17. Dam break: wave shape at t = 30 s and t = 65 s .The dam break wave shape is shown at different times in Figures 16 and 17. The initial data breaksand generates dispersive shock waves. A close-up on the profiles at t = 20 s , t = 30 s and t = 65 s shows two dispersive shock waves counter-propagating on both sides of the “dam”, and two rarefactionwaves moving in the direction of the center. This shows that our high-order numerical scheme was ableto capture accurately this phenomenon. Dispersive shock waves in a large class of dispersive shallow ater models were carried out in several works [3, 16, 28, 37, 36] and show a good agreement with ournumerical simulation. References [1]
J. L. Bona, T. Colin, and D. Lannes , Long wave approximations for water waves , Arch. Ration. Mech. Anal., 178(2005), pp. 373–410.[2]
P. Bonneton, F. Chazel, D. Lannes, F. Marche, and M. Tissier , A splitting approach for the fully nonlinear andweakly dispersive Green-Naghdi model , J. Comput. Phys., 230 (2011), pp. 1479–1498.[3]
C. Bourdarias, S. Gerbi, and R. Lteif , A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves , Comput. & Fluids, 156 (2017), pp. 283–304.[4] ,
A numerical scheme for the propagation of internal waves in an oceanographic model , in Finite volumesfor complex applications VIII—hyperbolic, elliptic and parabolic problems, vol. 200 of Springer Proc. Math. Stat.,Springer, Cham, 2017, pp. 101–108.[5]
J. Boussinesq , Th´eorie de l’intumescence liquide appel´ee onde solitaire ou de translation se propageant dans uncanal rectangulaire , C.R. Acad. Sci. Paris S´er. A-B, 72 (1871), pp. 755–759.[6]
J. Boussinesq , Th´eorie des ondes et des remous qui se propagent le long d’un canal rectangulaire horizontal, encommuniquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond. , J. Math.Pures Appl., 17 (1872), pp. 55–108.[7]
C. Burtea , New long time existence results for a class of Boussinesq-type systems , Journal de Math´ematiques Pureset Appliqu´ees, 106 (2016), pp. 203 – 236.[8]
F. Chazel , Influence of bottom topography on long water waves , M2AN Math. Model. Numer. Anal., 41 (2007),pp. 771–799.[9]
F. Chazel, D. Lannes, and F. Marche , Numerical simulation of strongly nonlinear and dispersive waves using aGreen-Naghdi model , J. Sci. Comput., 48 (2011), pp. 105–116.[10]
R. Cienfuegos, E. Barth´elemy, and P. Bonneton , A fourth-order compact finite volume scheme for fully nonlinearand weakly dispersive Boussinesq-type equations. II. Boundary conditions and validation , Internat. J. Numer. MethodsFluids, 53 (2007), pp. 1423–1455.[11]
W. Craig and C. Sulem , Numerical simulation of gravity waves , J. Comput. Phys., 108 (1993), pp. 73–83.[12]
W. Craig, C. Sulem, and P.-L. Sulem , Nonlinear modulation of gravity waves: a rigorous approach , Nonlinearity,5 (1992), pp. 497–522.[13]
V. Duchˆene, S. Israwi, and R. Talhouk , Shallow water asymptotic models for the propagation of internal waves ,Discrete Contin. Dyn. Syst. Ser. S, 7 (2014), pp. 239–269.[14]
V. Duchˆene, S. Israwi, and R. Talhouk , A new class of two-layer Green-Naghdi systems with improved frequencydispersion , Stud. Appl. Math., 137 (2016), pp. 356–415.[15]
A. Duran and F. Marche , Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations , Com-mun. Comput. Phys., 17 (2015), pp. 721–760.[16]
G. A. El, R. H. J. Grimshaw, and N. F. Smyth , Unsteady undular bores in fully nonlinear shallow-water theory ,Phys. Fluids, 18 (2006), pp. 027104, 17.[17]
K. S. Erduran, S. Ilic, and V. Kutija , Hybrid finite-volume finite-difference scheme for the solution of Boussinesqequations , Internat. J. Numer. Methods Fluids, 49 (2005), pp. 1213–1232.[18]
S. S. Frazao and Y. Zech , Undular bores and secondary waves -experiments and hybrid finite-volume modelling ,Journal of Hydraulic Research, 40 (2002), pp. 33–43.[19]
M. F. Gobbi, J. T. Kirby, and G. Wei , A fully nonlinear Boussinesq model for surface waves. II. Extension to O ( kh ) , J. Fluid Mech., 405 (2000), pp. 181–210.[20] A. E. Green and P. M. Naghdi , A derivation of equations for wave propagation in water of variable depth , J. FluidMech., 78 (1976), pp. 237–246.[21]
G.-S. Jiang and C.-W. Shu , Efficient implementation of weighted ENO schemes , J. Comput. Phys., 126 (1996),pp. 202–228.[22]
B. Khorbatly and S. Israwi , Full justification for the extended Green-Naghdi system for an uneven bottom withsurface tension . Preprint, hal-02994586, (2020).[23]
B. Khorbatly, R. Lteif, S. Israwi, and S. Gerbi , Mathematical modeling and numerical analysis for the higherorder Boussinesq system , Preprint, arXiv:2102.08045, (2021).[24]
B. Khorbatly, I. Zaiter, and S. Israwi , Derivation and well-posedness of the extended Green-Naghdi equations forflat bottoms with surface tension , Journal of Mathematical Physics, 59 (2018), p. 071501.[25]
J. Kirby , Nonlinear, dispersive long waves in water of variable depth , Advances in Fluid Mechanics, 10 (1996), p. 75.[26]
D. Lannes , The water waves problem , vol. 188 of Mathematical Surveys and Monographs, American MathematicalSociety, Providence, RI, 2013. Mathematical analysis and asymptotics.[27]
D. Lannes and F. Marche , A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient2D simulations , J. Comput. Phys., 282 (2015), pp. 238–268.[28]
O. Le M´etayer, S. Gavrilyuk, and S. Hank , A numerical scheme for the Green-Naghdi model , J. Comput. Phys.,229 (2010), pp. 2034–2045.[29]
F. Linares, D. Pilod, and J. Saut , Well-posedness of strongly dispersive two-dimensional surface wave Boussinesqsystems , SIAM J. Math. Anal., 44 (2012), pp. 4195–4221.[30]
P. A. Madsen, H. B. Bingham, and H. A. Sch¨affer , Boussinesq-type formulations for fully nonlinear and extremelydispersive water waves: derivation and analysis , Proceedings of the Royal Society of London. Series A: Mathematical,Physical and Engineering Sciences, 459 (2003), pp. 1075–1104. P. A. Madsen, R. Murray, and O. R. Sørensen , A new form of the Boussinesq equations with improved lineardispersion characteristics , Coastal Engineering, 15 (1991), pp. 371 – 388.[32]
P. A. Madsen and H. A. Sch¨affer , Higher-order Boussinesq-type equations for surface gravity waves: derivationand analysis , R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci., 356 (1998), pp. 3123–3184.[33]
Y. Matsuno , Hamiltonian formulation of the extended Green–Naghdi equations , Physica D: Nonlinear Phenomena,301-302 (2015), pp. 1 – 7.[34] ,
Hamiltonian structure for two-dimensional extended Green–Naghdi equations , Proc. R. Soc. A., 472 (2016).[35]
M. Ming, J. C. Saut, and P. Zhang , Long-time existence of solutions to Boussinesq systems , SIAM J. Math. Anal.,44 (2012), pp. 4078–4100.[36]
D. Mitsotakis, D. Dutykh, and J. Carter , On the nonlinear dynamics of the traveling-wave solutions of the serresystem , Wave Motion, 70 (2017), pp. 166 – 182. Recent Advances on Wave Motion in Fluids and Solids.[37]
D. Mitsotakis, B. Ilan, and D. Dutykh , On the Galerkin/finite-element method for the Serre equations , J. Sci.Comput., 61 (2014), pp. 166–195.[38]
J. Saut, C. C. Wang, and L. Xu , The cauchy problem on large time for surface-waves-type Boussinesq systems ii ,SIAM J. Math. Anal., 49 (2017), pp. 2321–2386.[39]
J.-C. Saut and L. Xu , The Cauchy problem on large time for surface waves Boussinesq systems , J. Math. PuresAppl. (9), 97 (2012), pp. 635–662.[40]
C.-W. Shu , Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservationlaws , in Advanced numerical approximation of nonlinear hyperbolic equations (Cetraro, 1997), vol. 1697 of LectureNotes in Math., Springer, Berlin, 1998, pp. 325–432.[41]
J. Stoker , Water Waves: The mathematical theory with applications , Interscience, New York, (1957).[42]
M. Tissier, P. Bonneton, F. Marche, F. Chazel, and D. Lannes , A new approach to handle wave breaking infully non-linear Boussinesq models , Coastal Engineering, 67 (2012), pp. 54 – 66.[43]
L. N. Trefethen , Spectral methods in MATLAB , vol. 10 of Software, Environments, and Tools, Society for Industrialand Applied Mathematics (SIAM), Philadelphia, PA, 2000.[44]
J. M. Witting , A unified model for the evolution nonlinear water waves , Journal of Computational Physics, 56(1984), pp. 203 – 236.[45]
V. E. Zakharov , Stability of periodic waves of finite amplitude on the surface of a deep fluid , Journal of AppliedMechanics and Technical Physics, 9 (1968), pp. 190–194.
Laboratoire de Math´ematiques UMR 5127 CNRS & Universit´e de Savoie Mont Blanc, Campus scientifique,73376 Le Bourget du Lac Cedex, France
Email address : [email protected] Lebanese American University (LAU), Graduate Studies and Research (GSR) office, School of Arts andSciences, Computer Science and Mathematics Department, Beirut, Lebanon