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.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
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,
Instead of considering total density, we may introduce the relative density
resulting in the following Euler equations,
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 \),
with
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.
Using this procedure, the following model is obtained from (2.3) (see [49, 64]):
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
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,
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,
Therefore, we are able to recover the particular case of the lake-at-rest stationary solutions with constant density,
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,
Equation (2.8) can be solved recursively by solving first the upper layer,
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
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,
with
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.
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,
where
and
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,
where \(\varvec{w}\) is the vector of the conserved variables,
The physical convective flux \(F_C(\varvec{w})\) is defined by,
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,
where
Finally, the term \(\varvec{T}(\varvec{w},\partial _x \varvec{w})\) corresponds to the mass, density, and momentum exchange between layers:
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\),
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\),
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,
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
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
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,
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:
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,
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,
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:
In particular, for a variable \(f_\alpha \) defined within the layer \(\alpha \), we shall write
The difference at the intercell \(i+\frac{1}{2}\) will be written as
We shall denote the average of a variable \(f_\alpha \) at the layer interface \(\alpha +\frac{1}{2}\) by
In addition, at the bottom or free surface interfaces, we shall assume
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
and
where \((f)_+\) is the positive part of f. Using (3.21), we define
where
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,
Once, the reconstructed states \({\mathbf {q}}_{i+\frac{1}{2}}^{HR,\pm }\) are defined, the Riemann solver can be expressed as,
and
where
with
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:
where
where, \(\widetilde{G_{\frac{1}{2}}}=\widetilde{G_{M+\frac{1}{2}}}=0\) and
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):
where
with
and
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,
where
with
where
and
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
where,
and
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
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,
where
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,
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}\),
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,
Since a subcell resolution \(K_{s}>N+1\) is admitted, this problem may be overdetermined. Hence, the following constraint must also be considered,
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,
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\)
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,
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,
Now, we can replace the numerical scheme (3.10) by the following equivalent well-balanced ADER-DG scheme,
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,
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,
Similarly,
where
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,
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}]\),
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,
if \(\left| \widehat{\varvec{w}}_{h,i}\right| <\delta _\epsilon \), then \(\varvec{w}_h(x,t^n)\) is redefined as follows,
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,
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:
and the bathymetry is given by,
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.
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,
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.
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,
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).
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,
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.
4.3 Smooth Distribution of Relative Density
In this example, we consider an initially smooth distribution of relative density given by
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.
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.
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
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
and the following discontinuous profile in relative density is imposed
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.
Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.
References
Adduce, C., Sciortino, G., Proietti, S.: Gravity currents produced by lock exchanges: experiments and simulations with a two-layer shallow-water model with entrainment. J. Hydraul. Eng. 138(2), 111–121 (2012). https://doi.org/10.1061/(ASCE)HY.1943-7900.0000484
Audusse, E., Bouchut, F., Bristeau, M.O., Klein, R., Perthame, B.: A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM J. Sci. Comput. 25(6), 2050–2065 (2004). https://doi.org/10.1137/S1064827503431090
Audusse, E., Bristeau, M.O.: A well-balanced positivity preserving “second-order’’ scheme for shallow water flows on unstructured meshes. J. Comput. Phys. 206(1), 311–333 (2005)
Audusse, E., Bristeau, M.O.: Finite-volume solvers for a multilayer saint-venant system. Appl. Math. Comput. Sci. 17, 311–320 (2007). https://doi.org/10.2478/v10006-007-0025-0
Audusse, E., Bristeau, M.O., Pelanti, M., Sainte-Marie, J.: Approximation of the hydrostatic Navier–Stokes system for density stratified flows by a multilayer model: kinetic interpretation and numerical solution. J. Comput. Phys. 230(9), 3453–3478 (2011)
Audusse, E., Bristeau, M.O., Perthame, B., Sainte-Marie, J.: A multilayer saint-venant system with mass exchanges for shallow water flows. Derivation and numerical validation. ESAIM Math. Model. Numer. Anal. 45(1), 169–200 (2011)
Bassi, C., Bonaventura, L., Busto, S., Dumbser, M.: A hyperbolic reformulation of the Serre–Green–Naghdi model for general bottom topographies. Comput. Fluids 212, 104716 (2020)
Bermúdez, A., Vázquez, M.: Upwind methods for hyperbolic conservation laws with source terms. Comput. Fluids 23(8), 1049–1071 (1994)
Bonaventura, L., Fernández-Nieto, E.D., Garres-Díaz, J., Narbona-Reina, G.: Multilayer shallow water models with locally variable number of layers and semi-implicit time discretization. J. Comput. Phys. 364, 209–234 (2018) https://doi.org/10.1016/j.jcp.2018.03.017. http://www.sciencedirect.com/science/article/pii/S0021999118301694
Bouchut, F., Zeitlin, V.: A robust well-balanced scheme for multi-layer shallow water equations. Discrete Continu. Dyn. Syst. Ser. B 13, 739–758 (2010). https://doi.org/10.3934/dcdsb.2010.13.739
Bristeau, M.O., Mangeney, A., Sainte-Marie, J., Seguin, N.: An energy-consistent depth-averaged euler system: derivation and properties. arXiv preprint arXiv:1406.6565 (2014)
Bürger, R., Fernández-Nieto, D., Andrés Osores, E.V.: A dynamic multilayer shallow water model for polydisperse sedimentation. ESAIM Math. Modell. Numer. Anal. (2019). https://doi.org/10.1051/m2an/2019032
Busto, S., Chiocchetti, S., Dumbser, M., Gaburro, E., Peshkov, I.: High order ADER schemes for continuum mechanics. Front. Phys. 8, 32 (2020)
Busto, S., Dumbser, M., Escalante, C., Gavrilyuk, S., Favrie, N.: On high order ADER discontinuous Galerkin schemes for first order hyperbolic reformulations of nonlinear dispersive systems. J. Sci. Comput. 87, 48 (2021)
Busto, S., Tavelli, M., Boscheri, W., Dumbser, M.: Efficient high order accurate staggered semi-implicit discontinuous Galerkin methods for natural convection problems. Comput. Fluids 198, 104399 (2020)
Busto, S., Toro, E.F., Vázquez-Cendón, M.E.: Design and analysis of ader-type schemes for model advection-diffusion-reaction equations. J. Comput. Phys. 327, 553–575 (2016)
Castro, M., Fernández-Nieto, E.: A class of computationally fast first order finite volume solvers: PVM methods. SIAM J. Sci. Comput. (2012). https://doi.org/10.1137/100795280
Castro, M., Gallardo, J.M., Parés, C.: High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. applications to shallow-water systems. Math. Comput. 75(255), 1103–1134 (2006)
Castro, M., Macías, J., Parés, C.: A q-scheme for a class of systems of coupled conservation laws with source term. Application to a two-layer 1-D shallow water system. ESAIM Math. Modell. Numer. Anal. 35(1), 107–127 (2001)
Castro, M., Pardo, A., Parés, C.: Well-balanced schemes based on a generalized hydrostatic reconstruction technique. Math. Models Methods Appl. Sci. 17(12), 2055–2113 (2007)
Castro, M.J., Parés, C.: Well-balanced high-order finite volume methods for systems of balance laws. J. Sci. Comput. 82(2), 48 (2020)
Casulli, V.: Semi-implicit finite difference methods for the two-dimensional shallow water equations. J. Comput. Phys. 86, 56–74 (1990)
Casulli, V.: A semi-implicit finite difference method for non-hydrostatic, free-surface flows. Int. J. Numer. Meth. Fluids 30(4), 425–440 (1999)
Casulli, V.: A high-resolution wetting and drying algorithm for free-surface hydrodynamics. Int. J. Numer. Meth. Fluids 60, 391–408 (2009)
Casulli, V.: A semi-implicit numerical method for the free-surface Navier–Stokes equations. Int. J. Numer. Meth. Fluids 74, 605–622 (2014)
Casulli, V., Cheng, R.: Semi-implicit finite difference methods for three-dimensional shallow water flow. Int. J. Numer. Methods Fluids 15, 629–648 (1992)
Charrier, D., Weinzierl, T.: Stop talking to me–a communication-avoiding ADER-DG realisation. SIAM J. Sci. Comput. (2018). Submitted to. arXiv:1801.08682
Cheng, Y., Dong, H., Li, M., Xian, W.: A high order central DG method of the two-layer shallow water equations. Commun. Comput. Phys. 28(4), 1437–1463 (2020). https://doi.org/10.4208/cicp.oa-2019-0155
Clain, S., Diot, S., Loubère, R.: A high-order finite volume method for systems of conservation laws-multi-dimensional optimal order detection (MOOD). J. Comput. Phys. 230(10), 4028–4050 (2011)
Cockburn, B., Hou, S., Shu, C.W.: The Runge–Kutta local projection discontinuous galerkin finite element method for conservation laws. IV: The multidimensional case. Math. Comput. 54(190), 545–581 (1990). http://www.jstor.org/stable/2008501
Cockburn, B., Lin, S.Y., Shu, C.W.: TVB Runge–Kutta local projection discontinuous galerkin finite element method for conservation laws iii: One-dimensional systems. J. Comput. Phys. 84(1), 90–113 (1989) https://doi.org/10.1016/0021-9991(89)90183-6. http://www.sciencedirect.com/science/article/pii/0021999189901836
Cockburn, B., Shu, C.: The Runge–Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems. J. Comput. Phys. 141(2), 199–224 (1998)
Cockburn, B., Shu, C.W.: Tvb Runge–Kutta local projection discontinuous galerkin finite element method for conservation laws ii: General framework. Math. Comput. 52(186), 411–435 (1989). http://www.jstor.org/stable/2008474
Cockburn, Bernardo, Shu, Chi-Wang.: The runge-kutta local projection \(p^1\)-discontinuous-galerkin finite element method for scalar conservation laws. ESAIM: M2AN 25(3), 337–361 (1991). https://doi.org/10.1051/m2an/1991250303371
De St. Venant, B.: Theorie du mouvement non-permanent des eaux avec application aux crues des rivers et a l’introduntion des marees dans leur lit. Acad. de Sci. Comptes Redus 73(99), 148–154 (1871)
Diot, S., Clain, S., Loubère, R.: Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials. Comput. Fluids 64, 43–63 (2012)
Dumbser, M., Balsara, D.S., Toro, E.F., Munz, C.D.: A unified framework for the construction of one-step finite volume and discontinuous galerkin schemes on unstructured meshes. J. Comput. Phys. 227(18), 8209–8253 (2008) https://doi.org/10.1016/j.jcp.2008.05.025. http://www.sciencedirect.com/science/article/pii/S0021999108002829
Dumbser, M., Castro, M., Parés, C., Toro, E.: ADER schemes on unstructured meshes for non-conservative hyperbolic systems: Applications to geophysical flows. Comput. Fluids 38, 1731–1748 (2009)
Dumbser, M., Castro, M., Parés, C., Toro, E.F.: ADER schemes on unstructured meshes for nonconservative hyperbolic systems: applications to geophysical flows. Comput. Fluids 38(9), 1731–1748 (2009). https://doi.org/10.1016/j.compfluid.2009.03.008
Dumbser, M., Casulli, V.: A staggered semi-implicit spectral discontinuous Galerkin scheme for the shallow water equations. Appl. Math. Comput. 219, 8057–8077 (2013)
Dumbser, M., Enaux, C., Toro, E.F.: Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws. J. Comput. Phys. 227(8), 3971–4001 (2008) https://doi.org/10.1016/j.jcp.2007.12.005. http://www.sciencedirect.com/science/article/pii/S0021999107005578
Dumbser, M., Facchini, M.: A local space-time discontinuous Galerkin method for Boussinesq-type equations. Appl. Math. Comput. 272, 336–346 (2016)
Dumbser, M., Hidalgo, A., Castro, M., Parés, C., Toro, E.: FORCE schemes on unstructured meshes II: Non-conservative hyperbolic systems. Comput. Methods Appl. Mech. Eng. 199, 625–647 (2010)
Dumbser, M., Zanotti, O., Loubère, R., Diot, S.: A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws. J. Comput. Phys. 278, 47–75 (2014). https://doi.org/10.1016/j.jcp.2014.08.009
Escalante, C., Dumbser, M., Castro, M.: An efficient hyperbolic relaxation system for dispersive non-hydrostatic water waves and its solution with high order discontinuous Galerkin schemes. J. Comput. Phys. 394, 385–416 (2019)
Escalante, C., de Luna, T.M., Castro, M.: Non-hydrostatic pressure shallow flows: Gpu implementation using finite volume and finite difference scheme. Appl. Math. Comput. 338, 631–659 (2018)
Fambri, F., Dumbser, M., Köppel, S., Rezzolla, L., Zanotti, O.: ADER discontinuous Galerkin schemes for general-relativistic ideal magnetohydrodynamics. Mon. Not. R. Astron. Soc. 477, 4543–4564 (2018)
Favrie, N., Gavrilyuk, S.: A rapid numerical method for solving Serre–Green–Naghdi equations describing long free surface gravity waves. Nonlinearity 30, 2718–2736 (2017)
Fernández-Nieto, E.H., Koné, E., Chacón-Rebollo, T.: A multilayer method for the hydrostatic navier-stokes equations: a particular weak solution. J. Sci. Comput. (2014). https://doi.org/10.1007/s10915-013-9802-0
Fernández-Nieto, E.D., Parisot, M., Penel, Y., Sainte-Marie, J.: A hierarchy of dispersive layer-averaged approximations of Euler equations for free surface flows. Commun. Math. Sci. 16(5), 1169–1202 (2018). https://doi.org/10.4310/CMS.2018.v16.n5.a1
Garres-Díaz, J., Bonaventura, L.: Flexible and efficient discretizations of multilayer models with variable density. Appl. Math. Comput. 402, 126097 (2021) https://doi.org/10.1016/j.amc.2021.126097. https://www.sciencedirect.com/science/article/pii/S0096300321001454
Gassner, G., Dumbser, M., Hindenlang, F., Munz, C.: Explicit one-step time discretizations for discontinuous Galerkin and finite volume schemes based on local predictors. J. Comput. Phys. 230(11), 4232–4247 (2011)
Gosse, L.: A well-balanced scheme using non-conservative products designed for hyperbolic systems of conservation laws with source terms. Math. Models Methods Appl. Sci. 11(02), 339–365 (2001)
Guerrero Fernández, E., Castro-Díaz, M..J., Luna, T..M..d: A second-order well-balanced finite volume scheme for the multilayer shallow water model with variable density. Mathematics 8(5), 848 (2020)
Hidalgo, A., Dumbser, M.: Ader schemes for nonlinear systems of stiff advection–diffusion–reaction equations. J. Sci. Comput. 48(1–3), 173–189 (2011)
Higdon, R.L.: Discontinuous galerkin methods for multi-layer ocean modeling: Viscosity and thin layers. J. Comput. Phys. (2020). https://doi.org/10.1016/j.jcp.2019.109018
Izem, N., Seaid, M., Wakrim, M.: A discontinuous Galerkin method for two-layer shallow water equations. Math. Comput. Simul. 120, 12–23 (2016). https://doi.org/10.1016/j.matcom.2015.04.009
Jackson, H.: On the eigenvalues of the Ader–Weno galerkin predictor. J. Comput. Phys. 333, 409–413 (2017)
Klaij, C.M., van der Vegt, J.J., van der Ven, H.: Space-time discontinuous Galerkin method for the compressible Navier–Stokes equations. J. Comput. Phys. 217(2), 589–611 (2006)
Kopriva, D.A., Gassner, G.: On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods. J. Sci. Comput. 44(2), 136–155 (2010)
Lane-Serff, G.F., Beal, L.M., Hadfield, T.D.: Gravity current flow over obstacles. J. Fluid Mech. 292, 39–53 (1995). https://doi.org/10.1017/S002211209500142X
Li, G., Li, J., Qian, S., Gao, J.: A well-balanced ADER discontinuous Galerkin method based on differential transformation procedure for shallow water equations. Appl. Math. Comput. (2021). https://doi.org/10.1016/j.amc.2020.125848
Li, G., Song, L., Gao, J.: High order well-balanced discontinuous galerkin methods based on hydrostatic reconstruction for shallow water equations. J. Comput. Appl. Math. 340, 546–560 (2018). https://doi.org/10.1016/j.cam.2017.10.027
de Luna, T.M., Fernández Nieto, E., Castro Díaz, M.J.: Derivation of a multilayer approach to model suspended sediment transport: application to hyperpycnal and hypopycnal plumes. Commun. Comput. Phys. 22(5), 1439–1485 (2017). https://doi.org/10.4208/cicp.OA-2016-0215
Morales de Luna, T., Castro Díaz, M., Parés, C.: Reliability of first order numerical schemes for solving shallow water system over abrupt topography. Appl. Math. Comput. 219(17), 9012–9032 (2013). https://doi.org/10.1016/j.amc.2013.03.033. https://www.sciencedirect.com/science/article/pii/S0096300313002865
Owren, B., Zennaro, M.: Derivation of efficient, continuous, explicit Runge–Kutta methods. SIAM J. Sci. Stat. Comput. 13, 1488–1501 (1992)
Parés, C.: Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM J. Numer. Anal. 44(1), 300–321 (2006). https://doi.org/10.1137/050628052
Reed, W.H., Hill, T.: Triangular mesh methods for the neutron transport equation. Tech. rep., Los Alamos Scientific Lab., N. Mex.(USA) (1973)
Rhebergen, S., Bokhove, O., van der Vegt, J.: Discontinuous Galerkin finite element methods for hyperbolic nonconservative partial differential equations. J. Comput. Phys. 227, 1887–1922 (2008)
Tavelli, M., Dumbser, M.: A high order semi-implicit discontinuous galerkin method for the two dimensional shallow water equations on staggered unstructured meshes. Appl. Math. Comput. 234, 623–644 (2014)
Titarev, V., Toro, E.: Ader: Arbitrary high order godunov approach. Journal of Scientific Computing 17(1–4), 609–618 (2002) https://doi.org/10.1023/A:1015126814947. https://www.scopus.com/inward/record.uri?eid=2-s2.0-0008392886&doi=10.1023%2fA%3a1015126814947&partnerID=40&md5=8bcd619a3597540cbc2f03446c333d4a
Titarev, V., Toro, E.: Ader schemes for three-dimensional non-linear hyperbolic systems. J. Comput. Phys. 204(2), 715–736 (2005) https://doi.org/10.1016/j.jcp.2004.10.028. http://www.sciencedirect.com/science/article/pii/S0021999104004358
Toro, E., Titarev, V.: Solution of the generalized riemann problem for advection-reaction equations. Proc. R Soc. A Math. Phys. Eng. Sci. 458(2018), 271–281 (2002) https://doi.org/10.1098/rspa.2001.0926. https://www.scopus.com/inward/record.uri?eid=2-s2.0-57249099681&doi=10.1098%2frspa.2001.0926&partnerID=40&md5=86bc8d2d1e3a77fe704c61aa299a53bb. Cited By 159
Toro, E., Titarev, V.: Derivative riemann solvers for systems of conservation laws and ader methods. J. Comput. Phys. 212(1), 150–165 (2006) https://doi.org/10.1016/j.jcp.2005.06.018. http://www.sciencedirect.com/science/article/pii/S0021999105003141
Toro, E.F., Millington, R., Nejad, L.: Towards very high order godunov schemes. In: Godunov methods, pp. 907–940. Springer (2001)
Tumolo, G., Bonaventura, L.: Simulations of Non-hydrostatic Flows by an Efficient and Accurate p-Adaptive DG Method, pp. 41–53. Springer International Publishing, Cham (2020). https://doi.org/10.1007/978-3-030-30705-9_5
Tumolo, G., Bonaventura, L., Restelli, M.: A semi-implicit, semi-Lagrangian, p-adaptive discontinuous Galerkin method for the shallow water equations. J. Comput. Phys. 232, 46–67 (2013)
van der Vegt, J., van der Ven, H.: Space-time discontinuous galerkin finite element method with dynamic grid motion for inviscid compressible flows: I. general formulation. J. Comput. Phys. 182(2), 546–585 (2002). https://doi.org/10.1006/jcph.2002.7185. http://www.sciencedirect.com/science/article/pii/S0021999102971858
van der Ven, H., van der Vegt, J.: Space-time discontinuous galerkin finite element method with dynamic grid motion for inviscid compressible flows: Ii. efficient flux quadrature. Comput. Methods Appl. Mech. Eng. 191(41), 4747–4780 (2002). https://doi.org/10.1016/S0045-7825(02)00403-6. http://www.sciencedirect.com/science/article/pii/S0045782502004036
Wang, Z., Liu, Y.: Extension of the spectral volume method to high-order boundary representation. J. Comput. Phys. 211(1), 154–178 (2006) https://doi.org/10.1016/j.jcp.2005.05.022. http://www.sciencedirect.com/science/article/pii/S0021999105002664
Wintermeyer, N., Winters, A.R., Gassner, G.J., Warburton, T.: An entropy stable discontinuous galerkin method for the shallow water equations on curvilinear meshes with wet/dry fronts accelerated by GPUs. J. Comput. Phys. 375, 447–480 (2018). https://doi.org/10.1016/j.jcp.2018.08.038
Wu, X., Kubatko, E.J., Chan, J.: High-order entropy stable discontinuous galerkin methods for the shallow water equations: curved triangular meshes and GPU acceleration. Comput. Math. Appl. 82, 179–199 (2021). https://doi.org/10.1016/j.camwa.2020.11.006
Zanotti, O., Dumbser, M.: Efficient conservative ader schemes based on weno reconstruction and space-time predictor in primitive variables. Comput. Astrophys. Cosmol. 3(1), 1 (2016)
Acknowledgements
This research has been partially supported by the Spanish Government and FEDER through the coordinated Research project RTI2018-096064-B-C1 and RTI2018-096064-B-C2, The Junta de Andalucía research project P18-RT-3163 and the Junta de Andalucia-FEDER-University of Málaga Research project UMA18-FEDERJA-161. The University of Malaga has also supported this research through the “Contrato Puente”.
Funding
Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors declare no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Fern�ndez, E.G., D�az, M.J.C., Dumbser, M. et al. An Arbitrary High Order Well-Balanced ADER-DG Numerical Scheme for the Multilayer Shallow-Water Model with Variable Density. J Sci Comput 90, 52 (2022). https://doi.org/10.1007/s10915-021-01734-2
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s10915-021-01734-2