Numerical Simulations of a Dispersive Model Approximating Free-Surface Euler Equations

In some configurations, dispersion effects must be taken into account to improve the simulation of complex fluid flows. A family of free-surface dispersive models has been derived in Fernández-Nieto et al. (Commun Math Sci 16(05):1169–1202, 2018). The hierarchy of models is based on a Galerkin approach and parameterised by the number of discrete layers along the vertical axis. In this paper we propose some numerical schemes designed for these models in a 1D open channel. The cornerstone of this family of models is the Serre – Green-Naghdi model which has been extensively studied in the literature from both theoretical and numerical points of view. More precisely, the goal is to propose a numerical method for the LDNH2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$LDNH_2$$\end{document} model that is based on a projection method extended from the one-layer case to any number of layers. To do so, the one-layer case is addressed by means of a projection-correction method applied to a non-standard differential operator. A special attention is paid to boundary conditions. This case is extended to several layers thanks to an original relabelling of the unknowns. In the numerical tests we show the convergence of the method and its accuracy compared to the LDNH0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$LDNH_0$$\end{document} model.


Introduction
Water waves and more generally water flows are of great interest in several scientific fields with applications to society issues such as protection of populations (tsunamis, floods, …) or energy production (water-turbines, …). Depending on the accuracy that is required in the applications, more or less complex models are used, from fully resolved to averaged equations.
In order to simulate the behaviour of free-surface fluid flows, let us consider the 2D Euler equations for an incompressible free-surface flow under gravity: ∇ ·û = 0, (1a) ∂ tû + (û · ∇)û + ∇p total = (0, −g), set in the moving domain (see Fig. 1) Here, I = (x , x r ) is a bounded interval of R and z b is the given topography (independent from time). The unknowns are the velocity fieldû = (û,ŵ), the pressure fieldp total and the water elevationη. The water height is deduced fromĥ(t, x) =η(t, x) − z b (x). Viscosity effects are not taken into account in this work but we refer to [14] for instance to deal with such terms. The model is supplemented with kinematic boundary conditions: ∂ tη +û s ∂ xη −ŵ s = 0, as well asp total t, x,η(t, x) = p atm (t, x),  (1) for some given atmospheric pressure p atm (t, x) > 0. Classically, the pressure field is decomposed into three parts:p where the hydrostatic part is given bŷ and whereq is the hydrodynamic pressure field -or commonly referred to as the nonhydrostatic component. Hydrostatic models such as the nonlinear shallow water equations [26] or the hydrostatic Navier-Stokes equations [3] are based on the assumptionq ≡ 0.
In addition to its complex mathematical structure (see e.g. [12,38,49]), Model (1) coupled to (2)(3)(4) is a real challenge owing to the fact the domain is moving: the water elevation is an unknown in itself. That is why, in spite of the increase of computer performance, reducedcomplexity models have been introduced, analysed and discretised. There exists an extensive literature about models approximating the Euler equations under simplifying hypotheses. Depending on the physical phenomena at stake and thus the fluid regime under study, some models turn out to be more accurate than others.
For example, the nonlinear shallow water equations (NLSW) [26,36,37,53] provide relevant results for large wavelengths but seem restrictive in other regimes, in particular due to the absence of dispersive effects. To go further, models were derived in a given regime of magnitude for the nonlinearity (ε: wave amplitude/water depth ratio) and for the frequency dispersion (μ: water depth/wavelength ratio). The competition between the two phenomena is responsible for the shape of water waves. Dispersion is necessary for instance for stratified flows or close to coastal areas. Different regimes lead to weakly/fully nonlinearweakly/fully dispersive models. For example, the NLSW equations correspond to μ = 0 for any ε = O(1).
It is worth noticing that classic dispersive models, usually written under a Boussinesq form, introduce high order derivatives for two unknowns (water height and velocity). On the contrary, non-hydrostatic pressure models introduce such effects by means of the nonhydrostatic component of the pressure. In [21] a semi-implicit finite difference model for non-hydrostatic free-surface flows is presented based on a vertical discretisation of the 3D free-surface Navier-Stokes equations. A similar approach is used in [63] where again the pressure is decomposed into hydrostatic and non-hydrostatic counterparts and a vertical discretisation is proposed for the 3D system. A more recent work on the subject is [22], where semi-implicit methods are extended to complex free-surface flows that are governed by the full incompressible Navier-Stokes equations and are delimited by solid boundaries and arbitrarily shaped free-surfaces. These approaches are comparable to the multilayer or layer-averaged framework, see for example [5,50,60]. In [5] and [60] the multilayer model is deduced by supposing within each layer a constant profile of the horizontal velocity. The multilayer model proposed in [50] is based on the irrotational hypothesis and that the horizontal velocity is quadratic (within each layer) and continuous at the interfaces. In [6] a numerical study of the linear dispersive celerity shows convergence of the multilayer model. In [35], a hierarchy of layer-averaged models has been derived: it is based on hypotheses made for the vertical velocity and pressure profiles. An explicit linear dispersion relation for each model is provided and authors proved that when the number of layer increases, celerity converges to the exact one given by Euler equations in the Stokes linear theory.
The latter model is considered here. Classic dispersive systems, and in particular DAE and SGN, may be written as a non-hydrostatic system. Readers can refer to [32] for the link between the two approaches. Non-hydrostatic formulation has several advantages from the numerical point of view, especially regarding boundary conditions. Notice that they do not rely on any irrotational assumption.
In recent years, an effort has been made in using hyperbolic approximations of dispersive an non-hydrostatic systems. See for instance [4,40] where hyperbolic depth-averaged equations use a hyperbolic approximation of a Boussinesq-type equation, or [54] for hyperbolic approximations of the Korteweg-de Vries equations. Other examples, based on non-hydrostatic models are given by [7,17,34] where authors proposed different hyperbolic approximations for the SGN equations, or [31,32] where authors proposed hyperbolic approximations for some well-known non-hydrostatic pressure and Boussinesq-type systems.
In this paper, we will focus on non-hydrostatic models. First, the non-hydrostatic model derived in [1], named DAE, which is an approximation of Euler equations when the frequency dispersion is small, is considered: where (u, w) is an approximation of (û,ŵ) along the water column, q is the non-hydrostatic component of the pressure field and γ > 0 a parameter. The case γ = 2 was studied in [16] and simulated in [2]. For a flat topography ∂ x z b ≡ 0, γ = √ 3 corresponds to the SGN equations [39]. This is a well-known dispersive system that has been extensively studied, see for instance [10,23,28,45] The SGN model, also studied in [35], is given by Moreover, both (5) and (6) are proved in [35] to belong to the same hierarchy of models depending on some orders of approximation for the set of unknowns 1 . Indeed, the key point is the dependence of the solution on the vertical coordinate z.
Remark that System (6) is made of first-order PDEs which account for conservation of mass (6a) and of momentum (6b-6d). (6e-6f) are diagnostic equations deduced from the incompressibility constraint (1a). There are several equivalent ways to write down the diagnostic equations (6e) and (6f) but this very definition is the only choice that provides a duality relation between the modified "pressure gradient" and the modified "velocity divergence".
This duality property is crucial from both theoretical and numerical points of view as highlighted in the following. It also allows to specify the boundary conditions: ① From the weak formulation associated to the underlying elliptic equation, we deduce the boundary conditions for the pressure fields (see Remark 3 which shows that it is necessary to specify q(∂ x (hq) + q b ∂ x z b ) at each boundary); ② We then infer from the numerical point of view how boundary conditions should be naturally imposed by symmetry in the resulting linear system (see § 16.3.2).
For the sake of simplicity and without restrictions 2 , we will consider in this paper a subcritical flow (which corresponds to the hyperbolic regime 3 in the Euler system (1)). Then, Model (6) is supplemented with the following boundary conditions: • For the non-hydrostatic pressure, we impose: for some flux χ ∈ R. The choice of these boundary conditions relies on the weak formulation associated to the underlying elliptic problem (see Remark 3) and on the boundary condition for the topography (7e).
• For the velocity field and the water height, it depends on the underlying hyperbolic regime. As mentioned above, we consider a subcritical flow so that, following the investigations of § 16.3.2, we impose for some constants Q u > 0, Q w ≥ 0 and h r > 0. We also assume that the topography satisfies the following boundary condition: It is shown in § 6.1 that System (6) can be rewritten under three different other formulations (Props. 1, 2 and 3). In particular, Proposition 1 will show that (6) is equivalent to the SGN equations under Boussinesq form, as presented in [45]. Then, from now on, System (6) coupled to boundary conditions (7) will be referred to indifferently as SGN or L DN H 2 (1). Although Model L DN H 2 (1) has a larger number of unknowns than L DN H 0 (1) (which may provide a richer modelling), the first goal of this paper is to design a numerical method for (6) with a reduced additional computational complexity compared to (5). See [56] for another strategy. The second goal of this paper is to extend the aforementioned numerical method designed for L DN H 2 (1) to its multilayer counterpart L DN H 2 (L) (for some number of layers L > 1) as derived in [35]. Indeed, L DN H 2 (1) relies on the approximation that the horizontal component of the velocity field does not depend on z, i.e. it is constant along the water column. In some cases (for instance when the flow is not shallow), it is necessary to add more degrees of freedom. The Galerkin method used for the semi-discretisation in z leads to L DN H 2 (L) where L is the number of vertical cells.
The paper is organised as follows: in Sect. 6, we first show some properties of the L DN H 2 (1) model: equivalent formulations (including the relationship with the SGN equations) and associated energy. Secondly, a numerical strategy to solve the monolayer model L DN H 2 (1) is presented: it consists in an iterative algorithm taking into account the gradientdivergence duality. In Sect. 17, the layer-averaged extension L DN H 2 (L) is reformulated with similar differential operators. The previous numerical method is then extended to the multilayer case, still based on the gradient-divergence duality. In Sect. 18, some classic numerical tests are presented to assess these strategies for L DN H 2 (1) and L DN H 2 (L), in particular in terms of accuracy with respect to L DN H 0 (L).

Analysis and Numerical Method for the LDNH 2 (1) Model
In this section we first study some of the properties of Model (6), such as some reformulations and its associated energy. Secondly, the design of a numerical discretisation of the model is presented.

Properties of the Model at the Continuous level
Let us set Let us also define the operators and notice that the following duality relation holds Given these notations, System (6) reads in a compact form System (10) has the same mathematical structure as the incompressible Euler equations with variable density except that differential operators are non-standard. Hence, this remark will help to design a similar numerical strategy.
It is shown below that this system can be rewritten under three different other forms (Props. 1, 2 and 3).
First, let us remark that System (10) is nothing but the SGN equations as presented in [45]:

Proposition 1 System (10) can be rewritten under the Boussinesq formulation
where Proof Indeed, we deduce from (6) that Inserting the two last equalities into (6b) leads to the expected result.
The main difference between the non-hydrostatic formulation (6) and the Boussinesq formulation given in Proposition 1 lies in the order of spatial derivatives (1 vs. 3). The consequence is a larger number of unknowns in (6) (6 vs. 2) but with at most first order derivatives. One key advantage of (6) is that the treatment of boundary conditions is easier as it will be shown later. Moreover, a smaller stencil is needed when one has to approach numerically first order compared to higher order derivatives.
A numerical algorithm to simulate (11) is designed in [11] by introducing a second-order parameterized perturbation and based on the inversion of I d + T . The parameter is set so that the linear dispersion relation is optimised with respect to the Airy relation. The numerical technique consists of a splitting method between the hyperbolic Saint-Venant equations (solved with a finite-volume scheme) and the dispersive part (solved with a finite-difference scheme). An extension to dimension 2 is provided in [47] by replacing the differential operator I d + T by a time independent "diagonal" approximation. In [28], a Discontinuous Galerkin Finite-Element method is applied.
We also mention the following formulation with high order derivatives of h: Here we used the standard notationξ : Inserting the two latter equalities into (6b), we obtain the expected result.

Remark 1
For a flat topography, we recover the model studied in [34,48]. In [48], the model and solved using a finite-difference method. The numerical technique used in [34] is based on an augmented-Lagrangian approach.
In the framework of projection methods, we exhibit another formulation. The two constraints (10c) can be replaced by applying the divergence (8b) to the momentum equation (10b) (initially divided by h):

Proposition 3 System (10) can be rewritten under the reduced formulation
Well-posedness of this system is studied in [8].

Remark 2
Notice that only h and u are involved (not w nor σ ) in the right hand side of (13). Moreover, we also notice that q b can be expressed directly from (13b) Inserting the latter expression into (13a) provides a unique equation for q: However, the complexity of the operators in the latter equation made us prefer working with (13).

Remark 3 Whether it be from the theoretical or numerical points of view, the equation of interest is
for some right hand side f. The well-posedness of Eq. (15) is proved in [8] by means of the Lax-Milgram theorem applied under some smoothness hypotheses on the water height h. In particular, straightforward computations show that which explains the present choice of boundary conditions (7a-7b).
Finally, we show that thanks to the correct formulation of diagnostic equations (6e) and (6f), associated to the definition (8b) of ∇ sgn · X that provides the duality relation (9) we can deduce an associated energy for SGN model. It is crucial from both theoretical and numerical points of view as highlighted in the following section, corresponding to the definition of the numerical method. (10) satisfy the following energy equality provided p atm and z b do not depend on time:

Proposition 4 Smooth solutions of System
Proof Let us multiply (6b) by u, (6c) by w and (6d) by σ so that Using the duality relation (9) as well as Eq. (6a) leads to the expected result.

Splitting Strategy at the Semi-Discrete Level
Let us write a semi-discretisation of System (10) based on a classic splitting technique for some time step t > 0, like e.g. in [1,30,33,47,56]. We first consider the hyperbolic step 4 coupled to boundary conditions (7c) and (7d). Any classic numerical method dedicated to the Shallow Water equations (including well-balanced schemes) can be used to solve it. Then the dispersive step reads which reads as a mixed velocity-pressure Darcy problem. Boundary conditions (7c), (7a) and (7b) are considered. In this paper we investigate a projection-correction approach.
To preserve the order of the method, we can consider the incremental method (see e.g. [41]) which consists of the following modified first step (which requires an initialisation of the pressure unlike the non-incremental version) and the modified second step This approach does not raise any additional issue, therefore for the sake of simplicity we present the results by considering the splitting strategy (16)- (17).
To make the reading easier, we shall denote from now on: X = X n+1 and or equivalently, expanding the non-classic operators: We recover the same operators as in the continuous case -see (13), only the right hand side is modified due to the splitting method. Hence the well-posedness investigated in [8] still holds.

Remark 4 The right hand side in (18') is not of order
inserting the values from the hyperbolic step (16), we get since 2 √ 3σ n + h n ∂ x u n = 0 and w n − u n ∂ x z b − √ 3σ n = 0, provided initial conditions satisfy the divergence constraints, i.e. initial conditions are well-prepared.
Once Equation (18) is solved for Q, the velocity field is updated using (17b)

Numerical Method at the Discrete Level
We focus in this section on the discretisation in space of the non-incremental version (16)(17). The incremental version does not raise major additional difficulties.
Let us consider a homogeneous Cartesian grid of interval Fig. 2). Standard approaches for Stokes-like problems rely on a staggered grid with velocity fields at the centre of the cells and pressure fields at the interfaces. A similar approach is also followed in [33]. The specific expression of ∇ sgn Q in the present problem -see (8b)induced a different choice, namely a collocated approach with all variables (X and Q) at the centres of the cells. This choice makes a discrete energy estimate easier to derive.

Hyperbolic Step
For the sake of completeness, let us briefly describe the numerical scheme applied to system (16), which consists of the classic shallow water part and two transport equations for w and σ .
Here, we considered a high-order Polynomial Viscosity Matrix finite-volume method using a second-order MUSCL state reconstruction operator similar to those proposed in [30,33]. Moreover, a third-order CWENO - [25] reconstruction has been implemented as well.
The semi-discrete in space high-order path-conservative scheme for System (16) reads as follows Dropping the dependence on time, D ± i+1/2 is given by a Polynomial Viscosity Matrix (PVM) path-conservative scheme Here, F(h X) = hu X and S(h X) = (gh, 0, 0) t . In the expression above, Q i+1/2 is the numerical viscosity that depends on the choice of the PVM method and (h X) ± i+1/2 is computed by means of a reconstruction procedure on the variables to the left ( − ) and right ( + ) of the inter-cell x i+1/2 . This is done using a standard MUSCL or CWENO reconstruction operator P i,(h X) (x).
accounts for the path-integral of nonconservative terms. In this work, the simple straight-line segment path is chosen, and thus the path-integral can be computed exactly using the midpoint rule. In particular, all the components are zero except the first one that reduces to The viscosity term Q i+1/2 results from the evaluation of a Roe Matrix at x i+1/2 by a polynomial that interpolates the absolute value function (see [20]). Here we are using a first-degree polynomial that uses some estimation of the slowest (S L ) and fastest (S R ) wave speeds. With this choice, the numerical scheme coincides with an extension for non-conservative systems of the standard HLL Riemann-solver.
Finally, the volume integral is numerically approximated using any Gaussian quadrature rule of the same order as the scheme.
The considered explicit high-order methods are well-balanced for water at rest solutions The eigenvalues of the hyperbolic operator in (16) are u n ± √ gh n and u n (transport equations for w and σ ). Therefore, the scheme is linearly L ∞ stable provided the following condition holds Moreover, the scheme is positive-preserving for the total water depth in the sense that: if h n i > 0 then h n+1 i > 0 for all i.

Non-Hydrostatic Step: Numerical Scheme for the Velocity-Pressure Problem
We could have directly discretised Eqs. (18'). However, in order to recover a discrete energy, we first discretise the velocity-pressure problem (17) and we then deduce a discretisation for (18') mimicking the continuous level. Finite differences are applied in the present work. Let us first consider the mixed formulation (17). Denote by H ∈ M N ,N (R) the diagonal matrix with entries H i,i = h * i . We expect that the discretisation of System (17) has the form where is an rectangular matrix to be specified; • 0 ∈ R 3N and 0 ∈ R 2N account for boundary conditions as detailed below.
The matrix in (20) is symmetric due to the duality relation (9). However, a naive approach where Eqs. (17b) and (17c) are discretised independently from each other would lead to a non-symmetric matrix. Hence we first discretise Eq. (17b) as Once B is built, B T X = 0 should be a consistent discretisation of (17c) up to boundary terms. More precisely, Equation (17b) is discretised for inner cells: which leads to B 12 is diagonal with entries (∂ x z b ) i while B 11 is tridiagonal: Hence, for inner nodes x i : which is consistent with (17c).
Boundary conditions. Let us now specify how boundary conditions are incorporated by means of ghost cells. BC (7a) at x = x 1/2 = x is discretised by Hence, the momentum equation for u is discretised in the first cell by: This provides the first line of Matrix B 11 as well as 0 1 = − h * 1/2 2h * 0 χ . It remains to specify h * 0 . Given that no boundary condition is imposed on h at x = x , we set h * 0 = 2h * 1 − h * 2 . We then observe that We remark that if we discretise BC (7c) by Note that this result implies that imposing ∂ x (hq) = χ at some boundary is consistent with prescribing hu = Q u . For the boundary conditions at x = x N +1/2 = x r , BC (7b) is discretised by which means that the momentum equation for u in the last cell becomes This yields the last line of B 11 and 0 N = 0. It implies that and 0 N = 0. Let us comment the latter result. Imposing q = 0 at some boundary is consistent with imposing ∂ x u = 0 at the same location. Even if there is no BC required in the hyperbolic problem, it is necessary from the numerical point of view in order to compute the flux at the corresponding interface. We mention that the BC handled in the first step of the algorithm must not be damaged by the second step. Likewise, still for numerical purposes, we impose

Numerical Scheme for the Pressure Problem
For a resolution of the mixed velocity-pressure problem (20) by means of the Uzawa method, see [8]. In the present paper, we rather combine both equations in (20) to obtain which is nothing but a consistent discretisation of the velocity-correction approach (18).

Remark 5
It is crucial that the initial data satisfy at the discrete level B T X * − 0 = 0 in order to prevent the propagation of errors of order O( t −1 ) as explained in Remark 4 Notice that Matrix C = B T H −1 B has the following structure: Then BQ = 0 ⇒ Q = 0 as B is one-to-one -see (21). More precisely, (24) provides the following discrete energy which is consistent with To solve (23), one may apply different strategies. One of them is to apply a direct resolution inverting C. Another variant is to use the Uzawa method. Here, we shall propose to follow an iterative process of Gauss-Seidel type: . This method converges due to Lemma 1. This is the latter method which has been chosen in the present work. Matrix C 11 is more or less the same matrix as for the DAE model (5). The extra additional time then corresponds to the iterative process (C 11 is factorised once and C 22 is diagonal).

Remark 6
Notice that we could explicitly solve the second block of (23): which corresponds to a discretisation of (14).

Remark 7
It is easy to check that the numerical scheme is well-balanced for water at rest solutions The hyperbolic step for (10a)-(10b) is well-balanced. We refer the reader to [18,20] for the details.
Since the first step of the numerical scheme is well-balanced for the hydrostatic underlying system, then the proposed iterative process in Lemma 1 yields to the resolution of two homogeneous linear systems with non-singular matrices C 11 and C 22 Therefore q = q b = 0, and thus the result follows.

Numerical Treatment of Wet-Dry Fronts
In order to ensure the robustness of the algorithm, a wet/dry treatment is needed for small or null values of the water height h. In particular, we use the following strategy: 1. The hydrostatic pressure terms gh∂ x (z b + h) at the horizontal velocity equations are modified for emerging bottoms to avoid spurious pressure forces (see [19]).
for some threshold ε > 0. In the tests described in this paper ε is set to 10 −6 .
Moreover, let us remark in vacuum situations the iterative process described in Lemma 1 leads to the resolution of two homogeneous linear systems with non-singular matrices C 11 and C 22 • C 11 q k+1 = 0; Therefore q = q b = 0 when the water thickness vanishes.

Notations and Model
Without any assumption on the shallowness of the flow, we can also approximate the Euler equations (1) by means of a "multilayer" model. The flow is split into arbitrary layers of constant height h α = h L for some integer L > 0 (see Fig. 3). We consider the dispersive layer-averaged model derived in [35] under the following reformulation: and for α ∈ {1, . . . , L} and the boundary condition The mass transfer term is given by To close the model, we define for any γ α+1/2 ∈ [0, 1] such that [35,Prop. 1] Let us recall the properties of this system: If (H , u α , w α , q α ) are smooth solutions to the multilayer model, we have under (25) (26) Moreover, if we take γ α+1/2 = 1 2 , then (26) is an equality. • The monolayer case L = 1 reduces to the SGN equations (6).

Differential Operators
Let us introduce the operators We get a global duality relation It is worth noticing that and We mention that the divergence operator stated in the present work is not the same as in [35]. We rather provide an equivalent formulation that satisfies the duality relation (27). In particular, σ α is involved instead of ∂ x u α which is important for the order of derivatives.
Let us proceed as in the single layer case with a splitting strategy between hyperbolic terms and non-hydrostatic terms. The latter part reads which implies The latter equation is expanded as, for α ∈ {1, . . . , L} and finally Let us remark is that the differential operator for q α in (30a) is independent from index α. Moreover, up to a coefficient L 2 , this operator is the same as in the monolayer case -see (18').

Structure of the Fully Discretised Equations
Our goal is to produce an algorithm that only relies on tools designed for the monolayer case. Let us split the gradient matrix B = B 1/2 from (21) into: B 11 is defined in equation (22) and B The fact that B 1 is independent from layer α (and exactly the same as in the monolayer case) is crucial for what is following.
Hence ∇ α ldnh Q α is approximated according to (28a) The last two terms account for interactions between layers: the model under investigation does not reduce to a monolayer model in each layer. Hence System (29) is approximated by ⎛ It is easy to verify that B + R T X = 0 is consistent with (28b). For the treatment of boundary conditions which are incorporated in vectors 0 and 0, we refer to the monolayer case ( § 16.3.2). As previously, the symmetry of the global matrix is due to the duality relation (27). From the discrete velocity-pressure problem, we deduce as previously a discrete pressure problem which reads Let us set which is blockwise tridiagonal, symmetric positive-definite; which has exactly the same structure as C in (24) with The blockwise components of C are then: (31) is a consistent discretisation of (30).

Iterative Scheme
To solve the 2N L × 2N L linear system (31), we apply an alternating direction-type method. More precisely, it amounts to solving iteratively: x-direction For each layer α, (30a) in q α knowing q α+1/2 which corresponds to a linear system with matrix C 11 given in (24) (the same for all layers which requires a single factorisation for all iterates at each time step): z-direction For each node x i , (30b-30c) in q α+1/2 knowing q α . It is a tridiagonal system with matrix S (i) solved by means of the Thomas' algorithm [64].
Indeed, we check that This strategy is equivalent to solving the following system by means of a Gauss-Seidel iterative procedure: ⎛ where unknowns have been re-labelled as ⎛ Remark 8 It must be underlined that this algorithm can be easily parallelised insofar as each direction (x and z) involves a blockwise diagonal matrix.

Numerical Simulations
Let us recall that the stability condition is prescribed by the hyperbolic part of the splitting strategy -see (19) with a CFL number to be specified. The numerical schemes presented in this paper, namely the resolution of (23) for L DN H 2 (L = 1) and the resolution of (31) for L DN H 2 (L), L ≥ 1, are assessed by means of some classic test cases and compared to numerical results obtained with the DAE model (5) and its multilayer counterpart described in [35].
It is a well-known fact that dispersive non-hydrostatic models usually require an extra dissipative mechanism to properly treat breaking waves near the coast. Otherwise, artificial overshooting may appear when the wave approaches the shore, see for instance [33].
A breaking criterion as the ones presented in [33] for non-hydrostatic pressure system or in [30] for a two-layer non-hydrostatic pressure system can be easily adapted to the formulations shown in this paper. Nevertheless, we intend here to do a fair comparison of the models L DN H 0 and L DN H 2 by themselves, without any other extra treatment such as breaking. Note that such a breaking criteria would be needed in some cases for very fine meshes, in particular in Tests 18.3 and 18.5.

Convergence Test
Let us start by a convergence test for both models, L DN H 0 and L DN H 2 , in order to assess the numerical strategy and the code. To do so, we consider the propagation of a solitary wave in a rectangular channel with constant topography, that is, a travelling wave that propagates at constant speed, without change in its shape. Let us remark that the particular definition of solitary waves depend on the particular model used, see for instance [61] and the references therein for further details on this topic. Let us consider here the one layer models (5) and (6). Following [15], solitary waves for these models are given by Here we set H * = 1, A = 0.1, g = 1 and z b = 0. The propagation of a solitary wave over a long distance is a standard assessment of stability and conservative properties of numerical schemes for Boussinesq-type equations [58,59,65]. A solitary wave propagates at constant speed and without change of shape over a horizontal bottom.
The domain is [−20, 20]. We perform the simulation with different numbers of volume cells at time t = 0.2 with a second-order scheme. The results are compared to the reference solution and the errors are shown in Tables 1 and 2. These results show the convergence towards the reference solution at second-order.

Still Water Steady-State
In this test we consider the same domain and boundary conditions as in the previous test.
The topography z b is given by We consider the initial condition given by h 0 (x) = −z b (x), and the rest of the flow variables are set to zero. The exact solution is a steady-state with still water, coinciding at all times with the initial condition. We have performed the numerical test for the L DN H 2 (L) system for L = 1, 4 with increasing number of cells. Tables 3 and 4 show the error observed at time T = 10, with a C F L number set to 0.9. It clearly shows that the scheme is well-balanced since these are all of the order of the machine precision.

Solitary Wave Propagation Over Reefs
A test case propagating a solitary wave over an idealised fringing reef assesses the ability of the model to handle nonlinear dispersive waves, breaking waves and bore propagation. The test configuration includes a fore reef, a flat reef, and an optional reef crest to represent fringing reefs commonly found in a tropical environment. Figure 4 shows  A solitary wave of amplitude 0.5 m is placed at point x = 10 m. Finally C F L = 0.9 and g = 9.81 m · s −2 . Free outflow boundary conditions are imposed. Figure 5 shows snapshots at different times, t √ g/H * = t 0 where H * = 1 m. Comparisons between experimental and simulated data allow to validate the numerical approach presented in this paper. Results are shown for L DN H 0 (L = 1) and L DN H 2 (L = 1) models. The water rushes over the flat reef without producing a pronounced bore-shape. The simulation also captures the offshore component of the rarefaction falls, exposing the reef edge, below the initial water level. The simulations match with experimental data and the L DN H 2 provides slightly better results.

Wave Propagation Over a Submerged Bar
The Dingemans experiment [27] of plunging breaking periodic waves over a submerged bar is considered. This case allows to study frequency dispersion characteristics and non-linear interactions. As waves propagate over a submerged bar, multiple phenomena occur, like the appearance of higher harmonics. The 1D domain [0, 30] is discretised with x = 0.005 and the bathymetry is defined on Fig. 6. Locations of the measurement points are specified in Table 5. The CFL number is set to C C F L = 0.9 and the gravity field to g = 9.81 m · s −2 . We run the numerical test from the "lake at rest" steady state as an initial condition. Boundary conditions correspond to free outflow at x = 30 and a sinusoidal wave train for η generated at x = 0. This is done as in [30] imposing in a relaxation zone: where A = 0.01 and T = 2.02 denote resp. amplitude and period. This test produces, up to the front slope, waves with wavenumbers k ≈ 0.63/H * and k ≈ 1.58H * respectively, where H * = 0.4 is the typical depth. Fig. 7 shows numerical results of time series of the free surface for Model L DN H 0 (L), for L ∈ {1, 2, 4}, while Fig. 8 concerns L DN H 2 (L) and Fig. 9 shows comparisons between L DN H 0 (L = 4) and L DN H 2 (L = 4).
Good agreements with experimental data are observed for all models up to Gauge #4. Beyond the bar, higher harmonics are released which explains discrepancies. We recover observations from the literature, such as [51] where σ -coordinates are used, or [30] where an enhanced two-layer version of the non-hydrostatic pressure system L DN H 0 is used. The results in [23] with a three-parameter Green-Naghdi model optimised for uneven bottoms, show the same level of agreement. Here, we would like to stress the ability of the proposed models to deal with a wide range of dispersive waves. The main difference between L DN H 0 (4) and L DN H 2 (4), as pictured on Fig. 9, can be seen for gauges #6 to #8 where the L DN H 2 model is more accurate.

Shoaling of a Solitary Wave on a Plane Beach
We finally consider the shoaling of a solitary wave on a beach with a constant slope (1 : 30) as described by Guibourg in [42] and then investigated in [11,28].
A sketch of the geometry is described on Fig. 10. The initial condition is a solitary wave at location x = 10 with amplitude A = 0.298, as described in [24,29]. The computational domain = [0, 27.5] is divided into cells of length x = 0.01. Free-outflow boundary conditions are considered and the CFL number is set to C C F L = 0.9.
We compare Models L DN H 0 (L) and L DN H 2 (L) for L ∈ {1, 2, 4}. Some temporal series of the free-surface elevation are measured at various locations (see Table 6) and compared with the corresponding numerical results.
Results are shown on Figs. 11. Numerical outputs for A = 0.289 predict the shoaling phenomenon during the wave run-up satisfactorily. The output clearly shows better performance for Model L DN H 2 , which is expected according to the linear dispersion relation of the continuum models.

Conclusion
This paper deals with some numerical approaches to simulate Models L DN H 0 and L DN H 2 introduced in [35]. The main objective is to compare their accuracy when applied to different standard test case scenarios. The L DN H 2 model presented in [35] may be seen as a multilayer extension of the Serre -Green-Naghdi equations. This model was derived from Euler equations assuming linear and quadratic vertical profiles for the vertical velocity and pressure. The L DN H 0 model presented in [35] differs from L DN H 2 due to the assumption of a linear vertical profile for pressure. One of the most attractive properties of the models relies on the increasing accuracy of the linear dispersion relation as the number of layers increases, whereas the model L DN H 2 shows better accuracy than L DN H 0 for the same number of layers.
The L DN H 0 model is solved by using a projection technique similar to the one introduced in [30]. The extension of this technique to L DN H 2 is not straightforward. We proposed a numerical method based on this projection technique to approximate the solution of the L DN H 2 model.
The complexity of the model L DN H 2 requires the design of an efficient strategy to solve it numerically for an increasing number of layers. To this aim, we have exploited here a duality relation at the continuous level. This allows to design an algorithm relying on an iterative process that solves a one-layer case in each iteration.
In particular, the algorithm may be decomposed into two different problems: The first one corresponds to a discrete parabolic problem per layer and the second to a tridiagonal linear system at each point of the horizontal discretisation that couples each layer. Moreover, (i) the matrix of the linear system of the parabolic problem per layer is the same for all layers, and (ii) each tridiagonal linear system at each point of the horizontal discretisation is independent. Then, we can observe that (i) since we have the same matrix for all layers, we may consider, for example, an LU-factorisation to reduce the computational time in all layers and (ii) this implies that, although this is not done in the paper, the proposed technique can be easily parallelised, especially when the number of layers increases.
Moreover, the final numerical scheme proposed here is also high-order, well-balanced for the water at rest solution, and positive-preserving for the total water depth. That results in an efficient and robust numerical scheme, even in the presence of wet/dry transitions. It is worth mentioning that, while an extensive literature is dedicated to the numerical resolution of shallow water flows such as the Depth-Averaged Euler equations or the Serre -Green- Naghdi equations, this is the first attempt, up to our knowledge, to design a robust and efficient numerical strategy for the multilayer extension of the SGN equations.
The proposed numerical scheme has been carefully validated, showing the second order of accuracy and comparing it with available experimental data. The obtained results exhibit an excellent fitting with the experiments and show that the proposed strategy is well-suited for most coastal processes: wave propagation, shoaling of the waves, run-up of waves onto a beach, higher dispersive harmonic waves, among others.
In the numerical tests presented in the paper, we can also see that for a fixed number of layers, the L DN H 2 is in many situations more accurate than L DN H 0 . This result shows that the L DN H 2 is of interest since the accuracy of both models has been compared by using a generalisation of classic projection techniques.
One should recall that L DN H 2 has twice more pressure unknowns than L DN H 0 . Therefore, if we apply a projection method to approximate the solution of L DN H 2 , it results in a more expensive algorithm, from the computational point of view, than for L DN H 0 . Hence, a fair comparison of both models from the computational point of view is beyond the scope of this paper.
Nevertheless, it is worth mentioning that this paper was aimed to analyse and propose a projection technique as a starting point for designing efficient and robust numerical methods of high-order for multilayer non-hydrostatic systems. As further research, it will be interesting to investigate the development of other numerical strategies for the L DN H 2 model and compare them with the one proposed here. For example, it may be of interest to consider relaxation techniques or rapid numerical methods and pseudo-compressibility approximations, as the ones proposed for one-layer dispersive shallow models in [31], [34] and [9], respectively.