An all-regime and well-balanced Lagrange-projection type scheme for the shallow water equations on unstructured meshes

TL;DR

This paper introduces a novel numerical scheme for the 2D shallow water equations that is well-balanced, capable of handling all flow regimes including low-Froude, and maintains positive water height on unstructured meshes.

Contribution

It proposes a Lagrange-projection based scheme that ensures well-balancedness, all-regime applicability, and positivity, improving robustness over traditional methods.

Findings

Successfully preserves lake at rest equilibrium solutions.

Effective in all flow regimes, including low-Froude.

Maintains positive water height throughout simulations.

Abstract

In this work, we focus on the numerical approximation of the shallow water equations in two space dimensions. Our aim is to propose a well-balanced, all-regime and positive scheme. By well-balanced, it is meant that the scheme is able to preserve the so-called lake at rest smooth equilibrium solutions. By all-regime, we mean that the scheme is able to deal with all flow regimes, including the low-Froude regime which is known to be challenging when using usual Godunov-type finite volume schemes. At last, the scheme should be positive which means that the water height stays positive for all time. Our approach is based on a Lagrange-projection decomposition which allows to naturally decouple the acoustic and transport terms. Numerical experiments on unstructured meshes illustrate the good behaviour of the scheme.

Figures

Tables

Table 1. Table 1: Well-balanced property. Errors for EXEX and IMEX schemes.

	${\parallel \rho -0.5\parallel }_{\infty }$	${\parallel 𝐮\parallel }_{\infty }$
EXEX	${\mathrm{2.6\hspace{0.17em}10}}^{-16}$	${\mathrm{1.3\hspace{0.17em}10}}^{-13}$
IMEX	${\mathrm{2.6\hspace{0.17em}10}}^{-16}$	${\mathrm{3.9\hspace{0.17em}10}}^{-8}$

Table 2. Table 2: Traveling vortex test case with flat bottom. Numbers of iterations and CPU times with low-Froude correction.

	Number of time steps	CPU time
EXEX	60264	1930
IMEX	689	175

Table 3. Table 3: Traveling vortex test case with non-flat bottom. Numbers of iterations and CPU times with low-Froude correction.

	Nb time steps	CPU time
EXEX	60264	15748
IMEX	2733	921

Equations

$${\partial }_{t}h+\nabla \cdot (h\mathbf{u})=0,$$

$${\partial }_t(h\mathbf{u})+\nabla \cdot (h\mathbf{u}\otimes \mathbf{u})+\nabla \frac{g{h}^2}{2}=-gh\nabla z,$$

$$h+z=\mathrm{constant},\mathbf{u}=0.$$

$$\tilde{t}=\frac{t}{T},$$

$$u_0=\frac{L}{T}\mathrm{and}{z}_0=h_0.$$

$$\mathrm{Fr}=\frac{u_0}{c_0},$$

$${\partial }_{\tilde{t}}\tilde{h}+{\nabla }_{\tilde{\mathbf{x}}}\cdot \left(\tilde{h}\tilde{\mathbf{u}}\right)=0,$$

$${\partial }_{\tilde{t}}\left(\tilde{h}\tilde{\mathbf{u}}\right)+{\nabla }_{\tilde{\mathbf{x}}}\cdot \left(\tilde{h}\tilde{\mathbf{u}}\otimes \tilde{\mathbf{u}}\right)+\frac{1}{{\mathrm{Fr}}^2}{\nabla }_{\tilde{\mathbf{x}}}\frac{{\tilde{h}}^2}{2}=-\frac{1}{{\mathrm{Fr}}^2}\tilde{h}{\nabla }_{\tilde{\mathbf{x}}}\tilde{z}.$$

$${\partial }_{t}h+\nabla \cdot \left(h\mathbf{u}\right)=0,$$

$${\partial }_t\left(h\mathbf{u}\right)+\nabla \cdot \left(h\mathbf{u}\otimes \mathbf{u}\right)+\frac{1}{{\mathrm{Fr}}^2}\nabla p=-\frac{1}{{\mathrm{Fr}}^2}h\nabla z.$$

$$h=h^{(0)}+h^{(1)}\mathrm{Fr}+h^{(2)}{\mathrm{Fr}}^2+\mathcal{O}({\mathrm{Fr}}^3)$$

$$p=p^{(0)}+p^{(1)}\mathrm{Fr}+p^{(2)}{\mathrm{Fr}}^2+\mathcal{O}({\mathrm{Fr}}^3)=p(h^{(0)})+h^{(1)}{p}^{\prime }(h^{(0)})\mathrm{Fr}+\mathcal{O}({\mathrm{Fr}}^2).$$

$$\nabla p^{(0)}+h^{(0)}\nabla z=0$$

$$\nabla p^{(1)}+h^{(1)}\nabla z=0$$

$${\partial }_t{h}^{(0)}+\nabla \cdot (h^{(0)}{\mathbf{u}}^{(0)})$$

$${\partial }_t(h^{(0)}{\mathbf{u}}^{(0)})+\nabla \cdot (h^{(0)}{\mathbf{u}}^{(0)}\otimes {\mathbf{u}}^{(0)})+\nabla p^{(2)}$$

$$\left({\int }_{\Omega }\nabla \cdot \mathbf{u}\mathrm{d}\Omega =0\text{ and }{\int }_{\Omega }\nabla \cdot (z\mathbf{u})\mathrm{d}\Omega =0\right)$$

$$0$$

$$={\int }_{\Omega }{\partial }_t(H-z)\mathrm{d}\Omega +{\int }_{\Omega }\nabla \cdot (\left(H-z\right){\mathbf{u}}^{(0)})\mathrm{d}\Omega$$

$$={\int }_{\Omega }{\partial }_{t}H\mathrm{d}\Omega +H{\int }_{\Omega }\nabla \cdot {\mathbf{u}}^{(0)}\mathrm{d}\Omega -{\int }_{\Omega }\nabla \cdot (z{\mathbf{u}}^{(0)})\mathrm{d}\Omega$$

$$={\partial }_t{\int }_{\Omega }H\mathrm{d}\Omega =|\Omega |{\partial }_t(h^{(0)}+z)=|\Omega |{\partial }_t{h}^{(0)}$$

$$\nabla \cdot {\mathbf{u}}^{(0)}=\nabla \cdot (\frac{z}{H}{\mathbf{u}}^{(0)}),$$

$$\left(1-\frac{z}{H}\right){\partial }_t{\mathbf{u}}^{(0)}+\nabla \cdot ({\mathbf{u}}^{(0)}\otimes {\mathbf{u}}^{(0)})+\frac{1}{H}\nabla p^{(2)}=\nabla \cdot (\frac{z}{H}{\mathbf{u}}^{(0)}\otimes {\mathbf{u}}^{(0)})-h^{(2)}\nabla \frac{z}{H}.$$

$$\{\begin{array}{rl}& h^{(0)}+z=H=\text{cst}\\ & \nabla \cdot {\mathbf{u}}^{(0)}=\nabla \cdot (\frac{z}{H}{\mathbf{u}}^{(0)})\\ & \left(1-\frac{z}{H}\right){\partial }_t{\mathbf{u}}^{(0)}+\nabla \cdot ({\mathbf{u}}^{(0)}\otimes {\mathbf{u}}^{(0)})+\frac{1}{H}\nabla p^{(2)}=\nabla \cdot (\frac{z}{H}{\mathbf{u}}^{(0)}\otimes {\mathbf{u}}^{(0)})-h^{(2)}\nabla \frac{z}{H}\end{array}$$

$$\{\begin{array}{rl}& h^{(0)}=\mathrm{cste}\\ & \nabla \cdot {\mathbf{u}}^{(0)}=0\\ & {\partial }_t{\mathbf{u}}^{(0)}+\nabla \cdot ({\mathbf{u}}^{(0)}\otimes {\mathbf{u}}^{(0)})+\frac{1}{h^{(0)}}\nabla p^{(2)}=0.\end{array}$$

$${\partial }_{t}h+h\nabla \cdot (\mathbf{u})=0,{\partial }_t(h\mathbf{u})+h\mathbf{u}(\nabla \cdot \mathbf{u})+\nabla p=-gh\nabla z,$$

$${\partial }_{t}h+(\mathbf{u}\cdot \nabla )h=0,{\partial }_t(h\mathbf{u})+(\mathbf{u}\cdot \nabla )(h\mathbf{u})=0.$$

$${\partial }_t\tau -\tau (\mathbf{x},t)\nabla \cdot \mathbf{u}=0,{\partial }_t\mathbf{u}+\tau (\mathbf{x},t)\nabla p=-\tau (\mathbf{x},t)\frac{g}{\tau }\nabla z.$$

$$\{\begin{array}{rlrl}{\partial }_t\tau & -\tau (\mathbf{x},t)\nabla \cdot \mathbf{u}& & =0,\\ {\partial }_t\mathbf{u}& +\tau (\mathbf{x},t)\nabla \Pi & & =-\tau (\mathbf{x},t)\frac{g}{\tau }\nabla z,\\ {\partial }_t\Pi & +\tau (\mathbf{x},t)a^2\nabla \cdot \mathbf{u}& & =\lambda (p^{\mathrm{EOS}}(\tau )-\Pi ),\end{array}$$

$$\{\begin{array}{rlrl}{\partial }_t\tau & -\tau (\mathbf{x},t^n)\nabla \cdot \mathbf{u}& & =0,\\ {\partial }_t\mathbf{u}& +\tau (\mathbf{x},t^n)\nabla \Pi & & =-\tau (\mathbf{x},t^n)\frac{g}{\tau }\nabla z,\\ {\partial }_t\Pi & +\tau (\mathbf{x},t^n)a^2\nabla \cdot \mathbf{u}& & =0,\end{array}$$

$$\{\begin{array}{rl}& {\partial }_{t}h+{\partial }_x(h{u}_1)=0,\\ & {\partial }_t(h{u}_1)+{\partial }_x\left(h{u}_1^2+g\frac{h^2}{2}\right)=-gh{\partial }_{x}z,\\ & {\partial }_t(h{u}_2)+{\partial }_x(h{u}_1{u}_2)=0.\end{array}$$

$$\{\begin{array}{rl}& {\partial }_t\tau -\tau (x,t^n){\partial }_x{u}_1=0,\\ & {\partial }_t{u}_1+\tau (x,t^n){\partial }_x\Pi =-\tau (x,t^n)\frac{g}{\tau }{\partial }_{x}z,\\ & {\partial }_t{u}_2=0,\\ & {\partial }_t\Pi +\tau (x,t^n)a^2{\partial }_x{u}_1=0,\end{array}$$

$${\partial }_t\phi +u_1{\partial }_x\phi =0,\phi \in \{h,u_1,u_2\}.$$

$${\tau }_j^{n+1-}={\tau }_j^n-{\tau }_j^n\frac{\Delta t}{\Delta x_j}\left(u_{j+1/2}^{♯}-u_{j-1/2}^{♯}\right)$$

Full text

1 Introduction

We are interested in the numerical approximation of the shallow water equations (SWE)

[TABLE]

where $\mathbf{x}\in\mathbb{R}^{2}\mapsto z(\mathbf{x})$ denotes a given smooth topography and $g>0$ is the gravity constant. Both the water depth $h$ and the velocity $\mathbf{u}=(u_{1},u_{2})\in\mathbb{R}^{2}$ depend on the space and time variables, namely $\mathbf{x}\in\mathbb{R}^{2}$ and $t\in[0,\infty)$. We assume that the initial water depth $h(\mathbf{x},t=0)=h_{0}(\mathbf{x})$ and velocity $\mathbf{u}(\mathbf{x},t=0)=\mathbf{u}_{0}(\mathbf{x})$ are given.

Let us briefly properties of system (1) in the case $\nabla z=0$ : the system is strictly hyperbolic over the phase space $\Omega=\{(h,h\mathbf{u})\in\mathbb{R}^{3}\leavevmode\nobreak\ |\leavevmode\nobreak\ h>0\}$. Moreover, if $\mathbf{n}\in\mathbb{R}^{2}$ is an arbitrary unit vector, the eigenstructure of (1) is composed by two genuinely nonlinear characteristic fields associated with the eigenvalues $\{\mathbf{u}^{T}\mathbf{n}-c,\mathbf{u}^{T}\mathbf{n}+c\}$, where $c\coloneqq\sqrt{gh}$ is the sound speed, and a linearly degenerated field associated with the eigenvalue $\mathbf{u}^{T}\mathbf{n}$. We recall also that the regions where $(\mathbf{u}^{T}\mathbf{n})^{2}<c^{2}$ (resp. $(\mathbf{u}^{T}\mathbf{n})^{2}>c^{2}$) are called subcritical or subsonic (resp. supercritical or supersonic).

We are interested in this work in developing a numerical scheme that satisfies the well-balanced property. More specifically we want our scheme to strictly preserve the "lake at rest" steady solutions, that are the states satisfying

[TABLE]

For a review on numerical schemes that satisfy the so-called well-balanced property we refer the reader to the pioneering work [BV94], books [Bou04] and [Gos13]. We also refer to [CKKS17] where the authors focus on the 1D case and already propose a well-balanced Lagrange-projection strategy. At last, in [CDCdL18] the proposed Lagrange-projection scheme is exact for a full set of equilibrium solutions (and not only the lake at rest).

The Lagrange-projection methodology is especially well suited for subsonic or near low-Froude number flows. We use an implicit-explicit strategy that allows to keep a stable scheme under a CFL time step limitation which is driven only by (slow) material waves and not by (fast) acoustic waves. The implicit-explicit Lagrange-projection [GR96] scheme is designed following the pioneering work [CNPT10]. More recent works are concernet with the case of Euler systems in the large friction or low-Mach regimes [CGK13, CGK14, CGK16] for single or two-phase flow models. The treatment of the low Froude number is considered through the so called all-regime (or asymptotic-preserving) property and follows the anti-diffusive technique on the pressure numerical flux introduced in [Del10] and also used in [CGK16].

The SWE has been largely studied and one can find nice overviews and references in the books [Bou04] and [Gos13]. The scheme proposed in [CKKS17] in one dimension has been studied in the framework of SWE and more specifically its behaviour for low-Froude number flows in [Zak17]. A different implicit-explicit methodology in two dimension context has been proposed by [BALMN14].

In section 2, we study the dimensionless system associated to the SWE (1) and its asymptotic limit in low Froude regimes. In section 3, we present the Lagrange-projection like acoustic / transport decomposition associated to system (1). In section 4, we present the schemes, the finite volume scheme in 1D, the study of its truncation error in low Froude regimes and the proposed correction, and finally the extension towards 2D schemes on unstructured meshes. At last, we show some numerical results in 2D to verify the well-balanced property and illustrate the behaviour of the scheme in different regimes, especially in the low Froude one.

2 Low Froude limit for continuous equations

2.1 Dimensionless shallow water equations

In this section, we briefly introduce the dimensionless SWE. These equations will be useful to study the low-Froude asymptotic behaviour of the solutions of (1). With this in mind, we define the following dimensionless quantities

[TABLE]

where $T$, $L$, $h_{0}$, $u_{0}$ and $z_{0}$ are respectively reference time, length, water height, velocity and topography such that

[TABLE]

Defining the Froude number $\mathrm{Fr}$ by

[TABLE]

where $c_{0}=\sqrt{gh_{0}}$ is the reference sound speed, easy calculations then give the dimensionless SWE

[TABLE]

2.2 Asymptotic equations in low Froude limit

In this section, we give the asymptotic behaviour of the solutions of the SWE equations in the low Froude limit. If we omit the tilde notation for the sake of readability in system (2) and, if we introduce the dimensionless pressure function $p(h)=\frac{h^{2}}{2}$, we get

[TABLE]

Let us assume that $h$ and $z$ admit the following expansions in powers of the Froude number:

[TABLE]

which gives in particular

[TABLE]

The governing equations give at order -2 and -1 with respect to the Froude number that

[TABLE]

The asymptotic behavior is then given by

[TABLE]

Now if we impose one of the following velocity boundary conditions

[TABLE]

integrating (4a) with respect to the space variable gives

[TABLE]

thus $\partial_{t}h^{(0)}=0$ and $h^{(0)}+z=H$ is constant both in space and time. This leads to $\nabla\cdot(h^{(0)}\mathbf{u}^{(0)})=0$ and therefore

[TABLE]

while the evolution of $\mathbf{u}$ is given by

[TABLE]

Notice that when the topography is flat, i.e. $z=0$, the three equations

[TABLE]

degenerate towards the incompressible Euler equations

[TABLE]

3 An acoustic/transport operator decomposition

Let us first introduce notations related to our discretization. We suppose that the computational domain $\Omega\subset\mathbb{R}^{2}$ is covered by $N$ polygonal cells $\left(\Omega_{j}\right)_{1\leq j\leq N}$. We consider $\Gamma$, a face of the cell $j$, and we suppose the following admissibility assumptions are satisfied:

• •

either there exists a single $1\leq k\leq N$ such that $\Gamma=\overline{\Omega_{j}}\cap\overline{\Omega_{k}}\not=\emptyset$. In this case we note $\Gamma=\Gamma_{jk}$ and $\Gamma_{jk}$ can either be a vertex or a single face of the mesh,

• •

either $\Gamma\subset\partial\Omega$ and we suppose that there exists a single $k>N$ that will help to index ghost values for boundary conditions and we shall note $\Gamma=\Gamma_{jk}$.

For $1\leq j\leq N$, we note $\mathcal{N}(j)$ the set of indices $k$ such that $\Gamma_{jk}$ is a face of $\Omega_{j}$ and if $k\in\mathcal{N}(j)$ we set $\mathbf{n}_{jk}$ to be the unit normal vector to $\Gamma_{jk}$ pointing out of $\Omega_{j}$.

We can now turn to the acoustic / transport decomposition of the system (1). If we develop the spatial derivatives and isolate the transport terms $(\mathbf{u}\cdot\nabla)\varphi$, where $\varphi=h,h\mathbf{u}$, we can use a splitting operator with respect to time to obtain on one hand the acoustic step

[TABLE]

and on the other hand the transport step

[TABLE]

With these notations, the Lagrange-projection algorithm is defined as follows: for a given discrete state $(h,h\mathbf{u})^{n}_{j}$, $j\in\mathbb{Z}$, defining $(h,h\mathbf{u})^{n+1}_{j}$ is a two-step process defined as follows

Update $(h,h\mathbf{u})^{n}_{j}$ to $(h,h\mathbf{u})^{n+1-}_{j}$ by approximating the solution of system (5), 2. 2.

Update $(h,h\mathbf{u})^{n+1-}_{j}$ to $(h,h\mathbf{u})^{n+1}_{j}$ by approximating the solution of system (6).

Relaxation approximation of the acoustic system. Before entering the details of these two steps in the following section, let us note that if we denote $\tau=1/h$, by simple manipulations system (5) can be recast into:

[TABLE]

Following [CNPT10], we will choose to approximate the solution of system (5) thanks to a Suliciu-relaxation process. More precisely we will solve

[TABLE]

with $p=p^{\mathrm{EOS}}(\tau)=g/(2\tau^{2})$, in the regime $\lambda\rightarrow+\infty$. The parameter $a$ is a constant that is chosen in agreement with the subcharacteristic stability conditions that will be given later. Over the time interval $\left[t^{n},t^{n}+\Delta t\right)$, we can account for the limit $\lambda\rightarrow+\infty$ by setting $\Pi(\mathbf{x},t^{n})=p^{\mathrm{EOS}}(\tau(\mathbf{x},t^{n}))$, and then solving the relaxed system with $\lambda=0$. We add another approximation by supposing that over $\left[t^{n},t^{n}+\Delta t\right)$ it is reasonable to replace $\tau(\mathbf{x},t)\partial_{x_{r}}$ by $\tau(\mathbf{x},t^{n})\partial_{x_{r}}$, $r=1,2$. Finally, we will define our approximation of the acoustic system (5) by solving

[TABLE]

over $\left[t^{n},t^{n}+\Delta t\right)$, with $\Pi(\mathbf{x},t^{n})=p^{\mathrm{EOS}}(\tau(\mathbf{x},t^{n}))$.

Note that the system (8) is rotational invariant. This will allow us in the following to define every flux for the two dimensional problem in the reference frame associated to each face. In this last referential, the problem will be reduced to a quasi-one dimensional problem that we study in the beginning of next section. One can also notice that the eigenstructure of (8) in the phase space $\left\{(h,h\mathbf{u}^{T},\Pi,z)\in\mathbb{R}^{5},h>0,z>0\right\}$ is very simple since it has three eigenvalues $\left\{-a,0,a\right\}$ all associated to linearly degenerated characteristic fields.

4 Finite volume approximation

In this paragraph, we present in details the first-order finite volume scheme associated with the acoustic / transport decomposition of section 3.

4.1 A well-balanced Lagrange-projection finite volume scheme in 1D

We start by considering one-dimensional problems and briefly recall the method proposed in [CKKS17]. In this case the Saint-Venant equations read

[TABLE]

The system associated with the acoustic step reads

[TABLE]

and the system that accounts for transport boils down to

[TABLE]

We suppose given a strictly increasing sequence $x_{j+1/2}\in\mathbb{R}$, for $j\in\mathbb{Z}$ and we consider the set of cells $\Omega_{j}=[x_{j-1/2},x_{j+1/2})$. The local space step is defined by $\Delta x_{j}=x_{j+1/2}-x_{j-1/2}$. We note $\Delta t>0$ the time step and we set $t^{n}=n\Delta t$ for $n\in\mathbb{N}$.

The following discretization strategy was presented in [CKKS17]: the acoustic step (8) is approximated by

[TABLE]

where for all $j$, $\Pi_{j}^{n}=g\frac{{h_{j}^{n}}^{2}}{2}$ and the numerical fluxes $u^{\sharp}_{j+1/2}$, $\Pi^{L,\sharp}_{j+1/2}$ and $\Pi^{R,\sharp}_{j-1/2}$ are defined by

[TABLE]

where $\mathbf{U}$ is the state $\begin{pmatrix}h\\ h\mathbf{u}\\ \Pi\\ z\end{pmatrix}$ and with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

If one chooses $\sharp=n$ (resp. $\sharp=n+1^{-}$) the system (9) provides a time-explicit (resp. time-implicit) discretization of the acoustic system (8). The approximation of the transport step is performed thanks to a standard upwind scheme for $\varphi\in\{h,hu_{1},hu_{2}\}$

[TABLE]

where

[TABLE]

Note that $\{gh\Delta z\}_{\Delta}$ that accounts for the gravity source term is always evaluated at time $t^{n}$, even for the time implicit scheme.

In the above formulas, the parameter $a$ is an approximation of the Lagrangian sound speed $hc=h\sqrt{gh}$ and must satisfy the sub-characteristic condition $a>hc$ which ensures that the relaxed system (7) is a dissipative approximation of the acoustic step (5) (see [Bou04, CC08, CC05, Des10] and the references therein). In order to limit the numerical diffusion we take a local approximation of the Lagrangian sound speed at every interface, given by $a_{j+1/2}=\kappa\max(h_{j}\sqrt{gh_{j}},h_{j+1}\sqrt{gh_{j+1}})$, where $\kappa>1$.

For detailed properties of the numerical scheme (9)-(10) we refer the reader to [CKKS17], nevertheless let us recall that: the overall discretization is conservative in the usual sense of finite volumes methods with respect to $(h,hu_{1},hu_{2})$. Moreover, the scheme is also well-balanced for lake at rest conditions: if $\mathbf{u}^{n}_{j}=0$ and $h^{n}_{j}+z^{n}_{j}=h^{n}_{j+1}+z^{n}_{j+1}$ for all $j\in\mathbb{Z}$, then $h^{n+1}_{j}=h^{n}_{j}$ and $\mathbf{u}^{n+1}_{j}=\mathbf{u}^{n}_{j}$, $j\in\mathbb{Z}$. At last, the time-implicit scheme is stable under a condition which does not depend either on the acoustic system or the sound speed $c$, but which only depends on the transport step and its material velocity $\mathbf{u}$ which is of particular interest in the low-Froude regime.

4.2 Truncation error in the low-Froude regime

In this paragraph, we consider the dimensionless shallow-water equation and we motivate a correction of the above scheme in order to make it efficient in low-Froude regimes. The correction is similar to the one in [CGK16] for low-Mach regimes and we focus on the explicit case $\sharp=n$.

In the following we will say that the flow is in the low Froude regime if $\mathrm{Fr}\ll 1$ and $\partial_{x}p+h\partial_{x}z=\mathcal{O}(\mathrm{Fr}^{2})$. Regarding the dimensionless equations (2), we can observe that, in this regime, the variations of the discharge $hu$ remain of order $1$ as expected.

We can express the fluxes given in the previous section, using the dimensionless quantities, which leads to

[TABLE]

[TABLE]

[TABLE]

if one focuses on the time-explicit scheme for the sake of simplicity.

If we compute the truncation errors in the fluxes above, using the fact that

[TABLE]

we obtain:

[TABLE]

[TABLE]

[TABLE]

At this stage, it is clear that the consistence errors are not uniform with respect to the Froude number in the pressure fluxes. In order to avoid large errors in the numerical diffusion terms when the Froude tends to zero, we propose to correct the flux formula of $\Pi_{\Delta}(\mathbf{U}_{L}^{\sharp},\mathbf{U}_{R}^{\sharp})$ by:

[TABLE]

which amounts to reduce the numerical diffusion on the pressure gradient in the low Froude regime. Indeed, we now get

[TABLE]

[TABLE]

and as long as we take $\theta_{j+1/2}=\mathcal{O}(\mathrm{Fr})$, we recover the uniform consistency of the global scheme with respect to the Froude number. In practice, we will set $\theta_{j+1/2}=\min\left(\frac{\lvert u^{n}_{j+1/2}\rvert}{\max(c_{j},c_{j+1})},1\right)$.

4.3 The Lagrange-projection scheme on 2D unstructured meshes

We now extend the Lagrange-projection scheme in two dimensions. Let $\mathbf{n}\in\mathbb{R}^{2}$ be a unit vector and $\mathbf{U}^{T}=(h,h\mathbf{u}^{T})$, we define

[TABLE]

Following standard lines, we take advantage of the rotational invariance of the acoustic system (8) to define the two-dimensional fluxes numerical fluxes (see for example [GR96, Bou04]). This leads to

[TABLE]

where

[TABLE]

that is to say

[TABLE]

with

[TABLE]

The source term is accounted for by the terms $\Pi_{jk}^{\sharp,\theta}$ since the fluxes $\Pi_{jk}^{\sharp,\theta}\mathbf{n}_{jk}$ in Equation (11b) are not symmetric, indeed $\Pi_{jk}^{\sharp,\theta}\mathbf{n}_{jk}\not=-\Pi_{kj}^{\sharp,\theta}\mathbf{n}_{kj}$ (even if $a_{jk}=a_{kj}$ and $\theta_{jk}=\theta_{kj}$, which is the case in practice).

As far as the transport step is concerned and in order to discretize the system (6), we use an explicit scheme between times $t^{n+1-}$ and $t^{n+1-}+\Delta t$, where the fluxes are chosen upwind with respect to the sign of $u_{jk}^{\sharp}$. If $\varphi\in\{h,h\mathbf{u}\}$, the scheme for the transport step reads

[TABLE]

where

[TABLE]

Note that one can rewrite the transport step (12) as follows

[TABLE]

where $L^{\sharp}_{j}=1+\Delta t\sum_{k\in\mathcal{N}(j)}\sigma_{jk}\,u_{jk}^{\sharp}$. Therefore, replacing the quantities $\varphi_{j}^{n+1-}$ using system (11) gives the following update formulas which take into account the acoustic and transport steps together:

[TABLE]

where the quantities $u_{jk}^{\sharp}$ and $\Pi_{jk}^{\sharp,\theta}$ are computed with $\mathbf{u}^{\sharp}$ and $\Pi^{\sharp}$, and the quantities $h_{jk}^{n+1-}$ and $(h\mathbf{u})_{jk}^{n+1-}$ with $\tau^{n+1-}$ and $\mathbf{u}^{n+1-}$ from system (11).

4.4 Stability and well-balanced properties

Let us first notice that the scheme is conservative with respect to the water height $h$, and with respect to $h\mathbf{u}$ if the topography is flat ($z=\mathrm{cste}$). In particular, it degenerates towards the scheme proposed by [CGK16] adapted to the framework of barotropic Euler system when the bottom is flat. Next, recall that from section 4.2, if $\theta$ is chosen to be like $\mathcal{O}(\mathrm{Fr})$ when $\mathrm{Fr}$ goes to [math], the truncation erro of the numerical scheme is uniform with respect to $\mathrm{Fr}$. At last, assuming that the time step $\Delta t$ is such that the CFL conditions associated to the acoustic step

[TABLE]

and to the transport step

[TABLE]

hold true, the water height $h_{j}^{n}$ is positive for all $j$ and $n>0$ provided that $h_{j}^{0}$ is positive for all $j$, for the time-explicit scheme corresponding to $\sharp=n$. Indeed, notice that $L_{j}^{n}$ turns out to be positive while the transport step correspond to a convex combination of states at time $t^{n+1-}$.

As far as the mixed implicit-explicit scheme corresponding to $\sharp=n+1-$ is concerned, the same properties hold true under the transport CFL condition

[TABLE]

Notice that the acoustic step is implicit and therefore is free of CFL condition.

Now we show that the schemes are well-balanced. We begin with the explicit scheme $\sharp=n$.

Proposition 1**.**

The full explicit scheme ($\sharp=n$) is well-balanced on 2D unstructured mesh in the sense that ($h_{j}^{0}+z_{j}=H=\text{cste}$ and $\mathbf{u}_{j}^{0}=\mathbf{0}$) $\Longrightarrow$ ($h_{j}^{n}+z_{j}=\text{cste}$ and $\mathbf{u}_{j}^{n}=\mathbf{0}$) $\forall j$.

Proof.

Assume that $h_{j}^{0}+z_{j}=H=\text{cste}$ and $\mathbf{u}_{j}^{0}=\mathbf{0}$ $\forall j$. We have

[TABLE]

Injecting those values in the acoustic step gives:

[TABLE]

Next, since $u_{jk}^{0}=0$, $\forall j,k$, the transport step is trivial and the variables $h$ and $h\mathbf{u}$ are unchanged at time $t^{1}$. ∎

Proposition 2**.**

The implicit-explict scheme ($\sharp=n+1-$) is well-balanced on 2D unstructured mesh in the sense that ($h_{j}^{0}+z_{j}=H=\text{cste}$ and $\mathbf{u}_{j}^{0}=\mathbf{0}$) $\Longrightarrow$ ($h_{j}^{n}+z_{j}=\text{cste}$ and $\mathbf{u}_{j}^{n}=\mathbf{0}$) $\forall j$.

Proof.

With the same calculus as in the explicit case, we can verify that the vector $(\mathbf{u}_{j}^{1-},\Pi_{j}^{1-})=(\mathbf{0},\Pi_{j}^{0})$ is the only solution of the coupled system over $(\mathbf{u},\Pi)$. Indeed we have in this case

[TABLE]

and therefore

[TABLE]

Then we easily get $\tau^{1-}_{j}=\tau^{0}_{j}$ so that the Lagrangian step is well-balanced. As before, we can conclude that the transport step is also well-balanced since $u_{jk}^{1-}=0$, $\forall j,k$. ∎

5 Numerical experiments

We present several test cases that aim at testing our scheme against classical flow configurations on unstructured meshes and also in the low Froude regime.

In the following, the EXEX scheme will refer to the full explicit scheme ($\sharp=n$) with the time step $\Delta t$ defined by

[TABLE]

where $K_{\mathrm{CFL}}=0.9$, $v_{jk}^{\mathrm{Acou}}=\tau_{j}a_{jk}$, $a_{jk}=1.01\max(h_{j}c_{j},h_{k}c_{k})$, $c_{j}=\sqrt{gh_{j}}$ and $v_{jk}^{\mathrm{Trans}}=\lvert u_{jk}^{n}\rvert$.

The IMEX scheme will refer to the implicit-explicit scheme ($\sharp=n+1-$) with the time step $\Delta t$ defined by

[TABLE]

where $v_{jk}^{\mathrm{Trans}}=\lvert u_{jk}^{n+1-}\rvert$. Thus, the time step of the IMEX scheme is not constrained by the acoustic waves.

Except if otherwise stated, we will always use the corrected numerical fluxes with

[TABLE]

so that $\theta$ approximates a local Froude number on every edge. On the other hand, for the sake of comparison we take $\theta=1$ on every edge for the fluxes without correction.

5.1 Test of the well-balanced property

In order to test the well-balanced property of the scheme, we first consider the following lake at rest initial condition:

[TABLE]

where $H=0.5$ is constant and the topography $z$ is a smooth bump defined by

[TABLE]

The physical domain $\left[0,1\right]\times\left[0,1\right]$ is discretized over a 20 000-cell triangular mesh. We impose Neumann boundary conditions and we observe the solution at final time $T_{f}=0.1$.

For both EXEX and IMEX schemes, the errors between the numerical and the exact solution, which is also the initial stationary condition, are machine epsilon as we can observe in Table 1.

5.2 Planar dam break test problem

We are interested in the behaviour of our schemes with regard to the propagation of a rarefaction wave and a shock wave. We use the same triangular mesh, boundary conditions and $T_{f}$ value as in section 5.1. The topography is also kept identical to the one given in (14), the velocity initialized to zero and the initial total water height $H=h+z$ is defined as follows:

[TABLE]

In Figure 1 we present the results for both the EXEX and IMEX schemes. We have performed a cut of the solution along the $y=0.5$ axis and compared it with the one computed by a genuine 1D code with a 200-cell uniform grid. We can observe that for both EXEX and IMEX the results of the 2D simulations are in agreement with the 1D results although they were computed with an unstructured mesh. It is worth noting that the 2D simulation manages to fairly preserve the planar structure of the approximate solutions.

5.3 Traveling vortex with flat bottom

In order to challenge our schemes with low Froude regimes, we consider a traveling vortex as in [BALMN14]. The exact solution of this test is detailed in [RB09]. For this test case we consider a flat bottom and we use a regular cartesian mesh of $160\times 160$ cells that discretizes the physical domain $\left[0,1\right]\times\left[0,1\right]$. The boundary conditions imposed are periodic along the $x$-direction and absorbing boundaries along the $y$-direction. The initial conditions are given by:

[TABLE]

where

[TABLE]

and

[TABLE]

Due to the periodic boundary conditions, the exact solution is periodic with period $T=\frac{5}{3}$ and given at any time $t>0$ by:

[TABLE]

We present the results of both the EXEX and IMEX schemes, with ($\theta=\mathcal{O}(\mathrm{Fr})$) and without correction ($\theta=1$) using $\epsilon=0.05$. The mapping of the velocity magnitude is displayed in Figure 2 and we can observe that the accuracy of the solution is really improved by the low-Froude correction. Furthermore, the accuracy of the solution between the EXEX and the IMEX scheme with low-Froude correction is comparable whereas it took about 100 times less time steps and 10 times less CPU time computation to reach the final time with the IMEX than with the EXEX scheme as we can see in Table 2.

5.4 Traveling vortex with non-flat bottom

We extend the physical domain of the traveling vortex test above to the rectangle $\left[0,2\right]\times\left[0,1\right]$. The boundary conditions and initial conditions for $h$ and $u$ are the same as in section 5.3. However we consider here a topography defined by $z(x,y)=10\exp\left(-5(x-1)^{2}-50(y-0.5)^{2}\right)$ following the idea of [BALMN14].

We do not have exact analytical solution because of the non-flat bottom but we still can compare in figure 3 the results between EXEX and IMEX schemes, with or without low Froude correction. Here again, the vortex structure of the flow is completely destroyed by numerical diffusion without low-Froude corrections $\theta=\mathcal{O}(\mathrm{Fr})$, with both schemes. The mapping of the Froude number is not given here, but is similar to the one of the velocity magnitude, with a range of values from $1.6\cdot 10^{-3}$ to $1.1\cdot 10^{-2}$. Finally, we can remark that the EXEX scheme took about 20 times more iterations and 15 times more CPU times than the IMEX scheme, both with low-Froude correction, as we can see in Table 3.

6 Conclusion

We have proposed a large time step and well-balanced scheme for the shallow-water equations in two dimensions for unstructured meshes. We studied the truncation of the scheme with respect to the Froude number $\mathrm{Fr}$ and gave a correction in accordance to the source term. By studying the one-dimensional case, we obtained proposed of modification of the scheme that allows to obtain a uniform truncation error with respect to the Froude number for one-dimensional flows.

Moreover, we showed that the semi-implicit scheme yields good numerical results for flows from low to high Froude values since its CFL condition is based on (slow) material waves only.

Further developments shall include extensions to high-order methods in multiple-dimensions, following for example what has already been achieved with Finite Volume or discontinuous Galerkin methods in 1D. We also intend to adapt the method presented in this work to other compressible flows models involving non-conservative terms.