Simulation of Complex Free Surface Flows

. We study the Serre–Green-Naghdi system under a non-hydrostatic formulation, modelling incompressible free surface ﬂows in shallow water regimes. This system, unlike the well-known (non-linear) Saint-Venant equations, takes into account the eﬀects of the non-hydrostatic pressure term as well as dispersive phenomena. Two numerical schemes are designed, based on a ﬁnite volume – ﬁnite diﬀerence type splitting scheme and iterative correction algorithms. The methods are compared by means of simulations concerning the propagation of solitary wave solutions. The model is also assessed with experimental data concerning the Favre secondary wave experiments [12]. Résumé. L’objet de l’étude est le système de Serre–Green-Naghdi aﬁn de modéliser des écoulements incompressibles à surface libre en eaux peu profondes, tout en tenant compte les eﬀets de la pression non-hydrostatique et les phénomènes dispersifs qui jouent un rôle important en particulier en océanographie côtière. Plusieurs schémas numériques sont présentés, basés sur une méthode de splitting, en séparant la partie hyperbolique et la partie non-hydrostatique. Les algorithmes sont comparés à travers des simulations avec des ondes solitaires. Le modèle est également validé avec des données expérimentales concernant les ondes de Favre [12].


Introduction
The simulation of incompressible, homogeneous, inviscid free surface fluid flows as a tool to describe and examine the propagation of surface water waves has been a major field of research oriented towards oceanographic applications for the past several decades.In general, such a physical system is governed by the free surface Euler equations, which, to this day, remains an analytically and numerically challenging topic.
For coastal, or large scale applications, it is commonly approximated by reduced-complexity models based on the shallowness of the fluid domain.Diverse models based on the (nonlinear) shallow water [or Saint Venant] system [14] have been proposed and analysed in the past (see for instance [4,5,8,11,23,24,26]).Although this model provides a reasonable tool for capturing the nonlinear transformation of waves, it fails to represent dispersive effects because of the underlying hydrostatic assumption.
The model investigated in the present article is the Serre-Green-Naghdi system (originally appeared in [18,28]), also known as the fully nonlinear Boussinesq system and referred to as LDNH 2 (1) in [17].This model is derived from the free surface Euler equations and reads [16,25] where x = x, u = u if d = 1 or x = (x, y), u = (u, v) if d = 2.The equations are set in a smooth bounded domain Ω ⊂ R d .
Unknowns.Here h : R + × R d → R + is the fluid height; u : R + × R d → R d is the vertically averaged horizontal velocity of the fluid, w : R + × R d → R is the vertically averaged vertical velocity and σ : R + × R d → R is the vertical correction to w. Finally q : R + × R d → R is the vertically averaged hydrodynamic pressure and q b : R + × R d → R is the hydrodynamic pressure at the bottom.
Data. z b : Ω → R is the time-independent bottom topography such that where n is the outward normal unit vector and g = 9.81 is the standard acceleration due to gravity.The system is complemented with multiple possible boundary conditions, depending on the physical setting modeled.
Boundary conditions.The boundary conditions applied to velocity unknowns are detailed in the numerical part of the paper.As for the pressure, we impose for some partition of ∂Ω = Γ in ∪ Γ out to be specified later.
The main purpose of the paper is to analyze some of the properties of the system in any dimension and to examine some 1D numerical algorithms.Developing robust and efficient algorithms for the numerical resolution of the Serre-Green-Naghdi system has drawn much attention lately, partially due to the knowledge from theoretical analyses and the surge in computational power (see for instance [6,13,15,[20][21][22]25]).
Our approach relies on a finite volume based splitting scheme that has shown promising results [1,2] when applied to the non-hydrostatic shallow water model introduced in [10].
The outline of the document is as follows.In Section 1, we present some properties of the model (equivalent formulations, energy balance, . . . ) and we prove that the corresponding projection step is well-posed.
In Section 2, we elaborate the finite volume -finite difference type splitting scheme.It consists in solving the hyperbolic part first, then the dispersive part.The hyperbolic prediction step involves the numerical resolution of the classic shallow water equations with a non-flat bottom topography.As for the resolution of the system arising from the non-hydrostatic correction step -which is the core of this work -we present two iterative algorithms: the Uzawa algorithm applied to the mixed velocity-pressure problem and the Gauss-Seidel method for the projection elliptic step.The different choices for the boundary conditions depending on the physical situation at hand, as well as their consequences on the discretization procedures, are also presented.
In Section 3, we describe some numerical test cases in order to assess the aforementioned algorithms.Simulations concerning the propagation of exact solitary waves over a flat bottom topography are run in order to compare the two algorithms, mainly in terms of computational time.It is followed by simulations of a dispersive dam-break problem.Finally, the numerical validation against the experimental data obtained from the Favre secondary wave experiment [12] is presented and compared with previous results [9].

Specific operators
Given the vectors we define the dispersive divergence and gradient operators and the source term vector Hence System (1) reads under the compact non-conservative form Lemma 1.1.The dispersive operators defined by (5) satisfy the Green-type formula:

Reformulations of the system
Let us first mention why we refer to System (6) as the Serre-Green-Naghdi model.Indeed: Remark 1.2.System (6) can be rewritten under the Boussinesq formulation where we denote for d = 1 and for d = 2 Here we made use of the notation System (1) has a similar structure to pressure-based problems: Proposition 1.3.For smooth solutions, Equations (6) can be rewritten as Remark 1.4.Equation (8c) expands as

Energy
The d-dimensional Serre-Green-Naghdi equations ( 6) admit an energy balance identity.
Proposition 1.5.Smooth solutions to System (6) satisfy the following identity: Proof.Let us multiply (1b) by u, (1c) by w and (1d) by σ, and sum the results.This leads to By making use of Equation (1a), of the duality relation (7) and since X satisfies Equations (1e-1f), we obtain the expected result.

Variational formulation
Let us focus in this part on the well-posedness of "elliptic" problems of type where operators are defined in (5), f is a given source term assumed in L 2 (Ω) and h is given under the hypothesis System ( 9) is supplemented with boundary conditions (3)(4).This kind of problem occurs at the continuous level in (8) and in the numerical strategy (see (14) below).
Let us define the following functional space.
Definition 1.6.For h ∈ H 1 (Ω), let us define the following function space endowed with the inner product and the corresponding norm Due to the duality relation (7) and the boundary conditions (3)(4), the variational problem associated to (9) reads: More precisely, The Cauchy-Schwarz inequality directly shows that A is continuous over H 1 h (Ω) × L 2 (Ω).Moreover: Proposition 1.8.Under Hypothesis (H) and for z b ∈ C 1 (Ω) satisfying (2), the operator A is coercive over Consequently, the variational formulation (10) is well-posed.

The splitting scheme
Let us take n ∈ N and discretize the system in time.Given a series of timesteps ∆t k , k ∈ N (defined later by an appropriate CFL condition) giving rise to the times t n = k≤n ∆t k , one has that for the step n + 1, the system can be split into two parts at the semi-discrete (in time) level: (1) A prediction step for the hyperbolic part of the system: for variables h and hX.Notice that this system is essentially the shallow water equations with uneven bottom topography: its resolution is classic and well-known.In the numerical simulations to follow (see Section 3), we chose a first order finite volume approach, a Rusanov solver with hydrostatic reconstruction for the scheme to be well-balanced [3].(2) A projection / correction step involving the non-hydrostatic (dispersive) part of the system: for variables X and Q.The novelty of the problem lies in the efficient resolution of the latter nonlinear problem.In what follows, we detail two techniques of discretization, both of them make use of some particular (structural) properties of the system and both of them are iterative methods for computational efficiency reasons.Notice that Equations (13b-13c) imply In the collocated framework, all variables are defined at the cell centers x i (i ∈ {0, . . ., N + 1}) with x 0 and x N +1 corresponding to ghost cells at the boundaries (see Figure 1).The use of a collocated scheme (especially for the pressure variables (q, q b )) is justified by the special form of the dispersive gradient ∇ SGN Q.
Hyperbolic step: System (12) consists of the classic 1D shallow water equations for (h, hu) and two transport equations for (hw, hσ).As it was already mentioned, a first order finite volume scheme was implemented: where X * = (1, X T ) T .The numerical fluxes F n i+1/2 (i ∈ {0, . . ., N }) approximate the physical fluxes at the interfaces; for the simulations, a Rusanov solver was chosen as an approximation.
Such a numerical scheme is stable, provided the following Courant-Friedrichs-Lewy (CFL) condition holds: , where c F is a constant depending on the chosen numerical solver.For the treatment of boundary conditions, we refer to Section 2.5.
Velocity-pressure system: System ( 13) is a Stokes-type system for which the water height is given: Let us denote by H n+1 ∈ M N,N (R) the diagonal N × N matrix with the corresponding diagonal elements h n+1 i , and by H n+1 ∈ M 3N,3N (R) the 3N × 3N block-diagonal matrix with block entries H n+1 .
The dispersive gradient ∇ SGN Q n+1 is discretized as BQ n+1 , as part of a finite difference type discretization scheme on the mesh of the interval I. Here, and in what follows throughout this section, boldfaced vectors refer to the discretized values over all cells (i ∈ {1, . . ., N }), such as for instance The matrix B ∈ M 3N,2N (R) has a relatively simple structure: where the matrix ) is an appropriate discretization for the term ∂ x (hq).We consider two different finite difference schemes for this (1) Second order central scheme on mesh points: for any interior cell (i ∈ {2, . . ., N − 1}) we have that Here, cell boundaries are x i+1/2 (i ∈ {0, . . ., N }) and the corresponding indices refer to the averaged values For this discretization, the corresponding tridiagonal matrix as the i th row.(2) Second order central scheme on adjacent cells: for any interior cell (i ∈ {2, . . ., N − 1}) we have that For this discretization, the corresponding tridiagonal matrix as the i th row.
When not specified, matrices B 11 and B can refer to either of these discretization schemes.Notice that the first and last rows of B 11 may contain only the truncated discretized formulae; the missing terms in the discretization correspond to information on the two ghost cells.With the establishment of the boundary conditions, the values on the ghost cells call for corrections on the first and last elements of the diagonals, up to an error term 0. The exact corrections to be made depend on the boundary conditions chosen, they will be specified in Section 2.5.
With this at our disposal, Equation (13b) is discretized as where 0 is the error terms (independent from the variables X and Q) due to the correction associated to the boundary conditions.
Owing to the duality relation (5), Equation (13c) admits a natural discretization: B T X n+1 = 0 since B T is a suitable discretization of the dispersive divergence, up to a boundary term correction induced error term of 0 (independent from variables X and Q).
To sum up, we are left with the resolution of the following linear system: In the next part, we detail two iterative algorithms for the resolution of this system, both of them using the symmetric structure of (17).We chose the iterative approach to circumvent solving directly the system of 5N equations.

Iterative methods: the Uzawa algorithm
Given the symmetric nature of System ( 17), a natural alternative to directly solving the mixed velocitypressure problem is to consider the Uzawa algorithm.This consists in the following consecutive iterations: where α is an iterative parameter and p ∈ N is the current iteration index.We remark that the matrix H is diagonal, its inversion induces no additional computational cost.As the initialization of the iterative process, Q n+1,0 = Q n is taken (since Q is not involved in the prediction step of the algorithm).
One additional concern about the initialization has to be adressed, which is the very first time step.Contrary to the variables in X, which are natural, physically easily measurable quantities, a priori knowledge on the hydrodynamic pressure components Q is rarely available.Throughout the implementation of the algorithm, an (arbitrary) initial quess of 0 was considered for these variables.This can obviously result in a high number of iterations during the first timestep, which can be spread out on multiple timesteps by imposing a maximal number of iterations p ≤ p max (with p max being constant).
Apart from the upper bound on the iterations, the Uzawa algorithm is run until the error on the correction is reasonably small, namely an error threshold of .The error itself is measured as the discrete L ∞ -norm of the difference between two consecutive iterations, that is: Notice that System (18) can be simplified into an iterative algorithm on Q only: which gives direct information on the convergence of the Uzawa algorithm [19] as a function of the eigenvalues Lemma 2.1.The algorithm converges iff the parameter α satisfies where Λ max is the largest eigenvalue of the aforementioned matrix ∆t n+1 C. The optimal iterative parameter is given by Λ min being the smallest eigenvalue of ∆t n+1 C.
However, obtaining information on the eigenvalues of this matrix is highly nontrivial, except for very particular cases (for example a lake at rest situation with flat bottom topography).

Iterative methods: the Gauss-Seidel approach
Multiplying the first equation of ( 17) by ∆t n+1 B T (H) −1 and substituting the second equation in it yields an equation similar to (19), with f = (f T q , f T q b ) T as the right hand side.Remark 2.2.In fact, this corresponds to a direct discretization of System (14).
The symmetric matrix C has a particular block structure of the following form with C 11 and C 22 being symmetric and positive definite N × N matrices.
Direct resolution of (20), much like the full velocity pressure system (17), may be time consuming, even though the matrix possesses a pleasant "quasi-diagonal" blockwise structure.Since C 22 is, in fact, a diagonal matrix, one could once again implement an Uzawa type iterative algorithm for the resolution of the discretized elliptic problem.However, we present instead a block Gauss-Seidel type method.This basically entails in performing a simple block-LU decomposition to separate C 12 for the explicit part of the algorithm and keeping the rest as an implicit part.For the iterative step k ∈ N we have the following two subsequent iterations to perform: up to a maximal iteration of p ≤ p max , or sufficiently small error term, the error being measured by the norm of the relative difference between two subsequent iterations, exactly the same way as in the Uzawa algorithm.
Once again, we initialize by taking q n+1,0 b = q n b .Owing to the positiveness of matrices C 11 and C 22 [27], one can show that Lemma 2.3.The block Gauss-Seidel type algorithm converges (with respect to the norm induced by the matrix C).

Boundary conditions
As far as boundary conditions are concerned, the three important variables are h, u (or equivalently the flux hu) and q (or its conservative counterpart hq).In what follows, we analyze different boundary conditions: the solid wall boundaries (mixed Dirichlet and Neumann boundary conditions) and a fluvial type inflow-outflow setting.We remark that boundary conditions on the boundary pressure term q b have no actual influence on the algorithm since boundary terms for q b do not appear in the scheme at any given moment.
In what follows we specify not only the appropriate boundary conditions but the calculations concerning the related error terms for both discretization schemes: ( ) the central scheme on mesh points and ( ) the central scheme on adjacent cells.

Solid wall boundaries
We assume that the boundary is a mirror, meaning that the spatial derivative (normal derivative) of quantities vanishes, except for the horizontal velocity, which simply vanishes at the boundary.In particular, on the left hand side boundary, we have that (∂ x h) i=1/2 = 0, u i=1/2 = 0, (∂ x w) i=1/2 = 0, and (∂ x σ) i=1/2 = 0, that is, on the ghost cells i = 0 these quantities are defined through basic discretization as For the Dirichlet boundary condition on the velocity u, we impose Neumann boundary condition on the pressure q, as well as on q b , therefore Similarly, one has on the right hand side boundary: and q bN +1 = q bN .( ) With the boundary conditions at our disposal, we can write the discretization of the momentum equation for the first cell: (h 3/2 (q 1 + q 2 ) − h 1/2 (q 0 + q 1 )) + q b1 (∂ x z b ) 1 = 0.
Since q 0 = q 1 , this means that and that no correction is required, 0 1 = 0. So, given the matrix B 1 , we can compute the approximation for the divergence Given that u 1 = −u 0 , we obtain that which is exactly a first order discretization of −(h∂ x u + 2 √ 3σ).Therefore no corrections are necessary either; 0 1 = 0.
With the same reasoning, we have that and that corrections are not necessary neither for the gradient nor for the divergence: 0 N = 0, 0 N = 0.
( ) With the boundary conditions at our disposal, we can write the discretization of the momentum equation for the first cell: Since q 0 = q 1 , and that h 0 = h 1 , this means that and that no correction is required, 0 1 = 0. So, given the matrix B 2 , we can compute the approximation for the divergence Given that u 1 = −u 0 we obtain that which is exactly a first order discretization of −(h∂ x u + 2 √ 3σ).Therefore no corrections are necessary either; 0 1 = 0.
With the same reasoning, we have that and that corrections are not necessary neither for the gradient nor for the divergence: 0 N = 0, 0 N = 0.

Fluvial inflow-outflow
For this case, we suppose that we are in a fluvial regime, meaning that |u| < √ gh.In this case, it is enough to prescribe either the fluid height or the flow flux at each of the boundaries.In our case we will prescribe a constant inflow flux Q inc > 0 on the left hand side boundary and a constant height H out on the right hand side boundary.The missing variables (h or hu respectively) can be recovered due to the underlying hyperbolic problem and the conservation of the associated Riemann-invariant on the outgoing characteristic.

Incoming flux on the left hand side:
We prescribe a Dirichlet type condition on the flux, that is (hu Since we are in a fluvial regime, the outgoing characteristic u+ √ gh is always positive, therefore the corresponding Riemann-invariant is preserved.This yields a second equation on h 0 and u 0 : Expressing u 0 from one equation and substituting it in the other one leads to a third order polynomial in Given that it could be expensive to solve this nonlinear equation, one can also take an extrapolation for h 0 , namely Once h 0 is determined, one can express u 0 from the inflow flux condition, so we have that For the vertical compontents of the velocity, we impose Neumann boundary condition, just as before, hence Since we imposed a Dirichlet type boundary condition on the velocity u, we prescribe a Neumann type boundary condition on the conservative pressure term hq and q b : ( ) With the boundary conditions at our disposal, we can write the discretization of the momentum equation for the first cell: Since q 0 = h1 h0 q 1 , means that and that no correction is required, 0 1 = 0. So, given the matrix B 1 , we can compute the approximation for the divergence Given that u 1 = 2Qinc−h0u0 h1 , we obtain that which is a first order discretization of −(h∂ It is at this point that we need that h 0 is not a function of u 1 .( ) Here we only treat the special case when h 0 = h 1 .With the boundary conditions at our disposal, we can write the discretization of the momentum equation for the first cell: Since q 0 = q 1 , and our additional hypothesis on h 0 , this means that and that no correction is required, 0 1 = 0. So, given the matrix B 2 , we can compute the approximation for the divergence Given that u 1 = 2Qinc−h0u0 h1 and that h 0 = h 1 , we obtain that which is a first order discretization of −(h∂ x u + 2 √ 3σ) up to the error term 0 1 = h1Qinc ∆x .Remark 2.4.Here we consider a Neumann boundary condition for w and σ.If we rather impose a Dirichlet boudnary condition with a vertical component of the incoming flux Q vert , then one has to take into account the divergence constraint.This means that, while one imposes the vertical flux on hw, hσ on the boundary has to verify the corresponding equation, yielding: Outgoing flow with prescribed height: For the outgoing flow, we impose a Dirichlet type boundary condition for the fluid height: Once again, the outgoing characteristic of the hyperbolic problem preserves the corresponding Riemann-invariant, so This means that the velocity on the right hand side phantom cell is given by For the vertical compontents of the velocity, we impose Neumann boundary condition, just as before, hence Since we imposed a Neumann type boundary condition on the velocity u, we prescribe a Dirichlet type boundary condition on the pressure terms.However, we have to distinguish the two discretizations when defining whether we impose Dirichlet type boundary conditions on the conservative variable hq or the non-conservative variable q to ensure the compatibility for the transposed quantities.
( ) In this case, we impose the Dirichlet type boundary condition on the nonconservative variable q: With the boundary conditions at our disposal, we can write the discretization of the momentum equation for the first cell: Since q N +1 = −q N , this means that and that no correction is required, 0 N = 0. So, given the matrix B 1 , we can compute the approximation for the divergence which is a first order discretization of −(h∂ x u + 2 √ 3σ).The corresponding correction term is given by ( ) In this case, we impose the Dirichlet type boundary condition on the conservative variable hq: With the boundary conditions at our disposal, we can write the discretization of the momentum equation for the first cell: and that no correction is required, 0 N = 0. So, given the matrix B 2 , we can compute the approximation for the divergence which is a first order discretization of −(h∂ x u + 2 √ 3σ) up to the error term 0 N = −

Numerical results
To compare the numerical strategies, three cases are considered: • An exact solution (solitary wave) where the three numerical methods are compared; • A dispersive shock wave from a dam break simulated with the Uzawa algorithm; • Comparisons with experimental data (Favre waves) and the numerical results obtained with the Uzawa algorithm.
In all the following simulations, the threshold used for the stopping criterion is set to = 10 −3 .

Solitary wave solution
It is a well known fact that the system admits solitary wave solutions in the case of constant atmospheric pressure and flat bottom topography.
Let us fix p atm ≡ cst and z b ≡ cst.For a > 0 and H 0 > 0, we set c = g(H 0 + a) and are solutions to System (6) in dimension 1.
For this test case, we consider the propagation of a solitary wave from the left to the right, initially centered at x 0 = 10 m, of relative amplitude a = 0.2 m, over a constant water depth H 0 = 1 m.The computational domain is 100 m long and discretized with 1280 cells.
Since the solitary wave is initially far from boundaries, the boundary conditions do not affect the computation, so we have chosen to impose free boundary conditions at the downstream and upstream ends for the sake of simplicity.We compare the water surface profile of our numerical solutions to the exact profile given by ( 22) at several times using two iterative methods, namely Uzawa (18) and Gauss-Seidel (21), and a direct resolution of (20).
Figure 2 shows that all schemes manage to capture the correct wave velocity but the Uzawa method seems to be less diffusive than other methods.The explanation may be related to a faster convergence of the iterative procedure with Uzawa than with Gauss-Seidel.
Table 1 shows the computational time and the relative error between the numerical solutions and the exact one given by ( 22): To obtain these results, we set a time T = 5 s and compute the numerical solution by increasing the number of cells from N = 80 to N = 1280.
x(m) Figure 2. Propagation of solitary wave over a flat bottom: water surface profiles at t = 5, 10, 15 and 20 s.Comparison of the numerical and exact solutions for the Uzawa, Gauss-Seidel and a direct method Remark 3.1.Note that the Gauss-Seidel method and the direct method were implemented using Python, whereas the Uzawa method was implemented using Fortran 95.This may explain the delay of computational times (CPU) between these methods.
From now on, all tests cases are run with the Uzawa method.

Dam-break
We consider the dam-break problem in order to study the ability of our numerical scheme to capture dispersives shocks.It is well known that discontinuous initial data of this type generate dispersive shock waves due to dispersive effects [6].This issue for non-linear and dispersive shallow water equations was also investigated in [6,13,21].We consider the following initial data: The computational domain is the interval (−300, 300) and discretized with 8000 cells.
Figure 3 shows the dispersive effects generated by System (1) when using discontinuous initial data.Here, we chose more cells (8000) because it is claimed in [7,21] that dispersive effects are hard to capture.It is recommended that high-order numerical schemes be used but that is not the topic of this paper.A third-order method is applied in [16].

Favre waves
In the following section, we compare our numerical scheme based on the Uzawa algorithm with experimental data.A good example of a physical situation in which the non-hydrostatic pressure terms play a significant role is the so-called Favre wave experiment.Favre waves are secondary free surface waves, undular bores that appear typically in a channel following a sudden change in the flow discharge.Such a situation occurs in the case when a sluice gate is suddenly closed, realised during the physical experiments carried out at EDF (Électricité de France) [12].A long and thin wave tank is considered with a small constant slope and open uphill entry and downhill exit as well as vertical walls (see Figure 4(a)).A stationary water flow is established in the basin up until a time t 0 when a gate rapidly descends at the exit level, blocking the flow and generating an upstream moving series of undulating waves.
In the experiments, a wave tank of length 40 m was taken; the lower 30 m section plays a significant role.The tank is of width 0.4 m, the walls on the side are vertical and the bottom is of a constant slope of At the lower end of the basin, a vertical solid gate is placed perpendicular above the water flow.Additionally, three sensors were also placed over the basin, each one measuring the water height at its horizontal position.The first sensor is situated 2 m from the exit of the tank, the second sensor is 15 m from the lower end, and finally the third sensor is at a distance of 30 m from the gate (see Figure 4(b)).As it was mentioned before, the gate is closed at time t 0 (when the flow has reached a stationary regime with the above given physical parameters).The closure takes place rapidly, in a time frame of T = 0.07s, when the gate completely blocks the exit, forming a solid wall boundary for the gravity flow.Measurements on the water height are available in the timeframe of [t 0 + 1.31, t 0 + 5.35] for the first sensor, in the timeframe of [t 0 + 10, t 0 + 16.9] for the second sensor, as well as in the timeframe of [t 0 + 21, t 0 + 28.92] for the third sensor.
In the context of an extended Saint-Venant system, in which non-hydrostatic pressure terms enrich the classic shallow water equations, first and second order finite volume solvers have already been developed and validated against this particular experimental case [9].
In our case, the simulations were implemented with a one horizontal dimensional system corresponding to the experiments (see Figure 4(b)).The implemented physical parameters were the same, except for the the flux,  which had to be recomputed for its one dimensional equivalent: Close to the gate (2 m), Figure 5 shows that the first dispersive wave is correctly captured by the scheme while the following ones are damped.The standard shallow water equations only provide a shock wave with the correct velocity but nor the shape neither the other waves.As expected, the Serre-Green-Naghdi equations provide a solution with a more relevant structure.Further (15 m), the damping is still observed but the period seems to be preserved.Once again, a higher-order method should be considered.

Conclusion
The non-hydrostatic formulation of the well-known Serre-Green-Naghdi equations is investigated in this paper.Unlike its Boussinesq formulation, only first order derivatives are involved in the present one.The consequence is a larger number of unknowns.The physical content of this model is richer than the shallow water equations but the computational cost is higher: at each time step, in addition to the resolution of a shallow water system, an elliptic equation must be solved involving non-standard operators.The latter one is proved to be well-posed in any dimension d ∈ {1, 2}.We designed some 1D numerical schemes based on the algebraic structure of the underlying linear system which mimics the properties of the continuous equation.They provide reliable results whether it be for analytical solutions or experimental solutions.It is a preliminary work which is extended in [16] where more numerical tests are carried out with a third-order method; in addition an extension of the numerical strategy to the multilayer extension of the Serre-Green-Naghdi model [17] will be proposed.

Figure 3 .
Figure 3. Dam-break: water surface at several times (a) The wave tank of the experiment (b) A cross section of the wavetank

Figure 4 .Figure 5 .
Figure 4. Schematic representation of the Favre wave experiment

Figure 6 .
Figure 6.Comparison between the numerical solution (obtained with N = 3000 nodes) and the experimental data for Sensor 2 (15 m).

Table 1 .
Relative L 2 -error table for the variables h and u and and time execution