A Hybrid Method to Solve Shallow Water Flows with Horizontal Density Gradients

Abstract

In this work, a model for shallow water flows that accounts for the effects of horizontal density fluctuations is presented and derived. While the density is advected by the flow, a two-way feedback between the density gradients and the time evolution of the fluid is ensured through the pressure and source terms in the momentum equations. The model can be derived by vertically averaging the Euler equations while still allowing for density fluctuations in horizontal directions. The approach differs from multi-layer shallow water flows where two or more layers are considered, each of them having their own depth, velocity and constant density. A Roe-type upwind scheme is developed and the Roe matrices are computed systematically by going from the conservative to the quasi-linear form at a discrete level. Properties of the model are analyzed. The system is hyperbolic with two shock-wave families and a contact discontinuity associated to interfaces of regions with density jumps. This new field is degenerate with pressure and velocity as the corresponding Riemann invariants. We show that in some parameter regimes numerically recognizing such invariants across contact discontinuities is important to correctly compute the flow near those interfaces. We present a numerical algorithm that correctly captures all waves with a hybrid strategy. The method integrates the Riemann invariants near contact discontinuities and switches back to the conserved variables away from it to properly resolve shock waves. This strategy can be applied to any numerical scheme. Numerical solutions for a variety of tests in one and two dimensions are shown to illustrate the advantages of the strategy and the merits of the scheme.

Appendices

Appendix A: Details and Properties of the 2-D Hybrid Numerical Scheme

The hybrid strategy is now carried over to the two-dimensional space. Local linearizations of eigenvalues, eigenvectors, and amplitudes are needed for both the conservative and primitive formulations in 2-D in order to apply the numerical scheme (32), (33), (34).

Conservative Formulation

Let us first consider the 2-D conservative formulation (1) and its quasilinear form (2) for the variables $\mathbf{W }^c = (\rho h, \; \rho h u, \; \rho h v, \; h)^T$. Although not necessary in the description of the numerical scheme, we include the corresponding local linearizations of the Roe matrices (3) and (4) and source terms

$$\begin{aligned}&{\hat{A}} = \begin{pmatrix} 0 &{} 1 &{} 0 &{} 0 \\ \frac{g}{2} \overline{h}-{\hat{u}}^2 &{} 2 {\hat{u}} &{} 0 &{} \frac{g}{2} \overline{\rho h} \\ -{\hat{u}} {\hat{v}} &{} {\hat{v}} &{} {\hat{u}} &{} 0 \\ -{\hat{u}} \widehat{1/\rho } &{} \widehat{1/\rho } &{} 0 &{} {\hat{u}} \end{pmatrix}, \Delta x \mathbf{S }^A = \begin{pmatrix} 0 \\ - \bar{\rho }{\tilde{c}}^2 \Delta _x B \\ 0 \\ 0 \end{pmatrix},\nonumber \\&{\hat{B}} = \begin{pmatrix} 0 &{} 0 &{} 1 &{} 0 \\ -{\hat{v}} {\hat{u}} &{} {\hat{v}} &{} {\hat{u}} &{} 0 \\ \frac{g}{2} \overline{h} -{\hat{v}}^2 &{} 0 &{} 2 {\hat{v}} &{} \frac{g}{2} \overline{\rho h} \\ -{\hat{v}} \widehat{1/\rho } &{} 0 &{} \widehat{1/\rho } &{} {\hat{v}} \end{pmatrix}, \Delta x \mathbf{S }^B = \begin{pmatrix} 0 \\ 0 \\ -\bar{\rho }{\tilde{c}}^2 \Delta _y B \\ 0 \end{pmatrix}. \end{aligned}$$

(60)

The corresponding matrices of eigenvalues are $\Lambda ^A = \text {diag}( {\hat{u}}-{\tilde{c}}, {\hat{u}}, {\hat{u}}, {\hat{u}}+{\tilde{c}})$ and $\Lambda ^B = \text {diag}( {\hat{v}}-{\tilde{c}}, {\hat{v}}, {\hat{v}}, {\hat{v}}+{\tilde{c}})$. The matrices of eigenvectors needed in the numerical scheme are

$$\begin{aligned} \mathbf{R }^{c,A} = \begin{pmatrix} 1 &{} 1 &{} 0 &{} 1 \\ {\hat{u}} -{\tilde{c}} &{} {\hat{u}} &{} 0 &{} {\hat{u}}+{\tilde{c}} \\ {\hat{v}} &{} 0 &{} {\tilde{c}} &{} {\hat{v}} \\ \widehat{1/\rho } &{} -\frac{\overline{h}}{\overline{\rho h}} &{} 0 &{} \widehat{1/\rho }, \end{pmatrix} , \mathbf{R }^{c,B} = \begin{pmatrix} 1 &{} 0 &{} 1 &{} 1 \\ {\hat{u}} &{} {\tilde{c}} &{} 0 &{} {\hat{u}} \\ {\hat{v}} - {\tilde{c}} &{} 0 &{} {\hat{v}} &{} {\hat{v}}+{\tilde{c}} \\ \widehat{1/\rho } &{} 0 &{} -\frac{\overline{h}}{\overline{\rho h}} &{} \widehat{1/\rho } \end{pmatrix}. \end{aligned}$$

(61)

Finally, the amplitudes in the decomposition of both the flux gradient and source discretization in each direction are:

$$\begin{aligned} \begin{array}{lllllll} \alpha _1^{c,A} &{} = &{} \frac{({\hat{u}}+{\tilde{c}})\Delta _x (\rho h)-\Delta _x (\rho h u)}{2{\tilde{c}}}+\frac{g\overline{\rho h}}{4 {\tilde{c}}^2}(\Delta _x h-\widehat{1/\rho } \Delta _x (\rho h)), \\ \alpha _4^{c,A} &{} = &{} \frac{(-{\hat{u}}+{\tilde{c}})\Delta _x (\rho h)+\Delta _x (\rho hu)}{2{\tilde{c}}}+\frac{g\overline{\rho h}}{4{\tilde{c}}^2}(\Delta _x h-\widehat{1/\rho } \Delta _x (\rho h)), \\ \\ \alpha _1^{c,B} &{} = &{} \frac{({\hat{v}}+{\tilde{c}})\Delta _y (\rho h)-\Delta _y (\rho hv)}{2{\tilde{c}}}+\frac{g\overline{\rho h}}{4{\tilde{c}}^2}(\Delta _y h-\widehat{1/\rho } \Delta _y (\rho h)), \\ \alpha _4^{c,B} &{} = &{} \frac{(-{\hat{v}}+{\tilde{c}})\Delta _y (\rho h)+\Delta _y (\rho hv)}{2{\tilde{c}}}+\frac{g\overline{\rho h}}{4{\tilde{c}}^2} \left( \Delta _y h -\widehat{1/\rho } \Delta _y (\rho h) \right) , \\ \\ \alpha _2^{c,A} &{} = &{} - \frac{g\overline{\rho h}}{2{\tilde{c}}^2} \left( \Delta _x h - \widehat{1/\rho } \Delta _x (\rho h) \right) ,\\ \alpha _3^{c,A} &{} = &{} \frac{\Delta _x (\rho hv)-\frac{g{\hat{v}}}{2{\tilde{c}}^2} \Delta _x (\rho h^2)}{{\tilde{c}}}, \\ \\ \alpha _2^{c,B} &{} = &{} \frac{\Delta _y (\rho hu)-\frac{g{\hat{u}}}{2{\tilde{c}}^2}\Delta _y (\rho h^2)}{{\tilde{c}}}, \\ \alpha _3^{c,B} &{} = &{} -\frac{g\overline{\rho h}}{2{\tilde{c}}^2} \left( \Delta _y h -\widehat{1/\rho } \Delta _y (\rho h) \right) \end{array} \end{aligned}$$

(62)

$$\begin{aligned}&\beta _1^{c,A} = \frac{\bar{\rho }{\tilde{c}}}{2} \Delta _x B, \beta _2^{c,A} = \beta _3 ^{c,A} = 0, \beta _4^{c,A} = -\frac{\bar{\rho }{\tilde{c}}}{2}\Delta _x B , \nonumber \\&\beta _1^{c,B} = \frac{\bar{\rho }{\tilde{c}}}{2}\Delta _y B, \; \beta _2^{c,B} = \beta _3^{c,B} = 0, \; \beta _4^{c,B} = -\frac{\bar{\rho }{\tilde{c}}}{2}\Delta _y B.\nonumber \\ \end{aligned}$$

(63)

Here $\Delta _x(\cdot )$ and $\Delta _y(\cdot )$ refer to the horizontal (east-west) and vertical (north and south) differences respectively. We note that when $\rho$ is constant the second terms in $\alpha _1,\alpha _2$ vanish, as well as $\alpha _2^{c,A}, \alpha _3^{c,B}$, recovering the regular scheme for the 2-D shallow water equations.

Primitive Formulation

Using the primitive variables, $\mathbf{W }^{nc} = (u, \; v, \; p, \; \rho )^T$, the 2-D non-conservative primitive system reads

$$\begin{aligned} \begin{pmatrix} u \\ v \\ p \\ \rho \end{pmatrix}_t + \begin{pmatrix} u &{} 0 &{} (\rho h)^{-1} &{} 0 \\ 0 &{} u &{} 0 &{} 0 \\ \rho h c^2 &{} 0 &{} u &{} 0 \\ 0 &{} 0 &{} 0 &{} u \end{pmatrix} \begin{pmatrix} u \\ v \\ p \\ \rho \end{pmatrix}_x + \begin{pmatrix} v &{} 0 &{} 0 &{} 0 \\ 0 &{} v &{} (\rho h)^{-1} &{} 0 \\ 0 &{} \rho h c^2 &{} v &{} 0 \\ 0 &{} 0 &{} 0 &{} v \end{pmatrix} \begin{pmatrix} u \\ v \\ p \\ \rho \end{pmatrix}_y = \begin{pmatrix} -g B_x \\ -g B_y \\ 0 \\ 0 \end{pmatrix}. \end{aligned}$$

(64)

Any consistent discretization of the above system automatically preserves u, v and p across contact discontinuities (in a flat topography). For consistency, we choose the same averages here as in the conservative formulation ${\hat{u}}, \overline{\rho h}, \widehat{1/\rho }, {\tilde{c}}$. The matrices of eigenvectors are

$$\begin{aligned} \mathbf{R }^{nc,A} = \begin{pmatrix} {\tilde{c}} &{} 0 &{} 0 &{} {\tilde{c}} \\ 0 &{} {\tilde{c}} &{} 0 &{} 0 \\ -\overline{\rho h} {\tilde{c}}^2 &{} 0 &{} 0 &{} \overline{\rho h} {\tilde{c}}^2 \\ 0 &{} 0 &{} 1 &{} 0, \end{pmatrix} , \mathbf{R }^{nc,B} = \begin{pmatrix} 0 &{} {\tilde{c}} &{} 0 &{} 0 \\ {\tilde{c}} &{} 0 &{} 0 &{} {\tilde{c}} \\ -\overline{\rho h} {\tilde{c}}^2 &{} 0 &{} 0 &{} \overline{\rho h} {\tilde{c}}^2 \\ 0 &{} 0 &{} 1 &{} 0 \end{pmatrix}. \end{aligned}$$

(65)

Finally, the amplitudes in the decomposition of both the flux gradient and source discretization in each direction are:

$$\begin{aligned} \begin{array}{llllllllllllllllll} \alpha _1^{nc,A} &{} = &{} \frac{\overline{\rho h} {\tilde{c}} \Delta _x u - \Delta _x p}{2\overline{\rho h} {\tilde{c}}^2}, &{} &{} \beta _1^{nc,A} &{} = &{} -\frac{g \Delta _x B}{2 {\tilde{c}}} , &{}&{} \alpha _1^{nc,B} &{} = &{} \frac{\overline{\rho h} {\tilde{c}} \Delta _y v-\Delta _yp}{2\overline{\rho h} {\tilde{c}}^2}, &{} &{} \beta _1^{nc,B} &{} = &{} -\frac{g\Delta _y B }{2 {\tilde{c}}} \\ \alpha _2^{nc,A} &{} = &{} \frac{\Delta _x v}{{\tilde{c}}},&{} &{} \beta _2^{nc,A} &{} = &{} 0 , &{} &{} \alpha _2^{nc,B} &{} = &{} \frac{\Delta _y u}{{\tilde{c}}}, &{} &{} \beta _2^{nc,B} &{} = &{} 0 \\ \alpha _3^{nc,A} &{} = &{} \Delta _x \rho , &{} &{} \beta _3^{nc,A} &{} = &{} 0, &{} &{} \alpha _3^{nc,B} &{} = &{} \Delta _y \rho , &{} &{} \beta _3^{nc,B} &{} = &{} 0 \\ \alpha _4^{nc,A} &{} = &{} \frac{\overline{\rho h} {\tilde{c}} \Delta _x u +\Delta _x(\rho h)}{2 \overline{\rho h} {\tilde{c}}^2}, &{} &{} \beta _4^{nc,A} &{} = &{} -\frac{g\Delta _x B}{2{\tilde{c}}}, &{} &{} \alpha _4^{nc,B} &{} = &{} \frac{\overline{\rho h} {\tilde{c}}\Delta _y v + \Delta _y p}{2\overline{\rho h} {\tilde{c}}^2}, &{} &{} \beta _4^{nc,B} &{} = &{} -\frac{g\Delta _y B}{2{\tilde{c}}}. \end{array} \end{aligned}$$

(66)

1.1 Appendix A.1: Summary of the 2-D Numerical Scheme

The 2-D numerical scheme is summarized as follows. Each cel (i, j) is updated as

$$\begin{aligned} \begin{array}{lll} \mathbf{W }_{i,j}^{n+1} &{} = &{} \mathbf{W }_{i,j}^n-\dfrac{\Delta t}{\Delta x} \left\{ A^+_{i-\frac{1}{2},j}\left( \mathbf{W }_{i,j}^n-\mathbf{W }_{i-1,j}^n\right) - \Delta x\; \mathbf{S }_{A,i-\frac{1}{2},j}^+ + A^-_{i+\frac{1}{2},j}\left( \mathbf{W }_{i+1,j}^n-\mathbf{W }_{i,j}^n\right) - \Delta x\; \mathbf{S }_{A,i+\frac{1}{2},j}^- \right\} \\ &{} &{} -\dfrac{\Delta t}{\Delta y} \left\{ B^+_{i,j-\frac{1}{2}}\left( \mathbf{W }_{i,j}^n-\mathbf{W }_{i,j-1}^n\right) -\Delta y\; \mathbf{S }_{B,i,j-\frac{1}{2}}^+ + B^-_{i,j+\frac{1}{2}}\left( \mathbf{W }_{i,j+1}^n-\mathbf{W }_{i,j}^n\right) -\Delta y\; \mathbf{S }_{B,i,j+\frac{1}{2}}^- \right\} . \end{array} \end{aligned}$$

Here,

$$\begin{aligned} \begin{array}{lll} A^+\Delta _x \mathbf{W } -\Delta x\; \mathbf{S }_A^+= \displaystyle {\sum _{\lambda _k^A> 0}} (\alpha _k^A \lambda _k^A - \beta _k^A) \mathbf{r }_k^A ~, \quad A^-\Delta _x \mathbf{W } - \Delta x\; \mathbf{S }_A^- = \displaystyle {\sum _{\lambda _k^A \le 0} } (\alpha _k^A \lambda _k^A - \beta _k^A) \mathbf{r }_k^A, \\ B^+\Delta _y \mathbf{W } -\Delta y\; \mathbf{S }_B^+ = \displaystyle {\sum _{\lambda _k^B > 0}} (\alpha _k^B \lambda _k^B - \beta _k^B) \mathbf{r }_k^B ~, \quad B^-\Delta _y \mathbf{W } -\Delta y\; \mathbf{S }_B^-= \displaystyle {\sum _{\lambda _k^B \le 0} } (\alpha _k^B \lambda _k^B - \beta _k^B) \mathbf{r }_k^B. \end{array} \end{aligned}$$

If the jump in density near cell (i, j) does not exceed a threshold, i.e., if $\max _{i-1 \le i' \le i+1,}$ ${j-1 \le ij \le j+1}( |\rho _{i'+1,j'+1}-\rho _{i',j'}| ) < \rho _o$, then

$$\begin{aligned} W = \begin{pmatrix} \rho h \\ \rho h u \\ \rho h v\\ h \end{pmatrix}, \mathbf{r }_1^A = \begin{pmatrix} 1\\ {\hat{u}} -{\tilde{c}} \\ {\hat{v}} \\ \widehat{1/\rho } \end{pmatrix} \mathbf{r }_2^A = \begin{pmatrix} 1\\ {\hat{u}} \\ 0 \\ -\frac{\overline{h}}{\overline{\rho h}} \end{pmatrix} \mathbf{r }_3^A = \begin{pmatrix} 0 \\ 0 \\ {\tilde{c}} \\ 0 \end{pmatrix} \mathbf{r }_4^A = \begin{pmatrix} 1\\ {\hat{u}}+{\tilde{c}} \\ {\hat{v}} \\ \widehat{1/\rho } \end{pmatrix},\\ \lambda _1^A = {\hat{u}}-{\tilde{c}}, \lambda _2^A ={\hat{u}}, \lambda _3^A = {\hat{u}}, \lambda _4^A ={\hat{u}}+{\tilde{c}}, \end{aligned}$$

$$\begin{aligned} \mathbf{r }_1^B= \begin{pmatrix} 1\\ {\hat{u}} \\ {\hat{v}} - {\tilde{c}} \\ \widehat{1/\rho } \end{pmatrix}, \mathbf{r }_2^B= \begin{pmatrix} 0 \\ {\tilde{c}} \\ 0 \\ 0 \end{pmatrix}, \mathbf{r }_3^B= \begin{pmatrix} 1\\ 0 \\ {\hat{v}} \\ -\frac{\overline{h}}{\overline{\rho h}} \end{pmatrix}, \mathbf{r }_4^B= \begin{pmatrix} 1\\ {\hat{u}} \\ {\hat{v}}+{\tilde{c}} \\ \widehat{1/\rho } \end{pmatrix},\\ \lambda _1^B = {\hat{v}}-{\tilde{c}}, \lambda _2^B = {\hat{v}}, \lambda _3^B = {\hat{v}}, \lambda _4^B = {\hat{v}}+{\tilde{c}}, \end{aligned}$$

$$\begin{aligned} \begin{array}{lllllll} \alpha _1^{A} &{} = &{} \frac{({\hat{u}}+{\tilde{c}})\Delta _x (\rho h)-\Delta _x (\rho h u)}{2{\tilde{c}}}+\frac{g\overline{\rho h}}{4 {\tilde{c}}^2}(\Delta _x h-\widehat{1/\rho } \Delta _x (\rho h)), \\ \alpha _4^{A} &{} = &{} \frac{(-{\hat{u}}+{\tilde{c}})\Delta _x (\rho h)+\Delta _x (\rho hu)}{2{\tilde{c}}}+\frac{g\overline{\rho h}}{4{\tilde{c}}^2}(\Delta _x h-\widehat{1/\rho } \Delta _x (\rho h)), \\ \\ \alpha _1^{B} &{} = &{} \frac{({\hat{v}}+{\tilde{c}})\Delta _y (\rho h)-\Delta _y (\rho hv)}{2{\tilde{c}}}+\frac{g\overline{\rho h}}{4{\tilde{c}}^2}(\Delta _y h-\widehat{1/\rho } \Delta _y (\rho h)), \\ \alpha _4^{B} &{} = &{} \frac{(-{\hat{v}}+{\tilde{c}})\Delta _y (\rho h)+\Delta _y (\rho hv)}{2{\tilde{c}}}+\frac{g\overline{\rho h}}{4{\tilde{c}}^2} \left( \Delta _y h -\widehat{1/\rho } \Delta _y (\rho h) \right) , \\ \\ \alpha _2^{A} &{} = &{} - \frac{g\overline{\rho h}}{2{\tilde{c}}^2} \left( \Delta _x h - \widehat{1/\rho } \Delta _x (\rho h) \right) , \\ \alpha _3^{A} &{} = &{} \frac{\Delta _x (\rho hv)-\frac{g{\hat{v}}}{2{\tilde{c}}^2} \Delta _x (\rho h^2)}{{\tilde{c}}}, \\ \\ \alpha _2^{B} &{} = &{} \frac{\Delta _y (\rho hu)-\frac{g{\hat{u}}}{2{\tilde{c}}^2}\Delta _y (\rho h^2)}{{\tilde{c}}}, \\ \alpha _3^{B} &{} = &{} -\frac{g\overline{\rho h}}{2{\tilde{c}}^2} \left( \Delta _y h -\widehat{1/\rho } \Delta _y (\rho h) \right) \end{array} \end{aligned}$$

$$\begin{aligned} \beta _1^{A} = \frac{\bar{\rho }{\tilde{c}}}{2} \Delta _x B, \beta _2^{A} = \beta _3 ^{A} = 0, \beta _4^{A} = -\frac{\bar{\rho }{\tilde{c}}}{2}\Delta _x B , \\ \beta _1^{B} = \frac{\bar{\rho }{\tilde{c}}}{2}\Delta _y B, \; \beta _2^{B} = \beta _3^{B} = 0, \; \beta _4^{B} = -\frac{\bar{\rho }{\tilde{c}}}{2}\Delta _y B. \end{aligned}$$

If the jump in the density exceed the chosen threshold ( $\max _{i-1 \le i' \le i+1,j-1 \le ij \le j+1}( |\rho _{i'+1,}$ ${j'+1}-\rho _{i',j'}| ) \ge \rho _o$ ), then

$$\begin{aligned} W = \begin{pmatrix} u \\ v \\ p \\ \rho \end{pmatrix}, r_1^A= \begin{pmatrix} {\tilde{c}} \\ 0 \\ -\overline{\rho h} {\tilde{c}}^2 \\ 0 \end{pmatrix} r_2^A= \begin{pmatrix} 0 \\ {\tilde{c}} \\ 0 \\ 0 \end{pmatrix} r_3^A= \begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \end{pmatrix} r_4^A= \begin{pmatrix} {\tilde{c}} \\ 0 \\ \overline{\rho h} {\tilde{c}}^2 \\ 0 \end{pmatrix},\\ \lambda _1^A = {\hat{u}}-{\tilde{c}}, \lambda _2^A ={\hat{u}}, \lambda _3^A = {\hat{u}}, \lambda _4^A ={\hat{u}}+{\tilde{c}}, \end{aligned}$$

$$\begin{aligned} \mathbf{r }_1^B= \begin{pmatrix} 0 \\ {\tilde{c}} \\ -\overline{\rho h} {\tilde{c}}^2 \\ 0 \end{pmatrix}, \mathbf{r }_2^B= \begin{pmatrix} {\tilde{c}} \\ 0 \\ 0 \\ 0 \end{pmatrix}, \mathbf{r }_3^B= \begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \end{pmatrix}, \mathbf{r }_4^B= \begin{pmatrix} 0 \\ {\tilde{c}} \\ \overline{\rho h} {\tilde{c}}^2 \\ 0 \end{pmatrix},\\ \lambda _1^B = {\hat{v}}-{\tilde{c}}, \lambda _2^B = {\hat{v}}, \lambda _3^B = {\hat{v}}, \lambda _4^B = {\hat{v}}+{\tilde{c}}, \end{aligned}$$

and

$$\begin{aligned} \begin{array}{llllllllllllllllll} \alpha _1^{A} &{} = &{} \frac{\overline{\rho h} {\tilde{c}} \Delta _x u - \Delta _x p}{2\overline{\rho h} {\tilde{c}}^2}, &{} &{} \beta _1^{A} &{} = &{} -\frac{g \Delta _x B}{2 {\tilde{c}}} , &{}&{} \alpha _1^{B} &{} = &{} \frac{\overline{\rho h} {\tilde{c}} \Delta _y v-\Delta _yp}{2\overline{\rho h} {\tilde{c}}^2}, &{} &{} \beta _1^{B} &{} = &{} -\frac{g\Delta _y B }{2 {\tilde{c}}} \\ \\ \alpha _2^{A} &{} = &{} \frac{\Delta _x v}{{\tilde{c}}},&{} &{} \beta _2^{A} &{} = &{} 0 , &{} &{} \alpha _2^{B} &{} = &{} \frac{\Delta _y u}{{\tilde{c}}}, &{} &{} \beta _2^{B} &{} = &{} 0 \\ \\ \alpha _3^{A} &{} = &{} \Delta _x \rho , &{} &{} \beta _3^{A} &{} = &{} 0, &{} &{} \alpha _3^{B} &{} = &{} \Delta _y \rho , &{} &{} \beta _3^{B} &{} = &{} 0 \\ \\ \alpha _4^{A} &{} = &{} \frac{\overline{\rho h} {\tilde{c}} \Delta _x u +\Delta _x(\rho h)}{2 \overline{\rho h} {\tilde{c}}^2}, &{} &{} \beta _4^{A} &{} = &{} -\frac{g\Delta _x B}{2{\tilde{c}}}, &{} &{} \alpha _4^{B} &{} = &{} \frac{\overline{\rho h} {\tilde{c}}\Delta _y v + \Delta _y p}{2\overline{\rho h} {\tilde{c}}^2}, &{} &{} \beta _4^{B} &{} = &{} -\frac{g\Delta _y B}{2{\tilde{c}}}. \end{array} \end{aligned}$$

(67)

Here the threshold and averages are given by Eqs. (28), (40), (46). In addition, the second order extension and entropy fix is conducted as in [24].

Appendix B: Derivation of the Model

We start with the 3D Euler equations

$$\begin{aligned} \partial _t \begin{pmatrix} \rho \\ \rho u \\ \rho v\\ \rho w \end{pmatrix} + \partial _x \begin{pmatrix} \rho u\\ \rho u^2+p \\ \rho uv \\ \rho u w \end{pmatrix} + \partial _y \begin{pmatrix} \rho v\\ \rho uv\\ \rho v^2 +p\\ \rho vw \end{pmatrix} +\partial _z \begin{pmatrix} \rho w\\ \rho uw\\ \rho v w\\ \rho w^2+p \end{pmatrix} = \mathbf{0 }, \end{aligned}$$

(68)

and integrate in the vertical direction from the bottom topography B(x, y) to the surface $B(x,y)+h(x,y,t)$, where h is the depth of water, using the relation

$$\begin{aligned} \int _B^{B+h} \partial _r f dz= \partial _r \int _B^{B+h} f dz-f_{\big |_{B+h}}\partial _r (B+h)+f_{\big |_B}\partial _r B, \end{aligned}$$

where r is any of the spatial variables x or y, and f is any of the conserved variables $\rho ,\rho u,\rho v$ or $\rho w$. For any conserved variable f, we also define the vertically averaged quantity

$$\begin{aligned} {\bar{f}} (x,y)= \frac{1}{h} \int _B^{B+h} f(x,y,z)dz. \end{aligned}$$

In addition, we assume a hydrostatic balance so that the pressure is given as $p(x,y,z,t) = g \int _z^{B+h} \rho (x,y,z') dz'$. If the density is approximately uniform in the vertical direction, the pressure is approximated as

$$\begin{aligned} p \approx g \bar{\rho }(x,y) (h+B-z). \end{aligned}$$

(69)

The vertically integrated pressure is then approximated as $g \rho h^2/2$, and denoted by p. Vertically integrating the first three components of Eq. (68) we get

$$\begin{aligned} \begin{array}{llllll} \partial _t (h \bar{\rho }) &{} + &{} \partial _x (h \overline{\rho u})+\partial _y (h \overline{\rho v}) - \rho _{{\big |}_{B+h}} S_{{\big |}_{B+h}} + \rho _{\big |_B} S_{\big |_B} = 0, \\ \\ \partial _t (h \overline{\rho u}) &{} + &{} \partial _x(h\overline{\rho u^2}+h\bar{p}) + \partial _y (h \overline{\rho uv}) - (\rho u)_{\big |_{B+h}} S_{{\big |}_{B+h}} + (\rho u)_{\big |_{B}} S_{\big |_B} = -g p_{\big |_B} \partial _x B,\\ \\ \partial _t (h \overline{\rho v}) &{} + &{} \partial _x (h \overline{\rho u v}) + \partial _y (h \overline{\rho v^2}+h \bar{p}) - (\rho v)_{\big |_{B+h}} S_{\big |_{B+h}} + (\rho v)_{\big |_B} S_{\big |_B} = - g p_{\big |_B} \partial _y B. \end{array} \end{aligned}$$

(70)

where

$$\begin{aligned} S_{B+h} = \partial _t(B+h) + u_{\big |_{B+h}} \partial _x (B+h)+v_{\big |_{B+h}} \partial _y (B+h)-w_{\big |_{B+h}} , \\ S_B = \partial _t B + u_{\big |_{B}} \partial _x B+v_{\big |_{B}} \partial _y B-w_{\big |_{B}}. \end{aligned}$$

The bottom and surface are assumed to be streamlines so that $S_{\big |_{B+h}}=S_{\big |_B} =0$ vanish. Assuming that the flow is shallow, we replace $\overline{\rho u}\approx \bar{\rho }\bar{u}, \overline{\rho u^2} \approx \bar{\rho }\bar{u}^2$ and so on, and drop the bars.

The system can be closed by assuming that the density $\rho$ is advected by the flow $\partial _t \rho + u \partial _x \rho +v \partial _y \rho =0$, and combining it with the first equation we get system (1).