Non-hydrostatic layer-averaged approximation of Euler system with enhanced dispersion properties

A new family of non-hydrostatic layer-averaged models for the non-stationary Euler equations is presented in this work, with improved dispersion relations. They are a generalisation of the layer-averaged models introduced in Fernández-Nieto et al. (Commun Math Sci 16(05):1169–1202, 2018), named LDNH models, where the vertical profile of the horizontal velocity is layerwise constant. This assumption implies that solutions of LDNH can be seen as a first order Galerkin approximation of Euler system. Nevertheless, it is not a fully (x, z) Galerkin discretisation of Euler system, but just in the vertical direction (z). Thus, the resulting model only depends on the horizontal space variable (x), and therefore specific and efficient numerical methods can be applied (see Escalante-Sanchez et al. in J Sci Comput 89(55):1–35, 2021). This work focuses on particular weak solutions where the horizontal velocity is layerwise linear on z and possibly discontinuous across layer interfaces. This approach allows the system to be a second-order approximation in the vertical direction of Euler system. Several closure relations of the layer-averaged system with non-hydrostatic pressure are presented. The resulting models are named LIN-NHk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_k$$\end{document} models, with k=0,1,2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0,1,2$$\end{document}. Parameter k indicates the degree of the vertical velocity profile considered in the approximation of the vertical momentum equation. All the introduced models satisfy a dissipative energy balance. Finally, an analysis and a comparison of the dispersive properties of each model are carried out. We show that Models LIN-NH1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_1$$\end{document} and LIN-NH2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document} provide a better dispersion relation, group velocity and shoaling than LDNH models.


Introduction
Many efforts have been devoted to including dispersive effects in models in fluid mechanics, and in particular, they have been widely studied in the context of shallow flows in the literature. There are essentially two approaches to consider dispersive effects in this framework: Boussinesq type systems and non-hydrostatic systems.
Boussinesq type systems are mainly based on classic shallow water system, whose unknowns are the total height of the fluid and the horizontal velocity. The system is then extended by introducing high order derivatives of the variables. The two pioneering works on Boussinesq models were proposed by Boussinesq (1872) and Peregrine (1967). Some other models were proposed afterwards: Madsen et al. (1991); Madsen and Sørensen (1992), Nwogu (1993), Serre (1953), Green and Naghdi (1976), Schäffer and Madsen (1995), Lannes and Bonneton (2009), among many others. For a review on Boussinesq type dispersive models, see Lannes (2013), Kirby (2016).
Contrary to Boussinesq systems, non-hydrostatic models only consider first order derivatives of the variables but with a larger number of unknowns. Usually, in addition to the water depth and horizontal velocity, we find the vertical velocity, non-hydrostatic pressure, and other unknowns related to the vertical profile of the different variables. One of the essential points of these models is the presence of constraint equations related to the incompressibility condition, see for example Yamazaki et al. (2008) and Bristeau et al. (2015).
Although both approaches seem quite different, in many cases they are equivalent, as it was shown in Escalante and Morales de Luna (2020). They introduced a general formulation that relates some well-known classic Boussinesq systems and non-hydrostatic models. Concretely, it was shown that many of these classic Boussinesq systems might be written as non-hydrostatic models. We want to point out that the non-hydrostatic formulation has several advantages from the numerical point of view. In particular, avoiding high order derivatives makes the numerical approximation easier. Moreover, the treatment of boundary conditions is also more straightforward.
A strategy to improve the dispersion relation of such systems is to increase the accuracy of the vertical profile for both velocity and pressure components. For this purpose, several attempts have also been introduced in the literature. The so-called multilayer approach considers non-material vertical interfaces, dividing the fluid into virtual vertical layers. Different profiles for the velocity and pressure within each layer may be assumed. Notice that the term multilayer has been used sometimes for stratified flow models, where one considers a constant density within each layer. This has no relation with the approach we are considering in this paper. To avoid any possible confusion, in Fernández-Nieto et al. (2018) it was proposed to name the former (non-stratified flows) as layer-averaged models. We shall do the same in what follows.
In Sainte-Marie (2011), a non-hydrostatic layer-averaged model is proposed, which corresponds to the extension of the hydrostatic model in Audusse et al. (2011). The momentum equations are approximated by considering a constant vertical profile for horizontal and vertical velocity components. Then, velocity may be discontinuous at the interfaces, although the pressure profile is supposed to be continuous. However, a linear profile of the vertical velocity within each layer is considered to approximate the incompressibility condition, which is a compatible condition of strong solutions in each layer for this equation.
A similar model was proposed in Bai and Cheung (2013), which assumes a constant profile for the horizontal velocity and linear profiles for the vertical velocity and pressure. The common ground for both references is that the horizontal velocity has a constant vertical profile within each layer and may be discontinuous at the interfaces. A different assumption is found in Lynett and Liu (2004a), where a layer-averaged model is derived by using a continuous global profile of the horizontal velocity, which is quadratic within each layer.
In Fernández-Nieto et al. (2018), a hierarchy of models is presented, with an associated energy balance, where again a constant profile is assumed for the horizontal velocity inside each layer. At the same time, several degrees of freedom are introduced for the vertical velocity and pressure, accounting for the vertical approximation of such variables. In that work, the linear dispersion relation was studied for the different models proposed. It was shown that the dispersion relation converges to the exact dispersion relation for Euler equations when the number of layers goes to infinity. A numerical strategy to solve these models is proposed in Escalante-Sanchez et al. (2021). Let us recall the notation of the models introduced in that paper as it will be helpful in what follows. These models are named as LDNH k where LDNH stands for Layerwise Discretisation Non-Hydrostatic, and k is the order of approximation. More precisely, LDNH 0 may be seen as a generalisation of models Yamazaki et al. (2008) and Bristeau et al. (2015), while LDNH 2 is as a generalisation of the Serre-Green-Naghdi model. LDNH 1 is an in-between intermediate model. Similarly, in Cantero-Chinchilla et al. (2018) authors derived a weighted-averaged non-hydrostatic pressure model, under the assumption of linear horizontal velocity, whereas the vertical velocity and the non-hydrostatic pressure was assumed to have a quadratic profile. This paper focuses on the derivation of 2D layer-averaged models where the vertical profile of the horizontal velocity is supposed to be layerwise linear and discontinuous at the interfaces. The extension to 2D is straightforward. Concretely, a family (with respect to the degree of approximation of the vertical velocity and pressure) of novel models is introduced, named as LIN-NH k models for k = 0, 1, 2, where the index k corresponds to the degree of approximation of the vertical velocity in the vertical momentum equation. Therefore, three models are proposed, based on three different profiles for the pressure, which are related to the approximation of the vertical momentum equation and the incompressibility condition. In addition, the pressure profile is considered as a polynomial of degree k + 1. All the proposed non-hydrostatic models also satisfy a dissipative energy balance. As we show in this paper, since more unknowns are introduced, linear dispersion relation, group velocity, and linear shoaling are greatly improved for such models. In particular, the LIN-NH 1 and LIN-NH 2 models exhibit excellent results for these dispersive properties. Moreover, another contribution is a general procedure to compute the group velocity and shoaling gradient from the dispersion relation for the wave celerity.
Although a general formulation based on a vertical Galerkin approach, i.e. polynomials with an arbitrary degree describing the vertical profile of the variables, is possible from the theoretical point of view (see Garres-Díaz et al. 2023 for the hydrostatic case), we remains here in the case of first order polynomials for the horizontal velocity. Going to the general case in this non-hydrostatic framework would lead to excessively complicated models, where we would also have an arbitrary (high) number of unknowns associated to the pressure needed, and therefore also of constraints related to the incompressibility condition. Moreover, to our point of view, this model, that is expected to be second order accurate in the vertical direction, is a good compromise between simplicity and accuracy, being already an important improvement with respect to previous multilayer non-hydrostatic systems where the horizontal velocity is supposed to be always constant (and therefore first-order accurate in the vertical direction).
In addition, the models proposed in this work can be seen as second order Galerkin vertical discretisations of Euler system. The main difference with respect to a fully Galerkin discretisation is the fact that we obtain a model depending only on the horizontal space variable. As a consequence, more efficient numerical methods may be designed. Actually, a specific and efficient numerical method was proposed for the layer-averaged LDNH models in Escalante-Sanchez et al. (2021). Here, we only focus on the deduction of the LIN-NH models and its properties. Taking into account the key role of designing efficient methods for the proposed models, it is out of the scope of this paper and must be addressed in the future.
The paper is organised as follows. In Sect. 2, we state the notation and the layer-averaged procedure for the layerwise linear horizontal velocity case. Section 2.1 focuses on the discrete spaces and hypotheses over the vertical profile of each variable (horizontal and vertical velocities and pressure), as well as the normal jump condition at the interfaces associated with the weak formulation of the Euler system. In Sect. 2.2 we describe the vertical averaging of the mass and horizontal momentum equations without detailing the pressure terms. Section 3 is devoted to derive the cascade of LIN-NH k non-hydrostatic models, for k = 0, 1, 2. The dispersive properties for these models are analysed in Sect. 4, and some conclusions are presented in Sect. 5. Finally, Appendix A contains a summary of all the models LIN-NH introduced in this work and the models LDNH proposed in Fernández-Nieto et al. (2018), including the layerwise vertical profile of the unknowns, and the coefficients of the dispersion relation associated to each model.

Initial system and layer-averaged approach
This section introduces the general settings associated with the layer-averaged approach for the Euler equations. In this work we shall assume a piecewise linear approximation of the horizontal velocity.
This layer-averaged approach can be seen as a technique to approximate the solution of the Euler equations in the framework of the Discontinuous Galerkin methods. The final model that we obtain is a system of partial differential equations. The solution of this system may be seen as a particular weak solution of Euler equations, in the sense that it corresponds to an approximated piecewise smooth weak solution that may be discontinuous at the internal interfaces. We remark that these interfaces are not physical or material interfaces. They represent a virtual decomposition or partition of the domain vertically. The procedure is as follows: first, we consider a vertical discretisation. Second, a layer-averaging approach is considered assuming an appropriate structure of the weak solutions. In particular, here, we shall assume that the velocity vector has horizontal components that are linear in the vertical z-direction within each layer. Third, some closure relations (or constraints) are needed, and they will be presented later in Sect. 3.1.
For the sake of simplicity, we shall consider here a 2-dimensional space, where the space variables are denoted by (x, z) ∈ R 2 . Nevertheless, it could be extended easily to the 3D case. In what follows, ∇ = (∂ x , ∂ z ) represents the usual differential operator. Let us consider now an incompressible fluid with constant density ρ ∈ R. Let us denote g ∈ R the gravity acceleration, p ∈ R the pressure, U = (u, w) ∈ R 2 the velocity vector, and a moving domain where the topography z b is assumed to be independent from time and H the water height (see Fig. 1). The Euler system is written, for (x, z) ∈ (t), as where g = (0, −g) ∈ R 2 and I is the identity matrix. Notice that we write the time derivative of the density in previous system, despite of assuming the case of an incompressible fluid. It is done for the sake of clarity and generality when writing the normal jump conditions associated to the layer-averaging approach. The system is completed with initial and boundary conditions, and the following kinematic conditions: Now, to apply the layer-averaged approach, the domain is subdivided along the vertical direction into L ∈ N layers with thickness h α (t, x), which are denoted by α (t). The layers are separated by L +1 interfaces L α+1/2 (t), defined by z = z α+1/2 (t, x) for α = 0, 1, . . . , L. The interfaces are assumed sufficiently smooth, at least C 1 (see Fig. 1). More explicitly, Notice then that where z b = z 1/2 and h α = z α+1/2 − z α−1/2 . Moreover, the total height of the fluid is equal to Let us also denote by z α the midpoint of the layer α , i.e., Furthermore, let us introduce the following notations: for any function f (t, x, z), we define: • Its approximations at the interfaces L α+1/2 for α = 0, . . . , L In the case of continuous functions, we shall simply write f α+1/2 = f + α−1/2 = f − α+1/2 ; • The average value over the layer α • The linear average of f over the layer α • The average value through the interface L α+1/2 • The jump across the interface L α+1/2 Note that, in general, f α = f α , but the equality holds if f is constant or linear in z over the layer α .

Layerwise approximation of the variables and normal jump conditions
This subsection deals with the main hypotheses upon velocity and pressure fields in the layer-averaged framework. In particular, we detail here the vertical profile for the variables considered in this work and its notation. We also write the normal jump condition associated with the weak formulation.

Choice of discrete spaces for velocity and pressure fields
The layer-averaged models presented in Fernández-Nieto et al. (2013, 2016 assume that the horizontal velocity is constant in the z-direction within each layer. In practice, this will formally limit us to first order approximations of the solution for the Euler equations. In this work we shall consider a piecewise linear profile in z for the horizontal velocity. Concerning the other variables, a piecewise parabolic profile in z is supposed for the vertical velocity, and a piecewise third order polynomial in z is used for the non-hydrostatic pressure. Nevertheless, some other alternative simplified equations will be presented afterwards, resulting in simplified models. More explicitly, let us specify the discretisation assumptions for the velocity unknowns where u α and w α are the horizontal and vertical components of the velocity field in layer α . Then, the following profiles and notation are chosen for this model: Horizontal velocity: a linear profile on the vertical direction is assumed within each layer where u α and λ α are the averaged horizontal velocity and the slope respectively in the z direction. In what follows and for the sake of simplicity, we shall not write the dependence of the variables on (t, x) explicitly unless necessary. It follows then that the limit at the interfaces u − α+1/2 and u + α−1/2 are obtained by evaluating (6) at z α+1/2 and z α−1/2 , respectively, i.e., Vertical velocity: focusing on the incompressibility condition and given the piecewise linear profile (6) for the horizontal velocity, a compatibility criterion leads us to consider a piecewise parabolic polynomial in z, that is Hence, there are three parameters to set. On the one hand, the integration of (8) over the layer yields Note that w α (z α ) = w α in this case because of the contribution of the second-order term. On the other hand, a second-order Taylor expansion, which is exact for a quadratic function, yields As we want w α to approximate w | α , let us mimic the incompressibility constraint ∂ z w = −∂ x u using (6). Hence we impose Therefore, combining (8), (9) and (10), we get for α = 1, . . . , L, with Evaluating (11) at z α±1/2 , we obtain Pressure: the total pressure is decomposed into a hydrostatic part g (z b + H − z), and a non-hydrostatic counterpart q: where we assume that the pressure is known at the surface, which is usually set to zero.
As mentioned above, we choose a layerwise linear horizontal velocity and a layerwise parabolic vertical velocity (potentially discontinuous across the interfaces). Consequently, the pressure must be a layerwise cubic function due to the vertical momentum equation. We assume in the present work that the pressure is continuous across interfaces.

Jump conditions and evolution equation for the layer midpoints
Following Fernández-Nieto et al. (2013), the model is deduced by looking for a particular piecewise smooth weak solution (ρ, U, p) of system (1). More explicitly, we search for a solution that satisfies the weak formulation of the system for a particular set of tests functions. We refer to Fernández-Nieto et al. (2013) for further details. In particular, they correspond to a classic, smooth solution inside each layer α , while it should satisfy the normal flux jump conditions across the interfaces L α+1/2 , α = 0, . . . , L. The jump condition for the mass conservation law reads and for the momentum conservation law We are considering here a constant density for the fluid. Hence, following Fernández-Nieto et al. (2013), (15a) implies where α+1/2 is the mass transfer term through the interface L α+1/2 given by In particular, using previous equation for α±1/2 and making the average, we infer that which gives the evolution of the midpoint at each layer α .

Layer-averaged approximation: mass and horizontal momentum equations
In this section we obtain the mass and horizontal momentum equations of the target model by a layer-averaged process. In what follows, we use a general depth-integration process for mass and horizontal momentum equations.

Mass conservation
The incompressibility condition is integrated over each layer, which leads to Taking into account the definition of the mass transfer terms (16), we obtain the mass conservation laws where we recall that α±1/2 account for mass transfer across interfaces L α±1/2 . Actually, they can be expressed in terms of the velocities by combining previous equations, getting, for α defined in (4). Moreover, summing Eq. (18) over α we get Therefore, it is equivalent to considering Eq. (18) and the set of equations defined by (19) and (20). Note that 1/2 and L+1/2 correspond to the mass exchange at the free surface and the bottom, respectively, and they are set to zero unless specified. In any case, both of them should be provided as data.
Averaged horizontal momentum conservation. The approximation of the horizontal momentum equation is deduced in two steps. First, a layer-integration of the horizontal momentum equation provides the horizontal momentum conservation law equation related to the averaged velocity, u α .
see (5) for notations. Note that here we have used the jump condition (15b), otherwise a term u ∓ α±1/2 would appear instead of u α±1/2 . This is done in what follows when the limit values of the horizontal and vertical velocities at the interfaces appear.
Secondly, we evaluate the mean deviation of the horizontal velocity and use (6) to derive the conservation law associated with λ α . To do so, the horizontal momentum equation is multiplied by (z − z α ) and then integrated inside the layer, yielding to for α = 1, . . . , L.
Note that the integrals involving the pressure in (21) and (22) have not been specified. These terms will be computed later on depending on the assumption on the pressure profile. The simplest choice corresponds to the assumption of hydrostatic pressure, and this model is presented in the following subsection. Non-hydrostatic pressure models are deduced in Sect. 3.

LIN-H model: linear horizontal velocity and hydrostatic pressure
Let us mention the case of hydrostatic pressure, where the resulting model is a particular case of the multilayer-moment model presented in Garres-Díaz et al. (2023). In that case, the final system is defined by (18), (21), (22), where, as it is usual in the framework of hydrostatic shallow flows, the pressure is defined by Now, the integrals associated to the pressure in (21) and (22) are Therefore, the final system for hydrostatic pressure is defined by the following set of equations: This model verifies exactly a dissipative energy balance, stated by the following theorem.
Theorem 1 System (23) satisfies the dissipative energy balance The proof can be seen as a particular case of the result stated in Theorem 2.

Remark 1
The hyperbolicity of layer-averaged models is still an open question. For the case of layerwise constant horizontal velocity, the hyperbolicity of the two-layer model (L = 2) was proven in Aguillon et al. (2018). For L > 2 it is an open problem, up to our knowledge. Actually, following (Garres-Díaz et al. 2023, Remark 6), a loss of hyperbolicity could arise. Moreover, in Garres-Díaz and Bonaventura (2021) it is concluded that the hyperbolicity of these layer-averaged systems might be lost in the case of strong vertical profiles, that is, when large differences between the layer velocities (u α ) and the averaged one (u) appear. Notice that in this case the mass transfer terms ( α+1/2 ) may become large. Nevertheless, in these previous works (Garres-Díaz and Bonaventura 2021; Garres-Díaz et al. 2023), no loss of hyperbolicity was found in the numerical simulations performed. Focusing now on the polynomial expansion of the velocity within each layer, the onelayer case corresponds to the so-called moment approach, giving as results the Shallow Water Moment models (see Kowalski and Torrilhon 2018). It was proven that only the case of linear velocity (system (23) with L = 1) leads to a hyperbolic model, whereas the resulting system in the general case (polynomials u ∈ P k>1 [z]) is not hyperbolic already for the onelayer case (see Koellermeier and Rominger 2020). Therefore, the same would hold for the multilayer case (L > 1 and u α ∈ P k>1 [z]). This is also one of the reasons for not going to this general case in this work.
authors introduced a hierarchy of non-hydrostatic layer-averaged models for Euler equations for d = 0 and k ≤ 1.
In the following subsections, firstly, the incompressibility condition within each layer is deduced. Secondly, three models are proposed, corresponding to different approximations of the vertical velocity and non-hydrostatic pressure counterparts. All of them satisfy a dissipative energy balance. We will detail the more complex model, where the vertical velocity is a layerwise quadratic polynomial according to the incompressibility condition, and therefore the non-hydrostatic pressure is a layerwise cubic function. Later, we see that if the vertical structure of the vertical velocity is simplified just for the vertical momentum equation, then the model degenerates to simpler models, but keeping still an analogous energy balance. This family of non-hydrostatic models will be denoted by Multilayer Horizontal Linear discretisation Non-Hydrostatic k model (from now on LIN-NH k ), where k is the degree of the polynomial approximating the vertical velocity in the vertical momentum conservation equation, that is, k = 0, 1, 2.

Averaged incompressibility condition
In the non-hydrostatic framework, we need some extra information or constraints to solve the resulting system, which now includes the degrees of freedom related to the non-hydrostatic pressure. Note that for each layer, there are three unknowns that define the non-hydrostatic pressure profile, namely q α+1/2 , q α and π α . Therefore, one would need three constraints at each layer.
On the one hand, two restrictions are associated to the relationship between the first and second order vertical derivative of the vertical and horizontal velocities. Concretely, in Sect. 2.1 the profile for the vertical velocity was obtained, where from the incompressibility inside each layer we got (12).

LIN-NH 2 model: non-hydrostatic model with w˛∈ P 2 [z]
In what follows we shall detail the first non-hydrostatic model that we propose here. It will be denoted by LIN-NH 2 model. Its final formulation is (31) below with the subsequent notations. Later on, in Sect. 3.3 we shall introduce as well two extra simplified models, named LIN-NH 1 and LIN-NH 0 models. Let us recall that, up to this point, we have evolution Eqs. (21) and (22) together with the restrictions (25) and the jump conditions (19) and (20).
In order to close the system, some extra equations are needed, which are obtained by successive integration on the vertical momentum conservation equation, as shown in what follows.
Let us first focus on the horizontal momentum Eqs. (21) and (22). In these equations, we compute, in the non-hydrostatic framework, the explicit pressure terms Let us focus now on the equations for the variables involved in the vertical velocity. To do so, we consider a basis in P 2 [z] of test functions and compute the averaged equations For the sake of simplicity, we omit here the details of the computations. Let us introduce the following notations: obtaining hence the LIN-NH 2 model, which reads for α = 1, . . . , L, combined with the following constraints, where In system (28) the velocities at the interfaces are given by and, therefore, The system is thus composed of 8 L + 1 equations for 8 L + 1 unknowns, Note that thanks to the vertical partition (4), using Eq. (18) for h α an explicit expression for α+1/2 is achieved, and then neither h α nor α+1/2 are unknowns anymore, unlike the total height H .
Finally, let us rewrite system (28) and constraints (29) in a compact way which will have the same structure for all models LIN-NH k : , where we introduced some notations, namely the vectors of unknowns and the following differential operators, for α ∈ {1, . . . , L} Then for α ∈ {2, . . . , L} we define These definitions of the NH-gradient and the NH-divergence operators satisfy the following duality relation, Remark 2 Notice that setting α = 0, α = 0 and π α = 0, system (31) reduces to the L DN H 2 (L) model described in Escalante-Sanchez et al. (2021) where the assumptions were: a layerwise-constant horizontal velocity, a layerwise-linear vertical velocity and a quadratic pressure.
An interesting property of such a system is that it satisfies a dissipative energy balance, as stated by the following theorem.

Theorem 2 The proposed LIN-NH 2 model, defined by (31) satisfies the dissipative energy balance
Proof The proof relies on algebraic relations. Let us multiply (31b) by X α . Notice that Given the expressions of interface velocities (30), we get Summing the last equality over all layers, we get which completes the proof.
Notice that thanks to the compact form (31) and the duality relation (32), the proof of Theorem 2 has been notably shortened. Let us also remark that in the case of smooth enough solutions, then the equality would hold in Theorem 2 instead of an inequality. This is also true for all the energy results in this work.

1-Layer model
In the 1-layer case, the model may be written in Boussinesq form by expressing all variables in terms of the averaged horizontal velocity and the fluid height as the only unknowns of the system. As we have said previously, this will result in a system of equations that contains third order derivatives in space and time. See for example, Fernández-Nieto et al. (2018), where it is shown that the LDNH 2 model with 1 layer coincides with the Serre-Green-Naghdi model. More explicitly, the model proposed in this work for 1-layer case may be written as follows, Note that by a subsequent substitution of the different unknowns of the system into the second and third equations, we easily obtain a system formed by three equations with the unknowns H , u and , which includes terms of third order derivatives in space and time.

Simplified non-hydrostatic models
We propose now to obtain two simplified models in order to reduce the complexity of the previous one. The main objective is to reduce the number of unknowns related to the pressure and consequently the number of constraints. From a numerical point of view, this would be important to reduce the computational cost. For example, if a projection method is considered, it will be needed to solve a linear system for the pressure where the number of equations is proportional to the number of pressure variables and layers. Then, for the simplified models, we consider a more straightforward structure for the vertical velocity and, therefore, the nonhydrostatic pressure. We remark that it is essential that, to obtain models with an associated exact energy balance, we must always consider the incompressibility equation without any simplification on the vertical velocity profile.
We mention that among the two simplified models presented below, named LIN-NH 1 and LIN-NH 0 , the one corresponding to a layerwise constant vertical velocity (LIN-NH 0 ) does not produce good results concerning the dispersive properties. However, we find interesting to show these results for the sake of completeness in the study, since it is natural to wonder whether it is a reasonable simplification or not. A possible reason for the poor results for that model could be that there should be a relation between the degrees of the polynomials approximating horizontal and vertical velocities. This fact would suggest that the degree of the polynomial approximation for the vertical velocity should be at least the one for the horizontal velocity, and it cannot be not be smaller.
The derivation of these simplified models can be obtained by following a similar approach as the one carried out in Sect. 3.2 for LIN-NH 2 model. This results in a similar compact form (31). Hence analogous energy balances as in Theorem 2 are obtained. For the sake of brevity, we do not give here these proofs, which can be easily carried out following the analogous steps described in the proof of Theorem 2.

LIN-NH 1 model: non-hydrostatic model with w˛∈ P 1 [z]
We consider first that the vertical velocity is a layerwise linear polynomial in z, w α ∈ P 1 [z] to approximate the vertical momentum equation, that is, Concerning the pressure, we consider a piecewise parabolic function, q α ∈ P 2 [z], given by Note that this pressure profile coincides with (14) by setting π α = (δq) α . Then, the unknowns related to the non-hydrostatic pressure are q α , and q α−1/2 for α = 1, . . . , L. Recall that q L+1/2 is a given data, usually assumed to vanish at the surface. The system thus comprises 6 L + 1 equations for 6 L + 1 unknowns, Analogously to what is done for the LIN-NH 2 model, the LIN-NH 1 model can be written in a compact form (31), where in this case with The following NH-gradient operator read in this case, for α ∈ {1, . . . , L}, and for α = 1, These differential operators verify in this case the following duality relation, This simplified system satisfies a dissipative energy balance, that coincides with the result stated at Theorem 2 by setting α = 0 and π α = (δq) α . We obtain in this case the following result.

Theorem 3 The proposed LIN-NH 1 model satisfies the dissipative energy balance
with X α is defined by (35).
The proof is analogous to LIN-NH 2 model.

LIN-NH 0 model: non-hydrostatic model with w˛∈ P 0 [z]
Finally, we consider even a simpler case where a layerwise constant vertical velocity, w α ∈ P 0 [z], is assumed in the vertical momentum equation. Therefore, we have the structure for the vertical velocity, and the pressure, which is linear q α ∈ P 1 [z], is given by Now, we have only one pressure unknown by layer, q α−1/2 . Making such assumptions is equivalent to considering q α = q α = q α+1/2 + q α−1/2 2 for α = 1, . . . , L.
LIN-NH 0 model can be written in a compact form (31), where in this case with with the following NH-gradient operator for α ∈ {1, . . . , L} and the NH-divergence operator, for α ∈ {2, . . . , L} and for α = 1, These differential operators verify in this case the following duality relation, This simplified system satisfies a dissipative energy balance, that coincides with the result stated at Theorem 2 by setting α = α = 0, π α = (δq) α , q α = q α . Furthermore, in Tables  3 and 4 in Appendix A, the NH-divergence and NH-gradient operators for LIN-NH and LDNH models are summarized. In these tables, we also observe how these operators, in the case of simplified models, are obtained through appropriate simplifications in the operators associated to LIN-NH 2 model.

Dispersion relations
In this section, we study the linear dispersion relations of the LIN-NH k systems, following the approach used in Escalante et al. (2019), Lynett and Liu (2004b). In particular, we focus on the linear dispersion relations for the wave celerity, the group velocity and the linear shoaling. In other words, a formal study is performed for the PDE systems introduced previously, focusing on the propagation of dispersive waves.
To do so, a flat bottom is assumed as usual (z b constant), and the governing equations are linearised around a steady-state solution. Then, a standard Stokes-type Fourier analysis is used to obtain the wave and group velocities. Finally, a shoaling analysis of the linearised equations is carried out.

Linear dispersion relation for the LIN-NH 2 model
Let us begin with the LIN-NH 2 model that is described by (28). Let us linearise the system around the steady state and let us consider η = z b + H the free-surface. We consider then the following asymptotic expansion and for any variable φ α ∈ {u α , λ α , w α , ϕ α , ψ α }, For the sake of simplicity, we shall neglect in this section the superindex (1) and bars ( · ) so that the notation is less cumbersome. Using this linearisation in (28), we shall neglect O( 2 ) terms and keep the system at first order: Remark that combining (42c), (42d), and (42i), we obtain Moreover, given that q L+1/2 = 0 and (δq) α = q α+1/2 − q α−1/2 we have Let us consider now a plane-wave solution of the linearised system in the form ( η, u α , λ α , w α , ϕ α , δq α , π α ) T e i(kx−ωt) , where here the hat notation ( · ) is new and is not related to (5b). Inserting the plane wave into (43) we get Doing so in (42e) gives The same may be done in the other equations and we obtain, using the relation Then (44) may be written as where ·, · represents the scalar product in R L and, performing a progressive substitution from bottom to top, we obtain This means that, provided A 2,L is invertible, In the particular case α = 1 L , A 2,L reduces to and assuming e, U = 0 we get Therefore, the following dispersion relation for the wave celerity c := ω k is obtained for LIN-NH 2 model

Linear dispersion relation for the LIN-NH 1 model
In the particular case of LIN-NH 1 , the equations that govern the system are those corresponding to (28) under the assumption of (δq) α = π α , as well as the restrictions (25). Doing a linearisation as in (40) and (41) for φ α ∈ {u α , λ α , w α , ϕ α }, keeping the first order terms and after inserting the plane wave solution ( η, u α , λ α , w α , ϕ α , δq α ) T e i(kx−ωt) , we obtain system (44) where, in this case, d α = 2 α H 2 0 k 2 12 . That leads us, for the particular case α = 1 L , to the following dispersion relation for the wave celerity c := ω k : where A 1,L is similar to (45), but with a change on the value of d, more precisely:

Linear dispersion relation for the LIN-NH 0 model
Remark that LIN-NH 0 is quite similar to LIN-NH 1 model, where some equations are removed and This means that the linearisation (40) and (41) for φ α ∈ {u α , λ α , w α }, keeping the first order terms, will result in a reduced version of (42). Then, following a similar procedure as in the previous subsections, we get that a plane wave solution of this reduced linearised system should satisfy that is the same as (44) except for the term −ωk 2 α H 2 0 12 u α , which is now missing in the last equation. As a consequence, for the particular case of α = 1/L, we will obtain a matrix A 0,L similar to A 1,L in (47), but with the only difference that the last term (1 + d) I is replaced by I. That is, the following dispersion relation for the wave celerity c := ω k for LIN-NH 0 model is satisfied:

Computation of the linear group velocity and the linear shoaling of dispersive systems
Let us introduce a general procedure that allows to compute the linear group velocity and shoaling for dispersive systems presented in this work. Moreover, we provide a very general methodology that can be applied to compute the group velocity and the linear shoaling of any dispersive PDE system once the wave celerity c 2 (k, H 0 ) of the system is known.

Group velocity
The group velocity c g (k, H 0 ) is essentially computed by taking the derivative of the wave celerity c, and is defined as Let us express the ratio of the wave celerity of a given dispersive system and the shallow water celerity as a function of k H 0 : For instance, for the LIN-NH 1 model, f is given in (46) as: Therefore, the group velocity can be written as follows where f means the derivative of the function f w.r.t. its argument.

Linear shoaling
The linear shoaling gradient γ is a non-dimensional quantity to measure the change in wave height in the presence of a bottom slope. It was firstly introduced in Madsen et al. (1991) and can be expressed as As usual, we will assume the dependency on the x−direction for H 0 , η and the wave number k. On the contrary, we suppose that the frequency ω does not depends on x, and therefore ∂ x ω = 0. The linear shoaling gradient can be determined for a given dispersive PDE system by assuming the constancy of the energy flux ∂ x η 2 c g = 0, and, therefore, Now, from (49), it yields where again ϒ stands for the derivative of the function ϒ w.r.t. its argument. Given the definition of the celerity c = ω k and the constancy of ω with respect to x, it follows that We can see again that both relative dispersion errors converge to 0 as the number of layers L increases. Nevertheless, the dispersion relation for the new model LIN-NH 2 is now much better than that of LDNH 0 . Even for the simple case of one layer, the relative dispersion error is relatively small. Therefore, we can state that LIN-NH 2 is an excellent choice in terms of dispersion accuracy. Despite the increased complexity of the model, it gives better results for large values of k H 0 even for a small number of layers.
In the same way, we compare in Figs. 5 and 6 the group velocity, and in Figs. 7 and 8 the shoaling gradient for the different models considered in this work. We remark the great performance of the LIN-NH 1 and LIN-NH 2 models over LIN-NH 0 , LDNH 0 and LDNH 2 .
In order to complete the dispersion study, we provide in Table 1  | f (s)| . Fig. 3 Comparison of the relative dispersion errors for the wave celerity for each model and different number of layers. Notice that the range for the y-axis in the subfigures corresponding with models LIN-NH 1,2 is smaller to better distinguish the results

Fig. 4
Comparison between the relative dispersion errors for the wave celerity of different models, for a fixed number of layers L = 1, 2, 3 Note that we consider relative errors for the wave celerity and group velocity. However, absolute errors for the shoaling gradient are considered to avoid singularities when dividing by γ Airy . It can be highlighted the great performance of the model LIN-NH 2 even for large ranges of the parameter k H 0 . Despite the complexity of the presented models, the analysis of the dispersive properties for the LIN-NH 2 model are outstanding. With no more than 2 layers, it can accurately simulate dispersive water waves of very high frequency.  Notice that the range for the y-axis in the subfigures corresponding with models LIN-NH 1,2 is smaller in order to better distinguish the results

Fig. 6
Comparison between the relative dispersion errors for the group velocity of different models, for a fixed number of layers L = 1, 2, 3

Conclusions
A family of layer-averaged non-hydrostatic models for Euler equations with layerwise linear horizontal velocities, named LIN-NH k with k = 0, 1, 2, has been presented. The derivation of these models follows from a layer-averaging procedure of the mass and momentum conservation equations, where normal jump conditions associated with the weak formulation of the problem are taken into account. In the notation LIN-NH k , the value k indicates the closure relation chosen for the vertical profile of the vertical velocity and pressure. Thus, k stands for the degree of the polynomial approximation of the vertical velocity, and k + 1 is the degree for the pressure.
Once the most complex model (LIN-NH 2 ) is derived, simplified LIN-NH 0,1 models are obtained from simple assumptions upon the vertical profile of the vertical velocity, and  therefore the pressure, just on the vertical momentum equation. It is important to remark that for the incompressibility constraints, the vertical profile of the velocity is not simplified for any model. For this very equation, the vertical velocity is always a layerwise parabolic polynomial. It is essential to prove a dissipative energy balance for all the proposed models. All the models are written in a compact form, which allows to prove the dissipative energy balance in a few lines. It is done by defining some non-hydrostatic differential operators, which satisfy a crucial duality relation. It is a remarkable issue, especially when designing numerical schemes for such models (see Escalante-Sanchez et al. 2021).
An analysis and a comparison of dispersion properties have been performed. In particular, we computed the linear dispersion, shoaling, and group velocity. Interestingly, models LIN-NH 1 and LIN-NH 2 provide accurate dispersion relations, although LIN-NH 1 is a simpler model. Note that the latter model is computationally cheaper since it has 2L equations and unknowns less than LIN-NH 2 model. LIN-NH 1,2 models exhibit excellent dispersion properties for a wide range of wave numbers even for few layers. Then, these models can accurately reproduce dispersive water waves of very high frequency. We have also shown that LIN-NH 1,2 models improve the results of the models introduced in Fernández-Nieto et al. (2018). On the other hand, the LIN-NH 0 model, which relies on a layerwise constant vertical velocity assumption, provides poor results in the linear dispersion analysis. That would suggest that there should be a restriction between the degrees of the polynomial approximation assumed for the different variables. In particular, we could induce that the degree of the vertical velocity approximation should be greater than for the horizontal velocity. The extension to the Navier-Stokes case for a general rheology, accounting for the contribution of the deviatoric stress tensor, is addressed in the second part of this paper. In the future, it would be also interesting to propose efficient numerical schemes for LIN-NH 1,2 models, following the strategy developed in Escalante-Sanchez et al. (2021) for LDNH models.
In this work, three layerwise models are presented. They are closely related with previous models LDNH 0,2 introduced in Fernández- Nieto et al. (2018). A comparison of all these models, including models LDNH 0,2 , has been performed in Sect. 4.5 in terms of dispersion properties.
For the sake of clarity, we include in what follows a summary of the details of these models by means of several tables. More precisely: • In Table 2 we show the selected approximation discrete spaces for each model, the vertical profiles of the velocity, non-hydrostatic pressure terms and dispersion relations.
Concerning the dispersion relations of these models, this is done by means of two variables (ζ 1 , ζ 2 ) which, following the notation introduced in Sect. 4, take the form ω 2 = gk 2 H 0 L A −1 e, e , with A = I + k 2 H 2 0 L 2 ζ 1 I + ζ 2 The dispersion relation for each model is then recovered by simply replacing the corresponding values of ζ 1 , ζ 2 given in Table 2.

Table 3
Summary of variables related to the velocity (X α ) field and its NH divergence (∇ NH · X α )  Table 3, the unknowns of the model related to the velocity, X α , are explicitly stated together with the non-hydrostatic divergence operator, ∇ NH · X α . • In Table 4, the unknowns related to the non-hydrostatic terms, Q α , and the non-hydrostatic gradient operator, ∇ NH Q α , are shown.