A Non-oscillatory Central Scheme for One-Dimensional Two-Layer Shallow Water Flows along Channels with Varying Width

Abstract

We present a new high-resolution, non-oscillatory semi-discrete central scheme for one-dimensional two-layer shallow-water flows along channels with non-uniform rectangular cross sections and bottom topography. The scheme extends existing central semi-discrete schemes for hyperbolic conservation laws and it enjoys two properties crucial for the accurate simulation of shallow-water flows: it preserves the positivity of the water height, and it is well balanced, i.e., the source terms arising from the geometry of the channel are discretized so as to balance the non-linear hyperbolic flux gradients. Along with a detailed description of the scheme and proofs of these two properties, we present several numerical experiments that demonstrate the robustness of the numerical algorithm.

Appendices

Appendix A: Proof of Theorem 1

Proof

(i) Starting with the initial conditions u _i (x)=0 (i=1,2), w ₁(x)=W ₁ and h ₂(x)=H ₂ for all x, we fix a space scale Δx and the corresponding partition of the solution domain, $\{ I_{j} \} := \{ [x_{j - \frac{1}{2}}, x_{j + \frac{1}{2}}] \}$. We then define the cell averages of the conserved quantities in the modified shallow-water model, (15a)–(15d), as

$$ \overline{A}^T_{1,j} := \overline{ \sigma}_j \overline{w}_{1,j} = \overline{\sigma }_j W_1 \quad\mbox{and}\quad A_{2,j} := \overline{\sigma}_j \overline{h}_{2,j}= \overline{\sigma}_j H_2 $$

(86)

and

(87a)

(87b)

The reconstructed point values of w ₁ and h ₂, clearly satisfy $w_{1, j \pm\frac{1}{2}}^{\pm} = W_{1}$, $h_{2,j\pm\frac {1}{2}} = H_{2}$, and those of Q _i , $Q_{j \pm\frac{1}{2}}^{\pm} = 0$ (i=1,2), thus the values of the water heights (38) (and the corresponding values of h ₂ recovered via the minmod reconstruction), and those of the bottom topography at the cell interfaces, $B_{j \pm \frac{1}{2}}^{\mp}$, satisfy

$$ h_{1,j + \frac{1}{2}}^{-} - h_{1,j - \frac{1}{2}}^{+} = - \bigl( B_{j + \frac{1}{2} }^{-} - B_{j - \frac{1}{2}}^{+} \bigr) \quad\mbox{and}\quad h_{2,j + \frac{1}{2}}^{-} - h_{2,j - \frac{1}{2}}^{+} = 0 $$

(88)

In view of this, the first and third components of the numerical fluxes $H_{j + \frac{1}{2}}$ in (22) read

(89a)

(89b)

That is,

(90a)

(90b)

which allows us to recover $\overline{w}_{1,j} (t + \Delta t)= W_{1}$ and h _2,j (t+Δt)=H ₂ exactly from (86).

Noting that, according to (38), $h_{1,j+\frac{1}{2}}^{+} = h_{1,j+\frac{1}{2}}^{-} =: h_{1,j+\frac{1}{2}}$ and that the minmod reconstruction guarantees $h_{2,j+\frac{1}{2}}^{+} = h_{2,j+\frac {1}{2}}^{-} =: h_{2,j+\frac{1}{2}}$, the second and fourth components of the numerical flux amount, respectively, to (40) and (41), and since (88) hold, they are balanced by (48)–(51), therefore

$$ \frac{d}{d t} \overline{Q}_{i,j} (t) = 0 \quad\Rightarrow\quad \overline{Q}_{i,j} (t + \Delta t) = \overline{Q}_{i,j} (t) = 0,\quad i = 1,2. $$

(91)

and u _i,j (t+Δt)≡0 (i=1,2) are also recovered exactly.

(ii) We begin by writing explicitly the cell average $A^{T}_{1,j}(t + \Delta t)$ when the system (22) is evolved with forward Euler’s ODE solver,

$$ \overline{A}^T_{1,j}(t + \Delta t) = \overline{A}^T_{1,j}(t) - \lambda \bigl[ H_{j + \frac{1}{2}}^{(1)}(t) - H_{j - \frac{1}{2}}^{(1)}(t) \bigr], $$

(92)

where λ=Δt/Δx. This amounts to

(93)

(where all the terms on the right hand side are understood to be evaluated at time t). Using $Q_{1,j \pm\frac{1}{2}}^{\pm} = \sigma_{j \pm\frac{1}{2}}^{\pm} h_{1,j \pm\frac{1}{2}}^{\pm} u_{1,j \pm \frac{1}{2}}^{\pm }$, we write

(94)

The terms involving $\sigma_{j \pm\frac{1}{2}} B_{j \pm\frac {1}{2}}^{\pm}$ and $\sigma_{j \pm\frac{1}{2}} B_{j \pm\frac{1}{2}}^{\mp}$ on the right hand side cancel, and since $a_{j \pm\frac{1}{2}}\sigma_{j \pm\frac {1}{2}} \geq |u_{1,j \pm\frac{1}{2}}^{\pm}| \sigma_{j \pm\frac{1}{2}}^{\pm }$, (58), and $h_{1,j \pm\frac{1}{2}}^{\mp} \geq0$, the CFL restriction (63) allows us to write

$$ \overline{A}^T_{1,j}(t + \Delta t) \geq \overline{A}^T_{1,j}(t) - \frac{1}{2} \overline {\sigma}_j \bigl(h_{1,j + \frac{1}{2}}^- + h_{1,j - \frac{1}{2}}^+ \bigr) = \overline{\sigma}_j \overline{B}_j, $$

(95)

from where (64) follows. □

Appendix B: On the Semi-Discrete Central Formulation for Hyperbolic Systems with Source Terms

The numerical scheme that we have presented in this work is based on the high-order semi-discrete central formulation for hyperbolic conservation laws first introduced in [17]. Here we outline the derivation of the equivalent semi-discrete formulation for hyperbolic systems with source terms,

$$ v_t + f(v)_x = S(v,x) $$

(96)

Before we proceed, we should clarify that this derivation consists of three steps: staggered evolution of cell averages, reprojection onto the original non-staggered grid, and the evaluation of the limit as Δt→0 of the finite difference $(\overline {v}^{n+1}_{j} - \overline{v}^{n}_{j}) /\Delta t$. A process that requires the formulation of several intermediate solutions presented below, however, for the actual implementation of the semi-discrete scheme one doesn’t need to calculate these solutions; the intermediate solutions are only introduced formally so as to form the finite difference and take the limit.

Starting with the cell averages $\{ \overline{v}^{n}_{j} \}_{j}$ over the partition $\{ I_{j} \}_{j} = \{ [x_{j -\frac{1}{2}}, x_{j + \frac{1}{2}}] \}_{j}$ at time t=t ⁿ, we use the estimates of the maximum speed of propagation, $a_{j + \frac{1}{2}}^{n}$, to define

$$ x_{j \pm\frac{1}{2},l}^n := x_{j \pm\frac{1}{2}} - a_{j \pm\frac {1}{2}}^n \Delta t \quad\mbox{and}\quad x_{j \pm\frac{1}{2},r}^n := x_{j \pm\frac{1}{2}} + a_{j \pm\frac{1}{2}}^n \Delta t, $$

(97)

and with these values we repartition the computational domain into two sets of cells: the first set,

$$\{ \tilde{I}_{j \pm\frac{1}{2}} \}_j = \bigl\{ \bigl[x_{j \pm\frac {1}{2},l}^n, x_{j \pm \frac{1}{2},r}^n\bigr] \bigr\}_j, $$

contains the neighborhood of the interfaces $x_{j \pm\frac{1}{2}}$ within which discontinuous solutions propagate, the second set,

$$\{ \tilde{I_{j}} \}_j = \bigl\{ \bigl[x_{j - \frac{1}{2},r}^n, x_{j + \frac {1}{2},l}^n\bigr] \bigr\}_j, $$

contains the portion of the original cells where the solution remains smooth over the interval [t ⁿ,t ⁿ⁺¹] (consult Fig. 12). We then formulate two sets of staggered cell averages following the recipe in [11] (a non-staggered extension of the original NT scheme proposed in [22]) including, in our case, the contribution from the source term. These two sets of solutions read:

(98)

(99)

where $\lambda= \frac{\Delta t}{\Delta x}$. The midpoint values $v^{n + \frac{1}{2}}(x)$ are approximated using Taylor’s theorem and the balance law (96),

$$ v^{n+\frac{1}{2}}(x) \approx v^n(x) + \frac{\Delta t}{2} v_t^n(x) \approx v^n(x) - \frac{\Delta t}{2} \bigl[ f\bigl(v^n(x)\bigr)_x - S \bigl(v^n(x), x\bigr) \bigr], $$

(100)

the numerical derivatives $(u_{x})_{j}^{n}$ via a minmod limiter,

$$ (v_x)_j^n = \frac{1}{\Delta x} \mathrm{minmod} \bigl( \overline {v}^n_{j + 1} - \overline{v}^n_j, \overline{v}^n_j - \overline{v}^n_{j - 1} \bigr), $$

(101)

and the flux derivatives can be calculated as f(v ⁿ(x)) _x =A(v)v _x , where $A(v) = \frac{df}{dv}$ is the Jacobian matrix of f, or component by component as

$$ f\bigl(v^n(x)\bigr)_x = \frac{1}{\Delta x} \mathrm{minmod} \bigl( f\bigl(v^n(x + \Delta x)\bigr) - f \bigl(v^n(x)\bigr), f\bigl(v^n(x)\bigr) - f \bigl(v^n(x - \Delta x)\bigr) \bigr). $$

(102)

A well-balance discretization of the integrals on the right hand side of (98) can be found in [24].

Fig. 12

Modified staggered central differencing and evolution

From these staggered averages, we form the interpolant

$$ R\bigl(x,t^{n+1}\bigr) = \sum _j \bigl\{ \bigl[\overline{\omega}_{j + \frac {1}{2}}^{n + 1} + (v_x)_{j + \frac{1}{2}}^{n + 1} (x - x_{j + \frac{1}{2}})\bigr] { \mathbf{1}}_{\bigl[x_{j+\frac{1}{2}, l}^n, x_{j+\frac{1}{2}, r}^n\bigr]} + \overline {\omega }_{j}^{n + 1} { \mathbf{1}}_{\bigl[x_{j-\frac{1}{2}, r}^n, x_{j+\frac{1}{2}, l}^n\bigr]} \bigr\} $$

(103)

where the derivatives $(v_{x})_{j \pm\frac{1}{2}}^{n + 1}$ are also approximated with a minmod limiter

$$ (v_x)_{j + \frac{1}{2}}^{n + 1} = \frac{2}{\Delta x} \mathrm {minmod} \biggl(\frac {\overline{\omega}_{j + 1}^{n+1} - \overline{\omega}_{j + \frac{1}{2} }^{n + 1}}{1 + \lambda(a_{j + \frac{1}{2}}^n - a_{j + \frac{3}{2}}^n)}, \frac{\overline{\omega}_{j + \frac{1}{2}}^{n+1} - \overline{\omega }_{j}^{n +1}}{1 + \lambda(a_{j - \frac{1}{2}}^n - a_{j + \frac {1}{2}}^n)} \biggr). $$

(104)

This interpolant is then reprojected onto the original (non-staggered) grid, obtaining the updated, non-staggered cell averages

(105)

In order to arrive at the semi-discrete formulation from (105), we consider the limit

$$ \frac{d}{dt} \overline{v}_j(t) = \lim_{\Delta t \rightarrow0} \frac {\overline{v}_j^{n+1} - \overline{v}_j^{n}}{\Delta t}. $$

(106)

To evaluate it, we write explicitly the finite difference using (105),

(107)

and expand it making use of (98) and (99)

(108)

In the limit Δt→0, the midpoint interface values approach

(109)

and we obtain (22)–(23)

(110)

with $H_{j \pm\frac{1}{2}}$ given by (19) as desired.