S Bellec

M Colin

M Ricchiuto

Discrete asymptotic equations for long wave propagation

Keywords: Asymptotic equations, discrete asymptotics, Boussinesq models, finite elements * Inria BSO, Team CARDAMOM Asymptotic equations, discrete asymptotics, Boussinesq models, finite elements

In this paper, we present a new systematic method to obtain some discrete numerical models for incompressible free-surface flows. The method consists in first discretizing the Euler equations with respect to one variable, keeping the other ones unchanged and then performing an asymptotic analysis on the resulting system. For the sake of simplicity, we choose to illustrate this method in the context of the Peregrine asymptotic regime, that is we propose an alternative numerical scheme for the so-called Peregrine equations. We then study the linear dispersion characteristics of our new scheme and present several numerical experiments to measure the relevance of the method.

Introduction

Wave transformation in near shore zone is well-described by the incompressible Euler equations. Due to their three-dimensional character, these equations are often too costly if one wants to perform numerical experiments, and often replaced by asymptotic depthaveraged models known as Boussinesq equations. A major characteristic of these models is their ability to describe the dispersive behaviour of wave propagation. Generally, the linear and nonlinear dispersion characteristics of the waves represented by Boussinesq models can be improved by including high order contributions in the double asymptotic expansion in terms of the ratios wave height over wave length (dispersion) and wave height over depth (nonlinearity) [START_REF] Lannes | Derivation of asymptotic two-dimensional timedependent equations for surface water wave propagation[END_REF]. Other techniques to improve the linear dispersion characteristics involve the inclusion of extra dispersive differential terms, derived either from a linear wave equation [START_REF] Beji | A formal derivation and numerical modeling of the improved Boussinesq equations for varying depth[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF], or by replacing depth-averaged values by point values at a properly chosen depth [START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF]. When numerically simulating the propagation of long waves, the physics represented by these continuous systems of Partial Differential Equations (PDE) is further filtered by the numerical scheme, and in particular by the form of the truncation error. For most of Boussinesq models, the task of designing an accurate numerical discretization is a nontrivial one, due to the presence of dispersion terms. Several approaches exist in literature, each with its own advantages and drawbacks. For details, the interested reader may refer to [START_REF] Bonneton | A splitting approach for the fully nonlinear and wakly dispersive Green-Naghdi model[END_REF][START_REF] Eskilsson | Spectral/HP discontinuous Galerkin methods for modelling 2D Boussinsq equations[END_REF][START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF][START_REF] Ricchiuto | Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF], to the review [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF], and references there in. The objective of this paper is to study the interaction scheme-PDE and to propose a framework to obtain new schemes with improved characteristics w.r.t. existing approaches. For this purpose, we introduce a new scheme reversing the model derivation procedure. More precisely, we propose to discretize partially the incompressible Euler equations with respect to one direction using a finite element method, and then follow Peregrine's derivation procedure. This new paradigm leads to a very promising scheme with nice dispersion properties.

The paper is organized as follows. In Section 2 we introduce some notation and recall our main result. The finite element discretization of a well known Boussinesq system, and most of the algebraic operators involved in our analysis are presented in Section 3. In Section 4, we detail the derivation of the new numerical scheme. The theoretical analysis of these discrete asymptotic models is presented in Section 5. Finally, Section 6 presents a numerical evaluation of the performances of the schemes, confirming our theoretical results. The paper is ended by an overlook of future developments related to the new approach proposed.

2 Setting, notations and main result.

Before going further, let us introduce some notations. For simplicity, in this article, we only deal with 2-D and 1-D problems. Denote by (x, z) respectively the horizontal and the vertical spatial dimension. Denote by d(x) the depth at still water and η(t, x) the surface elevation from its rest position. The total depth is then h(t, x) = d(x) + η(t, x) (see Figure 1). Let a be a typical wave amplitude, d 0 a reference water depth and λ a typical wavelength. In view of performing an asymptotic analysis, we introduce the nonlinearity parameter ε and the dispersion parameter σ defined by

ε = a d 0 , σ = d 0 λ .
Under the Boussinesq hypothesis ε = O(σ 2 ), Peregrine (see [START_REF] Peregrine | Long waves on a beach[END_REF]) first derived, from the Euler equations, the following standard system of Boussinesq's type

η t + [hū] x = 0, ūt + ūū x + gη x + ( d 2 6 ūtxx -d 2 [dū] txx ) = 0.
(

) 1
The model describes the evolution of the depth-averaged velocity ū and the surface elevation η within an accuracy of O(εσ 2 , σ 4 ) w.r.t. the Euler equations. The set of Equations ( 1) is now well-understood from the computation point of view and a classical numerical scheme can be obtained by using the finite element method in the following stetting. On an intervall [r, s], we introduce a set of nodes

r = x 0 < x 1 < • • • < x N = s,

Inria

where, for simplicity, we take a constant space step ∆ x = x j+1 -x j , ∀j ∈ {0, . . . , N }. We denote by E, Ū , D and H the vectors of the nodal values of η, ū, d and h. Similarly to what has been done in [START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF][START_REF] Walkley | A Numerical Method for Extended Boussinesq Shallow-Water Wave Equations[END_REF] (cf. also [START_REF] Ricchiuto | Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF] and references therein), we apply the P 1 Galerkin method to approximate the variational form of (1). In particular, we denote by {ϕ j } 0≤j≤N the standard piecewise linear continuous Lagrange basis, and introduce the discrete velocity, wave height and depth polynomials as follows

ū∆ (t, x) = N j=0 ūj (t)ϕ j (x), η ∆ (t, x) = N j=0 η j (t)ϕ j (x), d ∆ (x) = N j=0 d j ϕ j (x). (2)
The Galerkin approximation of (1), under the hypothesis of exact integration w.r.t. all the discrete polynomials involved, can be written in a compact matrix form

ME t + 1 3 2N (H Ū ) + H (N Ū ) + Ū (N H) = 0, M Ūt + 1 3 N ( Ū 2 ) + Ū (N Ū ) + gN E - 1 6 {D; Ūt } = 0, (3)
where the operation A B represents, for two given vectors A = (a i ) 0≤i≤N and B = (b i ) 0≤i≤N , the column vector (a i b i ) 0≤i≤N (note that, given a column vector A, notation A 2 simplifies A A). Remark that the matrices M, N and Q are the usual mass, derivation, and stiffness matrices arising in the Galerkin discretization and are detailed in the next section. Moreover, for given columns vectors A and B, we set

{A; B} = Q(A 2 B) + D (Q(A B) + 2(A B) (QA) -B (QA).
The aim of this paper is to propose a systematic method to obtain new numerical models describing free surface flows. It is based on the following idea : reverse the model derivation procedure and first discretize partially the Incompressible Euler equations and then derive fully discrete asymptotic equations by performing an asymptotic analysis. To illustrate the potential of this idea, we apply this method to the couple Euler-Peregrine equations by applying the Galerkin method to the variable x and then performing the asymptotic analysis of Peregrine's type to the resulting equations. Of course, when one deals with non-linear equations, this procedure does not commute with the classical one. In this paper, for simplicity, we deal with periodic boundary conditions. Note that our strategy can be adapted directly to the Dirichlet boundary conditions, since the crucial property (50) hold for this two kind of conditions. This strategy is similar to the one proposed for compressible multiphase flows in [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF]. As shown in the detailed derivation of the next sections, the new procedure leads to the following discrete equations approximating the discretized Euler system within an accuracy of O(εσ 2 , σ 4 )

ME t +M[H; Ū ] = 0, M Ūt + 1 3 N ( Ū 2 ) + Ū (N Ū ) + gN E + M D 2 6 (K 2 Ūt ) - D 2 (K[D; Ūt ]) = 0, (4)
having introduced the operator [• ; •] defined by

[A; B] = A (KB) + 1 3 K(A B) -M -1 (A (N B)) + 2M -1 (B (N A)) ,
with K = M -1 N . We see that, while involving similar algebraic operations, the new discretization is different from the classical ones, even for a simple case like Peregrine equations. The main differences are found in the treatment of the third order derivatives terms as well as in the nonlinear ones in the continuity (wave height) equation. We will show that scheme (3) also converges to an approximation of the Perigrine equations. However, both the linear phase relation, and the lienar shoaling gradient provided by (3) (see Section 5 and 6) are substantially closer to the exact ones than those given by [START_REF] Bonneton | A splitting approach for the fully nonlinear and wakly dispersive Green-Naghdi model[END_REF]. this gain in accuracy is confirmed by the numerical tests showing a reduction of the error of 5 to 10 times w.r.t. scheme (4). In the sequel, we show how to derive the scheme (3) and prove that not only they are consistent with system (1), but that they represent a substantial improvement w.r.t. the scheme obtained by discretizing directly the asymptotic equations (1). The paper is organized as follows. In Section 3, we recall the classic discretization of Peregrine system (1). In Section 4, we detail the derivation of the new numerical scheme. In Section 5, we present the study of the exact linear dispersion characteristics for the two schemes presented before. Then, we compare the results with the theoretical characteristics of Peregrine model. To conclude, in Section 6, we present two numerical computations in order to evaluate the accuracy of the new scheme.

Galerkin discretization of the Peregrine equations

We start by deriving a spatially discretized version of (1) using the Galerkin finite element scheme. Since we deal with varying bathymetry, we first rewrite (1) by developing the derivatives terms

η t + [(d + η)ū] x = 0, (5)
ūt + ūū x + gη x - d 2 3 ūtxx -dd x ūtx - dd xx 2 ūt = 0. (6)
According to the Galerkin method (see [START_REF] Ern | Theory and practice of finite elements[END_REF]), we multiply (5) by the Galerkin basis φ i and we integrate the resulting equation between x i-1 and x i+1 , for i ∈ {0, ..., N }, with the convention x -1 = x N and x N +1 = x 0 coming from the choice of the periodic boundary conditions to obtain

xi+1 xi-1 η ∆ (t, x)ϕ i (x)dx - xi+1 xi-1 (d ∆ (x) + η ∆ (t, x))ū ∆ (t, x)ϕ i (x)dx = 0. (7)
In the P 1 case, evaluating all the integrals exactly w.r.t. the linear spatial variation of d ∆ , η ∆ , and ū∆ , we derive from [START_REF] Ern | Theory and practice of finite elements[END_REF], using the decomposition (2) and evaluating all the integrals involving φ i , for all i ∈ {0, ..., N }

∆ x 6 ∂ ∂t (η i+1 + 4η i + η i-1 ) + 1 3 2 h i+1 ūi+1 -h i-1 ūi-1 2 + h i ūi+1 -ūi-1 2 + ūi h i+1 -h i-1 2 = 0. (8)
As already mentioned, for simplicity, the boundary conditions are assumed to be periodic, although other conditions can be embedded in the analysis. Applying the same procedure to Equation ( 6), one has, performing integration by parts, for all i ∈ {0, ..., N },

xi+1 xi-1 d dt ū∆ (t, x)ϕ i (x)dx - xi+1 xi-1 ū2 ∆ (t, x) 2 + gη ∆ (t, x))ϕ i (x)dx + xi+1 xi-1 -ϕ i (x)dd ∆ (ū ∆ ) tx + (d 2 ∆ (x)ϕ i (x)) x 3 (ū ∆ (t, x)) tx + (ϕ i (x)d ∆ (x)(ū ∆ (t, x)) t ) x 2 d ∆ dx = 0. (9)

Inria

Plugging relations (2) in [START_REF] Filippini | On the nonlinear behaviour of Boussinesq type models : amplitude-velocity vs amplitude-flux forms[END_REF], one derives after the computation of the different integrals, for all i ∈ {0, ..., N }

∆ x 6 d dt (ū i+1 +4ū i + ūi-1 ) + 1 3 ū2 i+1 -ū2 i-1 2 + ūi ūi+1 -ūi-1 2 + g η i+1 -η i-1 2 - 1 6∆ x d dt [2d i ūi (d i+1 -2d i + d i-1 ) + (d 2 i+1 ūi+1 -2d 2 i ūi + d 2 i-1 ūi-1 ) +d i (d i+1 ūi+1 -2d i ūi + d i-1 ūi-1 ) -ūi (d 2 i+1 -2d 2 i + d 2 i-1 )] = 0. (10)
In order to rewrite the discrete equations ( 9) in a simpler and more compact form, we define the mass matrix M, the first order derivative matrix N and the second order derivative matrix Q (which represents the laplacian operator with periodic boundary conditions) as follows

M = 1 6             4 1 0 • • • 0
1 0 • • • 0 1 4             , N = 1 2∆ x             0 1 0 • • • 0 -
. . . 1 1 0 • • • 0 -1 0             Q = 1 ∆ 2 x             -2 1 0 • • • 0 1 1 . . . . . . . . . 0
0 • • • 0 1 -2            
.

In addition, we introduce the following column vectors

Ū = (ū i ) 0≤i≤N , E = (η i ) 0≤i≤N , D = (d i ) 0≤i≤N , H = D + E.
Equation (8) can thus be written as

∆ x d dt ME + 1 3 2N (h i ūi ) i∈{0,..,N } + (h i (N Ū ) i ) i∈{0,..,N } + (ū i (N H) i ) i∈{0,..,N } = 0.
To further simplify the notation for given columns vectors A = (a i ) 0≤i≤N and B = (b i ) 0≤i≤N , we introduce the operator :

R N × R N → R N (A, B) → A B := (a i b i ) 0≤i≤N
In the sequel, for simplicity A 2 simplifies A A. As an example, the vector (h i (N Ū ) i ) i∈{1,..,n} can be rewritten as H (N Ū ). Owning this notation, the discrete equations ( 8)-( 10) can be compactly written as

M d dt E + 1 3 2N (H Ū ) + H (N Ū ) + Ū (N H) = 0, (11)
d dt M Ū + 1 3 N ( Ū 2 ) + Ū (N Ū ) + gN E - 1 6 {D; Ūt } = 0, (12)
RR n°8806

where we have introduced

{D; Ūt } = Q(D 2 Ūt ) + D (Q(D Ūt )) + 2(D Ūt ) (QD) -Ūt (QD 2 ). (13)
Equations ( 11)-( 12) will be taken in the sequel as the classical scheme for the Peregrine equations and be used in Sections 5 and 6 to make some comparison with the new scheme introduced in the next section.

4 A new setting for deriving discrete asymptotic models.

4.1 Semi-discretization of the 2D-Euler equations in non-dimensional form.

The aim of this section is to derive an alternative set of discrete equations, possibly having improved characteristics w.r.t. ( 11)-( 12), for example a better evaluation of the shoaling gradient phenomenon. For that purpose, we propose to discretize the 2D-Euler equations with respect to one direction, x for example, and then to perform an asymptotic analysis on the resulting equation, similar to the one used to obtain the Peregrine equations ( 1). The Euler equations written in terms of velocity (u, w), pressure p, constant density ρ and vertical gravity acceleration g reads :

u t + uu x + wu z + p x ρ = 0, (14)
w t + uw x + ww z + p z ρ + g = 0, (15)
u x + w z = 0, (16)
u z -w x = 0, (17)
where the last equation represents the irrotationality condition. In this paper, since our aim is to obtain a new scheme for the Peregrine system, we restrict ourselves to the 2D version of the Euler equation. We deal with periodic boundary conditions in the x direction, while on the free surface and sea-bed level we use the classical conditions :

• at the free surface

z = η w = η t + uη x , (18)
p = 0, (19)
• on the seafloor z

= -d w = -ud x . ( 20
)
Let d 0 be the averaged depth, a a typical wave amplitude, and λ a typical wavelength. The following usual non-dimensional variables are introduced

x = x λ , z = z d 0 , t = √ gd 0 λ t, η = η a , ũ = d 0 a √ gd 0 u, w = λ a 1 √ gd 0 w, p = p gd 0 ρ , ∆ x = ∆ x λ .
Inria Classically, we now consider the nonlinearity and dispersion parameters ε and σ defined by

ε = a d 0 , σ = d 0 λ .
Using the notation introduced above, the Euler equations and the irrotationality condition can be recast in a non-dimensional form as

εũ t + ε 2 ũũ x + ε 2 wũ z + px = 0, (21)
εσ 2 wt + ε 2 σ 2 ũ wx + ε 2 σ 2 w wz + pz + 1 = 0, ( 22
) ũx + wz = 0, ( 23
) ũz -σ 2 wx = 0 (so ũz = O(σ 2 )), (24)
The boundary conditions become :

• at the free surface z = εη

w = ηt + εũη x, (25)
p = 0, (26)
• at the bed z = -

d w = -dx ũ. ( 27
)
Our goal is to obtain a Boussinesq's type approximation of the Euler system ( 21)-( 27), under the assumptions ε << 1, σ << 1, and in the specific regime ε = O(σ 2 ), meaning that there exists constant C > 0 such that ε ≤ Cσ 2 . Using (24), we obtain from (21)

εũ t + ε 2 ũũ x + px = -ε 2 σ 2 w wx . ( 28
)
We now apply a Galerkin method on the variable x keeping t and z unchanged. It is assumed that the number of mesh points per wavelength is small, in particular ∆ x < 1, however, we assume

∆ x = O(σ) (it transpires that ∆ x = O(d 0 )).
In the sequel we drop the " ˜" and we introduce for all i ∈ {0, ..., N },

u i (t, z) = u(t, x i , z), w i (t, z) = w(t, x i , z), η i (t, z) = η(t, x i , z), p i (t, z) = p(t, x i , z).
In addition, the discrete horizontal velocity, wave height, depth, vertical velocity and pressure polynomials are written in the Galerkin basis as follows

u ∆ (t, x, z) = N j=0 u j (t, z)ϕ j (x), w ∆ (t, x, z) = N j=0 w j (t, z)ϕ j (x), η ∆ (t, x) = N j=0 η j (t)ϕ j (x), p ∆ (t, x, z) = N j=0 p j (t, z)ϕ j (x), d ∆ (x) = N j=0 d j ϕ j (x). (29)
As in Section 3, we focus on periodic boundary condition that is we introduce x -1 = x N and x N +1 = x 0 . The finite element discrete equations corresponding to (28)-( 22)-( 23)-( 24) can be written as, for all i ∈ {0, ..., N }

ε ∆ x 6 d dt (u i+1 + 4u i + u i-1 )+ ε 2 3 u 2 i+1 -u 2 i-1 2 + u i u i+1 -u i-1 2 + p i+1 -p i-1 2 = - ε 2 σ 2 3 w 2 i+1 -w 2 i-1 2 + w i w i+1 -w i-1 2 , ( 30
) RR n°8806 εσ 2 ∆ x 6 d dt (w i+1 + 4w i + w i-1 )+ ∆ x 6 d dz (p i+1 + 4p i + p i-1 ) + ∆ x = -ε 2 σ 2 u i w i+1 -w i-1 2 -w i u i+1 -u i-1 2 , (31)
u i+1 -u i-1 2 + ∆ x 6 d dz (w i+1 + 4w i + w i-1 ) = 0, (32)
∆ x 6 d dz (u i+1 + 4u i + u i-1 ) -σ 2 w i+1 -w i-1 2 = 0. ( 33
)
For the boundary conditions, we propose to integrate (25)-( 26) along the curve z = εη and equation ( 27) along the curve z = -d. For that purpose, we choose to introduce

ŵ∆ = N i=0 w t, εη(t, x i ) φ i (x), w∆ = N i=0 w t, -d(x i ) φ i (x), û∆ = N i=0 u t, εη(t, x i ) φ i (x), ǔ∆ = N i=0 u t, -d(x i ) φ i (x). (34)
The boundary conditions ( 25)-( 26)-( 27) then write

• at the free surface

∆ x 6 w i+1 (t,εη i+1 ) + 4w i (t, εη i ) + w i-1 (t, εη i-1 ) = ∆ x 6 d dt (η i+1 (t) + 4η i (t) + η i-1 (t) + 1 3 εη i+1 (t)u i+1 (t, εη i+1 ) -εη i-1 (t)u i-1 (t, εη i-1 ) 2
-εη i (t) u i+1 (t, εη i+1 ) -u i-1 (t, εη i-1 ) 2 + 2u i (t, εη i ) εη i+1 (t) -εη i-1 (t) 2 , (35)
∆ x 6 p i+1 (t, εη i+1 ) + 4p i (t, εη i ) + p i-1 (t, εη i-1 ) = 0, (36)
• at the bed

∆ x 6 w i+1 (t, -d i+1 ) + 4w i (t, -d i ) + w i-1 (t, -d i-1 ) = - 1 3
d i+1 u i+1 (t, -d i+1 ) -d i-1 u i-1 (t, -d i-1 ) 2 -d i u i+1 (t, -d i+1 ) -u i-1 (t, -d i-1 ) 2 + 2u i (t, -d i ) d i+1 -d i-1 2 . ( 37
)
Introducing the following column vector

W = (w i ) 0≤i≤N , U = (u i ) 0≤i≤N , E = (η i ) 0≤i≤N , P = (p i ) 0≤i≤N , I =    1 . . . 1    Ŵ = w i (εη i ) 0≤i≤N , Û = u i (εη i ) 0≤i≤N , W = w i (-d i ) 0≤i≤N , Ǔ = u i (-d i ) 0≤i≤N , Inria
we can rewrite Equations ( 30)-(37) into the following matrix-form :

ε d dt MU + ε 2 3 N (U 2 ) + U (N U ) + N P = - ε 2 σ 2 3 N (W 2 ) + W (N W ) , (38)
εσ 2 d dt MW + d dz MP + I = -ε 2 σ 2 (U (N W ) -W (N U )) , (39)
N U + M d dz W = 0, (40)
M d dz U -σ 2 N W = 0. ( 41
)
The boundary conditions become

• at the free surface

M Ŵ = d dt ME + ε 3 [N (E Û ) -E (N Û ) + 2 Û (N E)], (42)
M P = 0, (43)
• at the bottom

M W = - 1 3 N (D Ǔ ) -D (N Ǔ ) + 2 Ǔ (N D) . (44)
System ( 38)-( 44) represents the first step in our analysis. The next two sections are dedicated to the transformation of this system into an asymptotic set of equations.

4.2 Asymptotic expansions on the velocity U and the pressure p.

In this section, we derive an asymptotic expansion in terms of σ for the semi-discrete horizontal velocity U = U (t, z) following the procedure presented by Peregrine in [START_REF] Peregrine | Long waves on a beach[END_REF]. More precisely, we prove the following proposition.

Proposition 1. The pressure p and the velocity U satisfy expansion of the form

P = εE -zI + εσ 2 z 2 2 KU 0 + z[D; U 0 ] + O(ε 2 σ 2 , εσ 4 ) U = Ū + σ 2 D 2 6 (K 2 Ū ) - z 2 2 K 2 Ū -zK[D; Ū ] - D 2 (K[D; Ū ]) + O(εσ 2 , σ 4 ),
where the averaged velocity is defined in (56).

Proof. Since M is invertible, we obtain from the integration of (41) between 0 and an arbitrary depth z,

U (t, z) = U 0 (t) + O(σ 2 ), (45)
where U 0 (t) is a constant depending only on t and corresponds to the value of U at z = 0. We directly deduce that

Ǔ (t) = U 0 (t) + O(σ 2 ). (46)
Substituting relation (45) in equation ( 40) and setting K = M -1 N , we derive

d dz W = -KU 0 + O(σ 2 ). (47)
Integrating each line i ∈ {0, ..., N } of equation ( 47) with respect to z between -d i and an arbitrary depth z (-d i < z < η i ), using the boundary condition (44) and the estimates (45) on U and (46) on Ǔ , we obtain

W = -zKU 0 -D (KU 0 ) - 1 3 [K(D U 0 ) -M -1 (D (N U 0 )) + 2M -1 (U 0 (N D))] + O(σ 2 ). (48)
Plugging ( 48) in (41) and integrating the resulting equation between 0 and z, one derives the following expansion on

U U =U 0 -σ 2 z 2 2 K 2 U 0 + zK(D (KU 0 )) + z 3 K[K(D U 0 ) -M -1 (D (N U 0 )) + 2M -1 (U 0 (N D))] + O(σ 4 ). (49)
Looking for a similar expansion on the pressure array P , we substitute Equation (48) in Equation ( 39) to obtain

d dz MP = -I -εσ 2 d dt zKU 0 + M(D [KU 0 ]) + 1 3 N (D U 0 ) -D [N U 0 ] + 2U 0 [N D] + O(ε 2 σ 2 , εσ 4 ).
Using the fact that

MI = I, (50)
the previous relation becomes

d dz P = -I -εσ 2 d dt zKU 0 + D [KU 0 ] + 1 3 K(D U 0 ) -M -1 (D (N U 0 )) + 2M -1 (U 0 (N D)) + O(ε 2 σ 2 , εσ 4 ). (51)
In view of (51), it is natural to introduce the following bracket

[A; B] = A (KB) + 1 3 K(A B) -M -1 (A (N B)) + 2M -1 (B (N A)) . (52)
Remark that (49) then becomes

U = U 0 -σ 2 z 2 2 K 2 U 0 + z[D; U 0 ] + O(σ 4 ). (53)
Furthermore, integrating (51) with respect to z from an arbitrary depth to the free surface εη i , we can write

P = εE -zI + εσ 2 z 2 2 KU 0 + z[D; U 0 ] + O(ε 2 σ 2 , εσ 4 ) (54)
Substituting equations ( 54) and ( 49) in (38), we obtain an equation for the zero-th order velocity U 0 , equivalent to Equation 2.28 in [START_REF] Walkley | A Numerical Method for Extended Boussinesq Shallow-Water Wave Equations[END_REF], which reads :

d dt MU 0 + ε 3 N (U 0 U 0 ) + U 0 (N U 0 ) + N E = O(εσ 2 , σ 4 ). ( 55
)

Inria

Note that the choice of the constant of integration in (45) is not unique. However it transpires that the choice of U 0 (which is the value of the horizontal velocity U at z = 0) is not optimal as observed in [START_REF] Walkley | A Numerical Method for Extended Boussinesq Shallow-Water Wave Equations[END_REF]. This is why, in the sequel, we are going to get rid of it by introducing the averaged velocity matrix Ū = (ū i ) 0≤i≤N where

ūi = 1 d i + εη i εηi -di u i dz, (56)
and by looking for the equation satisfied by Ū . In this direction, we first derive the relation between U 0 and Ū . Equation (53) provides, for all i ∈ {0, ..., N },

u i = u 0 i -σ 2 z 2 2 K 2 U 0 i + z K[D; U 0 ] i + O(σ 4 ),
and by integration between -d i and εη i , we immediately get, using Taylor expansion, ūi = u 0 i -

1 εη i + d i σ 2 εηi -di z 2 2 dz K 2 U 0 i + εηi -di zdz K[D; U 0 ] i + O(σ 4 ), = u 0 i - σ 2 (d i + εη i ) d 3 i 6 K 2 U 0 i - d 2 i 2 K[D; U 0 ] i + O(εσ 2 , σ 4 ), = u 0 i -σ 2 d 2 i 6 K 2 U 0 i - d i 2 K[D; U 0 ] i + O(εσ 2 , σ 4 ).
This furnishes the desired relation

Ū = U 0 -σ 2 D 2 6 (K 2 U 0 ) - D 2 (K[D; U 0 ]) + O(εσ 2 , σ 4 ), (57)
or equivalently

U 0 = Ū + σ 2 D 2 6 (K 2 U 0 ) - D 2 (K[D; U 0 ]) + O(εσ 2 , σ 4 ). (58)
Then it transpires that

U 0 = Ū + O(ε, σ 2 ). (59)
Substituting (59) in (58), we derive

U 0 = Ū + σ 2 D 2 6 (K 2 Ū ) - D 2 (K[D; Ū ]) + O(εσ 2 , σ 4 ). (60)
Finally, plugging (60) into (53), one obtains the expansion of U as a function of the depth averaged velocity

Ū U = Ū + σ 2 D 2 6 (K 2 Ū ) - z 2 2 K 2 Ū -zK[D; Ū ] - D 2 (K[D; Ū ]) + O(εσ 2 , σ 4 ), (61)
which ends the proof of Proposition 1.

Depth-averaged equations

The aim of this section is to provide the final new discrete numerical model of Peregrine's type. In order to derive the equation on Ū (known as the momentum equation in the literature), we substitute (60) in (55) to obtain :

d dt M Ū + ε 3 N ( Ū 2 ) + Ū (N Ū ) + N E + σ 2 M d dt D 2 6 (K 2 Ū ) - D 2 K[D; Ū ] = O(εσ 2 , σ 4 ) (62)
In addition, to derive an equation on E (that is the continuity equation), we combine (42) and (44) to get

Ŵ -W = d dt E + ε M -1 3 [N (E Û ) -E (N Û ) + 2 Û (N E)] + M -1 3 [N (D Ǔ ) -D (N Ǔ ) + 2 Ǔ (N D)]. (63)
We integrate each lines of (63) between -d i and εη i , for all i ∈ {0, ..., N }, to obtain εηi -di

(KU ) i dz + Ŵi -Wi = 0,
which can be recast as

E t + [H; Ū ] + B = 0, (64)
where

B = εηi -di (KU ) i dz 0≤i≤N -[H; Ū ] + ε 3 [K(E Û ) -M -1 (E (N Û )) + 2M -1 ( Û (N E))] - 1 3 [K(-D Ǔ ) -M -1 (-D (N Ǔ )) + 2M -1 ( Ǔ (N (-D)))]. (65)
We can remark that the expression B is no more than a discretized version of the socalled Leibniz' Rule 1 . As a consequence, it transpires that B has the same accuracy of order O(εσ 2 , σ 4 ) than that of the equations and then can be neglected in the sequel (see Proposition 2 below for more details). In order to be more precise, we compute explicitly B by taking successively z = εη i and z = -d i in (61) to obtain the values of Û and Ǔ :

Û = Ū + σ 2 D 2 6 (K 2 Ū ) - ε 3 E 2 2 K 2 Ū -εE K[D; Ū ] - D 2 (K[D; Ū ]) + O(εσ 2 , σ 4 ) Ǔ = Ū + σ 2 - D 2 3 (K 2 Ū ) + D 2 (K[D; Ū ]) + O(εσ 2 , σ 4 ). (66)
By substituting Equations ( 61) and (66) into Equation (65), this provides the complete expression of

B B = σ 2 1 6 D (K(D 2 (K 2 Ū ))) - 1 6 D 3 (K 3 Ū ) - 1 9 K(D 3 (K 2 Ū )) + 1 9 M -1 D (N (D 2 (K 2 Ū ))) + 1 2 D 2 (K 2 [D; Ū ]) - 1 2 D (K(D (K[D; Ū ]))) + 1 6 K(D 2 (K[D; Ū ])) - 1 6 M -1 D (N D (K[D; Ū ])) - 2 9 M -1 (N D) (D 2 (K 2 Ū )) + 1 3 M -1 ((N D) (D (K[D; Ū ]))) + O(εσ 2 , σ 4
). 1 We recall the Leibniz' Rule : Given f (x, z), a(x) and b(x), where f and ∂f ∂x are continuous in x and z, and a and b are differentialble functions of x,

∂ ∂x b(x) a(x) f (x, z)dz = b(x) a(x) ∂ ∂x f (x, z)dz + f (x, b(x))b (x) -f (x, a(x))a (x).

Inria

Finally, our new non-dimensionalized system reads (note that we have multiply (64) by M)

d dt ME + M[H; Ū ] + MB = 0. ( 67
)
d dt M Ū + ε 3 N ( Ū 2 ) + Ū (N Ū ) + N E + σ 2 M d dt D 2 6 (K 2 Ū ) - D 2 K[D; Ū ] = O(εσ 2 , σ 4 ) (68)
To go further, we now investigate the behavior of the vector B by establishing the following proposition.

Proposition 2. For any bathymetry d, the additional term B in Equation (67) satisfies

B = O(εσ 2 , σ 4 ).
As a consequence, the numerical scheme (67)-(68) becomes

d dt ME + M[H; Ū ] = O(εσ 2 , σ 4 ), (69)
d dt M Ū + 1 3 N ( Ū 2 ) + Ū (N Ū ) + N E -σ 2 d dt d 2 0 3 (N K Ū ) = O(εσ 2 , σ 4 ), (70)
and is consistent with the Peregrine Equations (1).

Proof. For a better understanding, we first assume that the batymetry d = d 0 is constant.

In this setting, one has D = d 0 I and the operator D is no more than the multiplication by the real d 0 , that is, for example, D U = d 0 U . Hence B is equal to

B = σ 2 d 3 0 6 K 3 Ū - d 3 0 6 K 3 Ū - d 3 0 9 K 3 Ū + d 3 0 9 K 3 Ū + d 3 0 2 K 3 Ū - d 3 0 2 K 3 Ū + d 3 0 6 K 3 Ū - d 3 0 6 K 3 Ū + O(εσ 2 , σ 4 ). = 0 + O(εσ 2 , σ 4 ).
More generally, assume now that the bathymetry d is not constant. For any regular function u and its discrete version (u i ) 0≤i≤N , a Taylor expansion provides

u i+1 = u i + ∆ x u x (x i ) + ∆ 2 x 2 u xx (x i ) + ∆ 3 x 6 u xxx (x i ) + ..., (71)
and

u i-1 = u i -∆ x u x (x i ) + ∆ 2 x 2 u xx (x i ) - ∆ 3 x 6 u xxx (x i ) + ... (72)
Combining ( 71) and (72), we can prove that, for all i ∈ {0, ..., N }

(N U ) i = u x (x i ) + ∆ 2 x 6 u xx (x i ) + O(∆ 4 x ), (KU ) i = u x (x i ) + O(∆ 4 x ), (M -1 U ) i = u(x i ) - ∆ 2 x 6 u xx (x i ) + O(∆ 4 x ).
Plugging these expansions in Equations ( 67) and (68), we obtain, dividing by ∆ x at each point x i

η t + [hū] x - ∆ 2 x 6 η txx + [hū] xxx + h xx ūx + σ 2 ( d 2 d xx ūxxx 6 + 7 6 dd x d xx ūxx + ( 3 2 dd 2 xx d 2 x d xx + 5 6 dd x d xxx )ū x + dd xx d xxx ū) + O(∆ 4 x ) = O(εσ 2 , σ 4 ), ūt + ūū x + η x + σ 2 d 2 6 ūtxx - d 2 [dū] txx + ∆ 2 x 6 ūt + ūū x + η x + σ 2 d 2 6 ūtxx -σ 2 d 2 [dū] txx xx -ūx ūxx + σ 2 d 2 (d xx ūx ) x + O(∆ 4 x ) = O(εσ 2 , σ 4 ),
proving that our numerical scheme is consistent with the continuous Peregrine equations(1).

In addition, B is equal to

B = -σ 2 ∆ 2 x 6 d 2 d xx ūxxx 6 + 7 6 dd x d xx ūxx + ( 3 2 dd 2 xx d 2 x d xx + 5 6 dd x d xxx )ū x + dd xx d xxx ū + O(∆ 4 x ) + O(εσ 2 , σ 4 ),
from which it transpires that B contains only terms of order ∆ 2 x σ 2 , εσ 2 or σ 4 (actually, B is consistent with Leibniz' Rule). Recalling that ∆ x = O(σ), one has

B = σ 2 O(∆ 2 x) + O(εσ 2 , σ 4 ) = O(εσ 2 , σ 4 ),
which ends the proof of Proposition 2

To end this section, we return to the physical variables and neglect the contribution of B in (67)-(68) to obtain our new numerical scheme for the Peregrine Equations (1)

d dt ME + M[H; Ū ] = 0, (73)
d dt M Ū + 1 3 N ( Ū 2 ) + Ū (N Ū ) + gN E + M d dt D 2 6 (K 2 Ū ) - D 2 (K[D; Ū ]) = 0. ( 74
)
5 Study of the linear dispersion characteristics.

The aim of this section is to give some insights to measure the accuracy of the new method developped in the previous sections. For that purpose, we exhibit the dispersion relation as well as the shoaling coefficients of the linearized version of the scheme (73)-( 74). This study is widely inspired by the one proposed by Dingemans in [START_REF] Dingemans | Water wave propagation over uneven bottoms[END_REF] in the context of slowlyvarying water depth, that is we assume that d = d(βx) with β << 1. For the sake of completness, we also compare our computations with the ones performed on the linearized version of the classical scheme ( 11)-(12).

Linear characteristics of the new numerical model.

We first introduce the linearized version of the scheme (73)-(74) around the rest state which reads

d dt ME + M[D; Ū ] = 0, (75)
Inria d dt M Ū + gN E + M d dt D 2 6 (K 2 Ū ) - D 2 (K[D; Ū ]) = 0. (76)
As usual, when one deals with linear equations, a lot of computations can be performed explicitly. Indeed, differentiating (76) with respect to t, multiplying (75) by N and substituting the resulting equations, one obtains a decouple equation on the vector Ū :

M Ūtt -gN [D; Ū ] + M D 2 6 (K 2 Ū ) - D 2 (K[D; Ū ]) tt = 0. (77)
In order to exhibit the dispersion relation associated with (75)-(76), we then look for a plane-wave solution under the form

Ū = ūi 0≤i≤N where ūi = ū(t, x i ) with ū(t, x) = U(βx) exp -jωt + j β K(βx) , j 2 = -1. (78)
Owning the solution Ū , it is pertinent to introduce the wave number

k(βx) = ∂ ∂x 1 β K(βx) ,
and for all i = 0, ..., N , k i = k(βx i ). Then, we determine conditions on k and U so that Ū is a solution to the linear system (77). A Taylor expansion around the point x = x i provides directly

ūi+1 = 1 + β(j ∆ 2 x 2 k (βx i ) + ∆ x U (βx i ) U i ) ūi e jki∆x + O(β 2 ), (79) ūi
-1 = 1 + β(j ∆ 2 x 2 k (βx i ) -∆ x U (βx i ) U i ) ūi e -jki∆x + O(β 2 ). ( 80
)
In view of (77), we deduce that, ∀i ∈ {1, .., n}

(N Ū ) i = jk i sin c (k i ∆ x ) -β k 2 i ∆ 2 x 2 sin c (k i ∆ x ) k (βx i ) k i +β cos(k i ∆ x ) U (βx i ) U i ūi + O(β 2 ). ( 81
) (M Ū ) i = 1 3 (2 + cos(k i ∆ x )) + jβ k i ∆ 2 x 6 (cos(k i ∆ x ) k (βx i ) k i + 2 sin c (k i ∆ x ) U (βx i ) U i ) ūi + O(β 2 ). ( 82
)
Note that it is not possible to plug directly (81)-( 82) into (77), due to the presence of the vector (M -1 Ū ) in the bracket [D; Ū ]. Indeed, it is necessary to express each term (M -1 Ū ) with respect to ūi . To overcome this difficulty, the idea is to introduce the following new variables :

V = M -1 N D Ū , X = M -1 N Ū , Z = M -1 D (N Ū ) , W = M -1 ( Ū (N D)), Y = M -1 N X, T = D X + 1 3 (V -Z + 2W ) , S = M -1 N T.
and to perform asymptotic expansions of order β 2 on these variables. Using these new vectors, one can rewrite Equation (77) into

M Ūtt -gN T + M D 2 6 Y tt - D 2 S tt = 0 (83)
Note first that the vectors V , X, Z, W , Y and T depends only on Ū . It is then natural to introduce the following functions

v = V(βx) exp -jωt + j β K(βx) , x = X(βx) exp -jωt + j β K(βx) , z = Z(βx) exp -jωt + j β K(βx) , w = W(βx) exp -jωt + j β K(βx) , y = Y(βx) exp -jωt + j β K(βx) , t = T(βx) exp -jωt + j β K(βx) , s = S(βx) exp -jωt + j β K(βx) ,
where V(βx), X(βx), Z(βx), W(βx), Y(βx), T(βx) and S(βx) denote the amplitude of repectively v, x, z, w, y, t and s. Using Equations ( 81) and ( 82), one can rewite (83) for all i ∈ {1, .., n}

-ω 2 1 3 (2 + cos(k i ∆ x )) + jβ k i ∆ 2 x 6 (cos(k i ∆ x ) k (βx i ) k i + 2 sin c (k i ∆ x ) U (βx i ) ūi ) ūi -g jk i sin c (k i ∆ x ) -β k 2 i ∆ 2 x 2 sin c (k i ∆ x ) k (βx i ) k i + β cos(k i ∆ x ) T (βx i ) T(βx i ) T(βx i ) - ω 2 6 1 3 (2 + cos(k i ∆ x )) + jβ k i ∆ 2 x 6 (cos(k i ∆ x ) k (βx i ) k i + 2 sin c (k i ∆ x ) (d 2 (βx i )Y(βx i )) d 2 i Y(βx i ) ) d 2 i Y(βx i ) + (ω 2 ) 2 1 3 (2 + cos(k i ∆ x )) + jβ k i ∆ 2 x 6 (cos(k i ∆ x ) k (βx i ) k i + 2 sin c (k i ∆ x ) (d i (βx i )S(βx i )) dS(βx i ) ) d i S(βx i ) = 0.
It remains to estimate the first order derivatives U (βx i ), T (βx i ), Y (βx i ) and S (βx i ).

First observe that MX = N Ū which provides, using (81)-(82),

1 3 (2 + cos(k i ∆ x )) + jβ k i ∆ 2 x 6 cos(k i ∆ x ) k (βx i ) k i + 2 sin c (k i ∆ x ) X (βx i ) X(βx i ) X(βx i ) = jk i sin c (k i ∆ x ) -β k 2 i ∆ 2 x 2 sin c (k i ∆ x ) k (βx i ) k i + β cos(k i ∆ x ) U (βx i ) U(βx i ) U(βx i ) + O(β 2 ). ( 84
)
Collecting the terms of order β 0 in the last expression, one gets

X(βx i ) = jk i sin c (k i ∆ x ) 1 3 (2 + cos(k i ∆ x )) U(βx i ) + O(β). ( 85
)
Assuming that Equation ( 85) is valid for the continuous variable x, one derives by first applying the neperian logarithm and then by differentiating the resulting relation

X (βx i ) X(βx i ) = U (βx i ) U(βx i ) + cos(k i ∆ x ) sin c (k i ∆ x ) + k 2 i ∆ 2 x 3 sin c (k i ∆ x ) 1 3 (2 + cos(k i ∆ x )) k (βx i ) k i + O(β). ( 86
)
Inria Plugging (86) in Equation ( 84), one obtains

X(βx i ) = jk i sin c (k i ∆ x ) 1 3 (2 + cos(k i ∆ x )) U(βx i ) + β cos(k i ∆ x ) 1 3 (2 + cos(k i ∆ x )) + k 2 i ∆ 2 x 3 sin c (k i ∆ x ) 2 [ 1 3 (2 + cos(k i ∆ x ))] 2 U (βx i ) U(βx i ) + k 2 i ∆ 2 x 3 sin c (k i ∆ x )(cos(k i ∆ x ) -1) [ 1 3 (2 + cos(k i ∆ x ))] 2 + k 4 i ∆ 4 x 9 sin c (k i ∆ x ) 3 [ 1 3 (2 + cos(k i ∆ x ))] 3 k (βx i ) k i U(βx i )+O(β 2 ). ( 87
)
Following the same procedure, we express the variables Z and W in terms of U as follows

Z(βx i ) =jk i d i sin c (k i ∆ x ) 1 3 (2 + cos(k i ∆ x )) U(βx i ) + βd i βd i k 2 i ∆ 2 x 3 sin c (k i ∆ x ) 2 [ 1 3 (2 + cos(k i ∆ x ))] 2 d (βx i ) d i + cos(k i ∆ x ) 1 3 (2 + cos(k i ∆ x )) + k 2 i ∆ 2 x 3 sin c (k i ∆ x ) 2 [ 1 3 (2 + cos(k i ∆ x ))] 2 U (βx ) U(βx i ) βd i k 2 i ∆ 2 x 3 sin c (k i ∆ x )(cos(k i ∆ x ) -1) [ 1 3 (2 + cos(k i ∆ x ))] 2 + k 4 i ∆ 4 x 9 sin c (k i ∆ x ) 3 [ 1 3 (2 + cos(k i ∆ x ))] 3 k (βx i ) k i U(βx i ) + O(β 2 ), (88)
W(βx i ) = βd i 1 3 (2 + cos(k i ∆ x )) d (βx i ) d i U(βx i ). (89)
Recalling that MV = N (D Ū ), one can deduce the expression of V by substituting in (87) X by V and U by dU. Then T is computed directly by using the relation

T = D X + 1 3 (V -Z + 2W ) .
Furthermore, Y and S are obtained using the relations MY = N X and MS = N T and (87).

Coming back to (83), one obtains the expression of ω, using the expressions of the amplitudes of U , T , Y and S and collecting in the resulting equation the term of order β 0 , for all i = 0, ..., N ,

ω 2 gd i k 2 i = sin c (k i ∆ x ) 2 [ 1 3 (2 + cos(k i ∆ x ))] 2 + k 2 i d 2 i 3 sin c (k i ∆ x ) 2 . ( 90
)
Furthermore, collecting the term of order β, one obtains a relation between U, k, and the bathymetry d :

α (1) 1,i U (βx i ) U i + α (1) 2,i d (βx i ) d i + α (1) 3,i k (βx i ) k i = 0. (91)
Equation (91) describes the effect of linear shoaling since the numbers α

2,i and α

are known as the linear shoaling coefficients. Using (90), one can compute these three coefficients (we omitt the details for simplicity), for all i = 0, ..., N , α

1,i = 2 cos(k i ∆ x ) [ 1 3 (2 + cos(k i ∆ x ))] sin c (k i ∆ x ) + 2 k 2 i ∆ 2 x 3 sin c (k i ∆ x ) α (1)
,i = 2 1 + 2 cos(k i ∆ x ) 3 [ 1 3 (2 + cos(k i ∆ x ))] sin c (k i ∆ x ) - 1 + 2 cos(k i ∆ x ) 3 sin c (k i ∆ x ) [ 1 3 (2 + cos(k i ∆ x ))] k 2 i d 2 i α (1) 2
,i = cos(k i ∆ x ) 2 sin c (k i ∆ x ) 2 1 3 (2 + cos(k i ∆ x )) + 2 k 2 i ∆ 2 x 3 (cos(k i ∆ x ) -1) + k 4 i ∆ 4 x 9 sin c (k i ∆ x ) 2 [ 1 3 (2 + cos(k i ∆ x ))] + k 2 i d 2 i 3 cos(k i ∆ x ) 2 2 sin c (k i ∆ x ) [ 1 3 (2 + cos(k i ∆ x ))] 2 sin c (k i ∆ x ) -1 + 2 k 4 i ∆ 4 x 9 sin c (k i ∆ x ) 2 (1 -cos(k i ∆ x )) [ 1 3 (2 + cos(k i ∆ x ))] 2 (1) 3
.

Differentiating formally the dispersion relation (90) and assuming that ω is constant, we deduce that k i has to satisfy the following condition

k (βx i ) k i = -α (1) 4 d (βx i ) d i , (92)
where, for all i = 0, ..., N , α

k 2 i d 2 i 3 sinc(ki∆x) 2 [ 1 3 (2+cos(ki∆x))] 2 cos(ki∆x) sinc(ki∆x) + k 2 i ∆ 2 x 3 sinc(ki∆x) [ 1 3 (2+cos(ki∆x))] (1) 4,i = 1 2 1 -
.

This relation is used to compute formally k i for i = 0, ..., N and for a given bathymetry and therefore to obtain the coefficients α

2,i and α

3,i . Finally, following [START_REF] Dingemans | Water wave propagation over uneven bottoms[END_REF], we obtain an expression of the amplitude of the surface elevation A, for all i = 0, ..., N

A i A i = -α (1) s,i
d i i , (93) where α (1) s,i = α (1) 2,i -α (1) 3,i + α (1) 1,i α (1) 4,i α (1) 1,i . (94)
5.2 Linear characteristics of the classical Peregrine model.

We consider now the classical numerical model presented in Section 3. We reproduce the same procedure as in Section 5.1. The terms of order of β 0 gives the linear dispersion relation, for all i = 0...N ,

ω 2 gd i k 2 i = sin c (k i ∆ x ) 2 1 3 (2 + cos(k i ∆ x )) 1 3 (2 + cos(k i ∆ x )) + 1-cos(ki∆x) k 2 i ∆ 2 x 2 k 2 i d 2 i 3 , (95)
whereas the terms of order of β provides the linear shoaling coefficients as follows α

1,i = 2 1 3 (2 + cos(k i ∆ x )) sin c (k i ∆ x ) cos(k i ∆ x )+ 2 k 2 i d 2 i 3 sin c (k i ∆ x ) cos(k i ∆ x ) 1 -cos(k i ∆ x ) k 2 i ∆ 2 x 2 -sin c (k i ∆ x ) 3 + k 2 i ∆ 2 x 3 2 sin c (k i ∆ x ) 3 + k 2 i d 2 i 3 sin c (k i ∆ x ) 3 1 3 (2 + cos(k i ∆ x )) 1 -cos(k i ∆ x ) k 2 i ∆ 2 x 2 , α (2)
,i = 1 3 (2 + cos(k i ∆ x )) sin c (k i ∆ x ) 5 cos(k i ∆ x ) + 1 3 + k 2 i d 2 i 3 1 -cos(k i ∆ x ) k 2 i ∆ 2 x 2 sin c (k i ∆ x ) 5 cos(k i ∆ x ) + 1 3 -3 sin c (k i ∆ x ) + k 2 i ∆ 2 x 3 sin c (k i ∆ x ) 3 + k 2 i d 2 i 3 sin c (k i ∆ x ) 3 [ 1 3 (2 + cos(k i ∆ x ))] 1 -cos(k i ∆ x ) k 2 i ∆ 2 x 2 , Inria α (2) 2
,i = 1 3 (2 + cos(k i ∆ x )) cos(k i ∆ x ) 2 + k 2 i d 2 i 3 1 -cos(k i ∆ x ) k 2 i ∆ 2 x 2 cos(k i ∆ x ) 2 -sin c (k i ∆ x ) 2 cos(k i ∆ x ) + k 2 i ∆ 2 x 6 4 sin c (k i ∆ x ) 2 cos(k i ∆ x ) -3 sin c (k i ∆ x ) 2 -k 2 i d 2 i 1 -cos(k i ∆ x ) k 2 i ∆ 2 x 2 sin c (k i ∆ x ) 2 1 3 (2 + cos(k i ∆ x )) + k 4 i d 4 i 9 sin c (k i ∆ x ) 4 1 3 (2 + cos(k i ∆ x )) + k 2 i d 2 i 3 sin c (k i ∆ x ) 4 [ 1 3 (2 + cos(k i ∆ x ))] 2 1 -cos(k i ∆ x ) k 2 i ∆ 2 x 2 . (2) 3
Differentiating Equation (95), assuming that ω is constant, we obtain the relation that must satisfies each k i for all i = 0, ..., N

k i k i = -α (2) 4,i d i d i , (96)
where

α (2) 4,i = 1 2C 1 3 (2 + cos(k i ∆ x )) cos(k i ∆ x ) - k 2 i d 2 i 3 1 -cos(k i ∆ x ) cos(k i ∆ x ) k 2 i ∆ 2 x 2
, and

C = 1 3 (2 + cos(k i ∆ x )) sin c (k i ∆ x ) + k 2 i d 2 i 3 1 -cos(k i ∆ x ) sin c (k i ∆ x ) k 2 i ∆ 2 x 2 - sin c (k i ∆ x ) cos(k i ∆ x ) + k 2 i ∆ 2 x 3 cos(k i ∆ x ) sin c (k i ∆ x ) + k 2 i d 2 i 6 sin c (k i ∆ x )(1 -cos(k i ∆ x )) 1 3 (2 + cos(k i ∆ x )) k 2 i ∆ 2 x 2
.

Again, we obtain the relation between the amplitude A of the surface elevation and the bathymetry which reads, for all i = 1, ..., N ,

A i A i = -α (2) s,i d i d i , (97) where α (2)
s,i = α (2) 2,i -α (2) 3,i + α (2) 1,i α (2) 4,i α (2) 1,i
.

Analysis of the computations.

In this section, we study the linear dispersion relations derived in Section 5.1 and 5.2. More precisely, we draw the phase velocity and the amplitude of the wave with respect to the dispersion parameter σ in shoaling conditions for each scheme and we compare the results with the ones predicted by the linear theory associated to the Peregrine equations (1).

Phase velocity

For a given numerical scheme, the phase velocity is usually given by the relation

C = ω k .
As observe in literature, we can consider k and d as constant functions (k = k 0 , and d = d 0 see [START_REF] Walkley | A Numerical Method for Extended Boussinesq Shallow-Water Wave Equations[END_REF] for more details). Then, the phase velocity of our new numerical scheme (75)-( 76) is

C 2 (1) = gd sin c (k∆ x ) 2 [ 1 3 (2 + cos(k∆ x ))] 2 + k 2 d 2 3 sin c (k∆ x ) 2 , ( 98
)
while that derived for the classical scheme ( 11)-( 12) reads

C 2 (2) = gd sin c (k∆ x ) 2 1 3 (2 + cos(k∆ x )) 1 3 (2 + cos(k∆ x )) + 1-cos(k∆x) k 2 ∆ 2 x 2 k 2 d 2 3 . ( 99
)
Our aim is to plot the two curves (98)-(99) and compare with the one predicted by the linear theory :

C 2 P = gd 1 + k 2 d 2 3 . ( 100
)
We first fix the wavelength λ and we put ∆ x = λ N λ where N λ is the number of discretization points by wavelength. A direct computation gives :

k∆ x = 2π N λ .
We recall that σ = d λ , showing that kd = 2πσ. In Figure 2, we draw the relative errors between (98) and (99) and the phase velocity predicted by the linear theory. The error er is defined, for each scheme, by er = 100 C -C P C P . In Figure 2, one can observe that for N λ = 5, the error coming from the new scheme is acceptable (less than 1.6%) whereas the one of the classical scheme is greater than 5% for Inria depths bigger than 0, 3. For N λ = 10, although the error of the classical scheme is better than in the previous case (less than 2%), the one of the new scheme is much better and stay very close to 0. We conclude here that our new numerical scheme seems to reproduce much better the linear dispersive effects.

Linear shoaling test

We first recall the expression of the shoaling coefficients given by the linear theory associated with the Peregrine equations (1) (see [START_REF] Dingemans | Water wave propagation over uneven bottoms[END_REF]):

α 1 = 2, α 2 = 2 - k 2 d 2 3 , α 3 = 1, α 4 = 1 2 1 - k 2 d 2 3 .
Our aim is to compare, for a given situation, the evolution of the amplitude of the waves with respect to the space variable x given by the two relations ( 93) and ( 97) and the one derived from the linear theory. We consider a periodic wave with an initial amplitude a = 0.05 and a wavelength λ = 15 starting from the position x = 100. It propagates over an initial constant water depth d 0 = 13. The bottom is flat until x = 150 and it has a constant up-slope of 1 50 from x = 150 to x = 790. We compute the evolution of the wave amplitude with respect to x. For that, we propose the following procedure. Firstly, we integrate formally the relation between k and d ((92) for the new scheme and (96) for the classical one), given by the differential equation

k k = -α 4 d d .
We use a Strongly Stability-Preserving Runge-Kutta method (SSP-RK) to compute the solution k, using k 0 = 2π λ as initial condition. Then, we substitute this function k in the expression of α s , and we compute the amplitude of the wave A by integrating the relation ((93) for the new scheme and (97) for the classical one)

A A = -α s d d .
Again, we use a SSP-RK method and A(100) = 0.05 as initial condition. Fixing the value of ∆ x , it is possible to compute formally the surface elevation amplitude for each scheme.

The results are presented in Figure 3. We observe that when ∆ x is small, the curves of the two schemes matched the theoretical one, meaning that both schemes converges. However, when ∆ x becomes larger, one can see that the curve computed with the new scheme stay closed to the theoretical one while the one computed with the classical scheme underestimates drastically the amplitude of the wave.

Numerical experiments

This section is devoted to the investigation of the nonlinear behavioir of the two schemes (73)-( 74) and ( 11)- [START_REF] Lannes | Derivation of asymptotic two-dimensional timedependent equations for surface water wave propagation[END_REF]. To this end, we present two different test cases : the propagation of solitary waves over a flat bathymetry and a shoaling test for which some solitary wave profile is computed in a situation close to wave breaking (see [START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF] for details). These test cases bring to the fore major differences between the two numerical models and confirm the premiminary results of Section 5.

Inria

Soliton propagation

We first consider the propagation of an exact solitary wave solution to the Peregrine equations, with an amplitude equal to 0, 2 (details on the computation of this solution as well as mathematical conditions for the existence are given in [START_REF] Bellec | On the existence of solitary waves for Boussinesq type equations and a new conservative model[END_REF]) over a flat bathymetry d 0 = 1. The space intervall is equal to [0, 100]. In order to check our implementations, we have performed a grid convergence analysis. Numerical results have been compared with the initial profile (which is the profile of the exact solution). The meshes used contain respectively 1000, 2000, 4000 and 8000 points. In Figure 4, we have plotted the L 2 -norm of the rescaled error for each scheme. The slope obtained for the scheme (73)-(74) shows a convergence of order 2 while that corresponding to ( 11)-( 12) is equal to 1.7. We deduce that with the same initial finite elements method, the new procedure improves the order of convergence.

Shoaling test

This test case is dedicated to the wave shoaling characteristics of the two models considered here in conditions close to wave breaking. It has been initially proposed by Grilli et al. in [START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF] and consists in propagating a solitary wave of amplitude 0.2 m on a water depth d 0 = 1 m during a time 13 s. The shoaling takes place onto a constant slope of 1 35 (see Figure 5). Note that, to our knowledge, it is not possible to give an analytical solution in this configuration. We then decide to compute a reference solution using a very refined mesh of 8000 points and the scheme ( 11)-( 12). This solution is used as a standard in the sequel to make the comparisons.

In Figure 6, we have plotted the shoaling wave profile using succesively a mesh containing 250 points (left side) and 500 points (right side) for the two numerical schemes as well as the reference solution. Clearly, one can observe a major difference in the behaviour of the two schemes. For a mesh with 250 points (∆ x = 0, 3), the solution emanating from the new scheme (73)-(74) matches very well the reference curve, showing that the convergence as already occurs with a few numbers of points while it is obviously not the case for the scheme (11)- [START_REF] Lannes | Derivation of asymptotic two-dimensional timedependent equations for surface water wave propagation[END_REF]. Indeed, the green curve exhibit some amplitude and phase defects. To push on the analysis, we repeat the same computations with a grid mesh containing 500 points (∆ x = 0, 15). Again, we can see that the solution computed with the scheme (73)-(74) coincides perfectly with the reference curve. However, in this case, one can observe that the classical scheme provides better results without reaching the precision of the other scheme. Note that the conclusion doesn't changed if one computes the reference solution with the scheme (73)-(74). In addition, the Table 7 contains the relative errors in L 2 and L ∞ norms for each scheme, where the errors are defined by It transpires that the new scheme gives results very close to the exact solution, whereas the classical one need more points to give similar approximations. Finally, to compare the accuracy of the two models, we have performed a grid convergence analysis. In Figure 8, we have plotted the error defined in (101) in the L 2 -norm for each scheme, corresponding to succesively 500, 1000 and 2000 points. The slopes obtained for the error are similar for the two schemes. However, it is clear that the new model gives better results with 500 points that the classical scheme with 2000 points. We can deduce from this analysis that the new numerical scheme better reproduce physical effects and dispersion characteristics.

E 2 = ||E -E ref || 2 ||E ref || 2 , E ∞ = ||E|| ∞ -||E ref || ∞ ||E ref || ∞ , (101)

Conclusions and perspectives

We have presented a new systematic method to obtain discrete numerical model in the study of incompressible free surface flows. In order to evaluate the power of this method, we have considered the case of the so-called Peregrine equations and performed the computations in this academic situation. We have compared our new numerical scheme with the one obtained by performing directly a Galerkin method on the Peregrine equations.

Finally, by the use of several numerical experiments, we have shown the efficiency of our new scheme.

Moreover, we claim that it is possible to obtain other different numerical scheme, consistent with the Peregrine equations, by changing the initial discretization of the Euler equations. For example, by choosing the following discrete boundary conditions :

• At the free surface :

Ŵ = E t + K(E Û ) -E (K Û ).
• On the seafloor :

W = -K(D Ǔ ) -D (K Ǔ ) ,
we are able to derive a new numerical scheme in a conservative form

E t + K(H Ū ) = 0, Ūt + K Ū 2 -Ū (K Ū ) + gKE + D 2 6 (K 2 Ūt ) - D 2 K 2 (D Ūt ) = 0.
In the future, we plan to apply this new procedure to derive numerical schemes for Extended Boussinesq's models as well as for the Green-Naghdi equations. By this procedure, one of our goal is, for example, to enhance the linear dispersion characteristics of these models.

Figure 1 :

 1 Figure 1: Sketch of the free surface flow problem, main parameters description.

Figure 2 :

 2 Figure 2: Comparison of the phase velocity (N λ = 5 on the left, N λ = 10 on the right) of the classical and the new numerical scheme with the one given by the linear theory w.r.t σ.

Figure 3 :

 3 Figure 3: Evolution of the wave amplitude for the two numerical schemes and for Peregrine equations. Top left : ∆ x = 0, 5 (N λ = 30). Top right : ∆ x = 1, 5 (N λ = 10). Down : ∆ x = 3 (N λ = 5)

Figure 4 :

 4 Figure 4: Grid convergence results for the two numerical schemes.

Figure 5 :

 5 Figure 5: Shoaling of a solitary wave; computational configuration (Grilli and al ).

Figure 6 :

 6 Figure 6: Shoaling wave profiles of Peregrine schemes. Left : mesh with 250 points. Right : mesh with 500 points

Figure 7 :

 7 Figure 7: Relatives errors between the two numerical schemes and the reference solution for meshes with 250, 500,1000 and 2000 points

Figure 8 :

 8 Figure 8: Grid convergence results for the shoaling of the two numerical schemes.

  E ref being the reference solution and E the one of the two numerical schemes. The computations have been made with respectively 250, 500, 1000 and 2000 points.

	Relative error\ Points per meshes	250	500	1000	2000
	E 2 for the New scheme	1.55%	0.34%	0.09%	0.03%
	E 2 for the Classical scheme	31.25%	12.67%	3.70%	0.81%
	E ∞ for the New scheme	6.01%,	1.03%	0.07%	0.04%
	E ∞ for the Classical scheme	20.19%	9.67%	3.31%	0.85%
					Inria

RR n°8806

351, Cours de la Libération Bâtiment A 29

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgment

This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx Bordeaux -CPU (ANR-10-IDEX-03-02). Partial funding by the TANDEM contract (reference ANR-11-RSNR-0023-01 French Programme Investissements d'Avenir ) is also acknowledged.


Baptiste Mélès

Systems and Teaching in Stoic and Confucian Philosophies

published or not. The documents may come

L'archive ouverte

Introduction *

Systems have a lot of virtues. Among their abilities, the theoretical ones are often underlined. First, the system, defined as a rational organization of philosophical discourse, gives consistency: it assures that there is no contradiction between the propositions constituting the theory. It also offers unity: every part of the theory is bound to a single principle. It eventually allows for completeness: the author as well as the reader are sure that no proposition is missing, since each one stays exactly at its place, leaving no lacuna. Systematicity thus guarantees scientificity, which is the highest theoretical virtue for which one can hope, even when systematicity is not claimed and remains implicit.

But one is led to examine whether systems may also have some teaching virtues, especially in moral philosophy -thus not only theoretical virtues, but also practical ones. Are systems appropriate ways of teaching the path to wisdom? To consider this question, two philosophical schools will be compared: stoicism in the Western tradition and confucianism in the ancient Chinese tradition.

1 Systematicity as the royal way to wisdom It has often been said that the Stoics invented systematicity. Without, of course, being the first philosophers to try and develop a consistent theory, they are said to have paid particular attention to dividing philosophy in a few domains -ethics, physics and logic -and showing the links between those fields. Long and Sedley thus write that "Of all ancient philosophies, Stoicism makes the greatest claim to being utterly systematic. Arguably, the Stoics invented the notion of philosophy as "system 1 " [...]." I do not intend to examine whether this precedence is historically valid, but I will try and show which systematicity is engaged in Stoic philosophy.

Indeed the risk is to interpret retrospectively the ancient notion of system, keeping in mind the later fortune of this concept, and especially Kant's or Hegel's views about, and use of, systems. That Stoics have promoted systems for their theoretical virtues, above all for consistency and completeness, is beyond doubt, and -rightly -often emphasized. Victor Goldschmidt has shown how the Stoic theory of time not only belongs to physics, but is also correlated to ethics and even logic. For instance, physics shows that only the present exists, and that past and future do not exist in a strong sense; ethics thus claim that I have to keep my mind on the present, preventing myself from any remorse or hope2 . Ethics plays in this insightful analysis the role of an object for the system -time indeed has an ethical meaning -, but the system as educational method does not explicitly seem to possess any ethical virtue: only theoretical ones.

Less often mentioned are the educational virtues systems had for the Stoics, however important it may have been for them. The goal of Stoic philosophy is not exclusively theoretical. Wisdom is aimed beyond science, as Seneca writes: "I will state the difference between wisdom and philosophy. Wisdom is the human mind's good brought to perfection. Philosophy is the love and pursuit of wisdom; it strives for the goal which wisdom has achieved 3 ." And the fact that philosophy intends to lead to wisdom is corroborated by the ternary division of virtue (ἁρετή, arétè), i.e. excellence, reported by Aetius: "The Stoics said that wisdom is scientific knowledge of the divine and the human, and that philosophy is the practice of expertise (τέκνη, techne) in utility. Virtue singly and at its highest are triple -the physical one, the ethical one, and the logical one. For this reason philosophy also has three parts -physics, ethics and logic 4 ." Philosophy leads to wisdom because they share the same structure. There is more: wisdom preserves, at the highest level, the content of philosophy.

Even the order of the parts of philosophy may have an educational meaning. While it is an order of importance, it also refers to an order of teaching, showing what the last goals of philosophical learning are. Several different orders have been proposed: some Stoics "assign the first place to logic, the second to physics and the third to ethics; these include Zeno in his book On Discourse, Chrysippus, Archedemus and Eudromus. Diogenes of Ptolemais starts with ethics, Apollodores puts ethics second, while Panaetius and Posidonius start with physics 5 ." This order may also appear as an order of exhibition, if we presume that the most important parts of philosophy are treated at the end of the curriculum, insofar as they require the knowledge of the other parts.

Emphasizing the theoretical virtues of systems in Stoic philosophyconsistency, unity and completeness -as it has often, and rightly, been done, should therefore not hide some of its practical virtues recognized and promoted by the Stoics. Systematicity in itself does not directly lead to wisdom, but is the best way to express philosophy, which in turn leads to wisdom, once excellence is reached. It is not enough to say that systematicity does not obstruct the quest for wisdom: it is, in fact, the royal way that leads to it.

2 Systematicity as obstructing the Way to wisdom A glimpse at ancient Chinese philosophical texts suffices to establish that few are systematic. Why is that?

For theoretical reasons, one might think, concerning ancient Chinese cosmology, as it was developped since the Yìj īng (易經), which describes the world in an ever-changing state6 . For theoretical reasons too, if all phenomena are relative, as in Zhuāngzi (莊子, Tchouang-Tseu), or if language is not trustworthy, as in Lǎozi (老子, Lao-Tseu) 7 . With regard to these ancient Chinese philosophies, the lack of systematicity is not surprising.

More unexpected is the fact that even some relatively systematic philosophies can refuse the systematic way of teaching. It is mainly the case of Confucius' thought 8 . The comparaison may be all the more instructive as Confucian philosophy has often been compared with Stoicism 9 .

Of course, the word "system" is nowhere to be found in Confucius' work, which is not expressed in a systematic way, consisting in a set of words or acts of the Master related by his disciples, similarly to Epictetus' Discourses; but this did not prevent Confucius from expressing himself about notions which constitute the concept of system in the Western tradition.

"The way that can be spoken of Is not the constant way; The name that can be named Is not the constant name." See Laozi, Tao Te Ching (道得經, Dàodéjīng), Penguin Classics, translated by D. C. Lau, London: Penguin Books, 2006, 1: 《道可道非常道,名可名非常名》.

8 Kǒngfūzǐ (孔夫子). 9 Both link cosmology and ethics, cosmological and social order. The notion of duty (κατήκοντα or officia in Stoic philosophy) is not universel, but depends from the social position. This last point is exactly what distinguishes Mencius' (孟子, Mèngzi) orthodox confucianism from Mòzi's (墨子) thought.

Consistency is the minimal condition for a theory to be a system; it can without doubt be expected from Confucius. But there is more: even unitythe fact that several theories are correlated to a single principle -is claimed. Confucius strongly asserts this last criterion of "systematicity:" "There is a single thread binding (貫, guàn) my way (道, dào) together10 ." Further he also claims, in almost the same words: "I have a single thread binding it all together11 ." Thus Confucius does not deny, nor minimize, the consistency and unity of his thought: he emphasizes it.

Concerning completeness, some passages in the Analects, especially the sixteenth chapter, consist in a set of lists: three kinds of beneficial friendship, of harmful friendship, of beneficial pleasures, of harmful pleasures, etc. But these lists have no regularity, no rule of construction; they are very rare in Confucius' work; and above all, their authenticity is very doubtful. Confucius thus emphasizes unity but avoids completeness.

Confucius' relation to the form of his moral philosophy is itself regulated by educational interest. The way in which thoughts are exhibited is not neutral with respect to the content of the lessons. Confucian philosophy may be a system, but is not exposed as such. Keeping silent about what makes from the system a totality allows for the revelation of disciples' qualities 12 . The good student is the one who can develop by himself what follows from the thoughts of his master: "Ah, dear Zígòng (子貢)! I eventually can talk to you about the Book of Odes (詩經, Sh īj īng)! Tell such a man something and he can see its relevance to what he has not been told 13 ." This criterion, which characterizes good students, also enables the master to identify bad ones: "When I have pointed out one corner of a square to anyone and he does not come back with the other three, I will not point it out to him a second time 14 ." Thus the real wise is not the one who has reached excellence in all the fields of philosophy, as in the Stoic notion of wisdom; he is the one who can see much further than what the master says. As Confucius compares two of his disciples, Yán Huí (顏回) 15 and Zígòng, the latter tells him: "How dare I compare myself with Hui? When he is told one thing he understands ten. When I am told one thing I understand only two 16 !" Refusing any systematic discourse, any claim to sufficiency, is therefore a way to let the disciples become intellectually autonomous. The most important quality is not to know, be it at the highest level of excellence, but to study (學, xué) [START_REF]To see the importance of "study" (學, xué), first caracter of the Analects (I, 1)[END_REF] . Completeness is educationally dangerous: it prevents the students, but also the master, from studying.

Conclusion

Thus a difference emerges between two educational strategies. The Stoics choose systems in order to promote an idea of excellence, which implies an excellence of knowledge. Confucius, who avoids systems, reveals and shapes the understanding abilities of his disciples.

This conflict of strategies finds an analogon in mathematics, the theoretical science par excellence. While Euclides develops his mathematics in a systematic way, and tends to generality, the Nine Chapters, a classical book of Chinese mathematics, only exhibits particular cases. Not by inability, as Karine Chemla has shown [START_REF] Dunod | Les Neuf Chapitres. Le classique mathématique de la Chine ancienne et ses commentaires[END_REF] , but by choice: the reader has to find truths by himself. This confrontation reveals that systematicity does not only raise theoretical problems, but also educational ones.

A. A.Long and D. N. Sedley, The Hellenistic Philosophers, vol. 1, Cambridge: Cambridge University Press, 1987, p. 160.

Victor Goldschmidt, Le Système stoïcien et l'idée de temps, Collège philosophique, Paris:Vrin, 1953, ch.

3, iii, p. 73. 3 Seneca, Letters, 89,[4][5] op. cit., 26 G, p. 160. 4 Aetius I, Preface 2, SVF II, 35. See ibid., 26 A, p. 158. 5 Diogenes LaërtiusVII,[39][40][41] 26 B, p. 158.

Although we could have expected a systematic treatment from this book, which is structured by regular combinations of plain lines "-" (陰, yīn) and broken ones "--" (陽, yáng).

The Book of Way and Virtue begins with the words:

Confucius, The Analects ( 論語, Lúnyǔ), Penguin Classics, translated by D. C. Lau, London: Penguin Books, 1979, IV, 15: 《吾道一以貫之》.

Ibid., XV, 3: 《子一以貫之》.

A deep analysis of ancient Chinese teaching strategy, especially in Confucian philosophy, can be found in François Jullien, Le Détour et l'accès. Stratégies du sens en Chine, en Grèce, Paris: Grasset, 1995, ch. ix. In the same book, François Jullien highlights the theoretical reasons of this educational strategy and compares it, among other things, to military strategy.

 13 Confucius, op. cit., I, 15. See also III, 8 for an example of a good disciple, understanding and pursuing the master's thought, and therefore able to talk with him from the

Shījīng 14 Ibid

., VII, 8.15 Yán Huí, who died precociously, was Confucius' favorite disciple; see ibid., V, 25; VI, 2, 5 and 9; VII, 10;IX, 10, 19 and 20; XI, 2,[6][7][8][9][10] 22; XII, 1; and probably, following the tradition, VIII, 5.