An Arbitrary High Order Well-Balanced ADER-DG Numerical Scheme for the Multilayer Shallow-Water Model with Variable Density

Abstract

In this work, we present a novel numerical discretization of a variable pressure multilayer shallow water model. The model can be written as a hyperbolic PDE system and allows the simulation of density driven gravity currents in a shallow water framework. The proposed discretization consists in an unlimited arbitrary high order accurate (ADER) Discontinuous Galerkin (DG) method, which is then limited with the MOOD paradigm using an a posteriori subcell finite volume limiter. The resulting numerical scheme is arbitrary high order accurate in space and time for smooth solutions and does not destroy the natural subcell resolution inherent in the DG methods in the presence of strong gradients or discontinuities. A numerical strategy to preserve non-trivial stationary solutions is also discussed. The final method is very accurate in smooth regions even using coarse or very coarse meshes, as shown in the numerical simulations presented here. Finally, a comparison with a laboratory test, where empirical data are available, is also performed.

1 Introduction

When dealing with geophysical flows in shallow areas, it is a common practice to use vertically averaged models and, in particular, the so-called shallow water (or Saint-Venant) models (see [35]). In this framework, a vertically averaged horizontal velocity is considered, while the vertical component of the velocity is assumed small and it is neglected. Although this kind of approach has been largely used in many applications, this limits the amount of information that the model is able to provide. Although the lack on information in the vertical direction is a potential drawback of this technique, it is counterbalanced by much less computational effort compared to fully three-dimensional models, which is essential for efficient simulations of geophysical flows. Indeed, it reduces the dimension of the problem by one, making it more affordable to simulate large domains with a good horizontal resolution. However, there are geophysical flows where the vertical component may be important in order to understand the behavior of the full system. In recent years, multilayer shallow water models have been developed to address this issue, and are commonly considered as a good compromise between the traditional one layer shallow water model and the direct numerical simulation of the full 3D system. In the multilayer approach, a vertical division or discretization into a given number of layers is considered. For each layer, similar assumptions to the ones in standard shallow water models are imposed (see [4, 6, 53]). This allows, for instance, to recover a detailed vertical profile of the velocity. Let us remark that the layers are just a means of capturing complex vertical profiles and do not correspond necessarily to physical layers.

In the framework of multilayer shallow water systems, one may find different strategies concerning the interaction between the layers. For instance, [10, 19, 28, 57] assume immiscible layers. For the multilayer systems considered in [3, 49], the layers are not immiscible and a continuous exchange between the layers is assumed. This is done through the incorporation of the mass and momentum transfer terms that could be written as non-conservative terms appearing in the equations. Some recent works on multilayer shallow water systems even consider a variable number of layers, where the number of layers locally adapts depending on the presence of relevant vertical effects (see [9]).

Many recent works related to multilayer shallow water models focus in the incorporation and numerical treatment of non-hydrostatic effects, associated with an estimation of the vertical acceleration of the fluid derived from the Boussinesq-type wave equations (see [11, 46, 50] or [76]). By doing so, the dispersive properties of the model are enhanced and simulation of more physical effects becomes possible. Of course, the computational effort is bigger, but still more affordable than a fully three dimensional model. For recent first order hyperbolic reformulations of nonlinear dispersive shallow water systems, the reader is referred to [7, 14, 45, 48].

In this work, we employ a model described in [5, 54]. The considered model includes the effect of relative changes in density in the pressure terms of the system. There are many physical and practical situations where these density changes are relevant. One classical example can be found at the Strait of Gibraltar, where two different currents, varying in terms of temperature and salinity, meet. Indeed, the water from the Mediterranean sea is warmer but also has a higher salinity concentration than the water from the Atlantic ocean, resulting in a higher density for water coming from the Mediterranean sea. In order to accurately simulate this kind of situations, it is essential to take into account these changes in relative density.

Many multilayer models incorporating density variations may be found in the literature in recent years. In [12, 64], for instance, the differences in density are due to suspended sediments, where sediment species have their own sedimentary velocity. Another approach is given in [1], which considers a two layer model where the water entrainment is approximated by means of an heuristic-dependent source term. Furthermore, in [5] the authors propose a multilayer shallow water model that takes into account density variations due to a given tracer, such as the temperature. The model used in this paper follows a similar approach. Another example is [51], where the authors derive a semi-implicit time discretization approach for a variable density shallow water model with a variable number of layers to improve flexibility and efficiency. For an overview about the history of efficient semi-implicit discretizations of hydrostatic and non-hydrostatic shallow water flows, including flows with multi-valued free surface, the reader is referred to the work of Casulli, see e.g. [22,23,24,25,26].

When the vertical flow features are relevant and become the dominant effect, one may need a big number of layers in order to resolve them. As a consequence, the computational effort increases and we thus partially lose the advantage of using multilayer shallow water models when compared to the full three dimensional system. In this work we propose the use of the Discontinuous Galerkin (DG) method as a way to enhance the subcell resolution of the scheme, allowing coarser meshes while keeping a good accuracy in the vertical direction.

The Discontinuous Galerkin (DG) method itself dates back to the early work by Reed and Hill in [68]. It was later developed by Cockburn and Shu, who set a sound theoretical foundation for these numerical schemes in a series of well-known publications ([30, 31, 33, 34] among others), where the DG methods were extended to non-linear hyperbolic systems of conservation laws.

The DG framework provides the necessary tools that allow to reach arbitrary high order of accuracy in space and time. A natural approach regarding explicit DG methods is to use an explicit time discretization via high order Runge–Kutta schemes, leading to the family of Runge–Kutta discontinuous Galerkin schemes [30, 31, 33, 34]. Space-time DG methods where space-time basis and test functions are employed were first developed in [42, 59, 78, 79]. The resulting schemes are unconditionally stable. For semi-implicit DG schemes, the reader is referred to [15, 40, 70, 77]. An alternative approach to building high-order space-time DG approximations is to consider the Cauchy-Kovaleskaya procedure as the one used in ADER schemes first developed by Toro and Titarev in the finite volume method framework, see ( [16, 71,72,73,74,75]), using Taylor series expansions of the variables while substituting time derivatives by spatial derivatives, which are then treated by the DG scheme. It is well known that the Cauchy-Kovaleskaya (CK) technique often involves complex computations of time derivatives that can quickly become too cumbersome. To overcome these difficulties, one may use a local space-time DG predictor procedure, which allows to reach arbitrary high order in time while completely avoiding the CK procedure, see [37, 41]. The latter formulation of the ADER approach is based on the approximate solution of generalized Riemann problems by means of an element-local fixed point algorithm, the convergence of which was proven in [13], leading naturally to high order fully discrete one-step schemes. This shall be the approximation used in this paper.

DG and ADER-DG have been successfully applied to geophysical flows in the shallow water framework. See for instance the works for the one-layer shallow water model [62, 63, 81, 82] among others. Application of this technique to immiscible two-layer models can be found in [28, 39, 57]. Some works can be found also for multilayer flows in [56]. As far as we know, this is the first attempt to adapt this technique for variable density multilayer shallow water models, obtaining a well-balanced scheme for non-trivial steady states.

Notice that, no limiter is applied up to that point and therefore some way of coping with the Gibbs phenomena associated to non-linear high order schemes must also be provided. In this work, we use a multi-dimensional optimal order detection (MOOD) (see [29, 36]), which is a posteriori approach to the problem of limiting. In this way, the unlimited solution of the ADER-DG method is a candidate solution, which is then evaluated by the MOOD procedure to find its admissibility. If the solution is found inadequate in any given cell, then it will be discarded and recomputed by means of a robust finite volume solver, which would become the corrected solution in the cell. Notice that the use of a finite volume solver could destroy the subcell resolution of the DG scheme. In order to preserve it, an optimal subcell discretization is used, see [44].

One of the main advantages of the MOOD technique is that it allows to analyze different properties of the solutions, and therefore whether it is meaningful from the physical and numerical point of view. For instance, in this paper we consider that the solution must fulfill a relaxed discrete maximum principle to avoid spurious oscillations related to sharp gradients or discontinuities, but also some physical properties like positivity or a relative density equal or greater that a given threshold. The proper set up of these criteria is very important for the admissibility of the solutions, as it will be shown in this work.

Another major issue for the simulation of geophysical flows is the preservation of stationary solutions. Indeed, in practical situations, we often find flows that are small perturbations of an equilibrium state or they converge to one those for long simulation times. Thus, the numerical scheme must be able to exactly preserve this kind of solutions. In this work, we will present a technique to preserve stationary solutions in the framework of the ADER-DG schemes and Riemann solvers. As it will be later discussed, the considered model has several families of stationary solutions of interest, including stationary solutions with non-trivial density profiles.

This paper is organized as follows. In Sect. 2, the model and its key features are presented. The stationary solutions are also discussed in that section. In Sect. 3 the discretization of the model is detailed, including the limiting techniques considered here, alongside the strategy to preserve a certain class of stationary solutions exactly at the discrete level. Section 4 presents different numerical tests in order to analyze the behavior of the proposed numerical solver. Preservation of steady state solutions and a comparison with laboratory data is included. Finally, Sect. 5 gives some conclusions.

2 Model Description

In this section, we briefly recall the governing equations for the system studied here. It is not the aim of this paper to make an exhaustive derivation of the model. The interested reader may refer to [5] or [54] for further details. The model is derived from the incompressible Euler equations with free surface, where two continuity equations are considered: the one corresponding to an incompressible flow with constant density $\rho _0$ and another one corresponding to an advected flow with variable density $\rho _a$. The total density of the flow can be seen as a sum of the density $\rho _0$ plus the density fluctuation, $\rho _1$, with respect to $\rho _0$. The total density of the fluid is then,

$$\begin{aligned} \rho = \rho _0 + \rho _1. \end{aligned}$$

(2.1)

Instead of considering total density, we may introduce the relative density

$$\begin{aligned} \theta = \frac{\rho }{\rho _0} = 1+\frac{\rho _1}{\rho _0}, \end{aligned}$$

(2.2)

resulting in the following Euler equations,

$$\begin{aligned} \nabla \cdot \varvec{v}= & {} 0, \nonumber \\ \partial _t \theta + \nabla \cdot (\theta \, \varvec{v})= & {} 0, \nonumber \\ \partial _t( \theta \, \varvec{v})+ \nabla \cdot (\theta \varvec{v}\otimes \varvec{v})= & {} -g \, \theta \varvec{k}-\frac{1}{\rho _0}\nabla p, \end{aligned}$$

(2.3)

where $\varvec{v}$ is the velocity field, $\varvec{k}$ is the canonical vertical unitary vector, and g is the gravity. For the sake of simplicity, let us focus on the two dimensional (x, z) case.

A multilayer approach is considered, dividing the fluid into M layers, and the equations (2.3) are vertically averaged within each layer, under the usual hypotheses of shallow-water models. We then obtain the following density-dependent pressure expression in each layer $\alpha$,

$$\begin{aligned} p_\alpha (t,x,z) = p_{\alpha +\frac{1}{2}} + \theta _\alpha \,g\,(z_{\alpha +\frac{1}{2}}-z), \end{aligned}$$

(2.4)

with

$$\begin{aligned} p_{\alpha +\frac{1}{2}}(t,x) = p_S(t,x) + g\sum _{\beta =\alpha +1}^{M} \theta _{\beta }h_\beta (t,x). \end{aligned}$$

(2.5)

Here, $f_\alpha$ refers to the value of the variable f in the $\alpha$-layer, whereas $p_{\alpha +\frac{1}{2}}$ is the kinematic pressure at the interface between layers $\alpha$ and $\alpha +1$. Finally, $h_\alpha$ refers to the thickness of layer $\alpha$ and $p_S$ denotes the pressure at the free surface, which is usually set to zero. Figure 1 shows a sketch of the problem.

Fig. 1

Sketch of the multilayer approach

Using this procedure, the following model is obtained from (2.3) (see [49, 64]):

$$\begin{aligned} \left\{ \begin{array}{l} \partial _t h + \partial _x \left( h \sum _{\beta =1}^{M} { l_\beta } u_\beta \right) = 0, \\ \partial _t(h\theta _{\alpha })+\partial _x(h\theta _{\alpha }u_\alpha )= \frac{1}{l_\alpha }\left( \theta _{{\alpha +\frac{1}{2}}} G_{\alpha +\frac{1}{2}}-\theta _{{\alpha -\frac{1}{2}}} G_{\alpha -\frac{1}{2}}\right) , \\ \partial _t(h\theta _{\alpha }u_\alpha ) + \partial _x(h \theta _{\alpha }u_\alpha ^2)+gh\theta _{\alpha }\partial _x \eta +\frac{g\,l_\alpha }{2} \left( h \partial _x (h \theta _\alpha ) - h\theta _\alpha \partial _x h\right) \\ + g \sum _{\beta =\alpha +1}^M l_\beta \left( h \partial _{x} (h\theta _{\beta })-h\theta _\alpha \partial _x h\right) = \frac{1}{l_\alpha }\left( u_{{\alpha +\frac{1}{2}}} \theta _{{\alpha +\frac{1}{2}}} G_{{\alpha +\frac{1}{2}}} - u_{{\alpha -\frac{1}{2}}} \theta _{{\alpha -\frac{1}{2}}} G_{{\alpha -\frac{1}{2}}}\right) , \end{array} \right. \end{aligned}$$

(2.6)

where h is the total height of the water column, $\eta = h+z_b$ is the free surface, and $z_b$ is the bathymetry function. Finally, $u_\alpha$ refers to the horizontal velocity in the $\alpha$-layer and $G_{\alpha \pm \frac{1}{2}}$ $\alpha =1, \cdots M-1$, are the mass transfer terms between layers. Note that the full system (2.6) is obtained under the assumption of some closure relations, which will allow us to properly define the mass transfer terms. More explicitly, it is assumed that the layer thickness is proportional to the total height, $h_\alpha = l_\alpha h$, with $l_\alpha \in [0,1],$ $\alpha =1,\ldots ,M$ such that $\sum _{\alpha =1}^M l_\alpha = 1$. Under this assumption, the expression for the mass transfer terms is given by

$$\begin{aligned} G_{\alpha +\frac{1}{2}} = \sum _{\beta =1}^{\alpha }l_\beta \left( \partial _t h + \partial _x(h u_\beta ) \right) = \sum _{\beta =1}^{\alpha } l_\beta \left( \partial _x(h u_\beta ) - \partial _x \left( \sum _{\gamma =1}^{M} l_\gamma h u_\gamma \right) \right) . \end{aligned}$$

(2.7)

We will assume no mass transfer at the bottom or the free surface, that is $G_{1/2}=G_{M+1/2}=0$. Finally $\theta _{\alpha +1/2}$ and $u_{\alpha +1/2}$ are some approximations of u and $\theta$ at the layer interfaces. Typically,

$$\begin{aligned} u_{\alpha +1/2}=\frac{u_{\alpha +1} + u_{\alpha }}{2},\qquad \theta _{\alpha +1/2}=\frac{\theta _{\alpha +1}+\theta _{\alpha }}{2}, \qquad \alpha =1, \cdots , M-1, \end{aligned}$$

with $u_{1/2}=u_1$, $u_{M+1/2}=u_M$ and $\theta _{1/2}=\theta _1$, $\theta _{M+1/2}=\theta _M$.

Note that the full PDE system given by equations (2.6) has a number of different stationary solutions. The most intuitive ones are probably the steady states corresponding to homogeneous density profiles. For the case of constant density, the standard multilayer shallow-water equations are obtained,

$$\begin{aligned} \left\{ \begin{array}{l} \partial _t h + \partial _x \left( h\sum _{\beta =1}^{M}l_\beta u_\beta \right) = 0, \\ \partial _t(hu_\alpha ) + \partial _x(hu_\alpha ^2)+gh\partial _x \eta = \frac{1}{l_\alpha } \left( u_{\alpha +1/2} G_{\alpha +1/2} - u_{\alpha -1/2} G_{\alpha -1/2}\right) . \end{array} \right. \end{aligned}$$

Therefore, we are able to recover the particular case of the lake-at-rest stationary solutions with constant density,

$$\begin{aligned} \theta _\alpha = \text { const}, \quad u_\alpha = 0, \quad \text {for}\ \alpha =1,\ldots ,M, \quad \eta = \text { const}. \end{aligned}$$

However, system (2.6) also admits lake-at-rest stationary solutions corresponding to non-trivial density profiles. Indeed, stationary solutions with $u_\alpha =0, \, \alpha =1,\ldots ,M$ for the system (2.6) correspond to the solutions of the following ODE system,

$$\begin{aligned} P_\alpha :=gh\theta _{\alpha }\partial _x \eta +\frac{gl_\alpha }{2} \left( h \partial _x (h \theta _\alpha ) - h\theta _\alpha \partial _x h\right) + g \sum _{\beta =\alpha +1}^M l_{\beta }\left( h \partial _{x} (h\theta _{\beta })-h\theta _\alpha \partial _x h\right) =0.\nonumber \\ \end{aligned}$$

(2.8)

Equation (2.8) can be solved recursively by solving first the upper layer,

$$\begin{aligned} \theta _M\partial _{x}\eta + \frac{l_M}{2}h\partial _{x}\theta _M = 0, \end{aligned}$$

and then going downwards throughout the lower layers. Notice that there is an infinite number of solutions for $\eta (x) = h(x) + z_b(x)$ and $\theta _\alpha (x)$ for $\alpha =1,\ldots ,M$. However, we are interested in stationary solutions corresponding to a stratified fluid with constant free surface, that is

$$\begin{aligned} \eta (x) = h(x) + z_b(x) = \text { const}, \qquad \theta (z) = \theta _{surface} + \gamma (\eta -z). \end{aligned}$$

(2.9)

Unfortunately, the vertical discretization into layers of the model (2.6) does not allow a perfectly linear density profile as a stationary solution, unless the bathymetry is flat. Furthermore, it is easy to check that (2.9) is not a solution of (2.8) for a non-trivial bathymetry function. However, system (2.8) can be solved recursively, which results into a stratified density profile that could be seen as an approximation of (2.9) associated to the multilayer approach. In particular, those solutions are given by the following expression,

$$\begin{aligned} \begin{array}{c} u_\alpha = 0, \ \eta (x)=z_b(x)+h(x)=\text { const}, \\ \theta _M(x)= {{\bar{\theta }}}_M \ge 1, \\ \displaystyle {\theta _\alpha (x) = {{\bar{\theta }}}_\alpha \, h^{2(M-\alpha )}(x) + \sum _{\beta =\alpha +1}^{M}S_{2(M-\beta )}(M-\alpha +1) {{\bar{\theta }}}_\beta \, h^{2(M-\beta )}(x)}, \end{array} \end{aligned}$$

(2.10)

with

$$\begin{aligned} S_\beta (\alpha )= & {} (\beta +1)\cdot A_{\frac{\beta +2}{2}+1}(\alpha ),\\ A_p(k)= & {} {\left\{ \begin{array}{ll} 1 \qquad &{}\text {if}\ p\ge k,\\ (p-1)\prod _{\gamma = 2}^{k-p}(1+(p-2)C_{\gamma -1}) \qquad &{}\text {if}\ p<k,\\ \end{array}\right. }\\ C_\gamma= & {} C_{\gamma -1}-\frac{1}{Q_\gamma },\\ Q_\gamma= & {} Q_{\gamma -1}+\gamma +1,\\ C_0= & {} Q_0 = 1, \end{aligned}$$

where ${{\bar{\theta }}}_\alpha$ is a free choice constant parameter corresponding to the constant of integration of the ODE and set a particular profile for the stationary solution. Figures 2 and 3 depicts how the solution (2.10) behaves as we increase the number of layers.

Fig. 2

Solution of the ODE (2.8) for an stratified fluid with $M=5$ (left) and $M=10$ (right)

Fig. 3

Solution of the ODE (2.8) for an stratified fluid with $M=20$ (left) and $M=25$ (right)

Finally, we can give an approximation of the wave propagation speed of the PDE system (2.6). While a complete study of the hyperbolicity of the model (2.6) falls outside the scope of this work and the full spectral behavior of the model is unknown, empirically we observe that the model remains hyperbolic throughout all kind of situations. Furthermore, some previous work by Fernández et al. [12] allows us to consider a bound for the maximum and minimum wave propagation speed $\lambda _{\max }$, $\lambda _{\min }$, given by,

$$\begin{aligned} \lambda _{\min }\ge {{\bar{u}}} - \varPsi , \quad \lambda _{\max } \le {{\bar{u}}} + \varPsi \end{aligned}$$

(2.11)

where

$$\begin{aligned} \varPsi = \sqrt{\frac{2M-1}{2M} \bigg ( 2\sum _{\alpha =1}^{M} ({\bar{u}}-u_\alpha )^2 + gh \bigg ( 1 + \frac{1}{M}\sum _{\beta =1}^{M}(2\beta -1)\theta _\beta \bigg ) \bigg )}, \end{aligned}$$

(2.12)

and

$$\begin{aligned} {\bar{u}} = \frac{1}{M} \sum _{\alpha =1}^{M}u_\alpha . \end{aligned}$$

(2.13)

3 Numerical Scheme

In this section we describe an ADER-DG numerical scheme for (2.6) based on a discretization on uniform meshes with a posteriori subcell finite volume limiter of the MOOD type.

System (2.6) may be expressed in the following form,

$$\begin{aligned} \partial _t \varvec{w} + \partial _x \varvec{F}_C(\varvec{w}) +\varvec{P}(\varvec{w},\partial _x \varvec{w}, \partial _x \eta )-\varvec{T}(\varvec{w},\partial _x \varvec{w})=\varvec{0}, \end{aligned}$$

(3.1)

where $\varvec{w}$ is the vector of the conserved variables,

$$\begin{aligned} \varvec{w} = (h \, |\, h\theta _\alpha \,|\, h\theta _\alpha u_\alpha )^T \in {\mathbb {R}}^{2M+1}. \end{aligned}$$

(3.2)

The physical convective flux $F_C(\varvec{w})$ is defined by,

$$\begin{aligned} F_C(\varvec{w}) = \left( h u_\alpha \, | \, h \theta _\alpha u_\alpha \, | \, h \theta _\alpha \, u_\alpha ^2 \right) ^T \in {\mathbb {R}}^{2M+1}, \end{aligned}$$

(3.3)

and $\varvec{P}(\varvec{w}, \partial _x\varvec{w}, \partial _x \eta )$ corresponds to the pressure term, which depends on the relative density $\theta _\alpha$ and the free surface $\eta$, and has the following form,

$$\begin{aligned} \varvec{P}(\varvec{w},\partial _x {\varvec{w}}, \partial _x \eta )=\left( 0 \, | \, \varvec{0} \, | \, P_\alpha \right) \in {\mathbb {R}}^{2M+1}, \end{aligned}$$

(3.4)

where

$$\begin{aligned} P_\alpha = gh \theta _\alpha \partial _x \eta +\frac{g\,l_\alpha }{2} ( h \partial _x (h \theta _\alpha ) - h\theta _\alpha \partial _x h) + g \sum _{\beta =\alpha +1}^M l_\beta ( h\partial _x (h\theta _\beta ) -h \theta _\alpha \partial _x h). \end{aligned}$$

(3.5)

Finally, the term $\varvec{T}(\varvec{w},\partial _x \varvec{w})$ corresponds to the mass, density, and momentum exchange between layers:

$$\begin{aligned}&\varvec{T}(\varvec{w}, \partial _x \varvec{w})\nonumber \\&\quad = \left( 0 \, \Big | \, \frac{1}{l_\alpha }\left( \theta _{{\alpha +\frac{1}{2}}}G_{{\alpha +\frac{1}{2}}}-\theta _{{\alpha -\frac{1}{2}}}G_{{\alpha -\frac{1}{2}}}\right) \right. \nonumber \\&\left. \quad \Big | \, \frac{1}{l_\alpha }\left( u_{{\alpha +\frac{1}{2}}}\theta _{{\alpha +\frac{1}{2}}}G_{{\alpha +\frac{1}{2}}} - u_{{\alpha -\frac{1}{2}}}\theta _{{\alpha -\frac{1}{2}}} G_{{\alpha -\frac{1}{2}}} \right) \right) ^T \in {\mathbb {R}}^{2M+1}. \end{aligned}$$

(3.6)

We recall that $G_{\alpha +\frac{1}{2}}$ is described by (2.7).

The system of equations (3.1) is solved by applying the family of pure discontinuous Galerkin schemes ${\mathbb {P}}_N{\mathbb {P}}_N$, as described in [37], which provides high order of accuracy in space and time. The numerical scheme can be formulated as a predictor-corrector method: in the predictor step, a high order approximation of the solution of (3.1) is computed by solving a local Cauchy problem in the small, without interaction with the neighbor cells. Next, the corrector step will use this approximation with an explicit DG method that solves the Riemann problem taking into account the neighbor states.

We consider a computational domain $\varOmega$ discretized into a set of conforming elements $T_i=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$, $i=1,\ldots , N_s$, where $N_s$ is the total number of cells with a constant length $\varDelta x = x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$. As usual, the computational grid is the union of all the elements $T_i$,

$$\begin{aligned} \varOmega = \bigcup _{i=1}^{N} T_i. \end{aligned}$$

(3.7)

From now on, we shall use the following notation: for any variable f defined on a volume $T_i$, we denote by $f_{i\pm \frac{1}{2}}$ the values at the interface, depending on whether it is the right or left side of the cell. However, when the values correspond to projected or reconstructed variables into the interface, they will generally be denoted with the super index $f^\pm$, depending on whether they correspond to the left or to the right side of the intercell.

In the following, the discrete solution of the PDE system (3.1) at time $t^n$ is denoted by $\varvec{w}_h(x,t^n)$ and is defined in terms of piecewise polynomials of degree N on the spatial direction. We shall denote by ${\mathcal {U}}_h$ the space of piecewise polynomials up to degree N so that $\varvec{w}_h (\cdot ,t^n) \in {\mathcal {U}}_h$. In this work, a nodal basis defined by the Lagrange interpolation polynomials over the $(N+1)$ Gauss-Legendre quadrature nodes on the element $T_i$ is adopted. As usual in the discontinuous Galerkin (DG) approach, the discrete solution $\varvec{w}_h$ may be discontinuous across the intercells, as in finite volume methods. At each cell $T_i$, the discrete solution is written in terms of the nodal spatial basis functions $\varPhi _l(x)$ and some unknown degrees of freedom $\hat{\varvec{w}}_{i,l}^n$,

$$\begin{aligned} \varvec{w}_h(x,t^n) = \sum \limits _l \hat{\varvec{w}}^n_{i,l} \varPhi _l(x) := \hat{\varvec{w}}_{i,l}^n \varPhi _l(x), \quad \text { for } x \in T_i, \end{aligned}$$

(3.8)

where the Einstein summation convention over two repeated indices has been considered. The spatial basis functions are defined on the reference interval [0, 1] and the transformation from physical coordinates $x \in T_i$ to reference coordinates $\xi \in [0,1]$ is given by the linear mapping $x=x(\xi ) =x_i - \displaystyle \frac{\varDelta x}{2} + \xi \varDelta x$. With this choice, the spatial basis function is written in terms of the nodal basis function $\varphi _k(\xi )$, which satisfy the interpolation property $\varphi _k(\xi _j) = \delta _{kj}$, where $\delta _{kj}$ is the usual Kronecker symbol, $\xi _j$ are the nodal quadrature points, and the resulting basis is by construction orthogonal. Therefore, we write,

$$\begin{aligned} \varPhi _k(x) = \varphi _k(\xi ). \end{aligned}$$

Furthermore, due to this particular choice of a nodal basis, all integral operators can be decomposed into a sequence of one-dimensional operators acting only on the $N+1$ degrees of freedom in each dimension.

In order to derive the ADER-DG method, we first multiply the governing PDE system (3.1) with a test function $\varPhi _k \in {\mathcal {U}}_h$ and integrate over the space-time control volume $T_i \times [t^n,t^{n+1}]$. This leads to

$$\begin{aligned}&\int _{t^n}^{t^{n+1}} \int _{T_i} \varPhi _k\partial _t \varvec{w}\, dxdt \nonumber \\&\quad + \int _{t^n}^{t^{n+1}} \int _{T_i} \varPhi _k \left( \partial _x \varvec{F}_C(\varvec{w}) \,dxdt +\varvec{P}(\varvec{w},\partial _x \varvec{w}, \partial _x \eta )-\varvec{T}(\varvec{w},\partial _x \varvec{w}) \right) \,dx dt=\varvec{0}. \end{aligned}$$

(3.9)

As already mentioned, the discrete solution is allowed to jump across element interfaces, which means that the resulting jump terms have to be properly taken into account. In our scheme this is achieved via numerical flux functions in the form of approximate Riemann solvers that follows the path-conservative approach that was developed by Parés, Castro and collaborators in the finite volume framework [18, 67] and which has later been extended to the discontinuous Galerkin finite element framework in [38, 43, 69].

In classical Runge–Kutta DG schemes [32], only a weak form in space of the PDE is obtained, while time is still kept continuous, thus reducing the problem to a nonlinear system of differential equations, which is subsequently integrated with classical Runge–Kutta methods in time. In the ADER-DG framework, a completely different strategy is used. Here, higher order in time is achieved with the use of an element-local space-time predictor, denoted by ${\mathbf {q}}_h(x,t)$ in the following, and which will be discussed in more detail later. Using (3.8), integrating the first term by parts in time and integrating the flux derivative term by parts in space, taking into account the jumps between elements and making use of this local space-time predictor solution ${\mathbf {q}}_h$ instead of $\varvec{w}_h$, the weak formulation (3.9) can be rewritten as

$$\begin{aligned}&\left( \int _{T_i} \varPhi _k \varPhi _l \, dx \right) \left( \hat{\varvec{w}}^{n+1}_{i,l} - \hat{\varvec{w}}^{n}_{i,l} \right) - \int _{t^n}^{t^{n+1}} \! \! \int _{T_i^\circ } \left( \partial _x \varPhi _k \cdot \varvec{F}_C({\mathbf {q}}_h) \right) \, dx dt \nonumber \\&\quad + \int _{t^n}^{t^{n+1}} \! \! \varPhi _{k,i+\frac{1}{2}} {\mathcal {D}}_{i+\frac{1}{2}}^- \left( {\mathbf {q}}_{h,i+\frac{1}{2}}^-, {\mathbf {q}}_{h,i+\frac{1}{2}}^+,{z_b}_{h,i+\frac{1}{2}}^-,{z_b}_{h,i+\frac{1}{2}}^+ \right) \, dt \nonumber \\&\quad + \int _{t^n}^{t^{n+1}} \! \! \varPhi _{k,i-\frac{1}{2}}{\mathcal {D}}_{i-\frac{1}{2}}^+ \left( {\mathbf {q}}_{h,i-\frac{1}{2}}^-, {\mathbf {q}}_{h,i-\frac{1}{2}}^+,{z_b}_{h,i-\frac{1}{2}}^-,{z_b}_{h,i-\frac{1}{2}}^+ \right) \, dt \nonumber \\&\quad + \int _{t^n}^{t^{n+1}} \! \! \int _{T_i^\circ } \varPhi _k \left( \varvec{P}({\mathbf {q}}_h,\partial _x {\mathbf {q}}_h, \partial _x \eta _h)-\varvec{T}({\mathbf {q}}_h,\partial _x {\mathbf {q}}_h) \right) \, dx dt = \varvec{0}, \end{aligned}$$

(3.10)

where $T_i^\circ$ corresponds to the interior of $T_i$ and $\eta _h$ stands for the projection of $\eta$ onto the space ${\mathcal {U}}_h$. Moreover, ${z_b}_{h,i\pm \frac{1}{2}}^\pm$ are the extrapolated values of the bathymetry at the intercells. Remark that, due to the polynomial representation of $z_b$ at each cell, ${z_b}_{h,i+\frac{1}{2}}^- \ne {z_b}_{h,i+\frac{1}{2}}^+$ in general.

The first integral in (3.10) leads to the element mass matrix, which is diagonal since our basis is orthogonal. Note that the jumps across elements interfaces will be approximated with a Riemann solver, described in Sect. 3.2.

3.1 ADER-DG Space-Time Predictor

As already mentioned, the element-local space-time predictor is an important feature of ADER-DG schemes. The computation of the predictor solution ${\mathbf {q}}_h(x,t)$ is based on a weak formulation of the governing PDE system in space-time starting from the known solution $\varvec{w}_h(x,t^n)$. The PDE system is approximated with a so-called Cauchy problem in the small, i.e. without considering the interaction with the neighbour elements. For an element $T_i$, the predictor solution ${\mathbf {q}}_h$ is now expanded in terms of a local space-time basis,

$$\begin{aligned} {\mathbf {q}}_h(x,t) = \sum \limits _l \theta _l(x,t) \hat{{\mathbf {q}}}^i_l := \theta _l(x,t) \hat{{\mathbf {q}}}^i_l\,, \end{aligned}$$

(3.11)

with the multi-index $l=(l_0,l_1)$ and where the space-time basis functions $\theta _l(x,t) = \varphi _{l_0}(\tau ) \varphi _{l_1}(\xi )$ are again generated from the same one-dimensional nodal basis functions $\varphi _{k}(\xi )$ as before, i.e. the Lagrange interpolation polynomials of degree N passing through $N+1$ Gauss-Legendre quadrature nodes. The spatial mapping $x = x(\xi )$ is also the same as before and the physical time is mapped to the reference time $\tau \in [0,1]$ via $t = t^n + \tau \varDelta t$. Multiplication of the PDE system (3.1) with a test function $\theta _k$ and integration over the space-time control volume $T_i \times [t^n,t^{n+1}]$ yields the following weak form of the governing PDE, which we remark that it is different from (3.9), since now the test and basis functions are both space and time dependent:

$$\begin{aligned}&\int _{t^n}^{t^{n+1}} \int _{T_i} \theta _k(x,t)\partial _t {\mathbf {q}}_h\, dxdt \nonumber \\&\quad + \int _{t^n}^{t^{n+1}} \int _{T_i} \theta _k(x,t) \left( \partial _x \varvec{F}_C({\mathbf {q}}_h) \,dxdt +\varvec{P}({\mathbf {q}}_h,\partial _x {\mathbf {q}}_h, \partial _x \eta _h)\right. \nonumber \\&\left. \quad -\varvec{T}({\mathbf {q}}_h,\partial _x {\mathbf {q}}_h) \right) \,dx dt=\varvec{0}, \end{aligned}$$

(3.12)

where $\eta _h$ represents now the projection of $\eta$ onto the local space-time polynomials, that is, $\eta _h (x,t)= \varPi _1 {\mathbf {q}}_h (x,t)+{z_b}_h(x)$ with $\varPi _1 {\mathbf {q}}_h$ the first component of ${\mathbf {q}}_h$ and ${z_b}_h(x)$ the projection of the bottom topography, which only depends on space, onto the space ${\mathcal {U}}_h$.

Since we are only interested in an element local predictor solution, without considering interactions with the neighbor elements, the jump terms across interfaces are not taken into account at this stage. Instead, it will be done in the final corrector step of the ADER-DG scheme (3.10). This leads to,

$$\begin{aligned}&\int _{T_i} \theta _k(x,t^{n+1}) {\mathbf {q}}_h(x,t^{n+1}) \, dx \nonumber \\&\quad - \int _{T_i} \theta _k(x,t^{n}) {\mathbf {q}}^0_h(x,t^{n}) \,dx - \int _{t^n}^{t^{n+1}} \! \! \int _{T_i} \partial _t \theta _k(x,t) {\mathbf {q}}_h(x,t) \, dx dt \nonumber \\&\quad = - \int _{t^n}^{t^{n+1}} \int _{T_i} \theta _k(x,t) \left( \partial _x \varvec{F}_C({\mathbf {q}}_h) +\varvec{P}({\mathbf {q}}_h,\partial _x {\mathbf {q}}_h, \partial _x \eta _h)\right. \nonumber \\&\left. \quad -\varvec{T}({\mathbf {q}}_h,\partial _x {\mathbf {q}}_h) \right) \,dx dt. \end{aligned}$$

(3.13)

Using the local space-time ansatz (3.11), Eq. (3.13) becomes a local nonlinear system for the unknown degrees of freedom $\hat{{\mathbf {q}}}^i_{l}$ of the space-time polynomials ${\mathbf {q}}_h$. The solution of (3.13) can be easily found via a simple and fast converging fixed point iteration detailed e.g. in [41, 55] and the convergence of which was proven in [13]. For linear homogeneous systems, the iteration converges in a finite number of, at most, $N+1$ steps since the associated iteration matrix is nilpotent [58].

We stress that the choice of an appropriate initial guess ${\mathbf {q}}_h^0(x,t)$ for ${\mathbf {q}}_h(x,t)$ is of great importance in the convergence rate and thus the computational efficiency of the scheme. Several strategies exist to speed up the algorithm: in [83] it is suggested an extrapolation of ${\mathbf {q}}_h$ from the previous time interval $[t^{n-1},t^n]$ while in [55] the authors propose a second-order accurate MUSCL-Hancock-type approach based on discrete derivatives computed at time $t^n$. As an alternative, one can also use a Taylor series expansion of the solution ${\mathbf {q}}_h(x,t)$ about time $t^n$ and then use a continuous extension Runge–Kutta scheme (CERK) in order to generate the initial guess for the space-time predictor, as pointed out in [47]. For details, see [47] and [52, 66].

If an initial guess with polynomial degree $N-1$ in time is chosen, it is sufficient to use a single Picard iteration to solve (3.13) to the desired accuracy, (see [37]). For an efficient task-based formalism of ADER-DG schemes, see [27]. However, for the simulations considered in this work, we consider the initial guess to be $\varvec{q}_h^0(x,t) = \varvec{w}_h(x,t^n)$. This completes the description of the unlimited high order accurate and fully discrete ADER-DG schemes.

3.2 A Hydrostatic Reconstruction Riemann Solver

As in [54], here we use a path-conservative solver based on the hydrostatic reconstruction technique (see [2, 20, 65]). Let us briefly recall its definition,

$$\begin{aligned} {\mathcal {D}}_{i+\frac{1}{2}}^\pm \left( {\mathbf {q}}_{h,i+\frac{1}{2}}^-, {\mathbf {q}}_{h,i+\frac{1}{2}}^+,{z_b}_{h,i+\frac{1}{2}}^-,{z_b}_{h,i+\frac{1}{2}}^+ \right) = \varvec{D}_{i+\frac{1}{2}}^\pm \left( {\mathbf {q}}_{h,i+\frac{1}{2}}^{HR,-}, {\mathbf {q}}_{h,i+\frac{1}{2}}^{HR,+},{z_b}^{HR}_{h,i+\frac{1}{2}},\right) +\varvec{S}^\pm _{i+\frac{1}{2}}, \nonumber \\ \end{aligned}$$

(3.14)

where ${\mathbf {q}}_{h,i+\frac{1}{2}}^{HR,\pm }$ and ${z_b}^{HR}_{h,i+\frac{1}{2}}$ are the hydrostatic reconstructions of the extrapolated values ${\mathbf {q}}_{h,i+\frac{1}{2}}^\pm$ and ${z_b}_{h,i+\frac{1}{2}}^\pm$ while $\varvec{S}_{i+\frac{1}{2}}^\pm$ is the correction introduced in the solver due to the hydrostatic reconstruction procedure.

In order to make the notation easier, we shall neglect in what follows the subindex h, which accounts for the element-local space-time predictor ${\mathbf {q}}_h$.

Now, for any variable f, we shall define its average at the intercell $i+\frac{1}{2}$ as:

$$\begin{aligned} {\overline{f}}\equiv {\overline{f}}_{i+\frac{1}{2}}=\frac{1}{2}\left( f_{i+\frac{1}{2}}^-+f_{i+\frac{1}{2}}^+ \right) . \end{aligned}$$

(3.15)

In particular, for a variable $f_\alpha$ defined within the layer $\alpha$, we shall write

$$\begin{aligned} \overline{f_{\alpha ,}}_{i+\frac{1}{2}} = \frac{1}{2}\left( f_{\alpha ,i+\frac{1}{2}}^-+f_{\alpha ,i+\frac{1}{2}}^+ \right) . \end{aligned}$$

(3.16)

The difference at the intercell $i+\frac{1}{2}$ will be written as

$$\begin{aligned} \varDelta f \equiv (\varDelta f)_{i+\frac{1}{2}} = f_{i+\frac{1}{2}}^+-f_{i+\frac{1}{2}}^-. \end{aligned}$$

(3.17)

We shall denote the average of a variable $f_\alpha$ at the layer interface $\alpha +\frac{1}{2}$ by

$$\begin{aligned} \left\langle {{\overline{f}}} \right\rangle _{{\alpha +\frac{1}{2}},i+\frac{1}{2}}= & {} \frac{1}{2} \left( \overline{ f_{\alpha +1,}}_{i+\frac{1}{2}}+\overline{f_{\alpha ,}}_{i+\frac{1}{2}} \right) \nonumber \\= & {} \frac{1}{4}\left( f_{\alpha ,i+\frac{1}{2}}^-+f_{\alpha ,i+\frac{1}{2}}^+ + f_{\alpha +1,i+\frac{1}{2}}^-+f_{\alpha +1,i+\frac{1}{2}}^+\right) \end{aligned}$$

(3.18)

In addition, at the bottom or free surface interfaces, we shall assume

$$\begin{aligned} \left\langle {{\overline{f}}} \right\rangle _{\frac{1}{2},i+\frac{1}{2}} = \overline{f_{1,}}_{i+\frac{1}{2}}, \qquad \left\langle {{\overline{f}}} \right\rangle _{M+\frac{1}{2},i+\frac{1}{2}}=\overline{f_{M,}}_{i+\frac{1}{2}}. \end{aligned}$$

(3.19)

According to [54], ${\mathbf {q}}_{i+\frac{1}{2}}^{HR,\pm }$ and ${z_b}^{HR}_{i+\frac{1}{2}}$ are defined as follows: given the states ${\mathbf {q}}_{i+\frac{1}{2}}^\pm$ and ${z_b}_{i+\frac{1}{2}}^\pm$, we consider first

$$\begin{aligned} {z_b}^{HR}_{i+\frac{1}{2}}=\max \left( {z_b}_{i+\frac{1}{2}}^-, {z_b}_{i+\frac{1}{2}}^+\right) , \end{aligned}$$

(3.20)

and

$$\begin{aligned} h_{i+\frac{1}{2}}^{HR,\pm }=\left( h_{i+\frac{1}{2}}^\pm + {z_b}_{i+\frac{1}{2}}^\pm -{z_b}^{HR}_{i+\frac{1}{2}}\right) _+, \end{aligned}$$

(3.21)

where $(f)_+$ is the positive part of f. Using (3.21), we define

$$\begin{aligned} {\mathbf {q}}_{i+\frac{1}{2}}^{HR,\pm }=\left( h_{i+\frac{1}{2}}^{HR,\pm } \, \Big | \, h_{i+\frac{1}{2}}^{HR,\pm } \theta _{\alpha ,i+\frac{1}{2}}^\pm \, \Big | \, h_{i+\frac{1}{2}}^{HR,\pm } \theta _{\alpha ,i+\frac{1}{2}}^\pm u_{\alpha ,i+\frac{1}{2}}^\pm \right) ^T \in {\mathbb {R}}^{2M+1}, \end{aligned}$$

(3.22)

where

$$\begin{aligned} \theta _{\alpha ,i+\frac{1}{2}}^\pm =\frac{h\theta _{\alpha , i+\frac{1}{2}}^\pm }{h_{i+\frac{1}{2}}^\pm }, \qquad u_{\alpha ,i+\frac{1}{2}}^\pm =\frac{h\theta u_{\alpha , i+\frac{1}{2}}^\pm }{h\theta _{\alpha ,i+\frac{1}{2}}^\pm }. \end{aligned}$$

Remark that the notation (3.15)–(3.19) are now defined in terms of the reconstructed variables. For instance, (3.15) reads now as follows,

$$\begin{aligned} {\overline{f}}=\frac{1}{2}\left( f_{{i+\frac{1}{2}}}^{HR,+} + f_{{i+\frac{1}{2}}}^{HR,-}\right) . \end{aligned}$$

(3.23)

Once, the reconstructed states ${\mathbf {q}}_{i+\frac{1}{2}}^{HR,\pm }$ are defined, the Riemann solver can be expressed as,

$$\begin{aligned}&\varvec{D}^-_{i+\frac{1}{2}}\left( {\mathbf {q}}_{i+\frac{1}{2}}^{HR,-}, {\mathbf {q}}_{i+\frac{1}{2}}^{HR,+},{z_b}^{HR}_{i+\frac{1}{2}} \right) \nonumber \\&\quad = \frac{1}{2} \left( (1-\alpha _{1,{i+\frac{1}{2}}}) \varvec{E}_{{i+\frac{1}{2}}}-\alpha _{0,{i+\frac{1}{2}}}({\mathbf {q}}_{i+\frac{1}{2}}^{HR,+}-{\mathbf {q}}_{i+\frac{1}{2}}^{HR,-})\right) + \varvec{F}_C({\mathbf {q}}_{i+\frac{1}{2}}^{HR,-}), \end{aligned}$$

(3.24)

and

$$\begin{aligned}&\varvec{D}^+_{i+\frac{1}{2}}\left( {\mathbf {q}}_{i+\frac{1}{2}}^{HR,-}, {\mathbf {q}}_{i+\frac{1}{2}}^{HR,+},{z_b}^{HR}_{i+\frac{1}{2}} \right) \nonumber \\&\quad = \frac{1}{2} \left( (1+\alpha _{1,{i+\frac{1}{2}}}) \varvec{E}_{{i+\frac{1}{2}}}+\alpha _{0,{i+\frac{1}{2}}}({\mathbf {q}}_{i+\frac{1}{2}}^{HR,+}-{\mathbf {q}}_{i+\frac{1}{2}}^{HR,-})\right) - \varvec{F}_C({\mathbf {q}}_{i+\frac{1}{2}}^{HR,+}), \end{aligned}$$

(3.25)

where

$$\begin{aligned} \varvec{E}_{{i+\frac{1}{2}}}=\varvec{F}_C({\mathbf {q}}_{i+\frac{1}{2}}^{HR,+}) - \varvec{F}_{C}({\mathbf {q}}_{i+\frac{1}{2}}^{HR,-})+\varvec{P}_{{i+\frac{1}{2}}} - \varvec{T}_{{i+\frac{1}{2}}}, \end{aligned}$$

(3.26)

with

$$\begin{aligned} \varvec{P}_{{i+\frac{1}{2}}}=\left( \begin{array}{c} 0 \\ \hline \\ \varvec{0}\\ \hline \\ \displaystyle {g \overline{h\theta _\alpha } \varDelta \eta +\frac{g l_\alpha }{2} ( {\overline{h}} \varDelta (h\theta _\alpha ) - \overline{h\theta _\alpha } \varDelta h ) + g \sum _{\beta =\alpha +1}^M l_\beta ({\overline{h}}\varDelta (h\theta _\beta ) - \overline{h\theta _\alpha }\varDelta h)} \\ \end{array} \right) . \end{aligned}$$

Note that in order to make the notation less cumbersome, we have neglected the subindex ${i+\frac{1}{2}}$, so that here ${\overline{h}}$ stands for ${\overline{h}}_{i+\frac{1}{2}}$, $\varDelta h$ stands for $(\varDelta h)_{i+\frac{1}{2}}$ and so on.

The term $\varvec{T}_{{i+\frac{1}{2}}}$ corresponds to the discretization of the transfer terms. As pointed in [54], those terms should be discretized in a vertical upwind manner as follows:

$$\begin{aligned} \varvec{T}_{{i+\frac{1}{2}}} =\left( \begin{array}{c} 0 \\ \hline \\ \displaystyle {\frac{1}{l_\alpha } \left( \left( \theta G \right) ^{UP}_{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}} -\left( \theta G \right) ^{UP}_{{\alpha -\frac{1}{2}},{i+\frac{1}{2}}} \right) } \\ \hline \\ \displaystyle \frac{1}{l_\alpha } \left( \left( u \theta G \right) ^{UP}_{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}} -\left( u\theta G \right) ^{UP}_{{\alpha -\frac{1}{2}},{i+\frac{1}{2}}} \right) \end{array} \right) , \end{aligned}$$

where

$$\begin{aligned} \left( \theta G \right) ^{UP}_{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}}= & {} \left\langle {\overline{\theta }} \right\rangle _{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}}\widetilde{G_{{\alpha +\frac{1}{2}}}} - \frac{1}{2}|\widetilde{G_{{\alpha +\frac{1}{2}}}}|( \overline{\theta _{\alpha +1,}}_{{i+\frac{1}{2}}}-\overline{\theta _{\alpha ,}}_{{i+\frac{1}{2}}}), \\ \left( u\theta G \right) ^{UP}_{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}}= & {} \left\langle \overline{u\theta } \right\rangle _{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}}\widetilde{G_{{\alpha +\frac{1}{2}}}} - \frac{1}{2}|\widetilde{G_{{\alpha +\frac{1}{2}}}}|( \overline{(u\theta )_{\alpha +1,}}_{{i+\frac{1}{2}}}-\overline{(u\theta )_{\alpha ,}}_{{i+\frac{1}{2}}}), \end{aligned}$$

where, $\widetilde{G_{\frac{1}{2}}}=\widetilde{G_{M+\frac{1}{2}}}=0$ and

$$\begin{aligned} \widetilde{G_{{\alpha +\frac{1}{2}}}}=\widetilde{(G_{{\alpha +\frac{1}{2}}})}_{{i+\frac{1}{2}}}= \sum _{\beta =1}^\alpha l_\beta \left( \varDelta (h u_\beta ) -\varDelta \left( \sum _{\gamma =1}^M l_\gamma hu_\gamma \right) \right) \ 1 \le \alpha < M. \end{aligned}$$

The coefficients $\alpha _{0,{i+\frac{1}{2}}}$ and $\alpha _{1,{i+\frac{1}{2}}}$ are related with the numerical viscosity of the scheme and are defined in terms of the upper and lower bounds of the maximum and minimum of the waves speeds (see [17] for more details):

$$\begin{aligned} \alpha _{0,{i+\frac{1}{2}}}= \frac{\lambda _{{i+\frac{1}{2}}}^+ |\lambda _{{i+\frac{1}{2}}}^-|-\lambda _{{i+\frac{1}{2}}}^-|\lambda _{{i+\frac{1}{2}}}^+|}{\lambda ^+_{{i+\frac{1}{2}}}- \lambda ^-_{{i+\frac{1}{2}}}}, \qquad \alpha _{1,{i+\frac{1}{2}}}=\frac{|\lambda ^+_{{i+\frac{1}{2}}}|-|\lambda ^-_{{i+\frac{1}{2}}}|}{\lambda ^+_{{i+\frac{1}{2}}}-\lambda ^-_{{i+\frac{1}{2}}}}, \end{aligned}$$

where

$$\begin{aligned} \lambda _{{i+\frac{1}{2}}}^\pm ={{\bar{u}}}_{{i+\frac{1}{2}}}\pm \varPsi _{{i+\frac{1}{2}}}, \end{aligned}$$

with

$$\begin{aligned} {{\bar{u}}}_{{i+\frac{1}{2}}}=\frac{1}{M}\sum _{\alpha =1}^M \overline{u_{\alpha ,}}_{{i+\frac{1}{2}}}, \end{aligned}$$

and

$$\begin{aligned} \varPsi _{{i+\frac{1}{2}}} = \sqrt{\frac{2M-1}{2M} \bigg ( 2\sum _{\alpha =1}^{M} ({\bar{u}}_{{i+\frac{1}{2}}}-\overline{u_{\alpha ,}}_{{i+\frac{1}{2}}})^2 + g{\overline{h}}_{{i+\frac{1}{2}}} \bigg ( 1 + \frac{1}{M}\sum _{\beta =1}^{M}(2\beta -1)\overline{\theta _{\beta ,}}_{{i+\frac{1}{2}}} \bigg ) \bigg )}. \end{aligned}$$

Finally, the terms $\varvec{S}_{{i+\frac{1}{2}}}^\pm$ correspond to the corrections due to the use of the hydrostatic reconstruction and guarantee the consistency of the scheme as well as the well-balanced property [8] for the lake-at-rest with constant density steady states. $\varvec{S}_{{i+\frac{1}{2}}}^\pm$ are defined as,

$$\begin{aligned} \varvec{S}_{{i+\frac{1}{2}}}^\pm =\varvec{P}_{{i+\frac{1}{2}}}^\pm -\varvec{T}_{{i+\frac{1}{2}}}^\pm , \end{aligned}$$

where

$$\begin{aligned} \varvec{P}_{{i+\frac{1}{2}}}^\pm =\left( 0 \, \Big | \, \varvec{0} \, \Big | \, P_{\alpha ,{i+\frac{1}{2}}}^\pm \right) ^T \in {\mathbb {R}}^{2M+1}, \end{aligned}$$

(3.27)

with

$$\begin{aligned} P_{\alpha ,{i+\frac{1}{2}}}^\pm = \displaystyle {g \sum _{\beta =\alpha +1}^M l_\beta \overline{h^\pm }\varDelta h^\pm (\theta ^\pm _{\beta ,i+\frac{1}{2}} -\theta ^\pm _{\alpha ,i+\frac{1}{2}} )}, \end{aligned}$$

where

$$\begin{aligned} \overline{f^+}=\displaystyle \frac{f_{{i+\frac{1}{2}}}^{HR,+}+f_{i+\frac{1}{2}}^+}{2}, \qquad \overline{f^-}=\displaystyle \frac{f_{i+\frac{1}{2}}^-+f_{{i+\frac{1}{2}}}^{HR,-}}{2}, \end{aligned}$$

(3.28)

and

$$\begin{aligned} \varDelta f^+ = f_{i+\frac{1}{2}}^+ - f_{{i+\frac{1}{2}}}^{HR,+}, \qquad \varDelta f^- = f_{{i+\frac{1}{2}}}^{HR,-} - f^-_{i+\frac{1}{2}} . \end{aligned}$$

(3.29)

Note that the simplified expression (3.27) derives from the fact that, in the hydrostatic reconstruction, the primitive variables $u_\alpha ,$ $\theta _\alpha$, and $\eta$ remain constant.

The term $\varvec{T}^\pm _{{i+\frac{1}{2}}}$ is defined by

$$\begin{aligned} \varvec{T}^\pm _{{i+\frac{1}{2}}} =\left( \begin{array}{c} 0 \\ \hline \\ \displaystyle {\frac{1}{l_\alpha } \left( \left( \theta G \right) ^{UP,\pm }_{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}} -\left( \theta G \right) ^{UP,\pm }_{{\alpha -\frac{1}{2}},{i+\frac{1}{2}}} \right) } \\ \hline \\ \displaystyle \frac{1}{l_\alpha } \left( \left( u \theta G \right) ^{UP,\pm }_{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}} -\left( u\theta G \right) ^{UP,\pm }_{{\alpha -\frac{1}{2}},{i+\frac{1}{2}}} \right) \end{array} \right) , \end{aligned}$$

where,

$$\begin{aligned} \left( \theta G \right) ^{UP,\pm }_{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}}= & {} \left\langle \theta \right\rangle ^\pm _{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}}\widetilde{G_{{\alpha +\frac{1}{2}}}}^\pm - \frac{1}{2}|\widetilde{G_{{\alpha +\frac{1}{2}}}}^\pm |( \theta _{\alpha +1,{i+\frac{1}{2}}}^\pm -\theta _{\alpha ,{i+\frac{1}{2}}}^\pm ), \\ \left( u\theta G \right) ^{UP,\pm }_{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}}= & {} \left\langle u\theta \right\rangle ^\pm _{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}}\widetilde{G_{{\alpha +\frac{1}{2}}}}^\pm - \frac{1}{2}|\widetilde{G_{{\alpha +\frac{1}{2}}}}^\pm |\left( (u\theta )_{\alpha +1,{i+\frac{1}{2}}}^\pm -(u\theta )_{\alpha ,{i+\frac{1}{2}}}^\pm \right) , \end{aligned}$$

and

$$\begin{aligned} \left\langle \theta \right\rangle ^\pm _{{\alpha +\frac{1}{2}},i+\frac{1}{2}}=\frac{1}{2}\left( \theta ^\pm _{\alpha ,i+\frac{1}{2}}+\theta ^\pm _{\alpha +1,i+\frac{1}{2}} \right) . \end{aligned}$$

(3.30)

Similar for $\left\langle u\theta \right\rangle ^\pm _{{\alpha +\frac{1}{2}},{i+\frac{1}{2}}}$. Finally, $\widetilde{G_{{\alpha +\frac{1}{2}}}}^\pm$ is defined as

$$\begin{aligned} \widetilde{G_{{\alpha +\frac{1}{2}}}}^\pm = \sum _{\beta =1}^\alpha l_\beta \varDelta h^\pm \left( u^\pm _{\beta ,i+\frac{1}{2}}-\sum _{\gamma =1}^M l_\gamma u^\pm _{\gamma ,i+\frac{1}{2}}\right) , \end{aligned}$$

and $\widetilde{G_{\frac{1}{2}}}^\pm =\widetilde{G_{M+\frac{1}{2}}}^\pm =0$.

It is then straightforward to check that the numerical scheme (3.10) is exactly well-balanced for the solutions corresponding to water at rest with constant density.

3.3 Upwind Transference Terms at the Volume Integrals

As discussed before, according to [54], a direct discretization of transference terms may yield non-physical results, like relative densities below the unity, which directly conflicts with (2.2). Indeed, the transference terms can be interpreted as a vertical flux between layers, but by writing them as non-conservaive products, the information relative to the flux direction is lost. In [5], it is suggested that this behavior can be avoided by an upwind approximation of the transference terms that incorporate this information into the discretization.

Therefore, we perform a similar discretization as the one proposed at the intercells, that is, we consider a vertical upwinding discretization of those terms at each degree of freedom l using the following discretization,

$$\begin{aligned} \varvec{T}_{l} =\left( \begin{array}{c} 0 \\ \hline \\ \displaystyle { \left( \theta G \right) ^{UP}_{{\alpha +\frac{1}{2}},l} -\left( \theta G \right) ^{UP}_{{\alpha -\frac{1}{2}},l} } \\ \hline \\ \displaystyle \left( u \theta G \right) ^{UP}_{{\alpha +\frac{1}{2}},l} -\left( u\theta G \right) ^{UP}_{{\alpha -\frac{1}{2}},l} \end{array} \right) , \end{aligned}$$

where

$$\begin{aligned} \left( \theta G \right) ^{UP}_{{\alpha +\frac{1}{2}},l}= & {} \left\langle \theta \right\rangle _{{\alpha +\frac{1}{2}},l} {G_{{\alpha +\frac{1}{2}}}} - \frac{1}{2}|{G_{{\alpha +\frac{1}{2}}}}|(\theta _{\alpha +1,l}-\theta _{\alpha ,l}), \\ \left( u\theta G \right) ^{UP}_{{\alpha +\frac{1}{2}},l}= & {} \left\langle u\theta \right\rangle _{{\alpha +\frac{1}{2}},l}G_{{\alpha +\frac{1}{2}}} - \frac{1}{2}|G_{{\alpha +\frac{1}{2}}}|((u\theta )_{\alpha +1,l}-(u\theta )_{\alpha ,l}), \end{aligned}$$

with $\left( \theta G \right) ^{UP}_{\frac{1}{2},l} = \left( \theta G \right) ^{UP}_{M+\frac{1}{2},l}=\left( u \theta G \right) ^{UP}_{\frac{1}{2},l} =\left( u \theta G \right) ^{UP}_{M+\frac{1}{2},l}=0,$ and we also denote,

$$\begin{aligned} \left\langle \theta \right\rangle _{{\alpha +\frac{1}{2}},l} = \frac{1}{2}\left( \theta _{\alpha +1,l} + \theta _{\alpha ,l} \right) \end{aligned}$$

and analogously for $\left\langle u\theta \right\rangle _{{\alpha +\frac{1}{2}},l}$.

The term $G_{\alpha +\frac{1}{2}}$, defined in (2.7), is approximated according to the the polynomial basis used in the spatial discretization.

3.4 A Posteriori Subcell Finite Volume Limiter

The numerical discretization proposed before is unlimited, in the sense that there is no mechanism that prevents the appearance of Gibbs oscillations associated with high order in the proximity of shocks or steep gradients. Certainly, there is no need for any measure when the solution is smooth, but providing a non-linear limiter in the presence of shock waves or discontinuities is essential. The limiter used in this work follows the ideas presented in [44, 60, 80]. The limiting procedure acts in two steps: first it detects which cells need limiting, and second it introduces some kind of numerical viscosity into the solution in these regions.

The unlimited solution of the governing equation given by (3.13) is now considered as a candidate solution, $\varvec{w}_h^*(x,t^{n+1})$. This candidate solution will remain unchanged if it is considered acceptable. However, the candidate solution will be overridden when it does not satisfy some selection criteria. The cells where this occurs are denominated troubled cells. In the cells that are deemed troubled, we recompute the solution with a fully discrete second order accurate MUSCL-Hancock finite volume scheme. In order to do this, we first project the solution $\varvec{w}_h(x,t^{n})$ into a subgrid of $K_{s}$ elements in $T_i$ and denoted by ${\mathcal {S}}_{i,j}$ and verifying $T_i=\bigcup _j {\mathcal {S}}_{i,j}$ for $j=1,\ldots ,K_{s}$. The projection can be seen as an alternative data representation, denoted by $\varvec{v}_h(x,t^{n})$, and consists of a set of piecewise constant subcell averages. These averages are a $L_2$ projection that preserves the mean of $\varvec{w}_h(x,t^{n})$ in ${\mathcal {S}}_{i,j}$,

$$\begin{aligned} \varvec{v}_h(x,t^{n}) = \frac{1}{|{\mathcal {S}}_{i,j}|} \int _{{\mathcal {S}}_{i,j}} \varvec{w}_h(x,t^{n})\, dx, \qquad \forall x \in {\mathcal {S}}_{i,j} \subset T_i. \end{aligned}$$

(3.31)

The subcell averages $\varvec{v}_h(x,t^{n})$ are evolved with an explicit second order finite volume solver. The solution of this solver, $\varvec{v}_h(x,t^{n+1})$, has to be reconstructed into a limited DG polynomial. This is achieved through a classical least square reconstruction operator that preserves the average of the projected solutions,

$$\begin{aligned} \int _{{\mathcal {S}}_{i,j}} \varvec{w}_h(x,t^{n+1}) \, dx = \int _{{\mathcal {S}}_{i,j}} \varvec{v}_h(x,t^{n+1}) dx, \qquad {\mathcal {S}}_{i,j}\subset T_i. \end{aligned}$$

(3.32)

Since a subcell resolution $K_{s}>N+1$ is admitted, this problem may be overdetermined. Hence, the following constraint must also be considered,

$$\begin{aligned} \int _{T_i} \varvec{w}_h(x,t^{n+1}) \, dx = \int _{T_i} \varvec{v}_h(x,t^{n+1}) dx. \end{aligned}$$

(3.33)

In this way, the final solution $\varvec{w}_h(x,t^{n+1})$ is the candidate solution in regions where the solution is deemed adequate and the reconstructed solution in the regions that it is not. Note that $\varvec{v}_h(x,t^{n+1})$ is kept in memory in case that its cell $T_i$ is again considered troubled in the next time step. In this case, the finite volume solver uses $\varvec{v}_h(x,t^{n+1})$ as initial data and not the projected DG polynomial (3.32).

The main advantage of this approach to limitation is that it allows to verify any number of physical and numerical properties for the numerical scheme. Indeed, the nature of the limiter, which evaluates the suitability of the candidate solution $\varvec{w}_h^*(x,t^{n+1})$ a posteriori, allows to consider the cell as troubled not only in the presence of discontinuities, but also if the solution does not fulfill some prescribed properties. In our particular case, we set a numerical criteria and two physical ones.

Physical admissibility criteria The physical admissibility of the solution must be tightly correlated with the system of governing PDE (3.1). We set that the candidate solution, $\varvec{w}_h^*(x,t^{n+1})$, has to fulfill positivity of water column height, h, and that the relative density, $\theta _\alpha$, must always be equal or greater that one (see 2.2).

Numerical admissibility criteria A relaxed discrete maximum principle, adapted to polynomials, is used to detect discontinuities (see [44]). The criteria is applied in a posteriori manner as follows,

$$\begin{aligned} \min _{y\in {\mathcal {V}}_i} (\varvec{v}_h(y,t^n)) - \varvec{\delta } \le \varvec{v}_h(x,t^{n+1}) \le \max _{y\in {\mathcal {V}}_i} (\varvec{v}_h(y,t^n)) + \varvec{\delta }, \qquad \forall x \in T_i \end{aligned}$$

(3.34)

where the projection $\varvec{v}_h(x,t^{n+1})$ is used as a discrete form of the polynomial $\varvec{w}_h(x,t^{n+1})$, ${\mathcal {V}}_i$ is a set containing $T_i$ and its neighbors cells and $\varvec{\delta }$ is a small value that relaxes the criteria to allow some very small overshoot or undershoot and avoid roundoff errors that would arise if (3.34) is applied strictly.

A suitable number of subcells, $K_{s}$, also has to be chosen. In particular, we choose an optimal subgrid satisfying $K_{s} = 2N+1$. This choice is optimal in the sense that it allows to keep the same time step calculated for the DG polynomial and also to have a CFL number close to the theoretical maximum for the finite volume numerical scheme (see [44]). Note that in the case of a troubled cell next to a non-troubled cell, there will be a non-consistent numerical flux caused by the two different methods used in the troubled and non-troubled cell. Thus, it is important to update the flux in the non-troubled cell to be consistent with the flux calculated in the troubled cell, and keep intact the conservation properties of the numerical scheme.

3.5 Preserving Stationary Solutions in the ADER-DG Framework

In [21], the authors have proposed a general procedure to construct well-balanced high-order finite volumes schemes, i.e. numerical methods that satisfy the so-called C-property found in [8]. In this section, we propose to adapt the strategy described in [21] to define a well-balanced DG numerical scheme that exactly preserves the stationary solutions given in (2.10), that corresponds to a stationary stratified fluid.

The first step in this algorithm is to determine a local stationary solution of the family (2.10) at each cell denoted by $\varvec{u}_{e,i}(x), \ x \in T_i$ at each time step. Let us remark that the stationary solution is computed in every cell $T_i$ at each time step $t^n$. Nevertheless, we do not write the dependence on $t^n$ in order to make the notation more readable. Notice that $\varvec{u}_{e,i}$ is fully determined by $h_{e,i}$ and $\theta _{1,e,i}, \ldots \theta _{M,e,i}$, since stationary solutions (2.10) assumes $u_{\alpha ,e,i}=0$, $1 \le \alpha \le M$. According to (2.10), $\varvec{u}_{e,i}$ are determined by setting $M+1$ locally computed constants, denoted by ${{\bar{\eta }}}_i$, ${{\bar{\theta }}}_{\alpha ,i}$, $1 \le \alpha \le M$. In particular, ${{\bar{\eta }}}_i$

$$\begin{aligned} {{\bar{\eta }}}_{i}=\frac{1}{\varDelta x} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \left( h_h(x,t^n)+{z_b}_{h}(x) \right) \, dx, \end{aligned}$$

where we have denoted by $f_h$ the discrete representation of f onto the polynomial space ${\mathcal {U}}_h$. Similarly, the constants $\theta _{1,e,i}, \ldots \theta _{M,e,i}$ are computed as follows,

$$\begin{aligned} \frac{1}{\varDelta x}\int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} (h\theta )_{\alpha ,e,i}(x,{{\bar{\eta }}}_i, {{\bar{\theta }}}_{\alpha ,i}, \cdots , {{\bar{\theta }}}_{1,i}) \, dx=\frac{1}{\varDelta x} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} (h\theta )_{h,\alpha }(x,t^n) \, dx, \end{aligned}$$

for $1 \le \alpha \le M$. Once these constants are computed, the local stationary solution $\varvec{u}_{e,i}(x)$ is determined at each cell. Note that the local stationary solutions satisfy the pressure terms (2.8) at each cell,

$$\begin{aligned} \varvec{P}(\varvec{u}_{e,i},\partial _x \varvec{u}_{e,i},\partial _x{{\bar{\eta }}}_i)=0. \end{aligned}$$

(3.35)

Now, we can replace the numerical scheme (3.10) by the following equivalent well-balanced ADER-DG scheme,

$$\begin{aligned}&\left( \int _{T_i} \varPhi _k \varPhi _l \, dx \right) \left( \hat{\varvec{w}}^{n+1}_{i,l} - \hat{\varvec{w}}^{n}_{i,l} \right) - \int _{t^n}^{t^{n+1}} \! \! \int _{T_i^\circ } \left( \partial _x \varPhi _k \cdot \varvec{F}_C({\mathbf {q}}_h) - \varPhi _k \varvec{T}({\mathbf {q}}_h,\partial _x {\mathbf {q}}_h) \right) \, dx dt \nonumber \\&\quad + \int _{t^n}^{t^{n+1}} \! \! \varPhi _{k,i+\frac{1}{2}} {\mathcal {D}}_{i+\frac{1}{2}}^- \left( {\mathbf {q}}_{h,i+\frac{1}{2}}^-, {\mathbf {q}}_{h,i+\frac{1}{2}}^+,{z_b}_{h,i+\frac{1}{2}}^-,{z_b}_{h,i+\frac{1}{2}}^+ \right) \, dt \nonumber \\&\quad + \int _{t^n}^{t^{n+1}} \! \! \varPhi _{k,i-\frac{1}{2}}{\mathcal {D}}_{i-\frac{1}{2}}^+ \left( {\mathbf {q}}_{h,i-\frac{1}{2}}^-, {\mathbf {q}}_{h,i-\frac{1}{2}}^+,{z_b}_{h,i-\frac{1}{2}}^-,{z_b}_{h,i-\frac{1}{2}}^+ \right) \, dt \nonumber \\&\quad + \int _{t^n}^{t^{n+1}} \! \! \int _{T_i^\circ } \varPhi _k \left( \varvec{P}({\mathbf {q}}_h,\partial _x {\mathbf {q}}_h, \partial _x \eta _h)\right. \nonumber \\&\left. \quad -\varvec{P}((\varvec{u}_{e,i})_h,\partial _x(\varvec{u}_{e,i})_h , \partial _x (\eta _{e,i})_h)\right) \, dx dt = \varvec{0}, \end{aligned}$$

(3.36)

where $(\varvec{u}_{e,i})_h$ and $\partial _x(\varvec{u}_{e,i})_h$ are respectively the projections of $\varvec{u}_{e,i}$ and $\partial _x \varvec{u}_{e,i}$ into ${\mathcal {U}}_h$.

Moreover, the extrapolated values at the cell interfaces, denoted by ${\mathbf {q}}_{h,i\pm \frac{1}{2}}^\pm$, are computed in the following way,

$$\begin{aligned} {\mathbf {q}}_{h,i+\frac{1}{2}}^-=\varvec{u}_{e,i}(x_{i+\frac{1}{2}})+\hat{\varvec{q}}_{h,i+\frac{1}{2}}^-, \end{aligned}$$

(3.37)

where $\hat{\varvec{q}}_{h,i+\frac{1}{2}}^-$ is the extrapolation on the cell interface of the fluctuation $(\varvec{q}_{h,i}-(\varvec{u}_{e,i})_h)$, that is,

$$\begin{aligned} \hat{\varvec{q}}_{h,i+\frac{1}{2}}^-=(\varvec{q}_{h,i}-(\varvec{u}_{e,i})_h)(x_{i+\frac{1}{2}}). \end{aligned}$$

Similarly,

$$\begin{aligned} {\mathbf {q}}_{h,i+\frac{1}{2}}^+=\varvec{u}_{e,i+1}(x_{i+\frac{1}{2}})+\hat{\varvec{q}}_{h,i+\frac{1}{2}}^-, \end{aligned}$$

(3.38)

where

$$\begin{aligned} \hat{\varvec{q}}_{h,i+\frac{1}{2}}^+=(\varvec{q}_{h,i+1}-(\varvec{u}_{e,i+1})_h)(x_{i+\frac{1}{2}}). \end{aligned}$$

Finally, ${z_b}_{h,i+\frac{1}{2}}^\pm =z_b(x_{i+\frac{1}{2}})$. Note that if the bathymetry is a smooth function, it remains continuous at the cell interfaces and the hydrostatic reconstruction is no longer needed.

The ADER space-time predictor step is also modified according to this previous idea: an iterative algorithm that computes the fluctuation with respect to the local stationary solution $\varvec{u}_{e,i}$ is proposed. In particular, the following well-balanced space-time predictor is considered,

$$\begin{aligned}&\int _{T_i} \theta _k(x,t^{n+1}) {\hat{{\mathbf {q}}}}_h(x,t^{n+1}) \, dx \nonumber \\&\quad - \int _{T_i} \theta _k(x,t^{n}) {\hat{{\mathbf {q}}}}^0_h(x,t^{n}) \,dx - \int _{t^n}^{t^{n+1}} \! \! \int _{T_i} \partial _t \theta _k(x,t) {\hat{{\mathbf {q}}}}_h(x,t) \, dx dt \nonumber \\&\quad = - \int _{t^n}^{t^{n+1}} \int _{T_i} \theta _k(x,t) \left( \partial _x \varvec{F}_C({\mathbf {q}}_h) -\varvec{T}({\mathbf {q}}_h,\partial _x {\mathbf {q}}_h) \right) \,dx dt \nonumber \\&\quad -\int _{t^n}^{t^{n+1}} \int _{T_i} \theta _k(x,t) \left( \varvec{P}({\mathbf {q}}_h,\partial _x {\mathbf {q}}_h, \partial _x \eta _h)-\varvec{P}((\varvec{u}_{e,i})_h,\partial _x (\varvec{u}_{e,i})_h, \partial _x (\eta _{e,i})_h) \right) \,dx dt,\nonumber \\ \end{aligned}$$

(3.39)

where ${\hat{{\mathbf {q}}}}_h(x,t)$ is the projection of the fluctuation about the stationary solution $\varvec{u}_{e,i}$, $t \in [t^n, t^{n+1}]$,

$$\begin{aligned} {\hat{{\mathbf {q}}}}_h(x,t)=\varvec{q}_h(x,t)-(\varvec{u}_{e,i})_h(x). \end{aligned}$$

Finally, to clean possible spurious oscillations due to the absence of numerical viscosity in stationary solutions, we apply the following procedure: first we compute the average of the fluctuation with respect to the local stationary solution,

$$\begin{aligned} \widehat{\varvec{w}}_{h,i}=\frac{1}{\varDelta x}\int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \varvec{w}_h(x,t^n)-\varvec{u}_{e,i}(x) \, dx. \end{aligned}$$

if $\left| \widehat{\varvec{w}}_{h,i}\right| <\delta _\epsilon$, then $\varvec{w}_h(x,t^n)$ is redefined as follows,

$$\begin{aligned} \varvec{w}_h(x,t^n)=\varvec{u}_{e,i}(x)+\widehat{\varvec{w}}_{h,i}. \end{aligned}$$

Some numerical tests applying the techniques presented in this section can be found in Sect. 4.

3.6 Time Step Restriction

The DG method has a rather severe Courant-Friedrichs-Lewy (CFL) number that decreases with the order of the approximation polynomial, N. The time step for the explicit DG method in a multiple dimension framework follows,

$$\begin{aligned} \varDelta t \le \frac{1}{d} \, \frac{\text {CFL}}{2N+1} \, \frac{\varDelta x}{|\lambda _{max}|}, \end{aligned}$$

(3.40)

with $\text {CFL}<1$ and where d is the dimension space (here $d=1$) and $|\lambda _{max}|$ is an approximation of the maximum wave speed, given by (2.11). This restrictive CFL condition has to be seen with perspective. While it is true that it imposes a small and restrictive time step, the enhanced subcell resolution of the DG scheme allows coarse or even very coarse meshes, compensating the severe $\frac{1}{2N+1}$ restriction.

4 Numerical Tests

In this section we present some numerical tests in order to show the efficiency of the new well-balanced ADER-DG solver for multi-layer shallow water systems described previously. In particular, we consider a standard test to measure the order of accuracy of the numerical scheme, several test cases where the well-balanced properties of the scheme presented in Sect. 3.5 are shown, including a challenging test corresponding to a density stratified fluid with non-trivial density profiles. Additionally, a comparison of the numerical solution of the model and laboratory data of a simulation corresponding to an evolution of an internal density wave is also included. Finally, two more simulations of an initially smooth density profile and a highly discontinuous lock-exchange test are also considered.

Concerning the boundary conditions, it is well-known that the imposition of wall-type and/or open boundary condition is difficult for DG solvers. In this work, the numerical tests where open-boundary conditions is considered, we simply mark the boundary cells as always troubled and impose the boundary conditions using the second-order FV scheme. Finally, throughout all the numerical test presented in this section, the CFL parameter is set to 0.5.

4.1 Accuracy Test

We now seek to check numerically the order of accuracy of the numerical scheme. Since the proposed ADER-DG method is arbitrary high order in space and time, a second, third and fourth order version have been chosen for this test. We consider a 10 m long channel on the interval $[-5,5]$ with 5 vertical layers, discretized with 25, 50, 100, 200 and 400 uniform cells. The initial condition is given by:

$$\begin{aligned} u_\alpha (x,0)= & {} 0, \quad \theta _\alpha (x,0) = 1 + \frac{1}{100}\,e^{-5 x^2}, \quad 1 \le \alpha \le 5, \\ \eta (x,0)= & {} 1+\frac{1}{10}\,e^{-10\,x^2}, \end{aligned}$$

and the bathymetry is given by,

$$\begin{aligned} z_b(x)=\frac{1}{2}\,e^{-x^2}. \end{aligned}$$

The final simulation time is $t=0.5$ s and periodic boundary conditions are considered. The free surface and the velocity are depicted in Fig. 4 at the final simulation time $t=0.5$ s for the fourth order scheme using the 200 cells mesh.

Tables 1, 2 and 3 show the errors and order of accuracy for the conserved variables h, $h\theta _1$, and $h\theta _1 u_1$. These variables have been chosen since they correspond to the bottom layer and thus the pressure and transfer terms are more relevant. Similar results are obtained for the other unknowns. The numerical solutions are compared to a reference solution, which has been computed with the same numerical scheme on a finer mesh of 2400 points. As expected, the numerical accuracy is achieved.

Fig. 4

Accuracy test for $t=0.5$ s. Free surface (left) and velocity profiles (right)

Table 1 Order of accuracy for the second order scheme

Table 2 Order of accuracy for the third order scheme

Table 3 Order of accuracy for the fourth order scheme

4.2 Well-Balanced Property

We now intend to validate the well-balanced property of the scheme thanks to the techniques introduced in Sect. 3.5. Three test cases are considered: a classical lake at rest solution with constant density, a lake at rest solution with a non-trivial density profile and a perturbation of the free surface of the preceding test.

For the first one, we consider again a 10m long channel within the interval $[-5,5]$ and 5 vertical layers of constant density $\theta _\alpha =1$ for $\alpha =1,\ldots ,M$. As initial condition, we fix a constant free surface $\eta =2$, and bottom topography given by,

$$\begin{aligned} z_b(x) = \frac{1}{2} e^{-x^2}. \end{aligned}$$

(4.1)

The simulation is run up to $t=500$ seconds on a uniform mesh composed of 50 cells and a fourth order ADER-DG scheme in space and time is considered. Figures 5 and 6 show the surface and density profiles, respectively. As we can see, the scheme preserves the steady state solution up to machine precision for long simulations times.

Fig. 5

Lake-at-rest steady state with constant density. Left: initial condition. Right: free surface at $t=500$ s

Fig. 6

Lake-at-rest steady state with constant density. Relative density (left) and velocities (right) at $t=500$ s

More sophisticated stationary solutions can also be preserved, as described in Sect. 2. As it was discussed, system (2.6) also admits stationary solutions corresponding to steady lake-at-rest states with non-constant density profile. In particular, let us consider the bathymetry function (4.1) and the lake-at-rest steady state with constant free surface $\eta =2$ and the following $\theta _\alpha$ profiles for the particular case of three layers,

$$\begin{aligned} \theta _1= & {} {\bar{\theta }}_1 + 3\,{\bar{\theta }}_2\, h(x)^2, \end{aligned}$$

(4.2)

$$\begin{aligned} \theta _2= & {} {\bar{\theta }}_1 + {\bar{\theta }}_2 \, h(x)^2, \end{aligned}$$

(4.3)

$$\begin{aligned} \theta _3= & {} {\bar{\theta }}_1, \end{aligned}$$

(4.4)

where ${\bar{\theta }}_1=1.01$ and ${\bar{\theta }}_2=0.02$. The ability of the scheme to preserve such a stationary solution is tested in the following. To do so, we consider a fourth order ADER-DG scheme with 100 uniform cells in the interval $[-5,5]$, and we compute the numerical solution up to time 500 seconds. Open boundary conditions are set. Figures 7 and 8 show that the steady state is indeed preserved up to machine precision (Fig. 9).

Fig. 7

Lake-at-rest steady state with non-constant density profile. Left: surface and bottom. Right: relative density for each layer

Fig. 8

Difference between computed solution at time $t=500$ seconds and the original steady state for a lake-at-rest steady state with non-constant density profile: Left: difference on the relative densities. Right: difference on the velocities

Fig. 9

Free surface at time $t=0$ s (left) and $t=500$ s (right) for a lake-at-rest steady state with non-constant density profile

Finally, we consider a small perturbation of the previously described stationary lake-at-rest solution with non-constant density profile. The same bathymetry function, relative density profiles and number of cells are considered, but the free surface is now given by,

$$\begin{aligned} \eta (x,0) = 2 + \frac{1}{10} e^{-5\,x^2}. \end{aligned}$$

Wall boundary conditions are set.

The results are shown in Figs. 10, 11, 12, 13. As it can be seen, the numerical scheme reaches another stationary solution, as wall boundary conditions are set. As expected, the stationary free surface achieved at time $t=500$ seconds has increased, due to the initial perturbation (see 12). The final velocity profiles could be seen in Fig. 12, and they are close to 0. Finally, a new stratified solution is achieved.

Fig. 10

Perturbation of a steady state with non-constant density profile at $t=0$ seconds. Left: surface and bottom. Right: zoom on the free surface.

Fig. 11

Left: initial relative density profiles. Right: difference of relative densities at $t=0$ and $t=500$ s

Fig. 12

Free surface and velocities at time $t=500$ seconds

Fig. 13

Velocity profiles at time $t=100$ s (left) and $t=200$ s (right)

Fig. 14

Initial condition for a smooth distribution of $\theta$

Fig. 15

Evolution of a smooth distribution of $\theta$ at time $t=5$ s

Fig. 16

Evolution of a smooth distribution of $\theta$ at time $t=10$ s

Fig. 17

Evolution of a smooth distribution of $\theta$ at time $t=20$ s

4.3 Smooth Distribution of Relative Density

In this example, we consider an initially smooth distribution of relative density given by

$$\begin{aligned} \theta (x) = 1 + \frac{1}{100} \, e^{-10\, x^2}, \ \text { for } x\in [-4,4]. \end{aligned}$$

(4.5)

The free surface is initially constant and equal to $\eta = 2$, while the bathymetry is continuous and equal to 0. The spatial domain $[-4,4]$ is divided into just 30 cells and we consider $M=10$ layers. The relative density is assumed to be the same across all layers. Figure 14 shows this initial density distribution.

Fig. 18

Evolution of a smooth distribution of $\theta$ at time $t=30$ s

Fig. 19

Evolution of a smooth distribution of $\theta$ at time $t=40$ s

We consider open boundary conditions and the fifth order numerical scheme will be used. In the sequence of Figs. 15, 16, 17, 18 and 19 we see the time evolution of the density profiles. In order to better describe the behavior of the relative density in all these figures, a dual representation has been chosen. On the left, the actual profile and values of the relative density is displayed for a few selected layers. This is done for the sake of clarity. On the right, the vertical relative density profile is shown through a heat map whiting the domain.

We observe that the density in each layer tends to move towards the boundaries at different speeds. At time $t=5$ we see that a shock is formed, even though the initial density profile is smooth. The appearance of this shock means that the limiter is then activated. To see this fact, we have included in the figures a parameter $\beta$ that accounts for the limiter activation. It is equal to 1 when it is not activated and greater than 1 otherwise.

Fig. 20

Lock-exchange experiment in a flat channel: initial condition

Fig. 21

Lock-exchange experiment in a flat channel: evolution of the relative density at time $t=2.5$ s

Fig. 22

Lock-exchange experiment in a flat channel: evolution of the relative density at time $t=5$ s

Fig. 23

Lock-exchange experiment in a flat channel: evolution of the relative density at time $t=10$ s

Fig. 24

Lock-exchange experiment in a flat channel: evolution of the relative density at time $t=15$ s

Fig. 25

Lock-exchange experiment in a flat channel: evolution of the relative density at time $t=25$ s

Fig. 26

Lock-exchange experiment in a flat channel: comparison on the evolution of the front position computed with the numerical scheme versus the laboratory data for 20 layers (left) and 25 layers (right)

Fig. 27

Lock-exchange experiment in a flat channel: comparison on the evolution of the front position computed with the numerical scheme versus the laboratory data for 30 layers (left) and 40 layers (right)

Fig. 28

Relative density dam break problem: initial condition

Fig. 29

Relative density dam break problem at time $t=2.5$ s for a first order numerical scheme (upper) and a fourth order numerical scheme (bottom)

Fig. 30

Relative density dam break problem at time $t=5$ s for a first order numerical scheme (upper) and a fourth order numerical scheme (bottom)

Fig. 31

Relative density dam break problem at time $t=10$ s for a first order numerical scheme (upper) and a fourth order numerical scheme (bottom)

4.4 Lock-Exchange Experiment in a Flat Channel

We now intend to validate the numerical approach introduced here by comparison with a laboratory test presented in [1], where experimental data is available. The experiment consists in a lock-exchange process between two fluids with different densities $\rho _0$ and $\rho _1$ in a flat channel. The channel is 3 m long and the height of the water at the initial time is 0.3 m. The fluid with density $\rho _1$ is in a gatebox of 0.1 m placed on the left of the channel, which is then released into the fluid of density $\rho _0$. This initial condition can be seen in Fig. 20. The density $\rho _0$ is 1000 kg/m$^3$ while $\rho _1$ is 1034 kg/m$^3$. Therefore we set

Fig. 32

Relative density dam break problem at time $t=40$ s for a first order numerical scheme (upper) and a fourth order numerical scheme (bottom)

Fig. 33

Relative density dam break problem at time $t=60$ s for a first order numerical scheme (upper) and a fourth order numerical scheme (bottom)

$$\begin{aligned} \theta (x) = {\left\{ \begin{array}{ll} 1.034 \qquad \text {if}\ x\le 0.1,\\ 1.0 \qquad \quad \text {if}\ x>0.1. \end{array}\right. } \end{aligned}$$

(4.6)

A mesh of 80 uniform cells are used and a fourth order in space and time numerical scheme is considered. To mimic the laboratory experiment in [1], wall-type boundary conditions are considered.

The evolution of the density profile is shown in Figs. 20, 21, 22, 23, 24 and 25 for the particular case of $M=40$ layers. As in the previous case, a dual representation is introduced, with the actual profile and values of the relative density for a few selected layers on the left, and a heat map of the vertical density distribution on the right. Again, $\beta$ stands for the activation of the limiter.

The immediate formation of a shock can be observed, due to the discontinuous nature of the initial condition. We see that the limiter is being activated throughout the jump in density and it tracks the traveling wave corresponding to the perturbation of the free surface due to the density change. Some high frequency oscillations associated with the high order can also be appreciated, but the additional viscosity provided by the limiter manages to smooth the solution.

Figures 26 and 27 show the evolution of the front position over time. This is compared with the experimental results presented in [1]. We observe a progressive convergence to the empirical data as the number of layers is increased from $M=20, 25, 30,$ and $M=40$, where the convergence is excellent.

The agreement with experimental data is excellent, which is noteworthy since non-hydrostatic effects present in situations such as this (see for example [46]) generally prevent an accurate prediction of the front position of hydrostatic shallow water models. Let us also remark that this is achieved with a coarse discretization used in both the horizontal and vertical directions.

4.5 Lock Exchange on a Non-Flat Bathymetry

We now consider a dam break in relative density problem with a continuous bathymetry function in a channel on the interval $[-5,5]$. Initially, a constant free surface $\eta =2$ m is set. The bathymetry function is given by

$$\begin{aligned} z_b(x) = \frac{1}{2}\, e^{-x^2}, \end{aligned}$$

(4.7)

and the following discontinuous profile in relative density is imposed

$$\begin{aligned} \theta (x) = {\left\{ \begin{array}{ll} 1 \qquad &{}\text {if} \ x\le 0, \\ 1.15 &{}\text {if} \ x > 0. \end{array}\right. } \end{aligned}$$

(4.8)

Figure 28 shows the initial condition. The channel is discretized with $M=25$ layers on the vertical direction and only 50 cells on the horizontal direction. Free-flow boundary conditions are considered.

We show the simulations at different times obtained with the first and fourth order schemes in Figs. 29, 30, 31, 32 and 33. The enhanced resolution of the high order scheme allows to appreciate the sharp transition between fluids, that are expected to evolve in gravity currents over obstacles (see [61]). The first order scheme is unable to accurately describe this solution, with this coarse discretization, and ends up with nearly uniform density solution due to the numerical viscosity. Conversely, the fourth order scheme accurately describes the experiment and the typical profile in such physical situations is obtained.

5 Conclusions

We have presented a novel discretization for a multi-layer shallow water model with a density dependent pressure function. The proposed numerical scheme falls into the category of ADER discontinuous Galerkin finite element schemes with a posteriori subcell finite volume limiter. It is arbitrarily high order accurate in space and time, allowing to accurately capture detailed solutions, even with coarse or very coarse meshes. Moreover, a suitable well-balancing strategy has been presented that allows to preserve a set of non-trivial stationary solutions corresponding to stationary density stratified lake-at-rest type solutions. The numerical results are promising, with excellent agreement with experimental data. Moreover, the computational results satisfy a discrete maximum principle for the relative density, thanks to its finite volume based a posteriori subcell limiter technique.