Mathematical modeling and numerical analysis for the higher order Boussinesq system

. This study deals with higher-ordered asymptotic equations for the water-waves problem. We considered the higher-order/extended Boussinesq equations over a ﬂat bottom topography in the well-known long wave regime. Providing an existence and uniqueness of solution on a relevant time scale of order 1 / √ ε and showing that the solution’s behavior is close to the solution of the water waves equations with a better precision corresponding to initial data, the asymptotic model is well-posed in the sense of Hadamard. Then we compared several water waves solitary solutions with respect to the numerical solution of our model. At last, we solve explicitly this model and validate the results numerically.

In this paper, we investigate the one-dimensional flow of the free surface of a homogeneous, immiscible fluid moving above a flat topography z = −h 0 . The horizontal and vertical variables are denoted respectively by x ∈ R and z ∈ R and t ≥ 0 stands for the time variable. The free surface is parametrized by the graph of the function ζ(t, x) denoting the variation with respect to its rest state z = 0 (see Figure 1). The fluid occupies the strictly connected (ζ(t, x) + h 0 > 0) domain Ω t at time t ≥ 0 denoted by: 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 is divergence and curl free. We denote by (ρ, V ) the constant density and velocity field of the fluid. The first boundary condition at the free surface expresses a balance of forces. Kinematic boundary conditions are considered assuming that both the surface and bottom are impenetrable, that is no particle of fluid can cross. The set of equations describing the flow is now complete and is commonly known as the full Euler equations: (1.1) in (x, z) ∈ Ω t , t ≥ 0, P | z=ζ(t,x) = 0 for t ≥ 0, x ∈ R, ∂ t ζ − 1 + |∂ x ζ| 2 n ζ · V | z=ζ(t,x) = 0 for t ≥ 0, x ∈ R, −V · − → e z = 0 at z = −h 0 , t ≥ 0, , where n ζ = 1 1 + |∂ x ζ| 2 (−∂ x ζ, 1) T denotes the upward normal vector to the free surface.
The theoretical study of the above system of equations is extremely difficult due to its large number of unknowns and its time-dependent moving domain Ω t . In fact, we have a free boundary problem, in other words the domain is itself one of the unknowns. Using the assumption of irrotational velocity field, one can express the latter as the gradient of a 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's equation. Although the system now is simpler, a free boundary problem still exists. To get over this obstacle, Craig and Sulem [10,11] had an interesting idea following Zakharov work [42], consisting of a reformulation 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 trace of the velocity potential at the free surface, ψ(t, x) = ϕ(t, x, ζ(t, x)) = ϕ |z=ζ , the Dirichlet-Neumann operator is introduced where ϕ is defined uniquely from (ζ, ψ) as a solution of the following Laplace problem (see [28] for a complete and accurate analysis): with ∂ n = n.∇ x,z the normal derivative in the direction of the concerned vector n. Thus, the evolution of only the two variables (ζ, ψ) located at the free surface characterize the flow. This system is known by the Zakharov/Craig-Sulem formulation of the water-waves equations giving : The above system of equations has a particularly rich structure, and depending on the physical properties of the flow, it is possible to obtain solutions to (1.2) with different qualitative properties. Nonlinear effects, for example, become more important as wave amplitude increases. Although Zakharov's reformulation resulted in a reduced system of equations, the description of these solutions from a qualitative and quantitative point of view remains very complex. A remedy for this situation requires the construction of simplified asymptotic models whose solutions are approximate solutions of the full system. These approximate models allow to describe in a fairly precise way the behavior of the complete system in a specific physical regime. This requires a rescaling 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. The order of magnitude of these parameters makes it possible to identify the considered physical regime. We start by introducing respectively the commonly known nonlinear and shallowness parameters: where 0 ≤ ε ≤ 1 is often called nonlinearity parameter, while 0 ≤ µ ≤ 1 is called the shallowness parameter. In this manner, the dimensionless formulation of (1.2) reads: Let us now identify the asymptotic geophysical shallow-water (µ 1) category (or sub-regime) associated with our work. An additional assumption is made on the nonlinearity parameter, from which a diverse set of asymptotic models can be derived. More precisely, it is possible to deduce from (1.3) a (much simpler) asymptotic model that is more amenable to numerical simulations and have more transparent properties. For instance, taking ε ∼ µ into account, the flow under consideration is said to be in a small amplitude regime.
1.2. Shallow-water, flat bottom, small amplitude variations (µ 1, ε ∼ µ). In this paper, we restrict our work on the well-known long waves regime with a flat topography for which the "original" or "standard" Boussinesq system can be derived. Defining the depth-averaged horizontal velocity by : under the extra assumption ε ∼ µ, we can neglect the terms which are of order O(µ 2 ) in the Green-Naghdi equations (we refer to [17,16] for formal derivation and to [20,19,23,13] for well-posedness); then the standard Boussinesq equations reads: Many strategies exist to study the water-wave problem especially by deriving equivalent models with better mathematical structure such as well-posedness, conservation of energy, solitary waves, or physical properties (see for instance [2,29,3,7,32,35,34,5,6,28,36,37,38]). It is worth noticing that the well posed results for such model exist on a time scale of order 1/ √ ε (methods based on dispersive estimate in [42]) and 1/ε (energy estimate method in [28] ). A better precision is obtained when the O(µ 2 ) terms are kept in the equations: only O(µ 3 ) terms are dropped. Following the work in a series of papers on the extended Green-Naghdi equations [30,31,25,24], one may write the extended Boussinesq equations by incorporating higher order dispersive effects as follows: where h = 1 + εζ is the non-dimensionalised height of the fluid and we denote the three operators : 1.3. Presentation of the results. As mentioned before, we will first derive an extended Boussinesq equations in the same way as the derivation of the extended Green-Naghdi equations: we will keep every terms up to the third order in ε. This is done in the next section, section 2. Section 3 is devoted to the full justification of the extended Boussinesq system. We will firstly, in subsection 3.2, write the extended Boussinesq system in a quasilinear form. The linear analysis, performed in subsection 3.3 will permit by the energy estimate method to state, in the subsection 3.4, the main results of well-posedness, stability and convergence of the proposed extended Boussinesq system. As for usual Green-Naghdi and Boussinesq model, we are interested in the construction of a solution as a solitary wave. We will prove in section 4 that the profile of this solitary wave is a solution of a 3rd order non linear ordinary differential equation, ODE. Thus, it seems impossible to find an explicit form of this profile. Therefore, we will compute, using Matlab ODE solver ode45, an approximate profile. We will compare the obtained solutions with the solutions of water-waves equations and find that this solution is a better approximation than the solution of the original Green-Naghdi equation.
Lastly, instead of finding an analytical exact solitary wave, we will find an explicit solution with correctors in section 5.
1.4. Comments on the results. In this section we try to highlight the potential need of higher-ordered models and their benefits over the classical asymptotic ones. Despite having a more complicated structure than classical models, higher ordered models may still be considered simpler than the original full Euler system (1.1). In fact, as opposed to the full Euler system, these high order models enjoy a reduced structure in terms of number of equations, unknown numbers and dimension space which make them more suitable for theoretical and numerical study. Moreover, higher order approximations may have similar well-posedness results as classcial ones on relevant time scales due to standard mathematical tools. Based on section 3 and previous works [25,24] this can be concluded at least in the one-dimensional case. However, the advantage is obvious in terms of controlling the convergence precision of the approximation error with respect to Euler equations (see in particular Theorem 3 of section 3). On the other hand, while the solitary wave profile cannot be derived explicitly for higher order approximations, the numerical solution fits the corresponding one of the original Euler system much better than classical models (as shown by figure 2). The numerical solution computation requires simple discretization of a third-order nonlinear ODE using Matlab ode45 solver. Furthermore, it is noteworthy that by removing the ε 2 extended-Boussiseq ODE terms, the Green-Naghdi's ODE can be recovered.
1.5. Notation. We denote by C(λ 1 , λ 2 , ...) a constant depending on the parameters λ 1 , λ 2 , ... and whose dependence on the λ j is always assumed to be nondecreasing. The notation a b means that a ≤ Cb, for some non-negative constant C whose exact expression is of no importance (in particular, it is independent of the small parameters involved ).
We denote the L 2 norm | · | L 2 simply by | · | 2 . The inner product of any functions f 1 and f 2 in the Hilbert For any real constant s, x ) s/2 . For any functions u = u(t, X) and v(t, X) defined on [0, T ) × R d with T > 0, we denote the inner product, the L p -norm and especially the L 2 -norm, as well as the Sobolev norm, with respect to the spatial variable, , g and f g belonging to the domain of T .

The higher-order/extended Boussinesq equations
When the surface elevation is of small amplitude, that is, when an assumption is made on the nonlinearity parameter, the extended Green-Naghdi equations [30,31,25,24] can be greatly simplified. Based on this, the extended Boussinesq with ε ∼ µ reads for one-dimensional small amplitude surfaces: where the right-hand side is of order ε 3 , and we see the dependence on ε 2 in the left-hand side. Here h = 1 + εζ and we denote by Remark 1. Some of the components in the second equation's left-most term are of the size O(ε 3 ). They were kept to preserve the operator's = h + εT [h] − ε 2 ∂ 4 x good properties; otherwise, these properties would have been disrupted (see section 3.1).
2.1. The modified system. First of all, let us factorize all higher order derivatives (third and fifth) in the left-most term of the above system (2.1). In fact, we only have to factorize third-order derivatives and this is possible by setting ±ε 2 T [h](vv x ) in the second equation. An inconvenient feature appears in this left-most term due to the positive sign in front of the elliptic forth-order linear operator T which ravel the way towards well-posedness using energy estimate method. This obviously affect the invertibility of the factorized operator as we will see in section 3.1. For this reason we proceed as in [25,24] by using a BBM trick represented in the following approximate equation At this stage, it is noteworthy that from [25,24] one may conclude directly the well-posedness results for such system but when the effect of surface tension is taken into consideration, the existence time scale is up to order 1/ε. This presence of the surface tension was essential for controlling higher order derivatives yielding from the BBM trick (see remarks in [24]). In our case, the surface tension is neglected and thus we have to do proceed differently. The idea is to replace the capillary terms by a vanishing term ±ε 2 ζ xxx which will play a similar role. The term with a negative sign is used for a convenient definition of the energy space (see Definition 1) in such a way that the other term can be controlled. As a consequence, the existence time will get smaller with respect to the case of surface tension presence, i.e. the time scale reached is up to order 1/ √ ε. In view of the above notes (we refer to remarks 4 and 3 for more details), the modified system reads: where U = (ζ, v), h(t, x) = 1 + εζ(t, x) and denote by

Remark 2.
An equivalent formulation of system (2.2) has been numerically studied recently in [15]. This formulation is obtained by dividing the second equation of system (2.2) by the water height function, h and removing time dependency from the left-most factorized operator while keeping the same precision of the model. During the numerical computations this operator has to be inverted at each time step so one can solve system (2.2). The time dependency has to be amended in order to reduce the computational time.
We state here that the solution of (1.3) is also a solution to the extended Boussinesq system (2.2) up to terms of order O(ε 3 ).
. Proof. Equation one of (2.2) exactly coincides with that of (1.3). It remains to check that the second equation is satisfied up to a remainder R such that (2.4) holds. For this sake, we need an asymptotic expansion of ψ in terms of v which can be deduced from the work done in [25] as follows : Now we proceed iusing the same arguments as the ones used in Lemmas 5.4 and 5.11 in [28] to give some control on R ε 3 as follows : . Then we take the derivative of the second equation of (1.3) and substitute G[εζ]ψ and ψ by −ε∂ x (hv) and (2.5) respectively. Therefore, taking advantage of the estimates (2.6) provides the control of all terms of order ε 3 as in (2.4) with N large enough (mainly greater than 8).
3. Full justification of the extended Boussinesq system (µ 3 < µ 2 < µ 1, ε ∼ µ) The two main issues regarding the validity of an asymptotic model are the following: • Are the Cauchy problems for both the full Euler system and the asymptotic model well-posed for a given class of initial data, and over the relevant time scale ? • Can the water waves solutions be compared to the solutions of the full Euler system when corresponding initial data are close? If yes, can we estimate how close they are? When an asymptotic model answer these two questions, it is said to be fully justified. In the sequel, after the linear analysis of our model, we refer to section 3.4 to state the answers of these questions. Existence and uniqueness of our solution on a time scale 1/ √ ε is given by Theorem 1, while a stability property is provided by Theorem 2. Finally, the convergence Theorem 3 is stated and therefore the full justification of our model is proved.
Let us firstly state some preliminary results in the section below.
3.1. Properties of the two operators and −1 . Assume the nonzero-depth condition that underline the fact that the height of the liquid is always confined, i.e. : Under the above condition, let us introduce the operator , where much of the modifications in the previous section hinges on it, such as: The following lemma states the invertibility results of the operator on well chosen functional spaces.
Lemma 1. Suppose that the depth condition (3.1) is satisfied by the scalar function ζ(t, ·) ∈ L ∞ (R). Then, the operator is well defined, one-to-one and onto .
Proof. We refer to the recent works of two of the authors, [25, Lemma 1] and [24, Lemma 1], for the proof of this lemma.
Some functional properties on the operator −1 are given by the Lemma below.
Proof. We refer to the recent works of two of the authors, [25, Lemma 2] and [24, Lemma 2], for the proof of this lemma.

Quasilinear form.
In order to rewrite the extended Boussinesq system in a condensed form and for the sake of clarity, let us introduce an elliptic forth-order operator T [h] as follows: The first equation of the system (2.2) can be written as follows: Then we apply −1 to both sides of the second equation of the system (2.2), to get: Hence the higher order Boussinesq system can be written under the form: where the operator A is denied by: 3.3. Linear analysis. We consider the following linearized system around a reference state U = (ζ, v) T : The energy estimate method needs to define a suitable energy space for the problem we are considering here. This will permit the convergence of an iterative scheme to construct a solution to the extended Boussinesq system (2.2) for the initial value problem (3.6).
Definition 1 (Energy space). For all s ≥ 0 and T > 0, we denote by X s the vector space H s+2 (R)×H s+2 (R) endowed with the norm: ; X s ) endowed with its canonical norm. Remark 3. It is worth noticing that in the presence of surface tension the second term of the energy norm, |ζ x | 2 H s , is controlled by ε in front of it and this is sufficiently enough to give an existence time scale of order 1/ε. In fact, the second term here in | · | X s is due to the consideration of the vanishing term that is important for Definition 1 itself and for controlling higher order terms (see Proposition 2). Now we remark that a good suggestion of a pseudo-symmetrizer for A[U ] requires firstly the introduction of a forth-order linear operator J[h] as follows: where h = 1 + εζ . Thus a pseudo-symmetrizer for A[U ] is given by: . This is clearly necessary for controlling A 2 + A 3 (see Proposition 2).
Also, a natural energy for the initial value problem (3.6) is suggested to be as follows: Lemma 3 (Equivalency of E s (U ) and the X s -norm). Let s ≥ 0 and suppose that ζ ∈ L ∞ (R) satisfies consition (3.1). Then norm | · | X s and the natural energy E s (U ) are uniformly equivalent with respect to ε ∈ (0, 1) such that: Proof. We refer to the recent work of two of the authors [25, Lemma 3] for the proof of this important property.
The well-posedness and a derivation of a first energy estimate for the linear system is given in the following proposition. it holds that: for some λ T depending only on h −1 min , sup 0≤ √ εt≤T E s (U (t)) and sup 0≤ Proof. For the proof of the existence and uniqueness of the solution, we refer to the proof found in [20, Appendix A] which can be directly adapted to the problem we are considering here. Thereafter, we will focus our attention on the proof of the energy estimate (3.9). First of all, fix λ ∈ R. The proof of the energy estimate is centered on bounding from above by zero the expression e √ ελt ∂ t (e − √ ελt E s (U ) 2 ). For this sake, we use the fact that and J[h] are symmetric to evaluate the expression under the form: Now it remains to control the r.h.s components of the above equation. To do so, we firstly recall the commutator estimate we shall use due to Kato-Ponce [22] and recently improved by Lannes [27]: in particular, for any s > 3/2, and q ∈ H s (R), p ∈ H s−1 (R), one has: Also we shall use intensively the classical product estimate (see [1,27,22]): in particular, for any p, q ∈ H s (R 2 ), s > 3/2, one has: • Estimation of (SA[U ]Λ s ∂ x U, Λ s U ). We have: then it holds that: To control A 1 , by integration by parts, we have: Clearly, it holds that: By integrating by parts, it holds that: . Now using the fact that: x N , for any differentiable functions M , N and by integration by parts, we have: Although A 131 can be controlled directly with √ ε in front of the constant, one may improve this by ε instead. Indeed by integration by parts one has: Remark that h x = εζ x , then A 1311 posses sufficient ε's, unlike A 1312 on which we have to work a little more. Indeed, in view of (3.1) we have that h −1 > 0, then it holds: Again by integration by parts, we get : . Therefore one may control A 1312 by εC(h −2 min , |ζ| W 1,∞ , µ|v xxx | ∞ )E s (U ) 2 . Consequently, it holds: Collecting the information provided above we get: To control A 2 + A 3 , by remarking firstly that J[h] and T [h] are symmetric, and then by integration by parts after having performing some algebraic calculations and using (3.12), we have: Unfortunately, an inconvenient term appears in This term won't be controlled without gaining √ ε taken from h xx = εζ xx and the other √ ε sits in front of the constant. Due to this fact, it follows that: To control A 4 , by integration by parts, it holds: To control A 5 , by integration by parts, we have: Therefore, it holds that: Finally, by integration by parts, A 6 is controlled by εC |v x | ∞ E s (U ) 2 . Therefore, it holds: • Estimation of Λ s , A[U ] ∂ x U, SΛ s U . Let us remark that: To control B 1 , we use the expression of J[h] to write: Then by using the fact that: The √ ε in front of the constant is due to the inconvenient term represented by To control B 2 , by the expression of J[h] and (3.13), we have: Then, clearly the following estimate holds: To control B 3 , we have that is symmetric and that: Therefore, one may write:

At this point, using the expressions of T [h] and J[h], it holds:
Therefore, it holds that: which implies that: Thanks to the fact that, for all k ∈ N, h k − 1 = O(εζ) and using the explicit expression of combined with the identities: then by integration by parts and (3.10), it holds that: Also, by (3.10) it holds: and For controlling B 35 , the explicit expression of T [h] and (3.14) gives that: Thus, as a conclusion, it holds that: To control B 4 , as for B 3 and using (3.14) one may write: To control B 5 , using the expression of , (3.10) and (3. To control B 6 , using the same arguments as the ones used to control B 3 , using expression of , (3.10) and (3.14), it follows that: Now, using the expression of Q with the help of Lemma 2, estimate (3.10), in addition to (3.14) and the fact that [Λ s , ∂ x (M ·)]N = ∂ x [Λ s , M ]N , it holds: Eventually, as a conclusion, one gets: It is worth noticing that √ ε in front of the constant is due to B 1 and B 31 . • Estimation of Λ s ζ, [∂ t , J[h]]Λ s ζ . Using the expression of J[h] and by integration by parts, it holds that: • Estimation of Λ s v, [∂ t , ]Λ s v . It holds that: then by integration by parts: Finally, combining the above estimates in addition to that fact that H s (R) is continuously embedded in W 1,∞ (R), it holds that: Taking λ = λ T large enough (how large depending on sup C(h −1 min , E s (U )) such that the right hand side of the inequality above is negative for all t ∈ [0, T √ ε ], then it holds that: Thanks to Grönwall's inequality so that it holds and hence the desired energy estimate is finally obtained.

Main results.
Well-posedness of the extended Boussinesq system. Theorem 1 represents the well-posedness of the extended Boussinesq system (2.2) which holds in X s = H s+2 (R) × H s+2 (R) as soon as s > 3/2 on a time interval of size 1/ √ ε.
Theorem 1 (Local existence). Suppose that U 0 = (ζ 0 , v 0 ) ∈ X s satisfying (3.1) for any t 0 > 1 2 , s ≥ t 0 + 1. Then there exists a maximal time T max = T (|U 0 | X s ) > 0 and a unique solution U = (ζ, v) T ∈ X s Tmax to the extended Boussinesq system (2.2) with initial condition (ζ 0 , v 0 ) such that the non-vanishing depth condition (3.1) is satisfied for any t ∈ [0, Tmax √ ε ). In particular if T max < ∞ one has Proof. The proof follows same line as [25, Theorem 1] using the energy estimate proved in Proposition 2. This is due to the fact that in [25] a most general case is considered (i.e. the extended Green-Naghdi equations). Remark that the proof itself is an adaptation of the proof of the well-posedness of hyperbolic systems (see [1] for general details).
A stability property. Theorem 1 is complemented by the following result that shows the stability of the solution with respect to perturbations, which is very useful for the justification of asymptotic approximations of the exact solution. (The solution U = (ζ, v) T and time T max that appear in the statement below are those furnished by Theorem 1).
Theorem 2 (Stability). Suppose that the assumption of Theorem 1 is satisfied and moreover assume that with h(t, x) = 1 + ε ζ(t, x) and F = ( Proof. The proof consists on the evaluation of 1 2 d dt U 2 X s (R) . Knowing that fact, by subtracting the equations satisfied by U = (ζ, v) T and U = ( ζ, v) T , we obtain: Consequently, a similar energy estimate evaluation as in Proposition 2 yields the desired result.
Convergence. As a conclusion, the following convergence result states that the solutions of the full Euler system, remain close to the ones of the system we are considering, namely system (2.2), with a better precision as ε 3 is smaller.
Theorem 3 (Convergence). Let ε ∈ (0, 1), s > 3/2, and U 0 = (ζ 0 , ψ 0 ) T ∈ H s+N (R) 2 satisfying condition (3.1) where N is large enough, uniformly with respect to ε ∈ (0, 1). Moreover, assume U euler = (ζ, ψ) T to be a unique solution to the full Euler system (1.3) that satisfies the assumption of Proposition 1. Then there exists C, T > 0, independent of ε, such that • Our new model (2.2) admits a unique solution U xB = (ζ xB , v xB ) T , defined on [0, T √ ε ] with corresponding initial data (ζ 0 , v 0 ) T ; • The error estimate below holds, at any time Proof. The first point is provided by the local existence result Theorem 1. Thanks to Proposition 1, then the solution of the water wave equations (ζ, v) T solve our model (2.2) up to a residual R of order ε 3 . The error estimation then follows from the stability Theorem 2.

Explicit Solitary
Wave Solution of the extended Boussinesq system. Solitary waves were initially discovered in shallow water by J.S. Russell during his experiments to design a more dynamic canal boat [12]. Many partial differential equations have been derived in the literature to model the solitary wave observed by Russell. Such models are commonly known as the Korteweg-de Vries (KdV) scalar equation for a unidirectional flow or the coupled Boussinesq and Green-Naghdi evolution equations. These famous nonlinear and dispersive models describe the shallow water waves and admit explicit families of solitary wave solutions [4,33,26,39,8]. The explicit solitary solutions of different nonlinear PDE's can be calculated using many methods. One of these methods is replacing the partial differential equation by an ordinary one (ODE) and thus one can look for explicit solutions in terms of particular functions. This replacement can be done by setting a reference traveling wave and hence one look for traveling-wave solutions. In this section, we seek the explicit solution of traveling waves for the extended Boussinesq system. Let us recall that the extended Boussinesq system that we are considering can be written as: where h(t, x) = 1 + εζ(t, x) and denote by (4.2) In order to find solitary wave solutions of the extended Boussinesq system (4.1), we seek solutions in the form of the traveling wave ζ(t, where the constant c ∈ R is the velocity of the solitary wave. Plugging the above Ansatz into eq. (4.1) yields: We may now integrate and, using the vanishing condition at infinity to set the integration constant, we deduce from the first equation: Using (4.4), one can deduce that v c = cζ c + O(ε).
One can also check the following identity ζ c ζ c = (ζ c ζ c ) − 1 2 (ζ c ) 2 is true. Using the latter identities into the second equation of (4.3), we may now integrate and, using the vanishing condition at infinity to set

Explicit solution with correctors of order O(ε 3 ) for the extended Boussinesq equations
Another approach of dealing with nonlinear PDE's when looking for analytical exact solution is finding instead an explicit solution with correctors. Explicit solutions with correctors for asymptotic water waves models have been obtained in [21,18]. Actually, H s -consistent solutions are obtained to the models in the variable topography case using the analytic solution of the model in the flat topography configuration.
In what follows, we find an explicit solution with correctors of order O(ε 3 ) for the extended Boussinesq model (4.1) and validate the result numerically.
We start by defining an H s -consistent solution or in other words explicit solution with correctors of order O(ε 3 ).
The standard Boussinesq system can be easily obtained form the extended Boussinesq system (4.1) by dropping all terms of order O(ε 2 ). Thus the standard Boussinesq system can be written as: then (ζ, v) = (ζ 1 , v 1 ) + ε 2 (ζ 2 , v 2 ) is H s -consistent with the extended Boussinesq system (4.1).
Proof. First, we would like to mention that we denote by O(ε) any family of functions (f ε ) 0<ε<1 such that ( 1 ε f ε ) 0<ε<1 remains bounded in L ∞ [0, T √ ε ], H r (R) , for possibly different values of r. We may now proceed in proving the stated result.

Numerical validation
In this section, we numerically validate the result of Theorem 5. In fact, we consider the equations given by system (4.1) and we compute explicitly the solutions given by (5.10) and (5.11). Then, we compute the residues for both equations after substituting (5.10) and (5.11) correspondingly. First we have to set the initial conditions ζ 0 2 = v 0 2 = exp − 3πx 10 2 . We also choose the constant α = 1. The residues R 1 (ε) and R 2 (ε) of the first and second equation of the system (4.1) respectively, are defined as follow: where p ∈ {2, ∞}. The residues R p 1 (ε) and R p 2 (ε) for p = 1 and p = ∞ are computed for several values of ε, namely ε = 10 −1 , 10 −2 , 10 −3 , 10 −4 and 10 −5 , at time t = 1. The results are summarized in Table 1 and Figures 4 and 5 where we plot in a log-log scale the residues R p 1 and R p 2 for p = 1 and p = ∞ in terms of ε.  Table 1. The residues R 1 (ε) and R 2 (ε) for p = 2 (left) and p = ∞ (right)  One clearly sees that the curves of the residues for both p = 1 and p = ∞ are both parallel to ε 3 . This shows that the residues convergence rate is O(ε 3 ), which is in total agreement with our theoretical result.