An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves

An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves C. Escalante1 · E. D. Fernández-Nieto2 · T. Morales de Luna3 · M. J. Castro1 Abstract In this paper, we propose a two-layer depth-integrated non-hydrostatic system with improved dispersion relations. This improvement is obtained through three free parameters: two of them related to the representation of the pressure at the interface and a third one that controls the relative position of the interface concerning the total height. These parameters are then optimized to improve the dispersive properties of the resulting system. The optimized model shows good linear wave characteristics up to kH ≈ 10, that can be improved for long waves. The system is solved using an efficient formally second-order well-balanced and positive preserving hybrid finite volume/difference numerical scheme. The scheme consists of a two-step algorithm based on a projection-correction type scheme. First, the hyperbolic part of the system is discretized using a Polynomial Viscosity Matrix path-conservative finite-volume method. Second, the dispersive terms are solved using finite differences. The method has been applied to idealized and challenging physical situations that involve nearshore breaking. Agreement with laboratory data is excellent. This technique results in an accurate and efficient method. Keywords Dispersive waves · Non-hydrostatic · Shallow-water · Finite-volume · Finite-difference · Breaking waves 1 Introduction When modelling and simulating geophysical flows, the non-linear shallow-water equations, from now on SWE, is often a good choice as an approximation of the Navier–Stokes equations. Nevertheless, SWE does not take into account effects associated with dispersive waves. In recent years, a great effort has been made in the derivation of relatively simple mathematical models for shallow water flows that include long non-linear water waves. As computational power increases, Boussinesq type models [1,15,40,52,58,59,67–69] become more accessible. BC. Escalante [email protected] 1Universidad de Málaga, Málaga, Spain 2Universidad de Sevilla, Seville, Spain 3Universidad de Córdoba, Córdoba, Spain This means that one can use more sophisticated models to accurately describe reality, despite the higher computational cost. Moreover, to improve the non-linear dispersive properties of the model, information on the vertical structure of the flow should be included. The Boussinesq-type wave equations have prevailed due to their computational efficiency. The main idea is to include non-hydrostatic effects due to the vertical acceleration of the fluid in the depth-averaging process of the equations. For instance, one can assume that both non-linearity and frequency dispersion are weak and of the same order of magnitude. Since the early works of Peregrine [59], several improved and enhanced Boussinesq models have been proposed over the years: Madsen and Sørensen [52], Nwogu [58], Serre Green–Naghdi equations [40], and nonlinear and non-hydrostatic higher order shallow-water type models [16,22,34,71] among many others. One may use different approaches to improve non-linear dispersive properties of the models: considering a Taylor expansion of the velocity potential in powers of the vertical coordinate and in terms of the depth-averaged velocity [52] or the particle velocity components (u,w) at a chosen level [58]; including two scalars representing the vertical profile of the non-hydrostatic pressure [42]; using a better flow resolution in the vertical direction with a multi-layer approach [46,47,70]. Weakly-dispersive weakly-nonlinear Boussinesq systems ([16,52,58,71] among others) may offer good dispersive properties such as phase, group velocity and shoaling for kH values up to 3, where k denotes the wave number and H the typical depth. In intermediate waters, this range could not be enough, and much more information on the vertical structure must be retained. High order Boussinesq equations can offer better linear and non-linear dispersive properties. The counterpart is that extremely complex systems with high order derivatives arise (fifth order derivatives in [38]for kH up to 6), requiring an equally complex numerical scheme. The complexity increases for a two-horizontal dimension problem. The development of non-hydrostatic models for coastal water waves has been the topic of many studies over the past 15 years. Non-hydrostatic models are capable of solving many relevant features of coastal water waves, such as dispersion, non-linearity, shoaling, refraction, diffraction, and run-up. For example, in [22], Casulli splits the pressure into a hydrostatic and a non-hydrostatic part. Then, the system is approximated with a fractional step method, wherein the first step the hydrostatic free-surface flow model is approximated and in a second step a Poisson equation is solved. The decomposition of the pressure into a hydrostatic and non-hydrostatic parts is also considered in many other works. This decomposition has been considered for example by Ma et al. [50] and Gallerano et al. [37], where the system is reformulated in σ -coordinates. Stansby and Zhou [63] also considered this decomposition, where the main difference is that the σ -coordinate is not considered for the horizontal non-hydrostatic pressure gradients, to avoid the numerical truncate errors introduced by this approximation. Zijlema and Stelling, in [24,64,72], propose a multi-layer method based in a layered depth-averaging of the system. A similar procedure has been considered by Ai and Jin (see [4,5]). The main differences reside in the horizontal gradient of the non-hydrostatic pressure component, which is assumed to be independent of the depth at each layer. Another multi-layer model has been proposed in [10] where a numerical study of the linear dispersive celerity shows convergence to the results given by the linear Stokes theory. The linear shoaling properties of this model are studied in [11]. In [34], a hierarchy of layer-averaged multi-layer models has been proposed based on the hypothesis made for the vertical velocity and pressure profiles. Moreover, an explicit linear dispersion relation for each one is done. This allows proving analytically that when the number of layer increases, celerity converges to the exact one given by Euler equations in the Stokes linear theory. Bai and Cheung [9] proposed a modification of the depth-integrated two-layer model introduced in [8]. They include a free-parameter that relates the pressures at each of the two layers. This parameter allows improving the linear dispersion relation. Cui et al. introduced in [27] a two-layer model based in a depth-averaging of the system written in σ-coordinates. They introduced two free-parameters that relate the height of the two layers and the nonhydrostatic pressures unknowns of the two-layer system. The choice of these parameters is also related to the accuracy of the linear dispersive properties of the system. The paper is organized as follows. In Sect. 2we present a new two-layer depth-averaging system depending on three different free-parameters that are optimized to improve the dispersive relation of the system. The main hypothesis for the derivation is the assumption of a constant vertical profile of the horizontal velocity within each layer. This corresponds to the assumption of a shallow water regime as it is usual when deriving non-hydrostatic multilayer systems. The proposed two-layer system can be seen as a modification of one of the models presented in [34] with a correction in the non-hydrostatic pressure profile. This will leads to a system that improves the main dispersive properties of other two-layer non-hydrostatic pressure systems proposed in the literature (see for instance [27]), and can even improve the dispersive properties of some multi-layer systems with 5 layers (see for instance [11]). We recall that the resulting system is a model which is less expensive from the computational point of view while keeping good dispersive relations. Thus, the final model can be applied to intermediate waters for a wide range of waves for a kH up to 10. This is discussed through a linear analysis in Sect. 3, as well as numerical tests. It is worth mentioning that this approach leads to a system that does not present those aforementioned high order derivatives. In Sect. 4a numerical scheme is introduced based on a two-step projection-correction type scheme. In the first step we solve the hyperbolic hydrostatic system and in the second one, we include the non-hydrostatic effects. The proposed numerical scheme present a novel treatment of wet–dry fronts, that is very robust, as it will be shown in Sect. 6. No special treatment is required for the non-hydrostatic pressure since in the presence of wet–dry fonts it vanishes automatically thanks to the rewriting of the incompressibility equations. This is an improvement on non-hydrostatic pressure numerical schemes, where usually non-hydrostatic pressure is truncated to zero up to a threshold value. The reader should keep in mind that detailed small-scale breaking driven physics are not described by the model. This means that one has to include some breaking mechanism in the depth-integrated equations. In Sect. 5, a new robust and straightforward breaking mechanism is presented. The breaking mechanism consists in the extension to the two-layer system of the viscous mechanism used in [33]. The breaking criteria used here to switch on/off the breaking mechanism is a simplification of the criteria presented in [61] and is one of the simplest that can be found in the literature. Nevertheless some other and more sophisticated breaking criteria could be considered and implemented similarly. Finally, in Sect. 6, some numerical tests including comparisons with analytical solutions from the Stokes linear theory and laboratory data are shown. 2 The Non-hydrostatic Two-Layer System In this section, a two layer-averaged system derived from Euler equations is presented. This model is obtained by setting a non-material interface that separates two layers with different velocities and pressures. In the deduction of the equations is assumed that the horizontal velocity has a piecewise constant vertical profile, by while the vertical velocity and the nonhydrostatic counterpart pressure are piecewise linear, being the three unknowns discontinuous at the interface. The 2D (x,z)Euler system reads ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∂xu+∂zw=0, ∂tu+∂xu2+pT+∂z(uw)=0, ∂tw+∂x(uw)+∂zw2+pT=−g, (1a) (1b) (1c) where xand zdenote the horizontal Ox and the vertical Oz axis respectively. We consider this system for t>t0,x∈R,−H(x,t)≤z≤η(x,t), where ηis the unknown water elevation, His the bathymetry that can vary in space and time, uand ware the horizontal and vertical velocities. The water height is h=η+H. The model is completed with boundary conditions at the free surface ∂tη(x,t)+u(x,η(x,t), t)∂xη(x,t)−w(x,η(x,t), t)=0, pT((x,η(x,t), t)) =patm, (2a) (2b) where patm is the atmospheric pressure. At the bottom, no-penetration boundary condition is imposed u(x,−H(x,t), t)∂xH+w(x,−H(x,t), t)+∂tH=0.(3) The total pressure pTis supposed to be decomposed into a sum of a hydrostatic and nonhydrostatic part: pT=patm +g(η −z)+p(4) where p(x,z,t)is the non-hydrostatic pressure. Hereinafter, the atmospheric pressure will be supposed to be zero and the non-hydrostatic pressure is assumed to vanish at the free surface p(x,η(x,t), t)=0.(5) Given a positive water height h(x,t),weconsiderl1,l2∈[0,1]such that l1+l2=1, and we decompose the fluid along the vertical axis into two virtual layers of height hα=lαh,α=1,2. In what follows, α=1 stands for the lower layer and α=2 for the upper one (see Fig. 1). Let us denote L1(x,t)={z:−H(x,t)≤z≤zI(x,t)},L2(x,t)={z:zI(x,t)≤z≤η(x,t)}, where zI(x,t)=h1−H.zαwill denote the level of the middle point of the layers: z1(x,t)=zI−H 2,z2(x,t)=zI+η 2. For a given function (x,z,t)→f(x,z,t), in the case that it is discontinuous at the interface zI, let us denote fI+(x,t)=lim z→zI(x,t) z>zI(x,t) f(x,z,t), fI−(x,t)=lim z→zI(x,t) z<zI(x,t) f(x,z,t), Fig. 1Bilayer setting as well as fb≡f(x,−H(x,t), t), fη(x,t)≡f(x,η(x,t), t). As we mentioned before, to derive the model we assume a given vertical profile on the unknowns of the problem. First, a piecewise constant profile of the horizontal velocity is considered. We denote by u1(x,t)and u2(x,t)the horizontal velocities at layers L1(x,t) and L2(x,t), respectively. This is the main hypothesis for the derivation of the model, which corresponds to a shallow water regime. In fact, if we take p=0in(4)andl1=1,we can derive in the same way the SWE. From previous assumption and the incompressibility condition (1a), we obtain that the vertical velocity has a piecewise linear profile. Let us recall that, using the notation introduced before, its limits at the interface level z=zI(x,t)are denoted by wI±. The vertical velocity at the bottom and the free surface are denoted by wband wη, respectively. Finally, we also consider a piecewise linear profile of the non-hydrostatic pressure, which can be determined from its atmospheric value, which is supposed to be zero, its limits at the interface and the pressure at the bottom. Thus we have pη=0, and denote its limits at the interface level z=zI(x,t)by pI±, and the pressure at the bottom, pb(x,t). We also denote by p1(x,t)and p2(x,t)the non-hydrostatic pressure part evaluated at the middle level of the corresponding layer, that is, p1(x,t)=p(x,z1(x,t), t), p2(x,t)=p(x,z2(x,t), t). From the incompressibility condition (1a), we get that a weak solution which is discontinuous at z=zI(x,t)must verify the following jump condition (see [34,35]) wI+−wI−=(u2−u1)∂xzI.(6) Let us denote by w1and w2the vertical velocity at the middle level of the corresponding layer, that is, w1(x,t)=w(x,z1(x,t), t), w2(x,t)=w(x,z2(x,t), t). Then, the incompressibility condition (1a) inside each layer reads, ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂xu1+2w1−wb h1=0, ∂xu2+2w2−wI+ h2=0. The integration over each layer of the incompressibility Eq. (1a), combined with conditions (2a), (3), gives the mass equations at each layer, ⎧ ⎪ ⎨ ⎪ ⎩ ∂th1+∂x(h1u1)=−(wI−−(∂tzI+u1∂xzI)) , ∂th2+∂x(h2u2)=(wI+−(∂tzI+u2∂xzI)) . (8a) (8b) Let us define ΓI−=wI−−∂tzI−u1∂xzI, ΓI+=wI+−∂tzI−u2∂xzI. Using the jump condition (6), we get ΓI+=ΓI−and we may define ΓI=wI++wI− 2−∂tzI−u1+u2 2∂xzI.(9) Then Eq. (8) becomes ⎧ ⎪ ⎨ ⎪ ⎩ ∂th1+∂x(h1u1)=−ΓI, ∂th2+∂x(h2u2)=ΓI. (10a) (10b) Thus, ΓIcan be interpreted as an approximation of the mass transfer across the interface zI. Moreover, by combining Eqs. (10a)and(10b) we obtain ΓI=l1∂x(h2u2)−l2∂x(h1u1).(11) Following the procedure described in [34], and taking into account the piecewise linear discontinuous profile of non-hydrostatic pressure, we obtain the following horizontal momentum equations ∂x(h1u1)+∂xh1u2 1+h1p1+u1+u2 2ΓI−pb∂xH−pI−∂xzI=−gh1∂xη. (12) ∂x(h2u2)+∂xh2u2 2+h2p2−u1+u2 2ΓI+pI+∂xzI=−gh2∂xη. (13) Note that the gradients of p1and p2appear because, inside each layer, we assume that the pressure is linear. The following vertical momentum equations are also deduced, ∂x(h1w1)+∂x(h1u1w1)+w1+w2 2ΓI+pI−−pb=0.(14) ∂x(h2w2)+∂x(h2u2w2)−w1+w2 2ΓI−pI+=0.(15) We shall consider in this work the unknowns pband pI:= pI−, for the non-hydrostatic pressure unknowns. To close the system, we will suppose that pI+=γ1pb+γ2pI,γ 1+γ2= 0,(16) and therefore, due to the assumption on the non-hydrostatic vertical profile: p1=pb+pI 2,p2=γ1pb+γ2pI 2. The underlying reason to assume this discontinuity for the non-hydrostatic pressure is to introduce an artificial set of free parameters, which will help to improve the dispersive properties of the system (see Sect. 3). Remark 1 Notice that for γ1=0,γ 2=1 we recover the original system derived in [34], and for γ1=1,γ 2=δ δ−2the approach pI+−pI−=δ(p2−p1). We will avoid the case γ1+γ2=0 (see Sect. 3.1 and “Appendix D.3”). To sum up, collecting all the equations described before, the system reads as: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂th+∂x(l1hu1)+∂x(l2hu2)=0, ∂t(l1hu1)+∂xl1hu2 1+l1hp1 +u1+u2 2ΓI−pb∂xH−pI∂xzI=−gl1h∂xη, ∂t(l1hu2)+∂xl2hu2 2+l2hp2 −u1+u2 2ΓI+(γ1pb+γ2pI)∂xzI=−gl2h∂xη, ∂t(l1hw1)+∂x(l1hu1w1)+w1+w2 2ΓI=pb−pI, ∂t(l2hw2)+∂x(l2hu2w2)−w1+w2 2ΓI=γ1pb+γ2pI, ∂xu1+2w1−wb l1h=0, ∂xu2+2w2−wI+ l2h=0, (17a) (17b) (17c) (17d) (17e) (17f) (17g) where ⎧ ⎪ ⎨ ⎪ ⎩ wI+=u2∂xzI−∂x(l1hu1)−∂tH, wb=−u1∂xH−∂tH. (17h) (17i) The system depends on the parameters (l1,γ 1,γ 2), which need to be chosen. This will be done following a criterion that improves the dispersive relations of the system. Note that the proposed system has the same number of variables as when considering a continuous pressure, that is pI+−pI−=0, that are h,uα,w α,pI,pb and therefore the computational cost of the improved model will be the same than the original one proposed in [34]. Finally, system (17) can be written in the compact form ⎧ ⎨ ⎩ ∂tU+∂xF(U)+B(U)∂xU=G(U)∂xH+T(h,∂ xh,H,∂ xH,P,∂ xP), B(U,∂ xU,H,∂ xH)=0, (18) where we introduce the notation U= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ h qu,1 qu,2 qw,1 qw,2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,F(U)= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ l1qu,1+l2qu,2 q2 u,1 h+1 2gh2 q2 u,2 h+1 2gh2 qu,1qw,1 h qu,2qw,2 h ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,G(U)= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 gh gh 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , being qu,α =huα,qw,α =hwα,α=1,2 the horizontal and vertical discharges. B(U)is a matrix function such that B(U)∂xUinvolves the non-conservative products related to the mass transfer across interfaces that appear in the momentum equations B(U)= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 00 000 0−l2 u1+u2 2l2 u1+u2 200 0l1 u1+u2 2−l1 u1+u2 200 0−l2 w1+w2 2l2 w1+w2 200 0l1 w1+w2 2−l1 w1+w2 200 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . Note that the corresponding equations for the horizontal and vertical momentum at the layer Lα, have been divided by lα. The limit case l1=1,l2=0 corresponds to the one-layer model presented in [16]. For the sake of simplicity, we will assume that lα∈(0,1). Finally, the non-hydrostatic terms are given by P=⎛ ⎝ pb pI⎞ ⎠,T(h,∂ xh,H,∂ xH,P,∂ xP)=− ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 ∂x(hp1)−pb∂xH/l1−pI∂xzI/l1 ∂x(hp2)+(γ1pb+γ2pI)∂xzI/l2 (pI−pb)/l1 −(γ1pb+γ2pI)/l2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , where p1=pb+pI 2,p2=γ1pb+γ2pI 2, and B(U,∂ xU,H,∂ xH)=⎛ ⎝ l1h∂xqu,1−2qu,1∂xz1+2qw,1+2h∂tH 2l1h∂xqu,1+l2h∂xqu,2−2qu,2∂xz2+2qw,2+2h∂tH⎞ ⎠.(19) Operator (19) contains the incompressibility condition for each layer, and is obtained multiplying Eqs. (17f)and(17g)byl1h2and l2h2respectively. Equations (17h)and(17i)have been used as well. This allows to write the full two-layer system in terms of discharges. Remark 2 Let us consider the underlying hydrostatic system given by ∂tU+∂xF(U)+B(U)∂xU=G(U)∂xH.(20) Then the system (20) is hyperbolic for every l1∈(0,1). A proof is given in “Appendix A.1”. In the case of l1=1 2,the eigenvalues can be computed explicitly and are given by: λ1=u1,λ 2=u2,λ 3=u1+u2 2,λ 4,5=u1+u2 2±gh +3 4(u1−u2)2. An approximation of the eigenvalues for the general case are given in “Appendix A.2”. Remark 3 Concerning the well-possessedness of the proposed system, if we restrict just to the hyperbolic problem (20), without non-hydrostatic pressure, some studies can be found although there is no final answer. For classical two-layer shallow water systems some results can be found in [54–57]. In the case of multi-layer shallow water systems, we may cite the work [7], where the authors prove that the multilayer model is hyperbolic for the case of two layers, but the question remains open for the general case. Another work related to this topic is [36]. In this paper, it is shown that the solution of the multilayer model is a particular weak solution of the Navier-Stokes system. Concerning multilayer problems with non hydrostatic pressure, to our knowledge, not much is known concerning well-possessedness. 3 Linear Analysis of the New Model In this section, dispersive relations of the integrated two-layer system (17) are studied as usual (see [25,45,46,52,62]). In particular, we focus on the linear dispersion relation, the group velocity, and the linear shoaling. This constitutes a standard study of PDE systems for dispersive water waves modelling. The two first properties are related to the propagation of dispersive wave trains, and the latter with shoaling processes. This can occur when waves arrive at the continental shelf, from intermediate to shallow waters. In Sect. 3.1,aflatbottomHis assumed and the governing Eq. (17) are linearised around a steady state solution. Then, we conduct a standard Stokes-type Fourier analysis to obtain the linear dispersion relation and group velocity. In Sect. 3.2, a shoaling analysis of the linearised equations is carried out. The linear dispersion relation, group velocity and linear shoaling are derived following the procedure described in [62]. Remark that the studied dispersive relations are dependent on the parameters (l1,γ 1,γ 2). In Sect. 3.3 the aforementioned (l1,γ 1,γ 2) parameters will be chosen through examination of the derived dispersive relations following the standard procedure on such topics (see [46,52,62]). where  Ui+1/2= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝  hi+1/2  hi+1/2u1,i+1/2  hi+1/2u2,i+1/2  hi+1/2w1,i+1/2  hi+1/2w2,i+1/2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,Gi+1/2= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 g hi+1/2 g hi+1/2 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , being  hi+1/2=h+ i+1/2+h− i+1/2 2,uα,i+1/2= u+ α,i+1/2h+ i+1/2+u− α,i+1/2h− i+1/2 h+ i+1/2+h− i+1/2 for α=1,2. wα,i+1/2is defined in the same way as uα,i+1/2. The vector U± i+1/2is defined by a reconstruction procedure on the variables to the to the left (−) and right (+)oftheinter-cellxi+1/2: h± i+1/2,η ± i+1/2,u± α,i+1/2,w ± α,i+1/2. This reconstruction procedure is done using a MUSCL reconstruction operator (see [66]), combined with a minmod limiter. The MUSCL reconstruction operator takes into account the positivity of the water height. Finally, the variable H± i+1/2is recovered from H± i+1/2=h± i+1/2−η± i+1/2. This procedure allows the scheme to deal with emerging topographies, since the variable ηis reconstructed instead oh Hin those situations, which means that the gradient of ηis set to zero to avoid spurious non-physical pressure forces (see [18]). Qi+1/2is the viscosity matrix associated to the numerical method. For PVM schemes, Qi+1/2is obtained by a polynomial evaluation of the Roe Matrix. In this work, Qi+1/2is defined as Qi+1/2=β0Id +β1Ai+1/2, being β0=SR|SL|−SL|SR| SR−SL ,β 1=|SR|−|SL| SR−SL , where SL,SRare estimates of the smallest and biggest wave speeds respectively at the interface xi+1/2, as it is usually considered for such methods. Here, we use SL=min ˆui+1/2−gˆ hi+1/2,¯ui−ghi, SR=max ˆui+1/2+gˆ hi+1/2,¯ui+1+ghi+1, where ˆui+1/2=¯ui√hi+¯ui+1√hi+1 √hi+√hi+1 ,¯ui=l1u1,i+l2u2,i. Ai+1/2denotes the matrix Ai+1/2=JF(  Ui+1/2)+B(  Ui+1/2), being JFthe Jacobian matrix of the flux F. The scheme is a segment path-conservative extension of HLL scheme [41] for nonconservative systems. Notice that the term Qi+1/2Ai+1/2−1Gi+1/2H+ i+1/2−H− i+1/2 in (35) can be interpreted as the up-winding part of the source term discretization. This term makes no sense if any of the eigenvalues of Ai+1/2vanishes. One way to deal with this kind of resonant problems has been proposed in [20]. Remark 5 It is clear that the concept of weak solution and the definition of path-conservative schemes strongly depends on the chosen family of paths, which is a priori arbitrary. The crucial question is thus how to choose the “good” family of paths. The answer to this question is not easy and it is not the scope of this article. We refer to [21] for a detailed discussion on the subject. In this paper, we consider the canonical path that is the segment path. In the reference [13] a negative result of convergence for path-conservative numerical methods was given together with several numerical examples. Later, a new example was given in [2]. These results have led in many cases to an over-simplified picture according to which these convergence difficulties are a specific drawback of path-conservative methods. Moreover, the experience proves that path-conservative numerical methods and in particular PVM schemes converge with the expected order of accuracy and, under the adequate CFL condition, with the same stability property as their conservative counterparts for geophysical problems. The difficulty comes from the fact that the limits of numerical solutions may differ from the correct ones in the case of non-conservative systems as weak solutions of non-conservative systems may be defined in infinitely many ways (paths). In fact, the key point to understand the problem of convergence in non-conservative systems is the following: the limits of the numerical solutions satisfy a jump condition which is related to the numerical viscosity of the method and not to the physically relevant one. Of course, this phenomenon affects any numerical method in which the small-scale effects (the vanishing diffusion and/or dispersion) are not taken into account, regardless of whether it is path-conservative or not. For instance, even the Godunov method based on the right weak solutions of the Riemann problems fails, in general, to converge to the right solution: this is due to the numerical viscosity introduced in the averaging step. Remark 6 The numerical scheme considered for the underlying hyperbolic system (32a and 34) is well-balanced for water at rest solutions and linearly L∞−stable under the usual CFL condition. Moreover, the scheme is positive preserving for the water height for a smooth bathymetry under 1 2CFL condition. 4.3 Finite Difference Discretization for the Non-hydrostatic Terms In this subsection we describe the discretization of the non-hydrostatic terms. From Eq. (32b), for every k∈{1,2}of the Runge–Kutta stage, we get ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ q(k) u,1=q( k) u,1−Δt∂x(h(k)p(k) 1)−1 l1p(k) b∂xH+p(k) I∂xzI(k), q(k) u,2=q( k) u,2−Δt∂x(h(k)p(k) 2)+1 l2γ1p(k) b+γ2p(k) I∂xzI(k), q(k) w,1=q( k) w,1−Δt l1p(k) b−p(k) I, q(k) w,2=q( k) w,2−Δt l2γ1p(k) I+γ2p(k) b. (36a) (36b) (36c) (36d) Replacing Eqs. (36a–36d)into(32c), we obtain the Poisson-like equations a1∂xx p(k) b+a2∂xp(k) b+a3p(k) b+a4∂xx p(k) I+a5∂xp(k) I+a6p(k) I=RHS 1,(37) b1∂xx p(k) b+b2∂xp(k) b+b3p(k) b+b4∂xx p(k) I+b5∂xp(k) I+b6p(k) I=RHS 2,(38) where the coefficients aj,bjas well as the Right-Hand-Sides RH Sj,j∈{1,2},are given in “Appendix D.1”. We point out the dependency of the aforementioned coefficients aj≡aj(U(k)), bj≡bj(U(k)), RHSj≡RHSj(U( k)) and remark that, due to the form of the numerical scheme and the operator T,then one has that h(k)=h( k),and thus aj≡aj(U( k)), bj≡bj(U( k)), RHSj≡RHSj(U( k)). Equations (37)and(38) are discretized using second order finite differences. In order to obtain point value approximations of the non-hydrostatic pressure variables, terms p(k) I,i,p(k) b,i and RHSj,i,j∈{1,2},i∈{1,...,N}will be approximated at every mid point xiof each cell Ii. In the same way, the corresponding space derivatives will be approximated using compact second order finite differences ∂xP(k) i= P(k) i+1−P(k) i−1 2Δx,∂ xx P(k) i= P(k) i+1−2P(k) i+P(k) i−1 Δx2.(39) To compute the coefficients aj,bj,RHSj,that contains terms depending on U( k)as well as first derivatives ∂xU( k),we will use the averaged values on the cell Iiof variables U( k) and Has second order point value approximations of U( k) iand Hiat the centre of the cell Ii. The space derivatives are computed from the averaged values using compact second order finite differences ∂xU( k) i= U( k) i+1−U( k) i−1 2Δx,∂ xHi=Hi+1−Hi−1 2Δx.(40) After replacing (39)and(40)in(37)and(38), one has to solve a linear system AP=RHS,(41) where RHS =⎛ ⎝ RHS1 RHS2⎞ ⎠,RHS j=⎛ ⎜ ⎝ RHSj,1 . . . RHSj,N ⎞ ⎟ ⎠,P= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ p(k) b,1 . . . p(k) b,N p(k) I,1 . . . p(k) I,N ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , and A=D+M,D=⎛ ⎝T(1)0 0T(2)⎞ ⎠,M=⎛ ⎝0C(1) C(2)0⎞ ⎠,(42) being T(j),C(j)tridiagonal and symmetrical matrices of dimension N×N(see details in “Appendix D.2”). To solve the linear system (41), we propose an iterative linear solver based on a block version of the Jacobi method: ⎧ ⎨ ⎩ P(0)=0, DP(s)=RHS −MP(s−1),s∈1,2,..., (43) where at each iteration s,a tridiagonal linear system is efficiently solved by using the Thomas algorithm. As stopping criteria for the iterative linear solver (43), we choose ||P(s+1)−P(s)|| < tol. In the subsequent numerical test performed in this work, tol is set to 10−8. Once the non-hydrostatic pressure terms P(k)have been computed, the velocities U(k) can be updated from Eqs. (36a–36d), wherein a similar way, a second order point value approximation in the centre of the cell will be used. Remark 7 Although T(j)and C(j)are symmetrical matrices, in general the same does not hold true for the matrix A. A particular analysis of the singularity of the matrix Ahas been made for the case of small water height in “Appendix D.3”, but the same question remains open for the general case. Nevertheless, the numerical experiments performed show that the system is indeed invertible. Furthermore, the fast convergence of the iterative method used to solve the linear system let us make the hypothesis that it has a good condition number. Remark 8 The proposed iterative linear solver has some advantages over other standard linear solvers such as Krylov gradient-based methods (e.g., biconjugate gradient stabilized method) or LU factorization. A short discussion is summarized: – The proposed linear solver reaches the convergence in less than 20 iterations for the numerical tests showed in this work. – It can be easily implemented, and it is matrix-free, since only tridiagonal systems have to be solved, which is done with the matrix-free Thomas algorithm. This makes this algorithm extremely convenient for an implementation on GPU architectures, where memory is an important issue. –ALU factorization method it is not feasible since the matrix of the linear system A depends on time step. This will demand a considerable computational effort. Moreover, this choice is more memory consuming than the iterative method proposed in this work in (43). – Let us denote by  Pa convenient reorganization of the variable P,such that  P=p(k) b,1,p(k) I,1,...p(k) b,N,p(k) I,Nt . Then the corresponding associated matrix  Ais symmetric and penta-diagonal. This organization is more convenient for an iterative Krylov space based method. Nevertheless, we have found that this methods are more memory and time consuming, even with preconditioners, when compared with the iterative solver proposed in this work. – Considering future works and the extension of the scheme to bidimensional domains, this will lead to a Poisson-like equation in two dimensions. The proposed algorithm can be combined with a Jacobi or Jacobi relaxation method (see [3]) to take into account this extension. 4.4 Final Numerical Scheme For the sake of clarity, a guideline of the final numerical scheme will be highlighted in what follows. For every kth stage k∈{1,2}of the Runge–Kutta method, the problem (18) is split into two parts. A two-step projection-correction method is used: –Finite volume step (solving the hydrostatic system): From (34), solve explicitly the hyperbolic system at the kth stage of the Runge–Kutta: U( k)=U(k−1)−Δt ΔxD(k−1),− i+1/2+D(k−1),+ i−1/2+Ii, by means of a PVM path-conservative finite volume scheme (34) combining a MUSCL reconstruction operator to obtain the intermediate value U( k)in the two-step projectioncorrection method that contains the numerical solution of the underlying hyperbolic system. –Finite Difference step (non-hydrostatic pressure correction): –Solve the Poisson-like Eqs. (37)and(38) to obtain the non-hydrostatic pressures P(k)=⎛ ⎝ p(k) b p(k) I ⎞ ⎠. To do so, compact centred second order finite differences are used for the discretization of the derivatives that appear in (37)–(38), and a linear system is solved to obtain P(k). –With the computed non-hydrostatic pressure terms P(k), the velocities can be updated from (36a)–(36d), where in a similar way, a second order point value approximation in the centre of the cell will be used to compute the non-hydrostatic contribution. Finally, Un+1=1 2 Un+1 2 U(2). 4.5 Boundary Conditions In this work, three types of Boundary Conditions (BC) have been considered: periodic, outflow and generating/absorbing BCs. – Periodic BCs: given a spatial domain subdivided into a set of Ncomputing cells, cell I1and IN,which are the extremes of the domain, are considered as the same cell. This means that its neighbour cells are IN−1on the left and I2on the right. – Outflow BCs: homogeneous Neumann conditions are applied on the left and right boundaries. Since we use a second order MUSCL scheme, the usage of one ghost cell I0,IN+1 at each boundary is required in order to determine the values of the closest nodes to the boundary. The values of the variables at the ghost cells are extrapolated from the adjacent cells U0=U1,P0=P1,UN+1=UN,PN+1=PN. Nevertheless, reflections at the boundaries might perturb the numerical solution at the inner domain. As in many other works (see [43,45,53,60] among others), this condition is sometimes supplemented here with an absorbing BC, described below. – Generating/absorbing BCs: Periodic wave generation as well as absorbing BCs are achieved by using a generation/relaxation zone method similar to the one proposed in [53]. Generation/absorption of waves is achieved by simply defining a relaxation coefficient 0 ≤m(x)≤1,and a target solution (U∗,P∗). Given a width LRel of the relaxation zone on each boundary, we define kRel as the first natural number that kRelΔx≥LRel. The solution within the relaxation zone is then redefined to be, ∀i∈{1,...,kRel,N−krel,...N}:  Ui=miUi+(1−mi)U∗ i  Pi=miPi+(1−mi)P∗ i, where miis defined as mi=1−di LRel 2 , where diis the distance between the centre of the cells Iiand I1(or Iiand IN−k), in the case of i∈{1,...,k}(or i∈{N−k,...,N}). In our numerical experiments, we set L≤LRel ≤1.5L, being Lthe typical wavelength of the outgoing wave. Absorbing BC is the special case U∗ i=P∗ i=0,that will dump all the waves passing through. Remark 9 In Sect. 6, some laboratory tests are studied and the effects of the friction of the water with the bottom Hmust be considered. In order to include friction with the bottom, a Manning empirical formula is added to the right hand side of (18) as a source term S(U)=0−τ000 t,τ=gl1hn2u1|u1| (l1h)4/3,(44) and it is discretized in a semi-implicit way at the end of the second step of the proposed numerical scheme, at each Runge–Kutta stage. Remark 10 Since non-hydrostatic and friction terms appear only in the momentum equations, the final numerical scheme is well-balanced for the water at rest solution and positive preserving for the water height. Remark 11 The extension of the proposed numerical scheme to the case of moving bottom is straightforward. Remark 12 It is well known the fact that higher order schemes should be used when discretizing PDEs of dispersive nature, such as Serre–Green–Naghdi [40], Nwogu [58] among others, since the leading numerical error of the numerical scheme is of dispersive nature (see [25,45,48,60] among many other works where Boussinesq-type systems are discretized). In this case, a linear stability Von-Neumann analysis shows that second order schemes, as well as first order splitting schemes, may destroy the continuum dispersive relation of the original system. Although we did not perform such an analysis in this work, we obtain positive results in the numerical tests. Concretely, Sect. 6.2.2 includes a convergence analysis. Results in Table 1 show that the numerical scheme is formally second-order in space and time. We have also measured the impact on the phase velocity for the propagation of solitary waves for long integration times with the second order scheme. In Sect. 6.2.3 we have measured the variations on amplitude and celerity of a solitary wave of big amplitude. Results in Table 2 show excellent agreement, obtaining a variation of less than 0.1128%. 5 Breaking Wave Modelling and Wet–Dry Treatment 5.1 Breaking Wave Modelling As pointed in [61], in shallow water complex events can be observed related to turbulent processes. One of these processes corresponds to the breaking of waves near the coast. As it will be seen in the numerical tests proposed in this work, the model presented here cannot describe this process without an additional term which allows the model to dissipate the required amount of energy on such situations. When breaking processes occur, mostly close to shallow areas, two different approaches are usually employed when dispersive Boussinesqtype models are considered. Close to the coast where breaking starts, the SWE propagates breaking bores at the correct speed, since kH is small, and dissipation of the breaking wave is also well reproduced. Due to that, the simplest way to deal with breaking waves, when considering dispersive systems, consists in neglecting the dispersive part of the equation. This means to force the nonhydrostatic pressure to be zero where breaking occurs. Due to that, this technique has the advantage that only a breaking criterion is needed to stop and start it. However, the main disadvantage is that the grid-convergence is not ensured when the mesh is refined, and global and costly breaking criteria should be taken into account (see [43]). The other strategy, that will be adopted in this work, consists in the dissipation of breaking bores with a diffusive term. Again, breaking criteria to switch on/off the dissipation is needed. Usually, an eddy viscosity approach (see [61]) solves the matter, where an empirical parameter is defined, based on a quasi-heuristic strategy to determine when the breaking occurs. The main difficulty that presents this mechanism is that usually the diffusive term must be discretized implicitly due to the high order derivatives from the diffusion. Otherwise, it will lead to a severe restriction on the CFL number. As a consequence, an extra linear system has to be solved, losing efficiency. In any case, this challenge is overcome in [33] for the one layer non-hydrostatic system derived in [71], and a natural extension of this procedure for the two-layer case will be described in what follows. Let us consider the vertical component of the stress-tensor τzz =2ν∂zw, where ν(x,z,t)is a positive function. We use the same process as in Sect. 2to depth-average the vertical component of the stress-tensor. By taking into account the incompressibility condition and that the vertical velocity has a linear profile within each layer we obtain Lα ∂zτzz dz =−ςαhα∂xuα, where ςα=Lα∂zν. Let us denote T(h,∂ xh,H,∂ xH,P,∂ xP)as T(P,∂ xP)for the sake of brevity. Adding the proposed integrated viscosity terms, system (18) reads ⎧ ⎨ ⎩ ∂tU+∂xF(U)+B(U)∂xU=G(U)∂xH+T(P,∂ xP)−τ+Rb(U,∂ xU,ς 1,ς 2), B(U,∂ xU,H,∂ xH)=0, (45) where Rb(U,∂ xU,ς 1,ς 2)=000−ς1h∂xu1−ς2h∂xu2t. Now we close the system defining ςα(x,t). The linear combination of the non-hydrostatic pressures, pband pI, at the fourth and fifth component of T(P,∂ xP), can be expressed in terms of uαand its derivatives ∂xxuα,∂ xtuα,∂ xuαand ∂tuα. The proposed election of ςαin this work is based on the idea of cancelling those aforementioned ∂xuαterms with −ςαhα∂xuα. The procedure is detailed in “Appendix C”. The following definition, for the case of immobile bed ∂tH=0, is proposed: ⎧ ⎪ ⎨ ⎪ ⎩ ς1=w1−3u1∂xH, ς2=w2+3u2∂xzI. (46) Note that fourth and fifth components of Rbare essentially first order derivatives of u1,u2, and can be discretized explicitly. Finally, a breaking criteria to switch on/off the dissipation is needed. Following a natural and simpler extension of the criteria proposed by [61], wave energy dissipation associated with breaking begins when |∂x(l1hu1+l2hu2)|≥U1and continues as long as |∂x(l1hu1+ l2hu2)|≥U2,where U1=B1gh,U2=B2gh, denote the flow speeds at the onset and termination of the wave-breaking process and B1,B2 are calibration coefficients that should be calibrated through laboratory experiments. In this work, as in [61], we use B1=0.5andB2=0.15 for al the test cases studied. Note that the breaking criteria to switch on/off the dissipation is a simplified version of the one proposed in [61], that includes some improvements such as: taking into account a residence time for the activation/deactivation of the criteria to accounting a continuous dissipation; or computing in the breaking criteria |∂x(hu)|as 1 2(|∂x(hu)|+∂x(hu)) which would automatically become zero on the back of the crest of the wave. The breaking mechanism proposed in this work can be considered with this improved breaking criteria given in [61], as well as the one proposed in [43], which are more sophisticated and expensive. Nevertheless, although we have chosen fast and straightforward breaking criterion, the numerical tests in Sect. 6will show that this technique performs adequately. Moreover, the simple breaking mechanism considered in this work: corrects the classical overshoot that dispersive models present for the run-up of waves (see Fig. 13); ensures the grid convergence even if breaking mechanism is switching on/off during the experiment. Although the simple breaking criteria do not detect hydraulic jumps, we have observed that the proposed system without the breaking dissipation can handle well with hydraulic jumps. In any case, a more sophisticated breaking detector that reveals hydraulic jumps can be considered as well. Remark 13 In [33], the one layer non-hydrostatic system, which corresponds to the election of l2=0 for the presented system in this work, is studied. There, the classical eddy viscosity model approach for breaking waves is considered. The eddy viscosity model, which consists in a diffusion term in the horizontal momentum equation, is reinterpreted as a term of the form ςh∂xuat the vertical momentum equation. The resulting terms in this work in Rbcan be seen as the natural extension of the idea presented in [33], and thus as a natural extension of the classic eddy viscosity models for small-scale turbulence. 5.2 Wet–Dry Treatment For the computation of U( k)in the finite volume discretization of the underlying hyperbolic system, a wet–dry treatment adapting the ideas described in [18] is applied. The key to the numerical treatment for wet–dry fronts with emerging bottom topographies relies on: – The hydrostatic pressure terms gl1h∂xη, gl2h∂xηat the horizontal velocity equations are modified for emerging bottoms to avoid spurious pressure forces (see [18]). – To compute velocities appearing in Sect. 4.2 from the discharges, very small values for the water height hmay appear. To overcome the difficulties due to large round-off errors in computing the velocities uα,w αappearing in the scheme, we define the velocities analogously as in [44] applying the desingularization formula uα=√2hhuα h4+max(h4,4) ,w α=√2hhwα h4+max(h4,4) ,α∈{1,2}, which gives the exact value of uαand wαfor h≥, and gives a smooth transition of uα and wαto zero when htends to zero without truncation. In this work we set =10−5 for the numerical tests. A more detailed discussion about the desingularization formula can be seen in [44]. In the second step of the projection method, no special treatment is required due to the rewriting of the two incompressibility equations, that have been multiplied by h2, in terms of discharges. In presence of wet–dry fonts, non-hydrostatic pressure P vanishes and no artificial truncation up to a threshold value is needed. In “Appendix D.3”, a particular analysis is detailed for the case of h = , uα = wα = 0,∂x H = m. We also assume that 2≈0,·m≈0, to simplify the analysis. A more general case is considered in “Appendix D.3”. In this case, the Right hand side of the linear system becomes zero, and the matrix Aas well as the right hand side RHS defined in (41) simplifies to A=4(m2+1)⎛ ⎝I−I γ2Iγ1I⎞ ⎠,RHS =⎛ ⎜ ⎜ ⎜ ⎝ RHS1 RHS2 ⎞ ⎟ ⎟ ⎟ ⎠=⎛ ⎜ ⎜ ⎜ ⎝ 0 0 ⎞ ⎟ ⎟ ⎟ ⎠. Note that Ais invertible, since γ1+γ2= 0,and thus the non-hydrostatic pressure Pbecomes zero. Remark 14 The definition of the coefficients of the matrix A,are given in terms of l1h,l2h and η, and the topography Hdoes not appear. Thus, in the case of emerging bottoms, the derivatives of ∂xηare treated as the pressure terms ghα∂xηin the horizontal velocity equations, and the numerical scheme can deal with a wide variety of complex situations involving wet–dry fronts and emerging topographies, as we will see in the numerical tests. 6 Numerical Tests In this section we present five numerical tests for the model presented in (45) with the breaking and friction terms introduced previously. The chosen tests have been widely studied in dispersive water waves modelling (see [23,31–33,43,46,49,50,60]). Thus, this constitutes a proper way to validate the presented dispersive model and the numerical scheme with analytical solutions and real laboratory tests. The first three tests aim to validate the mathematical model and the numerical scheme with a comparison to analytical or reference solutions. The other two tests show that the solved equations can simulate complex laboratory experiments, showing comparisons with real data. Comparisons with the one-layer system results, which corresponds to the choice l2=0, are presented for some of the considered tests. These comparisons show that the one-layer system does not give good results in some cases. This means that at least a two-layer model should be used. The quantities of the parameters concerning the following numerical simulations are expressed in units of measure of the International System of Units. 6.1 Standing Wave in Closed Basin First, we shall study the dispersive properties of the model. To do so, we consider the analytical solution from the Stokes linear theory of a standing wave in a closed basin with length Land flat bottom H=10. The initial surface elevation is given by η(x,0)=Acos 2πx L, Fig. 10 Head-on collision of two solitary waves at times T√g(A+H)=0,25,50,75,100 Fig. 11 Head-on collision of two solitary waves. Final state after the interaction. Comparison of the computed free-surface with the exact soliton profiles the same speed but in opposite directions. Figure 11 shows the final state after the collision compared with the exact profile, showing an excellent agreement. 6.4 Solitary Wave Run-Up on a Plane Beach Synolakis [65] carried out laboratory experiments for incident solitary waves, to study propagation, breaking and run-up over a planar beach with a slope 1:19.85. Many researchers have used this data to validate numerical models. With this test case, we assess the ability of the model to describe shoreline motions and wave breaking when it occurs. Experimental data are available in [65] for surface elevation at different times. The bathymetry of the problem is described in Fig. 12. A solitary wave of amplitude A = 0.3 is placed at point x = 20 as initial condition for the free-surface elevation following the same procedure as in Sect. (6.2). A Manning coefficient of nm = 0.01 is used in order to define the glass surface roughness used in the experiments. The computational domain is [− 10, 40] and the numerical parameters used are Δx = 0.02, CFL = 0.9and g = 9.81. Free-outflow boundary conditions are considered. Simulations are carried out with the improved parameters (27). Fig. 12 Sketch of the bathymetry Fig. 13 Comparison of experiments data (red) and simulated ones with breaking terms proposed in (46) (blue) at different times during the run-up. Between bars, regions where breaking mechanism is active are shown (Color figure online) Figures 13 and 14 show snapshots, at different times, t√g/H=t0where H=1, comparing experimental and simulated data. Here we use the proposed ςαfor the breaking mechanism described in Sect. 5. Breaking parameters are set to B1=0.5andB2=0.15. Fig. 13 also shows where the breaking mechanism is active (region between the bars), and demonstrates the efficacy of the criteria. As we discussed in Sect. 5, the breaking criteria does not take into account the presence of hydraulic jumps. However, Fig. 14 shows that the system without the extra breaking dissipation proposed in this work can handle this situation correctly. The breaking mechanism also works properly in terms of grid convergence. Figure 15 shows the snapshots at times t√g/H=15 (run-up) and t√g/H=55 (run-down) for different mesh sizes. Besides, good results are obtained at maximum run-up, where friction terms play an important role. Note that no additional wet–dry treatment on the second step of the scheme is necessary. Fig. 14 Comparison of experiments data (red) and simulated ones with breaking terms proposed in (46) (blue) at different times during the run-down (Color figure online) Fig. 15 Comparison of free-surface simulation at times t√g/H=15 and t√g/H=55 for different mesh sizes Since the propagated wave in this test is weakly-nonlinear weakly-dispersive, the results obtained with the two-layer system are quite similar to the results obtained with the one-layer system (see [33]). Nevertheless, this test shows that the breaking mechanism introduced in this work, as well as friction terms, perform adequately for the two-layer system. Moreover, the corresponding discretization is robust and can deal with the presence of wet–dry fronts correctly. Fig. 16 Periodic waves breaking over a submerged bar. Sketch of the bathymetry 6.5 Regular Wave Propagation Over a Submerged Bar The experiment of plunging breaking periodic waves over a submerged bar by Dingemans [28] is considered here. The experiment studies the frequency dispersion characteristics and nonlinear interaction. As waves propagate over a submerged bar, multiple phenomena occur, such as non-linear shoaling, amplification of higher harmonics and wave breaking. The one-dimensional domain [0,30]is discretized with Δx=0.01 and the bathymetry is defined in the Fig. 16.TheCFLissetto0.5andg=9.81. Horizontal and vertical velocities and ηare set initially to 0. Boundary conditions correspond to free-outflow at x=30 and a sinusoidal wave train for ηgenerated at x=0. This is done using the target initial condition at the relaxation zone: η∗(t)=Asin 2π Tt,u∗ α(t)=gH 0η∗(t), w∗ α(t)=p∗ b(t)=p∗ I(t)=0,α=1,2 where Aand Tdenotes amplitude and period and H0=0.4. We focus on test cases (a),(b) and (c)described in the table below Case A (m) T (s) (a) 0.01 2.02 Non-breaking waves (b) 0.0145 2.525 Breaking waves (c) 0.0205 1.01 Non-breaking waves Finally, the relaxation length for the initial generating/absorbing boundary conditions are set to LRel =2 for cases (a)-(c). Depth gauges, which measure the free surface elevation, are placed for cases along x= gi,g1=4,g2=5.7,g3=10.5,g4=12.5,g5=13.5,g6=14.5,g7=15.7,g8= 17.3,g9=19,g10 =21. Test (a)and (c)produce, up to the front slope, waves with kH ≈0.63 and kH ≈1.58 respectively. Figure 17 shows the time series of free surface for case (a). The numerical model maintains a good agreement with the experimental data at each station. Minor discrepancies can be observed behind the bar at stations g5−g10,where higher harmonics are released. The test case (a)coincides with one of the proposed in [12], and similar good agreement with other numerical models can be found in the literature. Let us remark that in [50]σ-coordinate Fig. 17 Comparison of experiment data (a)(red points) and simulated ones with improved dispersive parameters (l(5) 1,γ(5) 1,γ(5) 2)(green) and (l(2L) 1,γ(2L) 1,γ(2L) 2)(blue) (Color figure online) discretization of the free-surface Euler equations is used. While in that paper the authors need three layers for the test, we can achieve good results with the technique proposed here that only considers two layers. The results in [25] with a three-parameter Green-Naghdi model optimized for uneven bottoms, are also in good agreement. Comparing to that, let us remark that the Green-Naghdi system described there includes high order terms using third order derivatives. We would like to stress the ability of the proposed two-layer system to deal with a wide range of dispersive waves, which is achieved without including high order terms in the system of equations. Only first order derivatives are used here. In [9] the test is also studied with a depth-integrated free-surface two-layer non-hydrostatic system. The results from case (c)are given in Fig. 18. In this situation, larger values of kH give short waves and make this case more challenging. Due to this fact, waves propagate and typical cnoidal profile can be observed over the top of the bar. This is well reproduced at gauges g1,g2,g3,g4. Behind the bar, at g5,g6,g7, higher harmonics appear. As we see, a good fit for the amplitude is observed. This test case was studied in [46], where excellent results are obtained at g5,g6,g8,g10,with a two-layer system that it has been also optimized up to the second-order nonlinear dispersive analysis. Figures 17 and 18 also show comparison with (l(5) 1,γ(5) 1,γ(5) 2)and (l(2L) 1,γ(2L) 1,γ(2L) 2) parameters. One can see a better matching in favour of the former. Fig. 18 Comparison of experiment data (c)(red points) and simulated ones with improved dispersive parameters (l(5) 1,γ(5) 1,γ(5) 2)(green) and (l(2L) 1,γ(2L) 1,γ(2L) 2)(blue) (Color figure online) The next case presented is (b), where waves start to break on the top of the bar. Numerical time series of surface elevation are shown in Fig. 19 for (l(5) 1,γ(5) 1,γ(5) 2).Hereweusethe proposed ςαfor the breaking mechanism described in Sect. 5.Breaking parameters are set to B1=0.5andB2=0.15. To evince that a breaking mechanism is needed, the corresponding numerical simulation without the breaking mechanism is also shown in Fig. 19.Thetest shows the need to consider a breaking mechanism. Before reaching the bar, both simulations coincide as expected, since the wave breaking starts around x=12, revealing that the start/stopping criteria performs properly. In all gauges that follow, the wave shape, frequency and amplitude are well reproduced in favour of the simulations with the breaking mechanism. In [29], the test case (b)is performed with a Serre Green–Naghdi system and a different and well-validated breaking mechanism. The results obtained in this work, are in accordance with the ones presented in [29]. The comparison with experimental data emphasizes the need to consider a dispersive model to faithfully capture the shape of the waves near the continental shelf, where a wide range of dispersive waves can be released. Thus, the two-layer system presented in this work can simulate an extensive range of frequencies of dispersive waves, which typically arises in nature from intermediate to shallow waters. Fig. 19 Comparison of experiment data (b) (red points) and simulated ones with breaking mechanism (green) and without (blue) (Color figure online) 7 Conclusions In this work, a non-hydrostatic two-layer model has been proposed to incorporate dispersive effects in the propagation of waves in a homogeneous, inviscid and incompressible fluid. The presented model corresponds to the system derived by E.D.Fernández Nieto et al. in [34]for the two-layer case, with a slight difference in the definition of the non-hydrostatic pressure at the interface. This allows us to derive a two-layer non-hydrostatic model with improved dispersive properties. The advantage is that this is done while keeping the computational cost of the model low. The numerical scheme employed combines a finite volume segment path-conservative scheme for the underlying hyperbolic system and finite differences for discretization of nonhydrostatic terms. Furthermore, it is second order accurate and it is well-balanced for the water at rest solution, positive preserving and linearly L∞-stable under the usual CFL condition. The extension to bidimensional domains has been done. More detail about the validation of the model is carrying out. The model is also able to simulate correctly the run-up and run-down process in wet–dry areas. Moreover, no numerical truncation for the non-hydrostatic pressure is needed at wet– dry areas, where non-hydrostatic pressure vanishes. To the best of our knowledge, this is an improvement in non-hydrostatic numerical schemes, where usually non-hydrostatic pressure is truncated to zero up to a threshold value. As it is well known, non-hydrostatic models need some dissipative mechanism for breaking waves to accurately model waves near the coastal areas (see [33]). We have proposed a new, efficient and straightforward breaking mechanism that performs appropriately in the tests shown. Numerical simulations show that the approach presented here correctly solves the propagation of eulerian solitary waves of high amplitude, preserving their shape for large integration times accurately. Comparison with experimental data is also presented. Experimental data justifies the need to incorporate dispersive effects to capture faithfully waves in the vicinity of the continental shelf. Moreover, complex processes such as run-up, shoaling, the appearance of higher harmonics and wet–dry areas, that can not be well reproduced with a one-layer non-hydrostatic system, are simulated successfully for the proposed tests. We would like to stress the good performance when dealing with a wide range of non-linear dispersive waves, as it was said in Sect. 6.5, where usually a Green–Naghdi type model or multilayer nonhydrostatic model is needed to correctly simulate the tests. The system proposed in this work can deal with such dispersive waves with a two-layer system, without the need of increasing the order of the derivatives o more layers which would increase the complexity of the model and its computational cost. Numerical tests also evince that the proposed model introduced in this work with the jump of the non-hydrostatic pressure at the interface (16), improves the results with respect to the two-layer system derived in [34]. This has been shown by comparing with analytical solutions from the Stokes first order theory. This improvement is achieved with the same computational cost than the original two-layer system derived in [34]. The two-layer model and the numerical scheme presented in this work provides thus an efficient and accurate approach to model dispersive effects in the propagation of waves near coastal areas and intermediate waters. Acknowledgements This research has been partially supported by the Spanish Government and FEDER through Research Project MTM2015-70490-C2-1-R and MTM2015-70490-C2-2-R, and Andalusian Government Research Project P11-FQM-8179. Funding was provided by Ministerio de Economía y Competitividad and Consejería de Economía, Innovación, Ciencia y Empleo, Junta de Andalucía. A Hyperbolicity of the Underlying Hydrostatic System A.1 Hyperbolicity Let us consider the underlying hydrostatic system (20) written in quasi-linear form ∂tU+A∂xU=G(U)∂xH,(47) where A=JF+B,beingJFthe Jacobian matrix of the flux F. The characteristic polynomial of the matrix Ais given by P(λ, l1)=1 2(u1−λ)( u2−λ)Q(λ, l1), being Q(λ, l1)a third order polynomial on λgiven by Q(λ, l1)=f(λ) −R(l1), f(λ) =(3u1−u2−2λ)(u2−λ)2−gh, R(l1)=l1K, K=(u1−u2)(u1−u2)2−4gh. For the sake of simplicity on the notation, we do not write explicitly the dependence on U. Let us study the hyperbolicity of the hydrostatic system. It is easily to check that λ1=u1,λ 2=u2 are eigenvalues of the system for every l1∈(0,1). It remains to check if the cubic polynomial Q(λ) has three distinct roots. We will prove that the cubic polynomial it has always three different roots for every l1∈[0,1],and in particular for every l1∈(0,1)as we just request. Let us give a sketch of the proof: 1. Let us remark that Qis a cubic polynomial on λsatisfying Q(−∞,l1)=+∞,Q(∞,l1)=−∞. 2. Note that f(λ) is a cubic polynomial that does not depend on l1. Moreover, it has two local extrema given by the roots of f(λ): λ±=u1+u2 2±1 3+u1−u2 22 . 3. A sufficient and necessary condition for the existence of three real and distinct roots of the cubic polynomial Q(λ, l1)is that: f(λ−)<R(l1)and f(λ+)>R(l1). 4. Note that when l1=0, the polynomial Q(λ, 0)has three roots: λ3=3u1−u2 2,λ 4,5=u2±gh, and thus f(λ−)<R(0)and f(λ+)>R(0). Similarly, when l1=1, the polynomial Q(λ, 1)has three roots: λ3=−u1+3u2 2,λ 4,5=u1±gh, and therefore f(λ−)<R(1)and f(λ+)>R(1). Thus, if we assume that K≥0,then R(1)≥R(l1)≥R(0)and therefore f(λ+)>R(1)≥R(l1) f(λ−)<R(0)≤R(l1). If we assume that K≤0,then R(1)≤R(l1)≤R(0)and therefore f(λ+)>R(0)≥R(l1) f(λ−)<R(1)≤R(l1). This concludes the proof. A.2 A First Order Approximation for the Eigenvalues Let us denote the eigenvalues that depends on l1as λ1(l1)=u1,λ 2(l1)=u2, as the known eigenvalues for any l1∈(0,1),and λ3(l1), λ4(l1), λ5(l1) as the eigenvalues that are roots of the cubic polynomial Q(λ, l1). As a particular case, we have found an explicit form of the eigenvalues for l1=1/2, λ1=u1,λ 2=u2,λ 3=u1+u2 2,λ 4,5=u1+u2 2±gh +3 4(u1−u2)2. Let us consider λ(l1)∈{λ3(l1), λ4(l1), λ5(l1)}, an eigenvalue that depends onl1and is a root of the cubic polynomial Q(λ, l1). We propose the following approximation of the eigenvalues, that gives the exact roots of the cubic polynomial Q(λ, l1)for l1∈{0,1/2,1}. λ3(l1)≈ λ3(l1)=3 2−2l1u1+2l1−1 2u2, λ4,5(l1)≈ λ4,5(l1)=l1u1+l2u2±gh +3l1l2(u1−u2)2. Another approximation for the eigenvalues is proposed in the following. Since λ(l1)is a root of Q(λ, l1), then Q(λ(l1), l1)=0, and deriving with respect to l1it yields λ(l1)=(u1−u2)2−4gh(u1−u2) 2(gh −3u1u2+3(u1+u2−λ(l1))λ(l1)). Thus, we propose to approximate the eigenvalues that are roots of Q(λ, l1)with the first order approximation λi(l1)≈ λi=λi(1/2)+λ i(1/2)(l1−1/2), i∈{3,4,5}, that can be explicitly computed, since λi(1/2)are known: λ3(l1)≈ λ3=u1+u2 2+(u1−u2)( 1−2l1)gh −1 4(u1−u2)2 gh +3 4(u1−u2)2+O(l2 1), λ4,5(l1)≈ λ4,5=u1+u2 2±gh +3 4(u1−u2)2 +(u1−u2)l1−1 2gh −1 4(u1−u2)2 gh +3 4(u1−u2)2+O(l2 1). This procedure is more rigorous and lead to a more sophisticated expressions of the approximated eigenvalues. 61. Roeber, V., Cheung, K.F., Kobayashi, M.H.: Shock-capturing Boussinesq-type model for nearshore wave processes. Coast. Eng. 57, 407–423 (2010) 62. Schffer, H.A., Madsen, P.A.: Further enhancements of Boussinesq-type equations. Coast. Eng. 26(1), 1–14 (1995) 63. Stansby, P.K., Zhou, J.G.: Shallow-water flow solver with non-hydrostatic pressure: 2D vertical plane problems. Int. J. Numer. Methods Fluids 28(3), 541–563 (1998) 64. Stelling, G., Zijlema, M.: An accurate and efficient finite-difference algorithm for non-hydrostatic freesurface flow with application to wave propagation. Int. J. Numer. Methods Fluids 43(1), 1–23 (2003) 65. Synolakis, C.E.: The runup of solitary waves. Fluid Mech. 185, 523–545 (1987) 66. Van Leer, B.: Towards the ultimate conservative difference scheme. V. A second order sequel to Godunov’s method. Comput. Phys. 32, 101–136 (1979) 67. Wei, G., Kirby, J.T., Grilli, S.T., Subramanya, R.: A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves. J. Fluid Mech. 294(–1), 71 (1995) 68. Whitham, G.B., Wiley, : Linear and nonlinear waves. Earthq. Eng. Struct. Dyn. 4(5), 518–518 (1976) 69. Witting, J.M.: A unified model for the evolution nonlinear water waves. J. Comput. Phys. 56(2), 203–236 (1984) 70. Wu, C.H., Young, C.-C., Chen, Q., Lynett, P.J.: Efficient nonhydrostatic modeling of surface waves from deep to shallow water. J. Waterw. Port Coast. Ocean Eng. 136(2), 104–118 (2010) 71. Yamazaki, Y., Kowalik, Z., Cheung, K.F.: Depth-integrated, non-hydrostatic model for wave breaking and run-up. Numer. Methods Fluids 61, 473–497 (2008) 72. Zijlema, M., Stelling, G.S.: Further experiences with computing non-hydrostatic free-surface flows involving water waves. Int. J. Numer. Methods Fluids 48(2), 169–197 (2005)