Numerical simulations of a dispersive model approximating free-surface Euler equations

Journal of Scientific Computing (2021) 89:55 https://doi.org/10.1007/s10915-021-01552-6 Numerical Simulations of a Dispersive Model Approximating Free-Surface Euler Equations Cipriano Escalante Sánchez1·Enrique D. Fernández-Nieto2· Tomás Morales de Luna1·Yohan Penel3,4 ·Jacques Sainte-Marie3,4 Received: 17 December 2020 / Revised: 25 May 2021 / Accepted: 7 June 2021 / Published online: 18 October 2021 © The Author(s) 2021 Abstract 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 LDNH 2model 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 projectioncorrection 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 LDNH 0model. BEnrique D. Fernández-Nieto [email protected] Cipriano Escalante Sánchez [email protected] Tomás Morales de Luna [email protected] Yohan Penel [email protected] Jacques Sainte-Marie [email protected] 1Dpto. de Matemáticas. Universidad de Córdoba, Campus de Rabanales, 14071 Córdoba, Spain 2Dpto. Matemática Aplicada I. ETS Arquitectura - Universidad de Sevilla., Avda. Reina Mercedes N. 2., 41012 Sevilla, Spain 3Inria Paris, 2 rue Simone Iff, CS 42112, 75589 Paris, Cedex 12, France 4Sorbonne Université, Université de Paris, CNRS, Laboratoire Jacques-Louis Lions, LJLL, F-75005 Paris, France 123 55 Page 2 of 35 Journal of Scientific Computing (2021) 89 :55 Contents 1 Introduction ............................................... 2 6 Analysis and Numerical Method for the LDNH 2(1)Model ...................... 6 6.1 Properties of the Model at the Continuous level .......................... 6 12.2Splitting Strategy at the Semi-Discrete Level ............................ 10 16.3Numerical Method at the Discrete Level .............................. 11 16.3.1Hyperbolic Step ........................................ 12 16.3.2Non-Hydrostatic Step: Numerical Scheme for the Velocity-Pressure Problem ....... 13 16.3.3Numerical Scheme for the Pressure Problem ........................ 15 16.3.4Numerical Treatment of Wet-Dry Fronts ........................... 16 17 Numerical scheme for the LDNH 2(L)Model ............................. 18 17.1Notations and Model ........................................ 18 17.2Differential Operators ........................................ 19 17.3Structure of the Fully Discretised Equations ............................ 20 17.4Iterative Scheme ........................................... 22 18 Numerical Simulations .......................................... 23 18.1Convergence Test .......................................... 23 18.2Still Water Steady-State ....................................... 25 18.3Solitary Wave Propagation Over Reefs ............................... 25 18.4Wave Propagation Over a Submerged Bar ............................. 26 18.5Shoaling of a Solitary Wave on a Plane Beach ........................... 27 19 Conclusion ................................................ 31 References .................................................. 33 1 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: ∇·ˆ u=0,(1a) ∂tˆ u+(ˆ u·∇)ˆ u+∇ˆptotal =(0,−g), (1b) set in the moving domain (see Fig. 1) (t)=(x,z)∈R2x∈I,zb(x)≤z≤ˆη(t,x). Here, I=(x,xr)is a bounded interval of Rand zbis the given topography (independent from time). The unknowns are the velocity field ˆ u=(ˆu,ˆw), the pressure field ˆptotal and the water elevation ˆη. The water height is deduced from ˆ h(t,x)=ˆη(t,x)−zb(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ˆη+ˆus∂xˆη−ˆws=0, ˆub∂xzb−ˆwb=0,(2) as well as ˆptotalt,x,ˆη(t,x)=patm(t,x), (3) 123 Journal of Scientific Computing (2021) 89 :55 Page 3 of 35 55 Fig. 1 Fluid domain for the Euler equations (1) for some given atmospheric pressure patm(t,x)>0. Classically, the pressure field is decomposed into three parts: ˆptotal(t,x,z)=ˆphydro(t,x,z)+ˆq(t,x,z) where the hydrostatic part is given by ˆphydro(t,x,z)=g(ˆη(t,x)−z)+patm(t,x), and where ˆqis 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 assumption ˆq≡0. Given this decomposition, BC (3) is equivalent to ˆqt,x,ˆη(t,x)=0.(4) In addition to its complex mathematical structure (see e.g. [12,38,49]), Model (1) coupled to (2-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 nonlinear – weakly/fully dispersive models. For example, the NLSW equations correspond to μ=0for any ε=O(1). We may refer to pioneering works e.g. Boussinesq [13] and Peregrine [57], models obtained to improve linear dispersion properties e.g. Madsen & Sorensen [52]and Nwogu [55], models suitable for deep waters [46], …For reviews of dispersive models, 123 55 Page 4 of 35 Journal of Scientific Computing (2021) 89 :55 see [43,45]. A focus will be made below on two systems: the Depth-Averaged Euler (DAE) system described in [1]andtheSerre – Green-Naghdi (SGN) equations (see [39,62]). 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: ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∂th+∂x(hu)=0, ∂t(hu)+∂x(hu2+hq)=−γ2 2q∂xzb−gh∂xη−h∂xpatm, ∂t(hw) +∂x(huw) =γq, γw+h∂xu−γ2 2u∂xzb=0, (5) where (u,w)is an approximation of (ˆu,ˆw) along the water column, qis 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 ∂xzb≡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 123 Journal of Scientific Computing (2021) 89 :55 Page 5 of 35 55 ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂th+∂x(hu)=0,(6a) ∂t(hu)+∂x(hu2)+∂x(hq)+qb∂xzb=−gh∂x(zb+h)−h∂xpatm,(6b) ∂t(hw) +∂x(huw) −qb=0,(6c) ∂t(hσ)+∂x(hσu)−2√3q−qb 2=0,(6d) w−u∂xzb−√3σ=0,(6e) 2√3σ+h∂xu=0.(6f) 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 unknowns1. Indeed, the key point is the dependence of the solution on the vertical coordinate z. More precisely, Equations (5) are referred to as LDNH 0(1)as an approximation of the Euler equations where (u,w,q)is a (P0,P0,P1)Galerkin approximation of (ˆu,ˆw, ˆq).Likewise, Equations (6) are referred to as LDNH 2(1)and (u,w,q)is a (P0,P1,P2)Galerkin approximation of (ˆu,ˆw, ˆq). In particular, (6) has a more accurate linear dispersion relation than (5) for large wavelengths [35,Fig.2]. 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 3which shows that it is necessary to specify q(∂x(hq)+qb∂xzb)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 restrictions2, we will consider in this paper a subcritical flow (which corresponds to the hyperbolic regime3in the Euler system (1)). Then, Model (6) is supplemented with the following boundary conditions: •For the non-hydrostatic pressure, we impose: At x=x,∂ x(hq)=χ;(7a) At x=xr,q=0,(7b) 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). 1The models of the hierarchy are denoted LDNH k(L)where LDNH stands for Layerwise Discretisation Non-Hydrostatic,kis the order of approximation and Lthe number of layers. 2The location of boundary conditions depends on the one hand on the hyperbolic step (hence the hyperbolic regime) and on the other hand on the compatibility between boundary conditions for both steps. The boundary condition for the water height does not interfere with the elliptic step. 3The hyperbolic part is not the core of this paper which is rather the elliptic part. 123 55 Page 6 of 35 Journal of Scientific Computing (2021) 89 :55 •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 At x=x,hu =Qu,hw=Qw,hσ=Qw √3;(7c) At x=xr,h=hr;(7d) for some constants Qu>0, Qw≥0andhr>0. We also assume that the topography satisfies the following boundary condition: ∂xzb(x)=0.(7e) Itisshownin§6.1 that System (6) can be rewritten under three different other formulations (Props. 1,2and 3). In particular, Proposition 1will 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 LDNH 2(1). Although Model LDNH 2(1)has a larger number of unknowns than LDNH 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 LDNH 2(1)to its multilayer counterpart LDNH 2(L)(for some number of layers L>1) as derived in [35]. Indeed, LDNH 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 LDNH 2(L)where Lis the number of vertical cells. The paper is organised as follows: in Sect. 6, we first show some properties of the LDNH 2(1)model: equivalent formulations (including the relationship with the SGN equations) and associated energy. Secondly, a numerical strategy to solve the monolayer model LDNH 2(1)is presented: it consists in an iterative algorithm taking into account the gradientdivergence duality. In Sect. 17, the layer-averaged extension LDNH 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 LDNH 2(1)and LDNH 2(L),inparticular in terms of accuracy with respect to LDNH 0(L). 6 Analysis and Numerical Method for the LDNH2(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. 6.1 Properties of the Model at the Continuous level Let us set X=⎛ ⎝ u w σ⎞ ⎠,Q=q qb,and Sh=⎛ ⎝−gh∂x(zb+h)−h∂xpatm, 0 0⎞ ⎠.(8a) 123 Journal of Scientific Computing (2021) 89 :55 Page 7 of 35 55 Let us also define the operators ∇sgn Q=⎛ ⎝ ∂x(hq)+qb∂xzb −qb −2√3q−qb 2⎞ ⎠and ∇sgn·X=2√3σ+h∂xu w−u∂xzb−√3σ(8b) and notice that the following duality relation holds X·∇sgn Q=∂x(hqu)−Q·(∇sgn·X). (9) Given these notations, System (6) reads in a compact form ⎧ ⎪ ⎨ ⎪ ⎩ ∂th+∂x(hu)=0,(10a) ∂t(hX)+∂x(huX)+∇sgn Q=Sh,(10b) ∇sgn·X=0.(10c) 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,2and 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 ∂th+∂x(hu)=0, Id+T[h,zb](∂tu+u∂xu)+g∂x(h+zb)+Q[h,zb]u=−∂xpatm,(11) where T[h,zb]v=R1[h,zb](∂xv) +R2[h,zb](v∂xzb), Q[h,zb]v=−2R1[h,zb](∂xv)2+R2[h,zb](v2∂2 xxzb), R1[h,zb]w=−1 3h∂x(h3w) −h 2w∂xzb, R2[h,zb]w=1 2h∂x(h2w) +w∂xzb. Proof Indeed, we deduce from (6)that σ=−h∂xu 2√3,w=−h∂xu 2+u∂xzb, qb=−h2 2∂x(∂tu+u∂xu)−2(∂xu)2+h∂xzb(∂tu+u∂xu)+u2∂2 xxzb, q=−h2 3∂x(∂tu+u∂xu)−2(∂xu)2+h 2∂xzb(∂tu+u∂xu)+u2∂2 xxzb. 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 1lies in the order of spatial derivatives (1 vs. 3). The consequence is a larger number of unknowns in (6)(6vs. 2) but with at most first order derivatives. One key advantage of (6) is that the treatment of boundary conditions is easier 123 55 Page 8 of 35 Journal of Scientific Computing (2021) 89 :55 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)isdesignedin[11] by introducing a second-order parameterized perturbation and based on the inversion of Id+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 Id+Tby 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: Proposition 2 System (10) also reads ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∂th+∂x(hu)=0, ∂t(hu)+∂xhu2+gh2 2+h2¨ h 3+h2˙u 2∂xzb+h2u2 2∂2 xxzb +hg+¨ h 2+˙u∂xzb+u2∂2 xxzb∂xzb=−h∂xpatm. Here we used the standard notation ˙ ξ:=∂tξ+u∂xξ. Proof We deduce from (6) ˙ h=−h∂xu=2√3σ=2(w −u∂xzb), qb=h˙w=h¨ h 2+h[∂t(u∂xzb)+u∂x(u∂xzb)], q=qb 2+h¨ h 12 =h¨ h 3+h 2[∂t(u∂xzb)+u∂x(u∂xzb)].  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 is reformulated by means of the change of variable (h,hu)→ h,u+∂x(h2˙ h) 3h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 ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ∂th+∂x(hu)=0,(12a) ∂t(hX)+∂x(huX)+∇sgn Q=Sh,(12b) −∇sgn·∇sgn Q h=−∇ sgn·Sh h+∇sgn·(u∂xX). (12c) 123 Journal of Scientific Computing (2021) 89 :55 Page 9 of 35 55 Eq. (12c) reads ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 12 q h−h∂x∂x(hq) h−6qb h−h∂xqb h∂xzb=2h(∂xu)2+h∂xg∂x(zb+h)+∂xpatm, (13a) 4+(∂xzb)2qb h−6q h+∂xzb ∂x(hq) h=u2∂2 xxzb−g∂x(zb+h)+∂xpatm∂xzb.(13b) Well-posedness of this system is studied in [8]. Remark 2 Notice that only hand uare involved (not wnor σ) in the right hand side of (13). Moreover, we also notice that qbcan be expressed directly from (13b) qb h=1 4+(∂xzb)26q h−∂xzb ∂x(hq) h+u2∂2 xxzb−g∂x(zb+h)+∂xpatm∂xzb. (14) Inserting the latter expression into (13a) provides a unique equation for q: 121+(∂xzb)2 4+(∂xzb)2 q h−h∂x4 4+(∂xzb)2 ∂x(hq) h+6∂xzb 4+(∂xzb)2 ∂x(hq) h −h∂x6∂xzb 4+(∂xzb)2 q h =2h(∂xu)2+h∂xg∂x(zb+h)+∂xpatm +h∂x∂xzb 4+(∂xzb)2u2∂2 xxzb−g∂x(zb+h)+∂xpatm∂xzb. 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 −∇sgn·1 h∇sgn Q=f(15) for some right hand side f. The well-posedness of Eq. (15)isprovedin[8] by means of the Lax-Milgram theorem applied under some smoothness hypotheses on the water height h.In particular, straightforward computations show that −I∇sgn·1 h∇sgn Q· Qdx=−[q(∂x(hq)+qb∂xzb)]∂I+I 1 h∇sgn Q·∇sgn  Qdx, 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·Xthat 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. Proposition 4 Smooth solutions of System (10) satisfy the following energy equality provided patm and zbdo not depend on time: 123 55 Page 16 of 35 Journal of Scientific Computing (2021) 89 :55 since His a diagonal matrix with positive entries provided the hyperbolic scheme is positive. Then BQ=0⇒ Q=0as Bis one-to-one – see (21). More precisely, (24) provides the following discrete energy CQ,Q=#H−1qb,qb$+12 %H−1q−qb 2,q−qb 2& +#H−1(B11q+B12qb), B11q+B12qb$, which is consistent with I 1 h∗|∇sgn Q|2dx=I 1 h∗q2 b+12 q−qb 22+∂x(h∗q)+qb∂xzb2dx. 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: 1. C11qk+1=1 tf−C12qk b; 2. C22qk+1 b=1 tfb−CT 12qk+1. This method converges due to Lemma 1. This is the latter method which has been chosen in the present work. Matrix C11 is more or less the same matrix as for the DAE model (5). The extra additional time then corresponds to the iterative process (C11 is factorised once and C22 is diagonal). Remark 6 Notice that we could explicitly solve the second block of (23): qb=C−1 22 1 tfb−CT 12q⇒ (C11 −C12C−1 22 CT 12)q=1 tf−C12C−1 22 fb, 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 h+zb=cst.u=w=σ=q=qb=0. 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 1yields to the resolution of two homogeneous linear systems with non-singular matrices C11 and C22 •C11qk+1=0; •C22qk+1 b=0. Therefore q=qb=0,and thus the result follows. 16.3.4 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(zb+h)at the horizontal velocity equations are modified for emerging bottoms to avoid spurious pressure forces (see [19]). 123 Journal of Scientific Computing (2021) 89 :55 Page 17 of 35 55 (a) (b) Fig. 3 Multilayer description 2. For the computation of values including a term in the form 1 h,for very small values of the water height h,weusetheideapresentedin[44], replacing 1 hby √2h 'h4+max h4,ε , 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 C11 and C22 •C11qk+1=0; •C22qk+1 b=0. Therefore q=qb=0when the water thickness vanishes. 123 55 Page 18 of 35 Journal of Scientific Computing (2021) 89 :55 17 Numerical scheme for the LDNH2(L)Model 17.1 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 Lfor some integer L>0(seeFig.3). We consider the dispersive layer-averaged model derived in [35] under the following reformulation: ∂th+∂xhu=0,u= L ( α=1 αuα, and for α∈{1,...,L} ∂t(hαuα)+∂xhαu2 α+hαqα+uα+1/2α+1/2−∂xzα+1/2qα+1/2 −uα−1/2α−1/2+∂xzα−1/2qα−1/2=−ghα∂xη, ∂t(hαwα)+∂x(hαuαwα)+wα+1/2α+1/2+qα+1/2−wα−1/2α−1/2−qα−1/2=0, ∂t(hασα)+∂x(hασαuα)=2√3qα−qα+1/2+qα−1/2 2 −α+1/2hα∂xuα 12 +wα+1/2−wα 2+α−1/2hα∂xuα 12 +wα−wα−1/2 2, together with diagnostic equations 2√3σα+hα∂xuα=0,α∈{1,...,L}, wα−wα−1−(uα−uα−1)∂xzα−1/2−√3(σα−1+σα)=0,α∈{2,...,L}, w1−u1∂xzb−√3σ1=0, and the boundary condition qL+1/2=0. The mass transfer term is given by α+1/2= L ( β=α+1 ∂xhβuβ−u. To close the model, we define •uα+1/2=(1−γα+1/2)uα+1+γα+1/2uα; •wα+1/2=(1−γα+1/2)wα+1+hα+1 2∂xuα+1+γα+1/2wα−hα 2∂xuα; for any γα+1/2∈[0,1]such that [35, Prop. 1] γα+1/2−1 2α+1/2≥0.(25) Let us recall the properties of this system: 123 Journal of Scientific Computing (2021) 89 :55 Page 19 of 35 55 •Let us set Kα=u2 α+w2 α+σ2 α 2.If(H,uα, wα,qα)are smooth solutions to the multilayer model, we have under (25) ∂tL ( α=1 hαKα+gzα+patm(x)+∂xL ( α=1 hαuαKα+qα+gη+patm(x)≤0. (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). 17.2 Differential Operators Let us introduce the operators Xα=⎛ ⎝ uα wα σα⎞ ⎠,Qα=qα qα−1/2, and ∇α ldnh Qα=⎛ ⎜ ⎝ ∂x(h∗qα) L−qα+1/2∂xz∗ α+1/2+qα−1/2∂xz∗ α−1/2 qα+1/2−qα−1/2 −2√3qα−qα+1/2+qα−1/2 2⎞ ⎟ ⎠,α∈{1,...,L}, ∇α ldnh·Xα=2√3σα+h∗ L∂xuα wα−wα−1−(uα−uα−1)∂xz∗ α−1/2−√3(σα−1+σα),α∈{2,...,L}, ∇1 ldnh ·X1=2√3σ1+h∗ L∂xu1 w1−u1∂xzb−√3σ1. We get a global duality relation L ( α=1 Qα·(∇α ldnh·Xα)=1 L∂x(hqαuα)− L ( α=1 Xα·∇α ldnh Qα.(27) It is worth noticing that ∇α ldnh Qα=⎛ ⎝ ∂x(h∗qα) L 0 −2√3qα ⎞ ⎠+⎛ ⎝ qα−1/2∂xz∗ α−1/2 −qα+1/2 √3qα+1/2⎞ ⎠−⎛ ⎝ qα+1/2∂xz∗ α+1/2 −qα−1/2 √3qα−1/2⎞ ⎠+⎛ ⎝ 0 0 2√3qα+1/2⎞ ⎠, (28a) and ∇α ldnh·Xα=2√3σα+h∗ L∂xuα wα−uα∂xz∗ α−1/2−√3σα−0 wα−1−uα−1∂xz∗ α−1/2−√3σα−1 −0 2√3σα−1.(28b) 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 ∂xuαwhich is important for the order of derivatives. 123 55 Page 20 of 35 Journal of Scientific Computing (2021) 89 :55 Let us proceed as in the single layer case with a splitting strategy between hyperbolic terms and non-hydrostatic terms. The latter part reads ⎧ ⎨ ⎩ h∗Xα−h∗X∗ α Lt+∇α ldnh Qα=0, ∇α ldnh·Xα=0, (29) which implies −∇α ldnh·∇α ldnh Qα h∗=− 1 Lt∇α ldnh·X∗. The latter equation is expanded as, for α∈{1,...,L} −h∗∂x∂x(h∗qα) h∗+12L2qα h∗=−L t2√3σ∗ α+h∗ L∂xu∗ α −Lh∗∂xqα+1/2∂xz∗ α+1/2−qα−1/2∂xz∗ α−1/2 h∗+6L2qα+1/2+qα−1/2 h∗,(30a) then for α∈{2,...,L} 2−∂xz∗ α+1/2∂xz∗ α−1/2qα+1/2 h∗+24+(∂xz∗ α−1/2)2qα−1/2 h∗ +2−∂xz∗ α−1/2∂xz∗ α−3/2qα−3/2 h∗ =− 1 Ltw∗ α−w∗ α−1−(u∗ α−u∗ α−1)∂xz∗ α−1/2−√3(σ∗ α−1+σ∗ α) −∂xz∗ α−1/2 L ∂xh∗(qα−qα−1) h∗+6qα−1+qα h∗,(30b) and finally (2−∂xz∗ 3/2∂xz∗ 1/2)q3/2 h∗+4+(∂xz∗ 1/2)2q1/2 h∗ =− 1 Ltw∗ 1−u∗ 1∂xzb−√3σ∗ 1−∂xzb L ∂x(h∗q1) h∗+6q1 h∗.(30c) Let us remark is that the differential operator for qαin (30a) is independent from index α. Moreover, up to a coefficient L2, this operator is the same as in the monolayer case – see (18’). 17.3 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=B1/2from (21) into: Bα−1/2=B1+Bα−1/2 2,where B1=⎛ ⎜ ⎜ ⎝ B11 0 0 0 −2√3IN0 ⎞ ⎟ ⎟ ⎠and Bα−1/2 2=⎛ ⎜ ⎜ ⎝ 0Bα−1/2 12 0−IN 0√3IN ⎞ ⎟ ⎟ ⎠, B11 is defined in equation (22)andBα−1/2 12 is a N×Ndiagonal matrix with entries (∂xz∗ α−1/2)i. The fact that B1is independent from layer α(and exactly the same as in the monolayer case) is crucial for what is following. 123 Journal of Scientific Computing (2021) 89 :55 Page 21 of 35 55 Hence ∇α ldnh Qαis approximated according to (28a) B1Qα+Bα−1/2 2Qα−Bα+1/2 2Qα+1+RQα+1=Bα−1/2Qα−Bα+1/2 2Qα+1+RQα+1, where R=⎛ ⎜ ⎜ ⎝ 0 0 0 0 02√3IN ⎞ ⎟ ⎟ ⎠and Qα= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ qα,1 . . . qα,N qα−1/2,1 . . . qα−1/2,N ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . 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 ⎛ ⎝H/(Lt)B+R B+RT0⎞ ⎠X Q=HX∗/(Lt)− 0  0. where B= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ B1/2−B3/2 20 ...... Bα−1/2−Bα+1/2 2 ...... 0BL−1/2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ∈M3NL,2NL(R), R=⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 0R0 ...... ...R 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠∈M3NL,2NL(R), X= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ . . . Uα Wα α . . . ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,Uα=⎛ ⎜ ⎝ uα,1 . . . uα,N ⎞ ⎟ ⎠,Q=⎛ ⎜ ⎜ ⎝ . . . Qα . . . ⎞ ⎟ ⎟ ⎠. It is easy to verify that B+RTX= 0is consistent with (28b). For the treatment of boundary conditions which are incorporated in vectors 0and  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 B+RTH−1B+RQ=1 LtB+RTX∗− 0−B+RTH−1 0.(31) Let us set 123 55 Page 22 of 35 Journal of Scientific Computing (2021) 89 :55 •C=B+RTH−1B+R∈M2NL,2NL(R)which is blockwise tridiagonal, symmetric positive-definite; •Cα−1/2=Bα−1/2TH−1Bα−1/2which has exactly the same structure as Cin (24) with B12 replaced by Bα−1/2 12 . The blockwise components of Care then: •C1,1=C1/2; •For α∈{2,...,L}: Cα,α =Cα−1/2+R−Bα−1/2 2T H−1R−Bα−1/2 2T=⎛ ⎝C11 Cα−1/2 12 Cα−1/2 12 T 2Cα−1/2 22 ⎞ ⎠ and Cα−1,α =CT α,α−1=Bα−3/2TH−1R−Bα−1/2 2 =⎛ ⎝0−6H−1−BT 11H−1Bα−1/2 12 02H−1−Bα−3/2 12 T H−1Bα−1/2 12 ⎞ ⎠. (31) is a consistent discretisation of (30). 17.4 Iterative Scheme To solve the 2 NL×2NLlinear system (31), we apply an alternating direction-type method. More precisely, it amounts to solving iteratively: x-direction For each layer α,(30a)inqαknowing qα+1/2which corresponds to a linear system with matrix C11 given in (24) (the same for all layers which requires a single factorisation for all iterates at each time step): C11qp+1 α=1 LtBT 11U∗ α−2√3∗ α−Cα−1/2 12 qp α−1/2; z-direction For each node xi,(30b-30c)inqα+1/2knowing qα. It is a tridiagonal system with matrix S(i)solved by means of the Thomas’ algorithm [64]. Matrix S(i)is a L×L tridiagonal symmetric positive-definite matrix with •S(i) 1,1=4+(∂xz∗ 1/2)2 i; •For α≥2, S(i) α,α =24+(∂xz∗ α−1/2)2 iand S(i) α,α−1=S(i) α−1,α =2−(∂xz∗ α−1/2)i (∂xz∗ α−3/2)i. Indeed, we check that S(i)x,x=2x2 1+2+(∂xzL−1/2)2 ix2 L + L−1 ( α=24x2 α+2(xα+xα−1)2+(∂xz∗ α−1/2)ixα−(∂xz∗ α−3/2)ixα−12≥0. 123 Journal of Scientific Computing (2021) 89 :55 Page 23 of 35 55 This strategy is equivalent to solving the following system by means of a Gauss-Seidel iterative procedure: ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ C11 0 ...E 0C11 S(1)0 ET... 0S(N) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ where unknowns have been re-labelled as ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ˆ qα=1 . . . ˆ qα=L ˇ qi=1 . . . ˇ qi=N ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,with ˆ qα=⎛ ⎜ ⎝ qα,1 . . . qα,N ⎞ ⎟ ⎠and ˇ qi=⎛ ⎜ ⎝ q1/2,i . . . qL−1/2,i ⎞ ⎟ ⎠. Remark 8 It must be underlined that this algorithm can be easily parallelised insofar as each direction (xand z) involves a blockwise diagonal matrix. 18 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 LDNH 2(L=1)and the resolution of (31)forLDNH 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 LDNH 0and LDNH 2by 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. 18.1 Convergence Test Let us start by a convergence test for both models, LDNH 0and LDNH 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 123 55 Page 24 of 35 Journal of Scientific Computing (2021) 89 :55 Table 1 Test 18.1 –L1errors with the LDNH 0-soliton and numerical orders of accuracy No. of cells hhuhw L1error Order L1error Order L1error Order 50 3.04e-03 0.00 1.81e-03 0.00 2.17e-03 0.00 100 9.26e-04 1.72 5.17e-04 1.81 5.97e-04 1.86 200 2.57e-04 1.85 1.39e-04 1.89 1.88e-04 1.67 400 6.96e-05 1.88 4.22e-05 1.72 8.44e-05 1.15 Table 2 Test 18.1 –L1errors with the LDNH 2-soliton and numerical orders of accuracy No. of cells hhuhw L1Error Order L1Error Order L1Error Order 50 1.21e-02 0.00 1.72e-02 0.00 1.22e-02 0.00 100 4.51e-03 1.42 4.02e-03 2.10 3.28e-03 1.89 200 1.33e-03 1.76 1.12e-03 1.85 8.50e-04 1.95 400 3.80e-04 1.80 3.24e-04 1.79 2.39e-04 1.83 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 h(t,x)=H∗+Asech21 H∗)Aγ 2(A+H∗)(x−ct),u(t,x)=c1−H∗ h(t,x), w(t,x)=− cA h(t,x))Aγ 2(A+H∗)sech31 H∗)Aγ 2(A+H∗)(x−ct) sinh 1 H∗)Aγ 2(A+H∗)(x−ct) (32) where Aand H∗are constant fixed values, c=√g(A+H∗)and •γ=2 in the case of LDNH 0, •γ=3/2 in the case of LDNH 2. Here we set H∗=1,A=0.1,g=1andzb=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 1and 2. These results show the convergence towards the reference solution at second-order. 123 Journal of Scientific Computing (2021) 89 :55 Page 25 of 35 55 Table 3 Test 18.2 – Well-balancing L1errors for the still water steady states with the LDNH 2(1)system No. of cells ηuw 50 2.27e-18 1.07e-17 2.62e-18 100 9.99e-18 3.28e-17 1.26e-17 200 7.77e-18 2.99e-17 1.49e-17 400 2.62e-17 4.89e-17 4.52e-17 Table 4 Test 18.2 – Well-balancing L1errors for the still water steady states with the LDNH 2(4)system, and u=1 4 4 ( α=1 uα,w=1 4 4 ( α=1 wα No. of cells ηuw 50 5.66e-18 6.50e-18 4.91e-18 100 6.66e-18 1.07e-17 8.98e-18 200 1.58e-17 1.26e-17 1.01e-17 400 2.84e-17 4.08e-17 3.25e-17 18.2 Still Water Steady-State In this test we consider the same domain and boundary conditions as in the previous test. The topography zbis given by zb(x)=0.9e−x2−1. We consider the initial condition given by h0(x)=−zb(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 LDNH 2(L)system for L=1,4 with increasing number of cells. Tables 3and 4show the error observed at time T=10, with a CFL 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. 18.3 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 4shows a sketch of the laboratory experiments carried out at the O.H. Hinsdale Wave Research Laboratory of Oregon State University. See for instance [66] for more details. The 1D domain [0,45]is discretised with x=0.045 m. A solitary wave of amplitude 0.5 m is placed at point x=10 m. Finally CFL =0.9and g=9.81 m ·s−2. Free outflow boundary conditions are imposed. Figure 5shows snapshots at different times, t√g/H∗=t0where H∗=1 m. Comparisons between experimental and simulated data allow to validate the numerical approach presented in this paper. Results are shown for LDNH 0(L=1)and LDNH 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 LDNH 2provides slightly better results. 123 55 Page 32 of 35 Journal of Scientific Computing (2021) 89 :55 (a) (b) Fig. 11 Test 18.5 with A=0.289 – blue lines indicate the numerical solution with models LDNH 0aand LDNH 2bfor different layers; red points represent the experimental measurements at the probe positions 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 LDNH 2is in many situations more accurate than LDNH 0. This result shows that the LDNH 2is of interest since the accuracy of both models has been compared by using a generalisation of classic projection techniques. One should recall that LDNH 2has twice more pressure unknowns than LDNH 0.Therefore, if we apply a projection method to approximate the solution of LDNH 2, it results in a more expensive algorithm, from the computational point of view, than for LDNH 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 LDNH 2model and compare them with the one proposed here. For example, it may be of interest to consider 123 Journal of Scientific Computing (2021) 89 :55 Page 33 of 35 55 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. Acknowledgements The research of this paper was partially supported by the Spanish Government (SG), the European Regional Development Fund (ERDF), the Regional Government of Andalusia (RGA), and the University of Málaga (UMA) through the projects of reference RTI2018-096064-B-C21, RTI2018-096064- B-C22 (SG-ERDF), UMA18-Federja-161 (RGA-ERDF-UMA), and P18-RT-3163 (RGA-ERDF). Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. References 1. Aïssiouene, N., Bristeau, M.-O., Godlewski, E., Mangeney, A., Parés, C.: A two-dimensional method for a family of dispersive shallow water models. SMAI J. Comp. Math. 6, 187–226 (2020) 2. Aïssiouene, N., Bristeau, M.-O., Godlewski, E., Sainte-Marie, J.: A combined finite volume-finite element scheme for a dispersive shallow water system. Netw. Heterog. Media 11(1), 1–27 (2016) 3. Allgeyer, S., Bristeau, M.-O., Froger, D., Hamouda, R., Jauzein, V., Mangeney, A., Sainte-Marie, J., Souillé, F., Vallée, M.: Numerical approximation of the 3d hydrostatic Navier-Stokes system with free surface. ESAIM Math. Model. Numer. Anal. 53(6), 1981–2024 (2019) 4. Antuono, M., Liapidevskii, V., Brocchini, M.: Dispersive Nonlinear Shallow-Water Equations. Stud. Appl. Math. 122(1), 1–28 (2009) 5. Bai, Y., Cheung, K.F.: Dispersion and nonlinearity of multi-layer non-hydrostatic free-surface flow. J. Fluid Mech. 726, 226–260 (2013) 6. Bai, Yefei, Cheung, Kwok Fai: Dispersion and kinematics of multi-layer non-hydrostatic models. Ocean Model. 92, 11–27 (2015) 7. Bassi, C., Bonaventura, L., Busto, S., Dumbser, M.: A hyperbolic reformulation of the Serre-Green-Naghdi model for general bottom topographies. Comp. Fluids 212, 104716 (2020) 8. Benyo, K., Charhabil, A., Debyaoui, M.-A., Penel, Y.: Simulation of complex free surface flows. Proc. of the 2019 session of CEMRACS 9. Bonnet-Ben Dhia, A., Bristeau, M., Godlewski, E., Imperiale, S., Mangeney, A., Sainte-Marie, J.: Pseudocompressibility, dispersive model and acoustic waves in shallow water flows. Preprint. HAL Id: hal- 02493518, (2020) 10. Bonneton, P., Barthelemy, E., Chazel, F., Cienfuegos, R., Lannes, D., Marche, F., Tissier, M.: Recent advances in Serre-Green Naghdi modelling for wave transformation, breaking and runup processes. Eur. J. Mech. B. Fluids 30(6), 589–597 (2011) 11. Bonneton, P., Chazel, F., Lannes, D., Marche, F., Tissier, M.: A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model. J. Comput. Phys. 230(4), 1479–1498 (2011) 12. Bouchut, F.: Nonlinear stability of finite volume methods for hyperbolic conservation laws and wellbalanced schemes for sources. Birkhauser, (2004) 13. Boussinesq, J.: Théorie des ondes et des remous qui se propagent le long d’un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond. J. Math. Pures Appl., pages 55–108, (1872) 14. Bristeau, M.-O., Di Martino, B., Guichard, C., Sainte-Marie, J.: Layer-averaged Euler and Navier-Stokes equations. Commun. Math. Sci. 15(5), 1221–1246 (2017) 15. Bristeau, M.-O., Mangeney, A., Sainte-Marie, J., Seguin, N.: An energy-consistent depth-averaged Euler system: Derivation and properties. Discret. Contin. Dynam. Syst. - Series B 20(4), 961–988 (2015) 123 55 Page 34 of 35 Journal of Scientific Computing (2021) 89 :55 16. Bristeau, M.-O., Sainte-Marie, J.: Derivation of a non-hydrostatic shallow water model; Comparison with Saint-Venant and Boussinesq systems. Discret. Contin. Dynam. Syst. - B 10(4), 733–759 (2008) 17. Busto, S., Dumbser, M., Escalante, C., Favrie, N., Gavrilyuk, S.: On high order ADER discontinuous galerkin schemes for first order hyperbolic reformulations of nonlinear dispersive systems. J. Sci. Comput. 87(2), 48 (2021) 18. Castro, M.J., Morales de Luna, T., Parés, C.: Well-Balanced Schemes and Path-Conservative Numerical Methods. In Handbook of Numerical Analysis, volume 18, pages 131–175. Elsevier, (2017) 19. Castro, M.J., Ferreiro, A.M. Ferreiro, García-Rodríguez, J.A., González-Vida, J.M., Macías, J., Parés, C., Vázquez-Cendón, M. Elena: The numerical treatment of wet/dry fronts in shallow flows: application to one-layer and two-layer systems. Math. Comp. Model., 42(3):419 – 439, 2005 20. Castro Díaz, M.J., Fernández-Nieto, E.: A class of computationally fast first order finite volume solvers: PVM methods. SIAM J. Sci. Comput. 34(4), A2173–A2196 (2012) 21. Casulli, V.: A semi-implicit finite difference method for non-hydrostatic, free-surface flows. Int. J. Numer. Methods Fluids 30(4), 425–440 (1999) 22. Casulli, V.: A semi-implicit numerical method for the free-surface navier-stokes equations. Int. J. Numer. Meth. Fluids 74(8), 605–622 (2014) 23. Chazel, F., Lannes, D., Marche, F.: Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model. J. Sci. Comput. 48(1), 105–116 (2011) 24. Clamond, D., Dutykh, D.: https://es.mathworks.com/matlabcentral/fileexchange/39189-solitary-water- wave, (2012) 25. Cravero, I., Puppo, G., Semplice, M., Visconti, G.: CWENO: uniformly accurate reconstructions for balance laws. Math. Comp. 87(312), 1689–1719 (2018) 26. Barré de Saint-Venant, A.: Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et à l‘introduction des marées dans leurs lits. C. R. Acad. Sci. 73, 237–240 (1871) 27. Dingemans, M.W.: Comparison of computations with Boussinesq-like models and laboratory measurements. Report H-1684.12, 32, Delft Hydraulics, (1994) 28. Duran, A., Marche, F.: Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations. Commun. Comput. Phys. 17(3), 721–760 (2015) 29. Dutykh, D., Clamond, D.: Efficient computation of steady solitary gravity waves. Wave Motion 51(1), 86–99 (2014) 30. Escalante, C., Castro, M.J., Fernández-Nieto, E.D., Morales de Luna, T.: An efficient two-layer nonhydrostatic approach for dispersive water waves. J. Sci. Comput. 79, 273–320 (2019) 31. Escalante, C., Dumbser, M., Castro, M.J.: An efficient hyperbolic relaxation system for dispersive nonhydrostatic water waves and its solution with high order discontinuous Galerkin schemes. J. Comput. Phys. 394, 385–416 (2019) 32. Escalante, C., Morales de Luna, T.: A general non-hydrostatic hyperbolic formulation for Boussinesq dispersive shallow flows and its numerical approximation. J. Sci. Comput. 83(3), 62 (2020) 33. Escalante, C., Morales de Luna, T., Castro, M.J.: Non-hydrostatic pressure shallow flows: GPU implementation using finite volume and finite difference scheme. Appl. Math. Comput. 338, 631–659 (2018) 34. Favrie, N., Gavrilyuk, S.: A rapid numerical method for solving Serre- Green-Naghdi equations describing long free surface gravity waves. Nonlinearity 30(7), 2718–2736 (2017) 35. Fernández-Nieto, E.D., Parisot, M., Penel, Y., Sainte-Marie, J.: A hierarchy of dispersive layer-averaged approximations of Euler equations for free surface flows. Commun. Math. Sci. 16(05), 1169–1202 (2018) 36. Ferrari, S., Saleri, F.: A new two-dimensional Shallow Water model including pressure effects and slow varying bottom topography. ESAIM: Math. Model. Numer. Anal., 38(2):211–234, (2004) 37. Gerbeau, J.-F., Perthame, B.: Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation. Discrete Contin. Dyn. Syst. Ser. B 1(1), 89–102 (2001) 38. Godlewski, E., Raviart, P.-A.: Numerical approximation of hyperbolic systems of conservation laws. Number 118 in Applied Mathematical Sciences. Springer-Verlag, (1996) 39. Green, A.E., Naghdi, P.M.: A derivation of equations for wave propagation in water of variable depth. J. Fluid Mech. 78(2), 237–246 (1976) 40. Grosso, G., Antuono, M., Brocchini, M.: Dispersive nonlinear shallow-water equations: some preliminary numerical results. J. Eng. Math. 67(1), 71–84 (2010) 41. Guermond, J.-L., Minev, P., Shen, J.: An overview of projection methods for incompressible flows. Comput. Methods App. Mech. Eng. 195(44–47), 6011–6045 (2006) 42. Guibourg, S.: Modèlisation numèrique et expèrimentale des houles bidimensionnelles en zone cotière. PhD thesis, Universitè Joseph Fourier, (1994) 43. Kirby, J.T.: Boussinesq models and their application to coastal processes across a wide range of scales. J. Waterway, Port, Coastal, Ocean Eng. 142(6), 03116005 (2016) 123 Journal of Scientific Computing (2021) 89 :55 Page 35 of 35 55 44. Kurganov, Alexander, Petrova, Guergana: A second-order well-balanced positivity preserving centralupwind scheme for the saint-venant system. Commun. Math. Sci. 5(1), 133–160 (2007) 45. Lannes, D.: The Water Waves Problem: Mathematical Analysis and Asymptotics, volume 188. Amer. Math. Soc., (2013) 46. Lannes, D., Bonneton, P.: Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation. Phys. Fluids 21(1), 016601 (2009) 47. Lannes, D., Marche, F.: A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations. J. Comput. Phys. 282, 238–268 (2015) 48. Le Métayer, O., Gavrilyuk, S., Hank, S.: A numerical scheme for the Green-Naghdi model. J. Comput. Phys. 229, 2034–2045 (2010) 49. LeVeque, Randall J.: Finite Volume Methods for Hyperbolic Problems. Cambridge Texts in Applied Mathematics. Cambridge University Press, (2002) 50. Lynett, P.J., Liu, P.L.-F.: Linear analysis of the multi-layer model. Coastal Eng. 51(5), 439–454 (2004) 51. Ma, G., Shi, F., Kirby, J.T.: Shock-capturing non-hydrostatic model for fully dispersive surface wave processes. Ocean Model. 43, 22–35 (2012) 52. Madsen, P.A., Murray, R., Sørensen, O.R.: A new form of the Boussinesq equations with improved linear dispersion characteristics. Coastal Eng. 15(4), 371–388 (1991) 53. Marche, F.: Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects. Eur. J. Mech. B. Fluids 26(1), 49–63 (2007) 54. Mazaheri, A., Ricchiuto, M., Nishikawa, H.: A first-order hyperbolic system approach for dispersion. J. Comput. Phys. 321, 593–605 (2016) 55. Nwogu, O.: Alternative form of Boussinesq equations for nearshore wave propagation. J. Waterway, Port, Coastal, Ocean Eng. 119(6), 618–638 (1993) 56. Parisot, M.: Entropy-satisfying scheme for a hierarchy of dispersive reduced models of free surface flow. Int. J. Numer. Methods Fluids 91(10), 509–531 (2019) 57. Peregrine, D.H.: Long waves on a beach. J. Fluid Mech. 27(04), 815–827 (1967) 58. Ricchiuto, M., Filippini, A.G.: Upwind residual discretization of enhanced boussinesq equations for wave propagation over complex bathymetries. J. Comput. Phys. 271, 306–341 (2014) 59. Roeber, V., Cheung, K.F., Kobayashi, M.H.: Shock-capturing boussinesq-type model for nearshore wave processes. Coastal Eng. 57, 407–423 (2010) 60. Sainte-Marie, J.: Vertically averaged models for the free surface non-hydrostatic Euler system: derivation and kinetic interpretation. Math. Models Methods Appl. Sci 21(03), 459–490 (2011) 61. Sánchez, C.E.: Mathematical and numerical modelling of dispersive water waves. PhD thesis, Universidad de Málaga, (2018) 62. Serre, F.: Contribution à l‘ètude des ècoulements permanents et variables dans les canaux. La Houille Blanche 8(6–7), 830–872 (1953) 63. Stelling, G., Zijlema, M.: An accurate and efficient finite-difference algorithm for non-hydrostatic freesurface flow with application to wave propagation. Int. J. Numer. Meth. Fluids 43(1), 1–23 (2003) 64. Thomas, L.: Elliptic problems in linear differential equations over a network. Columbia Univ., NY, Technical report, Watson scientific computing laboratory (1949) 65. Yamazaki, Y., Kowalik, Z., Cheung, K.F.: Depth-integrated, non-hydrostatic model for wave breaking and run-up. Int. J. Numer. Methods Fluids 61(5), 473–497 (2009) 66. Yao, Y., He, T., Deng, Z., Chen, L., Guo, H.: Large eddy simulation modeling of tsunami-like solitary wave processes over fringing reefs. Nat. Hazards Earth Syst. Sci. 19, 1281–1295 (2019) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123