High Order Finite Difference WENO Methods for Shallow Water Equations on Curvilinear Meshes

Abstract

A high order finite difference numerical scheme is developed for the shallow water equations on curvilinear meshes based on an alternative flux formulation of the weighted essentially non-oscillatory (WENO) scheme. The exact C-property is investigated, and comparison with the standard finite difference WENO scheme is made. Theoretical derivation and numerical results show that the proposed finite difference WENO scheme can maintain the exact C-property on both stationarily and dynamically generalized coordinate systems. The Harten-Lax-van Leer type flux is developed on general curvilinear meshes in two dimensions and verified on a number of benchmark problems, indicating smaller errors compared with the Lax-Friedrichs solver. In addition, we propose a positivity-preserving limiter on stationary meshes such that the scheme can preserve the non-negativity of the water height without loss of mass conservation.

Appendices

Appendix A WENO Interpolation

Given the point values $\{\mathbf{Q }_{i}\}$, we will approximate the values at half points $\mathbf{Q} _{i+1/2}^{\pm }$ as follows.

• (i) Compute an average state $\mathbf{Q} _{i+1/2}$, using either the simple arithmetic mean or a Roe average. Construct the right and left eigenvectors of the Jacobian $\partial \mathbf{F} / \partial \mathbf{Q}$, and denote their matrices by

  $$\begin{aligned} \mathbf {R}_{i+1/2}=\mathbf {R}(\mathbf{Q} _{i+1/2}),\quad \mathbf {R}^{-1}_{i+1/2}=\mathbf {R}^{-1}(\mathbf{Q} _{i+1/2}). \end{aligned}$$

• (ii) Project the conserved quantities $\mathbf{Q}$ to the local characteristic variables $\mathbf{v}$,

  $$\begin{aligned} \mathbf{v} _{j} = \mathbf {R}^{-1}_{i+1/2} \mathbf{Q} _{j}, \quad \text {for } j = i-2, \cdots , i+3. \end{aligned}$$

• (iii) Perform a scalar WENO interpolation on each component of the characteristic variable $\mathbf{v} _{j}$ to obtain the corresponding component of $\mathbf{v} _{i+1/2}^{\pm }$. Here, the procedure of a fifth-order WENO interpolation to obtain the k-th component $v^{-}_{k,i+1/2}$ is described, and the process to obtain $\mathbf{v} ^{+}_{i+1/2}$ is mirror-symmetric to the procedure described above with respect to the target point $x_{i+1/2}$.

  • (a) Choose one big stencil as $S = \{x_{i-2},\cdots ,x_{i+2}\}$, and three small stencils as $S^{(r)} _{i}= \{x_{i-r},\cdots , x_{i-r+2}\}$, $r=0,1,2$. On each small stencil, the standard interpolation gives

    $$\begin{aligned}&v^{(0)}_{k,i+1/2} = \frac{3}{8}v_{k,i} +\frac{3}{4}v_{k,i+1} -\frac{1}{8}v_{k,i+2},\\&v^{(1)}_{k,i+1/2} = -\frac{1}{8}v_{k,i-1} +\frac{3}{4}v_{k,i} +\frac{3}{8}v_{k,i+1},\\&v^{(2)}_{k,i+1/2} = \frac{3}{8}v_{k,i-2} -\frac{5}{4}v_{k,i-1} +\frac{15}{8}v_{k,i}, \\&v^\text {big}_{k,i+1/2} = d_{0} v^{(0)}_{k,i+1/2} +d_{1} v^{(1)}_{k,i+1/2} +d_{2} v^{(2)}_{k,i+1/2}, \end{aligned}$$

    with the linear weights being $d_{0}={5}/{16}$, $d_{1}={5}/{8}$ and $d_{2}={1}/{16}$.

  • (b) Compute nonlinear weights $\omega _r$ from the linear weights $d_r$,

    $$\begin{aligned} \omega _{r}=\frac{\alpha _{r}}{\alpha _{0}+\alpha _{1}+\alpha _{2}}, \quad \alpha _{r}=\frac{d_{r}}{(\beta _{r}+\epsilon )^2}, \quad r=0,1,2, \end{aligned}$$

    where $\epsilon = 10^{-6}$ is used to avoid division by zero, and the smoothness indicators are given by

    $$\left\{ \begin{aligned} \begin{aligned}&\beta _{0} = \frac{13}{12}\left( v_{k,i}-2v_{k,i+1}+v_{k,i+2}\right) ^2 +\frac{1}{4}\left( 3v_{k,i}-4v_{k,i+1}+v_{k,i+2}\right) ^2,\\&\beta _{1} = \frac{13}{12}\left( v_{k,i-1}-2v_{k,i}+v_{k,i+1}\right) ^2 +\frac{1}{4}\left( v_{k,i-1}-v_{k,i+1}\right) ^2,\\&\beta _{2} = \frac{13}{12}\left( v_{k,i-2}-2v_{k,i-1}+v_{k,i}\right) ^2 +\frac{1}{4}\left( v_{k,i-2}-4v_{k,i-1}+3v_{k,i}\right) ^2. \end{aligned} \end{aligned}\right.$$

    (A1)

  • (c) The WENO interpolation is defined by

    $$\begin{aligned} v^{-}_{k,i+1/2}=\sum _{r=0}^{2}\omega _{r}v^{(r)}_{k,i+1/2}. \end{aligned}$$

• (iv) Finally, project $\mathbf{v} _{i+1/2}^{\pm }$ back to the conserved quantities,

  $$\begin{aligned} \mathbf{Q} _{i+1/2}^{\pm } = \mathbf {R}_{i+1/2} \mathbf{v} _{i+1/2}^{\pm }. \end{aligned}$$

Appendix B Construction of $\sigma _{i+1/2}$

Compute the coefficient $\sigma _{i+1/2}$, which is based the smoothness indicators of $\mathbf {Q}^M$ on the WENO stencils.

• (i) Consider each component of $Q^M_l$, $l=1, 2, 3$. On the left biased stencil $\{x_{i-2},\cdots ,x_{i+2}\}$, calculate the smooth indicators on each smooth stencils $\beta _{r,l}$. The formulations of the smooth indicators are the same as those in WENO procedure (A1). We can obtain a candidate for the coefficient $\sigma ^{-}$ by

  $$\begin{aligned} \sigma ^{-}_l = \frac{\sigma _{l,\min }}{\sigma _{l,\max }} \end{aligned}$$

  with

  $$\begin{aligned}&\sigma _{l,\max } = 1 + \frac{|\beta _{0,l} - \beta _{2,l}|}{\epsilon +\min (\beta _{0,l}, \beta _{2,l}) }, \\&\sigma _{l,\min } = 1+ \frac{|\beta _{0,l} - \beta _{2,l}|}{\epsilon +\max (\beta _{0,l}, \beta _{2,l}) }, \end{aligned}$$

  and $\epsilon$ is a small positive number (taken to be $10^{-6}$ in all the numerical examples) to avoid division by zero.

• (ii) A similar formula to another candidate $\sigma ^+$ based on the right biased stencil $\{x_{i-1}, \cdots , x_{i+3}\}$.

• (iii) Then, we take the minimum of those candidates for all components:

  $$\begin{aligned} \sigma _{i+1/2} = \min _{l} (\sigma _{l}^-, \sigma _{l}^+). \end{aligned}$$

Notice that for the still water steady-state solution, i.e., $\mathbf {Q}^M=\text {const}$, the smoothness indicators would be zeroes. Hence, we can easily prove that $\sigma _{i+1/2}=1$.

Appendix C The HLLC Flux

The HLLC flux for shallow water equations in curvilinear coordinates is given by

$$\begin{aligned} \tilde{\mathbf{F }}^{\text {HLLC}}_{i+\frac{1}{2},j}=\left\{ \begin{array}{ll} \tilde{\mathbf{F }}_{\text L},&{}\text {if}\quad 0\leqslant S_{\text L},\\ \tilde{\mathbf{F }}_{*\text L},&{}\text {if}\quad S_{\text L}\leqslant 0\leqslant S_{*},\\ \tilde{\mathbf{F }}_{*\text R},&{}\text {if}\quad S_{*}\leqslant 0\leqslant S_{\text R},\\ \tilde{\mathbf{F }}_{\text R},&{}\text {if}\quad 0\geqslant S_{\text R}, \end{array}\right. \end{aligned}$$

(C1)

where

$$\begin{aligned} \tilde{\mathbf{F }}_{*\text L}&=\frac{S_{\text R}\tilde{\mathbf{F }}_{\text L}-S_{\text L}\tilde{\mathbf{F }}_{\text R}+S_{\text L}S_{\text R}\left( \mathbf{Q} ^{M}_{\text R}-\mathbf{Q} ^{M}_{\text L}\right) }{S_{\text R}-S_{\text L}} \nonumber \\&\quad +\frac{S_{\text L}}{S_{\text R}-S_{\text L}}\begin{pmatrix} 0\\ \frac{\tilde{\xi }_{y}}{\tilde{\xi }^{2}_{x}+\tilde{\xi }^{2}_{y}}h_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text R}-S_{\text R}\right) \\ \frac{-\tilde{\xi }_{x}}{\tilde{\xi }^{2}_{x}+\tilde{\xi }^{2}_{y}}h_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text R}-S_{\text R}\right) \end{pmatrix}, \end{aligned}$$

(C2)

$$\begin{aligned} \tilde{\mathbf{F }}_{*\text R}&=\frac{S_{\text R}\tilde{\mathbf{F }}_{\text L}-S_{\text L}\tilde{\mathbf{F }}_{\text R}+S_{\text L}S_{\text R}\left( \mathbf{Q} ^{M}_{\text R}-\mathbf{Q} ^{M}_{\text L}\right) }{S_{\text R}-S_{\text L}}\nonumber \\&\quad +\frac{S_{\text R}}{S_{\text R}-S_{\text L}}\begin{pmatrix} 0\\ \frac{\tilde{\xi }_{y}}{\tilde{\xi }^{2}_{x}+\tilde{\xi }^{2}_{y}}h_{\text L}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text L}-S_{\text L}\right) \\ \frac{-\tilde{\xi }_{x}}{\tilde{\xi }^{2}_{x}+\tilde{\xi }^{2}_{y}}h_{\text L}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text L}-S_{\text L}\right) \end{pmatrix}. \end{aligned}$$

(C3)

Here $\tilde{u}=\tilde{\xi }_{t}+\tilde{\xi }_{x}u+\tilde{\xi }_{y}v,\tilde{v}=\tilde{\xi }_{y}u-\tilde{\xi }_{x}v$. The wave speed estimates $S_{\text L}$ and $S_{\text R}$ are computed by

$$\begin{aligned} S_{\text L}&=\min \left( \tilde{u}_{\text L}-\sqrt{\left( \tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}\right) gh_{\text L}}, \tilde{u}_{\text R}-\sqrt{\left( \tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}\right) gh_{\text R}} \right) , \end{aligned}$$

(C4)

$$\begin{aligned} S_{\text R}&=\max \left( \tilde{u}_{\text L}+\sqrt{\left( \tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}\right) gh_{\text L}}, \tilde{u}_{\text R}+\sqrt{\left( \tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}\right) gh_{\text R}} \right) . \end{aligned}$$

(C5)

Take $\tilde{u}=\tilde{\eta }_{t}+\tilde{\eta }_{x}u+\tilde{\eta }_{y}v, \tilde{v}=\tilde{\eta }_{y}u-\tilde{\eta }_{x}v$ instead, we have

$$\begin{aligned} \tilde{\mathbf{G }}^{\text {HLLC}}_{i,j+\frac{1}{2}}=\left\{ \begin{array}{ll} \tilde{\mathbf{G }}_{\text L},&{}\text {if}\quad 0\leqslant S_{\text L},\\ \tilde{\mathbf{G }}_{*\text L},&{}\text {if}\quad S_{\text L}\leqslant 0\leqslant S_{*},\\ \tilde{\mathbf{G }}_{*\text R},&{}\text {if}\quad S_{*}\leqslant 0\leqslant S_{\text R},\\ \tilde{\mathbf{G }}_{\text R},&{}\text {if}\quad 0\geqslant S_{\text R}, \end{array}\right. \end{aligned}$$

(C6)

where

$$\begin{aligned} \tilde{\mathbf{G }}_{*\text L}= &\frac{S_{\text R}\tilde{\mathbf{G }}_{\text L}-S_{\text L}\tilde{\mathbf{G }}_{\text R}+S_{\text L}S_{\text R}\left( \mathbf{Q} ^{M}_{\text R}-\mathbf{Q} ^{M}_{\text L}\right) }{S_{\text R}-S_{\text L}}\nonumber \\&\quad +\frac{S_{\text L}}{S_{\text R}-S_{\text L}}\begin{pmatrix} 0\\ \frac{\tilde{\eta }_{y}}{\tilde{\eta }^{2}_{x}+\tilde{\eta }^{2}_{y}}h_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text R}-S_{\text R}\right) \\ \frac{-\tilde{\eta }_{x}}{\tilde{\eta }^{2}_{x}+\tilde{\eta }^{2}_{y}}h_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text R}-S_{\text R}\right) \end{pmatrix}, \end{aligned}$$

(C7)

$$\begin{aligned} \tilde{\mathbf{G }}_{*\text R}= &\frac{S_{\text R}\tilde{\mathbf{G }}_{\text L}-S_{\text L}\tilde{\mathbf{G }}_{\text R}+S_{\text L}S_{\text R}\left( \mathbf{Q} ^{M}_{\text R}-\mathbf{Q} ^{M}_{\text L}\right) }{S_{\text R}-S_{\text L}}\nonumber \\&\quad +\frac{S_{\text R}}{S_{\text R}-S_{\text L}}\begin{pmatrix} 0\\ \frac{\tilde{\eta }_{y}}{\tilde{\eta }^{2}_{x}+\tilde{\eta }^{2}_{y}}h_{\text L}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text L}-S_{\text L}\right) \\ \frac{-\tilde{\eta }_{x}}{\tilde{\eta }^{2}_{x}+\tilde{\eta }^{2}_{y}}h_{\text L}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text L}-S_{\text L}\right) \end{pmatrix} \end{aligned}$$

(C8)

and

$$\begin{aligned} S_{\text L}&=\min \left( \tilde{u}_{\text L}-\sqrt{\left( \tilde{\eta }_{x}^{2}+\tilde{\eta }_{y}^{2}\right) gh_{\text L}}, \tilde{u}_{\text R}-\sqrt{\left( \tilde{\eta }_{x}^{2}+\tilde{\eta }_{y}^{2}\right) gh_{\text R}} \right) , \end{aligned}$$

(C9)

$$\begin{aligned} S_{\text R}&=\max \left( \tilde{u}_{\text L}+\sqrt{\left( \tilde{\eta }_{x}^{2}+\tilde{\eta }_{y}^{2}\right) gh_{\text L}}, \tilde{u}_{\text R}+\sqrt{\left( \tilde{\eta }_{x}^{2}+\tilde{\eta }_{y}^{2}\right) gh_{\text R}} \right) . \end{aligned}$$

(C10)

We give the proof of $\tilde{\mathbf{F }}_{*\text L}$ as an example, then $\tilde{\mathbf{F }}_{*\text R},\tilde{\mathbf{G }}_{*\text L},\tilde{\mathbf{G }}_{*\text R}$ do not cause any trouble. Recall that

$$\begin{aligned} \tilde{\mathbf{F }}=\tilde{\xi }_{t}{} \mathbf{Q} +\tilde{\xi }_{x}{} \mathbf{F} +\tilde{\xi }_{y}{} \mathbf{G} =\begin{pmatrix} h\tilde{u}\\ h\tilde{u}u+\frac{1}{2}gh^{2}\tilde{\xi }_{x}\\ h\tilde{u}v+\frac{1}{2}gh^{2}\tilde{\xi }_{y} \end{pmatrix}, \end{aligned}$$

where $\tilde{u}=\tilde{\xi }_{t}+\tilde{\xi }_{x}u+\tilde{\xi }_{y}v, \tilde{v}=\tilde{\xi }_{y}u-\tilde{\xi }_{x}v$. Let

$$\begin{aligned} \bar{\mathbf{F }}=\begin{pmatrix} h\tilde{u}\\ h\tilde{u}^{2}+\frac{1}{2}gh^{2}\left( \tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}\right) \\ h\tilde{u}\tilde{v} \end{pmatrix}, \end{aligned}$$

and consider the matrix

$$\begin{aligned} \mathbf{M} =\begin{pmatrix} 1&{}0&{}0\\ \tilde{\xi }_{t}&{}\tilde{\xi }_{x}&{}\tilde{\xi }_{y}\\ 0&{}\tilde{\xi }_{y}&{}-\tilde{\xi }_{x} \end{pmatrix},\quad \mathbf{M} ^{-1}=\begin{pmatrix} 1&{}0&{}0\\ \frac{-\tilde{\xi }_{t}\tilde{\xi }_{x}}{\tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}}&{}\frac{\tilde{\xi }_{x}}{\tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}}&{}\frac{\tilde{\xi }_{y}}{\tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}}\\ \frac{-\tilde{\xi }_{t}\tilde{\xi }_{y}}{\tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}}&{}\frac{\tilde{\xi }_{y}}{\tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}}&{}\frac{-\tilde{\xi }_{x}}{\tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}}\\ \end{pmatrix}. \end{aligned}$$

Then we have $\mathbf{M} \tilde{\mathbf{F }}=\bar{\mathbf{F }}$. We approximate the HLLC flux of $\tilde{\mathbf{F }}$ by

$$\begin{aligned} \tilde{\mathbf{F }}^{\text {HLLC}}\approx \mathbf{M} ^{-1}\bar{\mathbf{F }}^{\text {HLLC}}. \end{aligned}$$

Since the third component of $\bar{\mathbf{F }}$ can be expressed by the first component and variable $\tilde{v}$, we may consider the first two components of $\bar{\mathbf{F }}$ at first, which actually is the HLL flux for SWEs problem in one dimension. Denote $S_{\text L}$ and $S_{\text R}$ as the estimates of wave speed $\tilde{u}-\sqrt{\left( \tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}\right) gh}$ and $\tilde{u}+\sqrt{\left( \tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}\right) gh}$. When $S_{\text L}\leqslant 0\leqslant S_{\text R}$, applying Rankine-Hugoniot conditions across the left and right waves respectively, we obtain

$$\begin{aligned} \frac{\bar{\mathbf{F }}^{\text {HLL}}-\bar{\mathbf{F }}_{\text L}}{\bar{\mathbf{Q }}^{\text {HLL}}-\bar{\mathbf{Q }}_{\text L}}&=S_{\text L},\quad \quad \quad \frac{\bar{\mathbf{F }}^{\text {HLL}}-\bar{\mathbf{F }}_{\text R}}{\bar{\mathbf{Q }}^{\text {HLL}}-\bar{\mathbf{Q }}_{\text R}}=S_{\text R}, \end{aligned}$$

where $\bar{\mathbf{F }}=\left( h\tilde{u},h\tilde{u}^{2}+\frac{1}{2}gh^{2}\left( \tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}\right) \right) ^{\mathrm {T}}$ and $\bar{\mathbf{Q }}=\left( h, h\tilde{u}\right) ^{\mathrm {T}}$. Combine the two formulas and substitute $\bar{\mathbf{Q }}^{\text {HLL}}$, and we have

$$\begin{aligned} \bar{\mathbf{F }}^{\text {HLL}}=\frac{S_{\text R}\bar{\mathbf{F }}_{\text L}-S_{\text L}\bar{\mathbf{F }}_{\text R}+S_{\text L}S_{\text R}\left( \bar{\mathbf{Q }}_{\text R}-\bar{\mathbf{Q }}_{\text L}\right) }{S_{\text R}-S_{\text L}}. \end{aligned}$$

Therefore, the flux can be written as

$$\begin{aligned} \bar{\mathbf{F }}^{\text {HLL}}_{i+\frac{1}{2}}=\left\{ \begin{array}{ll} \bar{\mathbf{F }}_{\text L},&{}\text {if}\quad 0\leqslant S_{\text L},\\ \frac{S_{\text R}\bar{\mathbf{F }}_{\text L}-S_{\text L}\bar{\mathbf{F }}_{\text R}+S_{\text L}S_{\text R}\left( \bar{\mathbf{Q }}_{\text R}-\bar{\mathbf{Q }}_{\text L}\right) }{S_{\text R}-S_{\text L}},&{}\text {if}\quad S_{\text L}\leqslant 0\leqslant S_{\text R},\\ \bar{\mathbf{F }}_{\text R},&{}\text {if}\quad 0\geqslant S_{\text R}. \end{array}\right. \end{aligned}$$

Then we consider all components together, the eigenvalues of $\bar{\mathbf{F }}$ are

$$\begin{aligned} \tilde{u}-\sqrt{\left( \tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}\right) gh},\quad \tilde{u},\quad \tilde{u}+\sqrt{\left( \tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}\right) gh}, \end{aligned}$$

which respectively represents the wave speeds $S_{\text L}$, $S_{*}$, and $S_{\text R}$. The third component of the flux can be expressed in terms of the first component and the variable $\tilde{v}$, that is, $\bar{\mathbf{F }}^{(3)}=\bar{\mathbf{F }}^{(1)}\tilde{v}$. Note that it is only the third component of the flux changes across the middle wave, we retain the HLL flux for the first two components of the flux and combine $\tilde{v}$ with the first component of the HLL flux. One can obtain an expression for the third component as

$$\begin{aligned} \bar{\mathbf{F }}_{i+\frac{1}{2},j}^{(3)}=\left\{ \begin{array}{ll} \bar{\mathbf{F }}_{i+\frac{1}{2},j}^{(1)}\tilde{v}_{l},\quad \text {if}\quad S_{*}\geqslant 0,\\ \bar{\mathbf{F }}_{i+\frac{1}{2},j}^{(1)}\tilde{v}_{r},\quad \text {if}\quad S_{*}< 0, \end{array}\right. \end{aligned}$$

where

$$\begin{aligned} S_{*}=\frac{h_{\text L}S_{\text R}\left( \tilde{u}_{\text L}-S_{\text L}\right) -h_{\text R}S_{\text L}\left( \tilde{u}_{\text R}-S_{\text R}\right) }{h_{\text L}\left( \tilde{u}_{\text L}-S_{\text L}\right) -h_{\text R}\left( \tilde{u}_{\text R}-S_{\text R}\right) }. \end{aligned}$$

If $S_{\text L}\leqslant 0\leqslant S_{*}$,

$$\begin{aligned} \bar{\mathbf{F }}^{\text {HLLC}}_{*\text L}&=\frac{1}{S_{\text R}-S_{\text L}} \left[S_{\text R}\begin{pmatrix} h_{\text L}\tilde{u}_{\text L}\\ h_{\text L}\tilde{u}_{\text L}^{2}+\frac{1}{2}gh_{\text L}^{2}\tilde{\xi }^{2}\\ h_{\text L}\tilde{u}_{\text L}\tilde{v}_{\text L} \end{pmatrix}-S_{\text L}\begin{pmatrix} h_{\text R}\tilde{u}_{\text R}\\ h_{\text R}\tilde{u}_{\text R}^{2}+\frac{1}{2}gh_{\text R}^{2}\tilde{\xi }^{2}\\ h_{\text R}\tilde{u}_{\text R}\tilde{v}_{\text L} \end{pmatrix}\right. \\&\quad \left. +S_{\text L}S_{\text R}\begin{pmatrix} h_{\text R}-h_{\text L}\\ h_{\text R}\tilde{u}_{\text R}-h_{\text L}\tilde{u}_{\text L}\\ h_{\text R}\tilde{v}_{\text L}-h_{\text L}\tilde{v}_{\text L} \end{pmatrix} \right] \\&=\frac{1}{S_{\text R}-S_{\text L}} \left[ S_{\text R}\begin{pmatrix} h_{\text L}\tilde{u}_{\text L}\\ h_{\text L}\tilde{u}_{\text L}^{2}+\frac{1}{2}gh_{\text L}^{2}\tilde{\xi }^{2}\\ h_{\text L}\tilde{u}_{\text L}\tilde{v}_{\text L} \end{pmatrix}-S_{\text L}\begin{pmatrix} h_{\text R}\tilde{u}_{\text R}\\ h_{\text R}\tilde{u}_{\text R}^{2}+\frac{1}{2}gh_{\text R}^{2}\tilde{\xi }^{2}\\ h_{\text R}\tilde{u}_{\text R}\tilde{v}_{\text R} \end{pmatrix} \right. \\&\quad +S_{\text L}S_{\text R}\begin{pmatrix} h_{\text R}-h_{\text L}\\ h_{\text R}\tilde{u}_{\text R}-h_{\text L}\tilde{u}_{\text L}\\ h_{\text R}\tilde{v}_{\text R}-h_{\text L}\tilde{v}_{\text L} \end{pmatrix}\\&\quad \left. +S_{\text L}\begin{pmatrix} 0\\ 0\\ h_{\text R}\tilde{u}_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \end{pmatrix}-S_{\text L}S_{\text R}\begin{pmatrix} 0\\ 0\\ h_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \end{pmatrix} \right] \\&=\frac{1}{S_{\text R}-S_{\text L}} [ S_{\text R}\bar{\mathbf{F }}_{\text L}-S_{\text L}\bar{\mathbf{F }}_{\text R}+S_{\text L}S_{\text R}\left( \bar{\mathbf{Q }}_{\text R}-\bar{\mathbf{Q }}_{\text L}\right) ] \\&\quad +\frac{S_{\text L}}{S_{\text R}-S_{\text L}}\begin{pmatrix} 0\\ 0\\ h_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text R}-S_{\text R}\right) \end{pmatrix}, \end{aligned}$$

where $\tilde{\xi }^{2}=\tilde{\xi }_{x}^{2}+\tilde{\xi }_{y}^{2}$. Thus we obtain our approximation of flux

$$\begin{aligned} \tilde{\mathbf{F }}^{\text {HLLC}}_{*\text L}&=\mathbf{M} ^{-1}\bar{\mathbf{F }}^{\text {HLLC}}_{*\text L}\\&=\frac{1}{S_{\text R}-S_{\text L}} [ S_{\text R}{} \mathbf{M} ^{-1}\bar{\mathbf{F }}_{\text L}-S_{\text L}{} \mathbf{M} ^{-1}\bar{\mathbf{F }}_{\text R}+S_{\text L}S_{\text R}( \mathbf{M} ^{-1}\bar{\mathbf{Q }}_{\text R}-\mathbf{M} ^{-1}\bar{\mathbf{Q }}_{\text L}) ] \\&\quad +\frac{S_{\text L}}{S_{\text R}-S_{\text L}}{} \mathbf{M} ^{-1}\begin{pmatrix} 0\\ 0\\ h_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text R}-S_{\text R}\right) \end{pmatrix}\\&=\frac{1}{S_{\text R}-S_{\text L}}\left\{ S_{\text R}\tilde{\mathbf{F }}_{\text L}-S_{\text L}\tilde{\mathbf{F }}_{\text R}+S_{\text L}S_{\text R}\left( \mathbf{Q} _{\text R}-\mathbf{Q} _{\text L}\right) \right\} \\&\quad +\frac{S_{\text L}}{S_{\text R}-S_{\text L}}\begin{pmatrix} 1&{}0&{}0\\ \frac{-\tilde{\xi }_{t}\tilde{\xi }_{x}}{\tilde{\xi }^{2}}&{}\frac{\tilde{\xi }_{x}}{\tilde{\xi }^{2}}&{}\frac{\tilde{\xi }_{y}}{\tilde{\xi }^{2}}\\ \frac{-\tilde{\xi }_{t}\tilde{\xi }_{y}}{\tilde{\xi }^{2}}&{}\frac{\tilde{\xi }_{y}}{\tilde{\xi }^{2}}&{}\frac{-\tilde{\xi }_{x}}{\tilde{\xi }^{2}}\\ \end{pmatrix}\begin{pmatrix} 0\\ 0\\ h_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text R}-S_{\text R}\right) \end{pmatrix}\\&=\frac{S_{\text R}\tilde{\mathbf{F }}_{\text L}-S_{\text L}\tilde{\mathbf{F }}_{\text R}+S_{\text L}S_{\text R}\left( \mathbf{Q} _{\text R}-\mathbf{Q} _{\text L}\right) }{S_{\text R}-S_{\text L}}+\frac{S_{\text L}}{S_{\text R}-S_{\text L}}\begin{pmatrix} 0\\ \frac{\tilde{\xi }_{y}}{\tilde{\xi }^{2}_{x}+\tilde{\xi }^{2}_{y}}h_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text R}-S_{\text R}\right) \\ \frac{-\tilde{\xi }_{x}}{\tilde{\xi }^{2}_{x}+\tilde{\xi }^{2}_{y}}h_{\text R}\left( \tilde{v}_{\text R}-\tilde{v}_{\text L}\right) \left( \tilde{u}_{\text R}-S_{\text R}\right) \end{pmatrix}. \end{aligned}$$

Finally, we change $\mathbf {Q}_{\text L,\text R}$ to $\mathbf {Q}^M_{\text L,\text R}$, in order to satisfy the still water steady state condition, namely

$$\begin{aligned} u=v=0,\quad h+b= \text {const}. \end{aligned}$$

For the still water steady state condition, the variables $\tilde{u}$ and $\tilde{v}$ become

$$\begin{aligned} \tilde{u}&=\tilde{\xi }_{t}+\tilde{\xi }_{x}u+\tilde{\xi }_{y}v=\tilde{\xi }_{t},\nonumber \\ \tilde{v}&=\tilde{\xi }_{y}u-\tilde{\xi }_{x}v=0. \end{aligned}$$

Hence,

$$\begin{aligned} \mathbf{Q} ^{M}_{\text L}=\mathbf{Q} ^{M}_{\text R}, \end{aligned}$$

and

$$\begin{aligned} \tilde{v}_{\text R}=\tilde{v}_{\text L}=0. \end{aligned}$$

Finally, we obtain the desired convex combination expression

$$\begin{aligned} \tilde{\mathbf{F }}_{*\text L}^{\text {HLLC}}=\frac{S_{\text R}\tilde{\mathbf{F }}_{\text L}-S_{\text L}\tilde{\mathbf{F }}_{\text R}}{S_{\text R}-S_{\text L}}=\alpha \tilde{\mathbf{F }}_{\text L}+\beta \tilde{\mathbf{F }}_{\text R},\quad \alpha +\beta =1, \end{aligned}$$

which directly leads to the exact C-property.

Appendix D Design of $\Lambda$ in PP-Limiter

• (i) If $p_{i+1/2,j}\geqslant 0$, $p_{i-1/2,j}\geqslant 0$, $q_{i,j+1/2}\geqslant 0$, and $q_{i,j-1/2}\geqslant 0$, then

  $$\begin{aligned} \Lambda ^{-}_{i+1/2,j} = \Lambda ^{+}_{i-1/2,j} = \Lambda ^{-}_{i,j+1/2} = \Lambda ^{+}_{i,j-1/2} = 1. \end{aligned}$$

• (ii) If $p_{i+1/2,j}< 0$, $p_{i-1/2,j}\geqslant 0$, $q_{i,j+1/2}\geqslant 0$, and $q_{i,j-1/2}\geqslant 0$, then

  $$\begin{aligned} \Lambda ^{-}_{i+1/2,j} = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i+1/2,j}}\right) ,\quad \Lambda ^{+}_{i-1/2,j} = \Lambda ^{-}_{i,j+1/2} = \Lambda ^{+}_{i,j-1/2} = 1. \end{aligned}$$

• (iii) If $p_{i+1/2,j}\geqslant 0$, $p_{i-1/2,j}<0$, $q_{i,j+1/2}\geqslant 0$, and $q_{i,j-1/2}\geqslant 0$, then

  $$\begin{aligned} \Lambda ^{+}_{i-1/2,j} = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i-1/2,j}}\right) ,\quad \Lambda ^{-}_{i+1/2,j} = \Lambda ^{-}_{i,j+1/2} = \Lambda ^{+}_{i,j-1/2} = 1. \end{aligned}$$

• (iv) If $p_{i+1/2,j}\geqslant 0$, $p_{i-1/2,j}\geqslant 0$, $q_{i,j+1/2}<0$, and $q_{i,j-1/2}\geqslant 0$, then

  $$\begin{aligned} \Lambda ^{-}_{i,j+1/2} = \min \left( 1,\frac{\Gamma _{i,j}}{ q_{i,j+1/2}}\right) ,\quad \Lambda ^{+}_{i+1/2,j} = \Lambda ^{-}_{i-1/2,j} = \Lambda ^{+}_{i,j-1/2} = 1. \end{aligned}$$

• (v) If $p_{i+1/2,j}\geqslant 0$, $p_{i-1/2,j}\geqslant 0$, $q_{i,j+1/2}\geqslant 0$, and $q_{i,j-1/2}<0$, then

  $$\begin{aligned} \Lambda ^{+}_{i,j-1/2} = \min \left( 1,\frac{\Gamma _{i,j}}{ q_{i,j-1/2}}\right) ,\quad \Lambda ^{+}_{i+1/2,j} = \Lambda ^{-}_{i-1/2,j} = \Lambda ^{-}_{i,j+1/2} = 1. \end{aligned}$$

• (vi) If $p_{i+1/2,j}< 0$, $p_{i-1/2,j}<0$, $q_{i,j+1/2}\geqslant 0$, and $q_{i,j-1/2}\geqslant 0$, then

  $$\begin{aligned} \Lambda ^{-}_{i+1/2,j} =\Lambda ^{+}_{i-1/2,j} = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i+1/2,j}+p_{i-1/2,j}}\right) ,\quad \Lambda ^{-}_{i,j+1/2} = \Lambda ^{+}_{i,j-1/2} = 1. \end{aligned}$$

• (vii) If $p_{i+1/2,j}< 0$, $p_{i-1/2,j}\geqslant 0$, $q_{i,j+1/2}<0$, and $q_{i,j-1/2}\geqslant 0$, then

  $$\begin{aligned} \Lambda ^{-}_{i+1/2,j} =\Lambda ^{-}_{i,j+1/2} = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i+1/2,j}+q_{i,j+1/2}}\right) ,\quad \Lambda ^{+}_{i-1/2,j} = \Lambda ^{+}_{i,j-1/2} = 1. \end{aligned}$$

• (viii) If $p_{i+1/2,j}< 0$, $p_{i-1/2,j}\geqslant 0$, $q_{i,j+1/2}\geqslant 0$, and $q_{i,j-1/2}<0$, then

  $$\begin{aligned} \Lambda ^{-}_{i+1/2,j} =\Lambda ^{+}_{i,j-1/2} = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i+1/2,j}+q_{i,j-1/2}}\right) ,\quad \Lambda ^{+}_{i-1/2,j} = \Lambda ^{-}_{i,j+1/2} = 1. \end{aligned}$$

• (ix) If $p_{i+1/2,j}\geqslant 0$, $p_{i-1/2,j}<0$, $q_{i,j+1/2}<0$, and $q_{i,j-1/2}\geqslant 0$, then

  $$\begin{aligned} \Lambda ^{+}_{i-1/2,j} =\Lambda ^{-}_{i,j+1/2} = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i-1/2,j}+q_{i,j+1/2}}\right) ,\quad \Lambda ^{-}_{i+1/2,j} = \Lambda ^{+}_{i,j-1/2} = 1. \end{aligned}$$

• (x) If $p_{i+1/2,j}\geqslant 0$, $p_{i-1/2,j}<0$, $q_{i,j+1/2}\geqslant 0$, and $q_{i,j-1/2}<0$, then

  $$\begin{aligned} \Lambda ^{+}_{i-1/2,j} =\Lambda ^{+}_{i,j-1/2} = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i-1/2,j}+q_{i,j-1/2}}\right) ,\quad \Lambda ^{-}_{i+1/2,j} =\Lambda ^{-}_{i,j+1/2} = 1. \end{aligned}$$

• (xi) If $p_{i+1/2,j}\geqslant 0$, $p_{i-1/2,j}\geqslant 0$, $q_{i,j+1/2}<0$, and $q_{i,j-1/2}<0$, then

  $$\begin{aligned} \Lambda ^{-}_{i,j+1/2} =\Lambda ^{+}_{i,j-1/2} = \min \left( 1,\frac{\Gamma _{i,j}}{ q_{i,j+1/2}+q_{i,j-1/2}}\right) ,\quad \Lambda ^{-}_{i+1/2,j} =\Lambda ^{+}_{i-1/2,j} = 1. \end{aligned}$$

• (xii) If $p_{i+1/2,j}< 0$, $p_{i-1/2,j}< 0$, $q_{i,j+1/2}< 0$, and $q_{i,j-1/2}\geqslant 0$, then

  $$\begin{aligned}&\Lambda ^{-}_{i+1/2,j} =\Lambda ^{+}_{i-1/2,j} \\&\quad =\Lambda ^{-}_{i,j+1/2} = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i+1/2,j}+p_{i-1/2,j}+q_{i,j+1/2}}\right) ,\quad \Lambda ^{+}_{i,j-1/2} = 1. \end{aligned}$$

• (xiii) If $p_{i+1/2,j}< 0$, $p_{i-1/2,j}<0$, $q_{i,j+1/2}\geqslant 0$, and $q_{i,j-1/2}<0$, then

  $$\begin{aligned}&\Lambda ^{-}_{i+1/2,j} =\Lambda ^{+}_{i-1/2,j} =\Lambda ^{+}_{i,j-1/2} \\&\quad = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i+1/2,j}+p_{i-1/2,j}+q_{i,j-1/2}}\right) ,\quad \Lambda ^{-}_{i,j+1/2} = 1. \end{aligned}$$

• (xiv) If $p_{i+1/2,j}< 0$, $p_{i-1/2,j}\geqslant 0$, $q_{i,j+1/2}<0$, and $q_{i,j-1/2}<0$, then

  $$\begin{aligned}&\Lambda ^{-}_{i+1/2,j} =\Lambda ^{-}_{i,j+1/2} =\Lambda ^{+}_{i,j-1/2}\\&\quad = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i+1/2,j}+q_{i,j+1/2}+q_{i,j-1/2}}\right) ,\quad \Lambda ^{+}_{i-1/2,j} = 1. \end{aligned}$$

• (xv) If $p_{i+1/2,j}\geqslant 0$, $p_{i-1/2,j}<0$, $q_{i,j+1/2}<0$, and $q_{i,j-1/2}<0$, then

  $$\begin{aligned}&\Lambda ^{+}_{i-1/2,j} =\Lambda ^{-}_{i,j+1/2} =\Lambda ^{+}_{i,j-1/2} \\&\quad = \min \left( 1,\frac{\Gamma _{i,j}}{ p_{i-1/2,j}+q_{i,j+1/2}+q_{i,j-1/2}}\right) ,\quad \Lambda ^{-}_{i+1/2,j} = 1. \end{aligned}$$

• (xvi) If $p_{i+1/2,j}< 0$, $p_{i-1/2,j}<0$, $q_{i,j+1/2}<0$, and $q_{i,j-1/2}<0$,

  • if (48) is satisfied with $\theta _{i\pm 1/2,j} = \theta _{i,j\pm 1/2} =1$, then

    $$\begin{aligned} \Lambda ^{+}_{i,j-1/2} = \Lambda ^{-}_{i+1/2,j} = \Lambda ^{+}_{i-1/2,j} = \Lambda ^{-}_{i,j+1/2} = 1; \end{aligned}$$

  • else

    $$\begin{aligned}&\Lambda ^{-}_{i+1/2,j} = \Lambda ^{+}_{i-1/2,j} = \Lambda ^{-}_{i,j+1/2} = \Lambda ^{+}_{i,j-1/2} = \frac{\Gamma _{i,j}}{ p_{i+1/2,j} + p_{i-1/2,j} + q_{i,j+1/2} + q_{i,j-1/2} }. \end{aligned}$$