Open Water Flow in a Wet/Dry Multiply-Connected Channel Network: A Robust Numerical Modeling Algorithm

Abstract

Our goal was to develop a robust algorithm for numerical simulation of one-dimensional shallow water flow in a complex multiply-connected channel network with arbitrary geometry and variable topography. We apply a central-upwind scheme with a novel reconstruction of the open water surface in partially flooded cells that does not require additional correction. The proposed reconstruction and an exact integration of source terms for the momentum conservation equation provide positivity preserving and well-balanced features of the scheme for various wet/dry states. We use two models based on the continuity equation and mass and momentum conservation equations integrated for a control volume around the channel junction to its treatment. These junction models permit to simulate subcritical and supercritical flows in a channel network. Numerous numerical experiments demonstrate the robustness of the proposed numerical algorithm and a good agreement of numerical results with exact solutions, experimental data, and results of the previous numerical studies. The proposed new specialized test on inundation and drying of an initially dry channel network shows the merits of the new numerical algorithm to simulate the subcritical/supercritical open water flows in the networks.

1 Introduction

The Saint–Venant equations are the common basis for the 1D modeling of the open flow in river and channel networks. Among different numerical methods in the various computer codes used in hydraulic engineering [e.g., HEC-RAS (Brunner 2016), ISIS 1D (CH2M HILL 2016), FEQ (Franz and Melching 1997), CHARIMA (Holly et al. 1990), NETSTARS (Lee and Hsieh 2003), FLDWAV (Fread and Lewis 1998)], the Preissmann (1961) scheme is one of the most widely used for the solution of one-dimensional unsteady open-flow channel problem. The Preissmann scheme is a weighted four-point nonlinear implicit finite-difference scheme. The advantage of this scheme is the efficient double-sweep algorithm for the numerical solution of the Newton–Raphson linearization of the Saint–Venant equations in a multiply-connected channel network. One of the disadvantages is that the scheme does not provide positivity preserving of a numerical solution, and therefore it cannot be applied to simulate water flow in drying channels. The positivity of water flow depth is especially important for schemes used for simulations of mountain river flow. Steep bottom slopes can generate artificial “drying” and consequently crash the calculations due to negative values of depth in numerical solutions if the Preissmann scheme or other schemes without positivity preserving are used for modeling. The need to develop robust schemes for modeling mountain river flows was revealed last years during computational efforts to simulate radionuclide transport in river networks of mountain watersheds contaminated after the Fukushima Daiichi nuclear accident (Kurikami et al. 2014; Zheleznyak et al. 2016; Kivva et al. 2018). Moreover, a local convergence of the Newton–Raphson method requires appropriate initial conditions that may be problematic for a large river network (Yu et al. 2017).

Many different schemes for hyperbolic equations that preserve nonnegativity of water flow depth were developed after the pioneering publication of Godunov (1959). We refer the reader to Boris and Book (1973), van Leer (1979), Zalesak (1979), Rodionov (1987), Nessyahu and Tadmor (1990), Cockburn (1989), Kivva (2008), Zhang and Shu (2011), Xing and Shu (2011), Xing (2016), where different approaches are presented. Among them, central (Kurganov and Tadmor 2000) and central-upwind (Kurganov and Levy 2002) schemes are widely used in a wide range of applications due to their simplicity, efficiency, and robustness. For the shallow water flow in open channels, these schemes were applied in Balbás and Karni (2009), Balbás and Hernandez-Duenas (2014), Hernandez-Duenas and Beljadid (2016). Balbás and Karni (2009) presented a central finite-volume scheme of second-order accuracy to simulate one-dimensional shallow water flows along channels with non-uniform rectangular cross-sections and bottom topography. Balbás and Hernandez-Duenas (2014) extended this scheme for channels with cross-sections of arbitrary shape and bottom topography. A central-upwind scheme with artificial viscosity was proposed by Hernandez-Duenas and Beljadid (2016) for shallow water flows in channels with arbitrary geometry and variable topography. Hodges (2019) does a comprehensive review of the Preissmann scheme versus Godunov-like methods for solving the Saint–Venant equations. He notes that there is no consensus on the best method. There are numerous examples of successful implementation of codes based on the Preissmann schemes for large-scale river networks. However, in contrast, with the Preissmann scheme, Godunov type schemes and other finite volume methods for high-resolution modeling could correctly preserve shocks and transcritical flows. The limitation of the Preissmann scheme to simulate the transcritical flows is taken into account in the TELEMAC MASCARET modeling system which library includes the numerical kernel based on the Preissmann scheme for plane rivers without dams, and a finite volume Roe scheme for dam-break simulations (Goutal and Maurel 2002; Goutal et al. 2012). Among other numerical methods for solving the Saint–Venant equations, we mention finite-volume methods that do not use the Godunov approach, finite-difference methods other than the Preissmann scheme, finite element method, finite element/volume method, dual-finite volume method, discontinuous Galerkin methods, residual distribution methods, and Lagrangian particle methods. For a detailed review of these methods, we refer to Yoshioka et al. (2015), Lai and Khan (2012, 2018), Hodges (2019). Yoshioka et al. (2015) note that most of the conventional numerical methods cannot appropriately handle river networks with loops and flow transitions.

Taking into account the needs for a robust computational algorithm for modeling multiply-connected river networks, we implemented a central-upwind scheme (Kurganov and Petrova 2007) for open water flow in multiply-connected channel networks as a practical alternative to the Preissmann scheme. In the framework of second-order central-upwind schemes, we propose a novel reconstruction of water surface elevation that does not require additional correction in partially flooded cells. The basic piecewise linear reconstruction in a central-upwind scheme for partially flooded cells can produce a negative value of water depth on cell boundary. Additional correction of the negative water depths proposed by Kurganov and Petrova (2007) artificially distort the free water surface that leads to the appearance of fictitious water velocities. Under additional correction to eliminate fictitious water velocities Bollermann et al. (2013) suggested to take into account water surface elevation in neighboring flooded computational cells. However, on a coarse grid, a partially flooded cell may not have neighboring flooded cells. Our reconstruction algorithm is a generalization of the reconstruction (Bollermann et al. 2013) and it is applicable for any bottom topography and its flooding.

The semi-discretization form of a central-upwind scheme is obtained by approximation of integral equations for the mass conservation and momentum balance. The integral form of the momentum conservation equation includes terms that account for the forces due to changes in channel width and bed elevation (Cunge et al. 1980). Assuming that the water and bed surface elevations and channel cross-section depend linearly on the water and bed elevations, and that channel cross sections are defined at cell faces, we can exactly calculate the source term integrals on any internal cell interval. The exact values of the source term integrals with our proposed reconstruction will guarantee that the hydrostatic fluxes at a cell are balanced by the source terms that account for the sloping bed and variable channel width for various wetting/drying steady states at rest.

One of the challenges is the treatment of a channel junction in the modeling of open water flow in a multiply-connected channel network. In many used hydraulic engineering codes, it is assumed that two groups of hydraulic conditions are held at a channel junction (Cunge et al. 1980). One is that the sum of the discharges has to be zero at the junction, and a second is that the water levels or energy levels at the ends of linked channels (reaches) are equal at the junction. Such hydrodynamic models as CHARIMA, FEQ, MIKE-11 (DHI 2017), RIVICE (Carson and Sydor 2013) and ONE-D (Environment Canada and Environment 1995) use the equality of water surface elevations at the junction. Others such as HEC-RAS and FEQ are energy-based models. Usually, these hydraulic conditions are specified as interior boundary conditions. In the general case, the assumption of equality of inflow and outflow discharges at a junction does not, in principle, permit to simulate open water flow in a channel network under wetting/drying conditions. Yoshioka et al. (2015, 2016) use as the internal boundary conditions at the junction the continuity equation discretized on a dual mesh around it and the momentum flux at the junction calculated as a weighted linear combination of the momentum fluxes inflowing into it. Sanders et al. (2001) and Bellamoli et al. (2018) simulated open water flow in channel networks by nesting a 2D model at junctions. Contarino et al. (2016) developed the implicit solver for the junction-generalized Riemann problem which was applied to construct a high-order ADER scheme for stiff hyperbolic balance laws in networks. We consider a channel junction as a computational cell that may have more than two inflows and outflows. For such cell, we integrate the continuity and momentum conservation equations over a control volume around the channel junction. We use in simulations two kinds of junction models: (1) based only on continuity equation for subcritical flow, and (2) based on mass and momentum conservation equations for supercritical flow. The first junction model is an analog of the model of the equality of water surface elevations. The junction treatment with the continuity equation also applied in the staggered Abbott-Ionescu type schemes (Cunge et al. 1980).

Constraints on the time step under which the scheme is positivity preserving may be very restrictive since they depend on the ratio of wetted cross-sectional areas at cell faces. To overcome this, similarly to (Bollermann et al. 2013), (Bollermann et al. 2011), we limit the outflow fluxes from the draining cell by introducing the so-called local draining time that is smaller than the CFL time restriction. This approach provides positivity preserving of the scheme without a reduction of the CFL time restriction.

The schemes which can simulate wetting/drying processes are important for water flow computation in river channels, especially, for mountain rivers. The wetting/drying possibilities are key requirements for the schemes used for 1D and 2D modeling of long-wave run-up in coastal areas [see, e.g., (Shokin et al. 2016)]. The paper presents the scheme and the algorithms for modeling water flow in a wet/dry multiply-connected channel network.

The paper is organized as follows. The definition of a channel network and open-flow water equations in channels are presented in Sect. 2. In Sect. 3.1, we briefly describe the central-upwind scheme from (Kurganov and Petrova 2007) implemented for the Saint–Venant equations. In Sect. 3.2, we present relations for the exact discretization of source terms in the Saint–Venant equations. We describe the new reconstruction method in Sect. 3.3. The slope friction treatment and channel junction models are given in Sects. 3.4 and 3.5, respectively. Boundary conditions, which are usually applied in modeling a subcritical open flow in a channel, are presented in Sect. 3.6. The positivity preserving and well-balancing properties of the scheme are proven in Sects. 3.7 and 3.8, respectively. Section 4 contains different numerical tests illustrating the merits of the scheme. Section 5 provides some concluding remarks.

2 Channel Network Definition and Water Flow Equations

The definition sketch of a channel network is shown in Fig. 1 that is similar to CHARIMA (Holly et al. 1990). We will consider any channel network as a certain directed graph in which the edges are separate channels (links), and the vertices are points (nodes) in which channels may merge into a single channel or one channel may diverge into separate ones. The graph may also include loops that correspond to channels around islands.

Fig. 1

Definition sketch of the channel network

The link represents the flow path between two nodes, and each link l is divided into non-uniform computational cells, which are numbered from 1 at the upper end of the link to n_l at the lower end of the link. Note that “lower” and “upper” generally, although not necessarily, mean “downstream” and “upstream”, respectively. Each link should have at least one computational cell.

The natural channel cross-section is replaced by its piecewise linear approximation (Fig. 2), which is assumed to be known for each computational cell edge.

Fig. 2

Channel cross-section and its piecewise linear approximation

One-dimensional shallow water flow in a natural channel with irregular cross-sections is described by the Saint–Venant equations (Cunge et al. 1980):

$$\frac{\partial A}{\partial t} + \frac{\partial Q}{\partial x} = 0$$

(1)

$$\frac{\partial Q}{\partial t} + \frac{{\partial \left( {{{Q^{2} } \mathord{\left/ {\vphantom {{Q^{2} } A}} \right. \kern-0pt} A} + gI_{1} } \right)}}{\partial x} = gI_{2} - gA(S_{0} + S_{f} )$$

(2)

where t is the time; x is the distance along the longitudinal axis of the channel; Q is the discharge; A is the wetted cross-section area; g is the gravitational constant; S_f is the friction slope due to the bed resistance; S₀ is the bed slope; I₁ is the hydrostatic pressure force; I₂ is the wall pressure force.

The friction slope S_f is assumed to be given by the Manning equation

$$S_{f} = \frac{{n^{2} Q\left| Q \right|}}{{A^{2} R^{{{4 \mathord{\left/ {\vphantom {4 3}} \right. \kern-0pt} 3}}} }}$$

(3)

where n is the Manning’s roughness coefficient; R is the cross-sectional hydraulic radius.

The hydrostatic pressure force, wall pressure force, wetted cross-section area and bed slope are defined, respectively, as (Cunge et al. 1980)

$$\begin{aligned} I_{1} = & \int_{0}^{h(x,t)} {(h - y)\sigma (x,y)dy} ; \\ I_{2} = & \int_{0}^{h(x,t)} {(h - y)\frac{\partial \sigma (x,y)}{\partial x}dy} \end{aligned}$$

(4)

$$A = \int_{0}^{h(x,t)} {\sigma (x,y)dy} ;\;S_{0} = \frac{\partial B}{\partial x}$$

(5)

where h is the water depth; $\sigma (x,h)$ is the channel width of a cross-section at point x and water depth h; B is the bed surface elevation (talweg or thalweg in geomorphology).

3 Numerical Discretization

We consistently renumber all channels (links) and nodes in a network, that is, each channel and node has its own unique number. The numbering of channels does not depend on the numbering of nodes.

In general, the index of the computational cells consists of the channel number and the cell number. The edges of the channel connected to one node are indexed by the node number and the channel number. Further, to simplify formulas where there will be no confusion, we omit some of the indices.

For simplicity, we assume that a channel geometry and a bottom topography do not change with time and depend only on the spatial coordinate. The other variables depend both on time and space. Also, where it will be possible, we will omit the temporal and spatial arguments.

3.1 Semi-Discrete Central-Upwind Scheme

A central-upwind scheme proposed by Kurganov and Petrova (2007) is used for numerical discretization of the Saint–Venant Eqs. (1) and (2).

We divide the spatial domain into grid cells $\left[ {x_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} ,x_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} } \right]$ of length $\Delta x_{j}$, where $x_{j}$ is the center of a grid cell, and denote by $\bar{U}_{j} (t)$ the cell averages of the solution $U = \left( {A,Q} \right)^{T}$ of (1) and (2) computed at time t

$$\bar{U}_{j} = \frac{1}{{\Delta x_{j} }}\int\limits_{{x_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} }}^{{x_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} }} {U(x,t)dx}$$

(6)

Then, a semi-discretization of (1) and (2) can be written as the following system of ODEs

$$\frac{d}{dt}\bar{U}_{j} (t) = - \frac{{H_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} (t) - H_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} (t)}}{{\Delta x_{j} }} + \frac{1}{{\Delta x_{j} }}\int\limits_{{x_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} }}^{{x_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} }} {\left[ {S(U,x) + G(U,x)} \right]dx}$$

(7)

where $H_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}$ are numerical fluxes at the cell interfaces $x_{{j \pm {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}$, $S(U,x) = \left( {0,gI_{2} - gAS_{0} } \right)^{T}$ and $G(U,x) = \left( {0, - gAS_{f} } \right)^{T}$.

The central-upwind numerical fluxes $H_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}$ are given by

$$\begin{aligned} H_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} & & & & & = \;\frac{{a_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } F(U_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } ,B_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} ) - a_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } F(U_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } ,B_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )}}{{a_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } - a_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } }} \\ & \quad + \frac{{a_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } a_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } }}{{a_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } - a_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } }}\left[ {U_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } - U_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } } \right] \\ \end{aligned}$$

(8)

where $U_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm } = \left( {A_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm } ,Q_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm } } \right)^{T}$ are the right/left-sided values of the piecewise linear reconstruction U at $x = x_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}$; $F(U,B) = \left( {Q,{{Q^{2} } \mathord{\left/ {\vphantom {{Q^{2} } A}} \right. \kern-0pt} A} + gI_{1} } \right)^{T}$ is the flux at the cell interface, and the one-sided local speeds are obtained using the eigenvalues of the Jacobian

$$\begin{aligned} a_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } = \hbox{max} \left\{ {0,u_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } + c_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } ,u_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } + c_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } } \right\} \hfill \\ a_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } = \hbox{min} \left\{ {0,u_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } - c_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ + } ,u_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } - c_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ - } } \right\} \hfill \\ \end{aligned}$$

(9)

Here $u = {Q \mathord{\left/ {\vphantom {Q A}} \right. \kern-0pt} A}$ and $c = \sqrt {g{A \mathord{\left/ {\vphantom {A {\sigma_{T} }}} \right. \kern-0pt} {\sigma_{T} }}}$, and $\sigma_{T}$ is the width of the channel at the water surface. The wetted water area at the cell interface may be very small and may lead to large values of flow velocity. In order to prevent this, we use the regularization technique suggested in (Kurganov and Petrova 2007)

$$u^{ \pm } = \frac{{\sqrt 2 A^{ \pm } Q^{ \pm } }}{{\sqrt {(A^{ \pm } )^{4} + \hbox{max} \left( {(A^{ \pm } )^{4} ,\varepsilon )} \right)} }}$$

(10)

where $\varepsilon$ is a small a priori chosen positive number. For consistency, the values of the discharge at cell interfaces are recomputed using $Q_{{j \pm {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm } = A_{{j \pm {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm } u_{{j \pm {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm }$, where $u_{{j \pm {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm }$ are given by (10).

The interface values $U_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm }$ are obtained from the cell averages $\bar{U}_{j}$ by a piecewise linear reconstruction which will be described below.

3.2 Auxiliary Relationships

In this section, we obtain relationships for computing the water volume between two specified channel cross-sections for an arbitrary free surface of the water. We also obtain formulas for calculating the hydrostatic pressure force I₁ and the wall pressure force I₂.

We replace the bed function B and wetted area A on the interval $\left[ {x_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} ,x_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} } \right]$ by their piecewise linear approximation

$$B(x) = B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} + (B_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )\frac{{x - x_{j} }}{{\Delta x_{j} }},$$

(11)

$$A(x,t) = \frac{1}{{\Delta x_{j} }}\int_{0}^{h(x,t)} {\left[ {\sigma_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} (y)(x - x_{{i - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} ) + \sigma_{{i - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} (y)(x_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - x)} \right]dy}$$

(12)

Assuming that the open water surface w is linear in a computational cell $\left[ {x_{L} ,x_{R} } \right]$, below we will use the following relations. Water volume V(x₁, x₂; w_L, w_R) between points x₁ and x₂ (Fig. 3) equals to

$$V(x_{1} ,x_{2} ;w_{L} ,w_{R} ) = \int_{{x_{1} }}^{{x_{2} }} {A(x,h)dx} = \int_{{x_{1} }}^{{x_{2} }} {\int_{0}^{h(x,t)} {\sigma (x,y)} dydx}$$

Fig. 3

Water surface w in a partially flooded cell $\left[ {x_{L} ,x_{R} } \right]$. Calculation of the water volume between x₁ and x₂

Changing the order of integration, we obtain

$$\begin{aligned} V(x_{1} ,x_{2} ;w_{L} ,w_{R} ) & = \frac{{(x_{2} - x_{1} )(2x_{R} - x_{1} - x_{2} )}}{2\Delta x}\int_{0}^{{h_{1} }} {\sigma_{L} dy - \frac{{(x_{R} - x_{2} )^{2} }}{2\Delta x}} \int_{{h_{1} }}^{{h_{2} }} {\sigma_{L} dy} + \frac{\Delta x}{{2\Delta h^{2} }}\int_{{h_{1} }}^{{h_{2} }} {\sigma_{L} \, (h_{R} - y)^{2} \, dy} \\ & \quad + \;\frac{{(x_{2} - x_{1} )(x_{1} + x_{2} - 2x_{L} )}}{2\Delta x}\int_{0}^{{h_{1} }} {\sigma_{R} dy + \frac{{(x_{2} - x_{L} )^{2} }}{2\Delta x}} \int_{{h_{1} }}^{{h_{2} }} {\sigma_{R} dy} - \frac{\Delta x}{{2\Delta h^{2} }}\int_{{h_{1} }}^{{h_{2} }} {\sigma_{R} \, (y - h_{L} )^{2} \, dy} \\ \end{aligned}$$

(13)

where $h_{i} = w_{i} - B_{i}$, $\Delta x = x_{R} - x_{L}$ and $\Delta h = h_{R} - h_{L}$. Note that for $h_{1} = h_{2}$, the second, third, fifth and sixth terms on the right-hand side of (13) are zero.

Similarly, the hydrostatic pressure force I₁ at a point x₁ can be calculated as

$$I_{1} (x_{1} ,h_{1} ) = \frac{{(x_{R} - x_{1} )}}{\Delta x}\int_{0}^{{h_{1} }} {(h_{1} - y) \, \sigma_{L} (y) \, dy} + \frac{{(x_{1} - x_{L} )}}{\Delta x}\int_{0}^{{h_{1} }} {(h_{1} - y) \, \sigma_{R} (y) \, dy}$$

(14)

Integration of the wall pressure force I₂ and $AS_{0}$ over an interval $\left[ {x_{1} ,x_{2} } \right]$ yields the following expressions

$$I_{2} (x_{1} ,x_{2} ) = \int_{{x_{1} }}^{{x_{2} }} {I_{2} dx} = \frac{{x_{2} - x_{1} }}{\Delta x}\int_{0}^{{h_{1} }} {(\sigma_{R} - \sigma_{L} )\left( {\frac{{h_{1} + h_{2} }}{2} - y} \right)dy} + \frac{1}{2\Delta h}\int_{{h_{1} }}^{{h_{2} }} {(\sigma_{R} - \sigma_{L} )(h_{2} - y)^{2} dy}$$

(15)

$$B_{x} (x_{1} ,x_{2} ) = \int_{{x_{1} }}^{{x_{2} }} {A\frac{\partial B}{\partial x}dx} = \frac{{B_{R} - B_{L} }}{2\Delta x}V(x_{1} ,x_{2} ;w_{L} ,w_{R} )$$

(16)

Remark

• 1. Note that for linear approximation of the water surface and bottom on an interval $\left[ {x_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} ,x_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} } \right]$, if the values of $w_{{j \mp {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm }$ and $B_{{j \pm {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}$ are known at the cell faces $x_{{j \mp {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}$, then it is easy to calculate $x^{ * }$ (Fig. 3) and integrals (13)–(16) for any internal interval $\left[ {x_{1} ,x_{2} } \right]$ and nonnegative $h(x,t)$.

• 2. The integrals (13)–(16) can be also computed for any piecewise-linear approximation of the water surface on an interval $\left[ {x_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} ,x_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} } \right]$.

3.3 Reconstruction of the Interface Values $U_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm }$

We present an algorithm of cell linear approximation of the $U_{j} (x)$ that takes into account all possible combinations of bottom topography and cell flooding. The slope of a cell linear reconstruction is defined by $U_{x,j}$ that is calculated as $U_{x,j} = \hbox{min} \bmod (\Delta^{ - } U_{j} ,\Delta^{ + } U_{j} )$ where $\Delta^{ - } U_{j}$ and $\Delta^{ + } U_{j}$ are the backward and forward difference derivatives. We denote by $h_{av,j}$, w_j and $w_{st,j}$ the depth of water which the surface is parallel to the cell bed, water surface elevation and the still water surface elevation in the cell j, respectively. In addition, we denote $h^{ \ge 0} = \hbox{max} (0,h)$.

Substituting $w_{1} = w_{2} = w_{st,j}$ into Eq. (13) and taking into account that $h_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - h_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} = B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - B_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} = - \Delta B$, we obtain the following equation for finding $w_{st,j}$

$$\begin{aligned} \frac{1}{2} \int_{0}^{{(w_{st,j} - B_{{j - 1/2}} )^{ \ge 0} }} {\sigma_{{j - 1/2}} dy} + \frac{1}{2}\int_{0}^{{(w_{st,j} - B_{{j + {1 / 2}}} )^{ \ge 0} }} {\sigma_{{j + {1 / 2}}} dy} \hfill \\ + \frac{1}{{2\Delta B^{2} }}\int_{{(w_{st,j} - B_{{j - {1 / 2}}} )^{ \ge 0} }}^{{(w_{st,j} - B_{{j + {1 / 2}}} )^{ \ge 0} }} {(w_{st,j} - B_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - y)^{2} \sigma_{{j - {1 / 2}}} dy} \hfill \\ - \frac{1}{{2\Delta B^{2} }}\int_{{(w_{st,j} - B_{{j - {1 / 2}}} )^{ \ge 0} }}^{{(w_{st,j} - B_{{j + {1 / 2}}} )^{ \ge 0} }} {(y - w_{st,j} + B_{{j - {1 / 2}}} )^{2} \sigma_{{j + {1 / 2}}} dy = \bar{A}_{j} } \end{aligned}$$

(17)

Equation (17) is a quartic equation for the piecewise linear approximation of the channel cross-section, and we can use the Newton–Raphson method or Ferrari’s method for its solution. We calculate $w_{st,j}$ only for cells in which still water may occur.

Similarly, substituting $h_{1} = h_{2} = h_{av,j}$ into Eq. (13), we obtain the following equation for finding $h_{av,j}$

$$\frac{1}{2}\int_{0}^{{h_{av,j} }} {\sigma_{{j - 1/2}} \; dy} + \frac{1}{2}\int_{0}^{{h_{av,j} }} {\sigma_{j+1 / 2} \; dy} = \bar{A}_{j}$$

(18)

which is a quadratic equation since

$$\begin{aligned} \int_{0}^{h} {\sigma_{j +1/2} \, dy} & \; & = \;\sum\limits_{k = 0}^{p - 1} {(h_{j+1/2,k + 1} - h_{j+1/2,k} )\frac{{\sigma_{j+1/2,k} + \sigma_{j+1/2,k + 1} }}{2}} \\ & + (h - h_{j+1/2,p} )\left[ {\sigma_{j+1/2,p} + \frac{1}{2}\frac{{\sigma_{j+1/2,p + 1} - \sigma_{j+1/2,p} }}{{h_{j+1/2,p + 1} - h_{j+1/2,p} }}(h - h_{j+1/2,p} )} \right] \quad {\text{for}}\;h_{j+1/2,p} < h \le h_{j+1/2,p + 1} \\ \end{aligned}$$

(19)

We divide conditionally all computational cells into two groups: “wet” and “dry”. A cell j will be called “wet” if $w_{st,j} \ge \hbox{max} (B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} ,B_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )$, and “dry” otherwise. We assume $w_{j} = w_{st,j}$ for the “wet” cells and $w_{j} = h_{av,j} + {{(B_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} + B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )} \mathord{\left/ {\vphantom {{(B_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} + B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )} 2}} \right. \kern-0pt} 2}$ for the “dry” cell. For “dry” cells in which still water may occur, we assume $w_{j} = w_{st,j}$ and $l_{j} = {{(w_{st,j} - B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )} \mathord{\left/ {\vphantom {{(w_{st,j} - B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )} {\left| {\Delta B_{j} } \right|}}} \right. \kern-0pt} {\left| {\Delta B_{j} } \right|}}$. For all other cells $l_{j} = 1$.

Depending on the combination of two neighboring cells, we use different formulas for approximating the forward and backward difference derivatives for w and Q. We approximate the backward difference derivative $\Delta^{ - }$ by the following way (Fig. 4):

Fig. 4

Approximation of the backward difference derivative $\Delta^{ - } w$

• 1. j − 1 and j cells are both “wet”. Then

  $$\Delta^{ - } w_{j} = \frac{{w_{st,j} - w_{st,j - 1} }}{{x_{j} - x_{j - 1} }};\;\Delta^{ - } Q_{j} = \frac{{\bar{Q}_{j} - \bar{Q}_{j - 1} }}{{x_{j} - x_{j - 1} }}$$

  (20)

• 2. j − 1 cell is “wet” and j cell is “dry”

  • if $B_{j+1/2} > B_{j-1/2}$

    $$\Delta^{ - } w_{j} = 2\frac{{w_{st,j} - w_{st,j - 1} }}{{\Delta x_{j - 1} + l_{j} \Delta x_{j} }};\;\Delta^{ - } Q_{j} = 2\frac{{\bar{Q}_{j} - \bar{Q}_{j - 1} }}{{\Delta x_{j - 1} + l_{j} \Delta x_{j} }}$$

    (21)

  • otherwise

    $$\begin{aligned}\Delta^{ - } w_{j} = &2\frac{{B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} + h_{av,j} - w_{st,j - 1} }}{{\Delta x_{j} }};\\ \Delta^{ - } Q_{j} = &2\frac{{\bar{Q}_{j} - \bar{Q}_{j - 1} }}{{\Delta x_{j} }}\end{aligned}$$

    (22)

• 3. j − 1 cell is “dry” and j cell is “wet”

  • if $B_{j -3/2} > B_{j -1/2}$

    $$\begin{aligned}\Delta^{ - } w_{j} =& 2\frac{{w_{st,j} - w_{st,j - 1} }}{{(1 + l_{j - 1} )\Delta x_{j - 1} + \Delta x_{j} }};\\\Delta^{ - } Q_{j} = &2\frac{{\bar{Q}_{j} - \bar{Q}_{j - 1} }}{{(1 + l_{j - 1} )\Delta x_{j - 1} + \Delta x_{j} }}\end{aligned}$$

    (23)

  • otherwise

    $$\begin{aligned}\Delta^{ - } w_{j} = &2\frac{{w_{st,j} - B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - h_{av,j - 1} }}{{\Delta x_{j} }};\\\Delta^{ - } Q_{j} = &2\frac{{\bar{Q}_{j} - \bar{Q}_{j - 1} }}{{\Delta x_{j} }}\end{aligned}$$

    (24)

• 4. j − 1 and j cells are both “dry”

  • if $B_{j - 3/2} > B_{j - 1/2}$ and $B_{j + 1/2} > B_{j - 1/2}$

    $$\begin{aligned}\Delta^{ - } w_{j} =& 2\frac{{w_{st,j} - w_{st,j - 1} }}{{(1 + l_{j - 1} )\Delta x_{j - 1} + l_{j} \Delta x_{j} }};\\\Delta^{ - } Q_{j} =& 2\frac{{\bar{Q}_{j} - \bar{Q}_{j - 1} }}{{(1 + l_{j - 1} )\Delta x_{j - 1} + l_{j} \Delta x_{j} }}\end{aligned}$$

    (25)

  • otherwise

    $$\Delta^{ - } w_{j} = \frac{{B_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - B_{{j - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} }}{{\Delta x_{j} }};\;\Delta^{ - } Q_{j} = 0$$

    (26)

The forward derivatives $\Delta^{ + }$ for w and Q are calculated in the similar way.

The interface point-values $w_{j \mp 1 / 2}^{ \pm }$ are obtained by a piecewise linear reconstruction

$$w_{j -1/2}^{ + } = \left\{ {\begin{array}{*{20}l} {w_{j} - 0.5 \, l_{j} \, w_{x,j} \,\Delta x_{j} ,\qquad \qquad {\text{if }} \, l_{j} \ge 0} \\ {w_{j} - 0.5 \, (1 - l_{j}) \, w_{x,j} \, \Delta x_{j} ,\quad \; {\text{otherwise}}} \\ \end{array} } \right.$$

(27)

and

$$w_{j +1/2}^{ - } = \left\{ {\begin{array}{*{20}l} {w_{j} + (1 - 0.5 l_{j}) \, w_{x,j} \, \Delta x_{j} ,\quad {\text{if }}l_{j} \ge 0} \\ {w_{j} + 0.5(1 + l_{j}) \, w_{x,j} \, \Delta x_{j} ,\quad {\text{otherwise}}} \\ \end{array} } \right.$$

(28)

According to (Kurganov and Tadmor 2000), this reconstruction will be second-order accurate if the approximate values of the derivatives $w_{x,j}$ are at least first-order approximations of corresponding exact derivatives. To ensure the non-oscillatory property in our numerical scheme, we evaluate $w_{x,j}$ using the minmod limiter (Nessyahu and Tadmor 1990), (Sweby 1984), (van Leer 1974)

$$w_{x,j} = \hbox{min} \bmod (\Delta^{ - } w_{j} ,\Delta^{ + } w_{j} )$$

(29)

where $\hbox{min} \bmod (a,b) = 0.5(\text{sgn} (a) + \text{sgn} (b)) \, \hbox{min} (\left| a \right|,\left| b \right|)$.

The interface values $Q_{{j \mp {1 / 2}}}^{ \pm }$ we calculate in the similar way. In the next step, water depth $h_{{j \pm {1 / 2}}}^{ \pm }$ at the cell faces is computed from the water surface elevation. Then, the area $A_{{j \pm {1 / 2}}}^{ \pm }$ is calculated from water depth based on channel geometry. Finally, we determine the numerical fluxes $H_{{j \pm {1 / 2}}}$ and $\int_{{x_{{j - {1 / 2}}} }}^{{x_{{j + {1 / 2}}} }} {S(U,x)} dx$ from (7), (8) using formulas (14)–(16) to compute I₁, $I_{2} (x_{{j - {1 / 2}}} ,x_{{j + {1 / 2}}} )$, and $B_{x} (x_{{j - {1 / 2}}} ,x_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )$.

3.4 Friction Slope Treatment and Time Evolution

We calculate the second component $G_{j}^{(2)}$ of $\int_{{x_{{j - {1 / 2}}} }}^{{x_{{j + {1 / 2}}} }} {G(U,x)} dx$ using the midpoint rule and de-singularization procedure

$$\frac{1}{{\Delta x_{j} }}\int_{{x_{{j - {1 / 2}}} }}^{{x_{{j + {1 / 2}}} }} {G(U,x)} dx = gn^{2} \frac{{\left| {\bar{Q}_{j} } \right|}}{{\hbox{max} (\bar{A}\,R^{{{4 \mathord{\left/ {\vphantom {4 3}} \right. \kern-0pt} 3}}} ,\varepsilon_{1} )}}\bar{Q}_{j} = G_{j}^{(2)} \bar{Q}_{j}$$

(30)

where $\varepsilon_{1}$ is a priori chosen small positive number.

In the general case, the friction slope term in (7) can be a source of stiffness for the ODE system when the depth of water is very small (Chertock et al. 2015). The explicit treatment of the friction term imposes a severe time step restriction, the result of which is that the time step should be several times less than a typical time step under the CFL conditions. To overcome this difficulty, we use the forward Euler method of time integration of the ODE system (7) with an implicit treatment of only a part of the friction slope term. As a result, we obtained the following fully discrete central-upwind scheme

$$\bar{A}_{j} (t + \Delta t) = \bar{A}_{j} (t) - \frac{\Delta t}{{\Delta x_{j} }}\left( {H_{j+1/2}^{(1)} (t) - H_{j -1/2}^{(1)} (t)} \right)$$

(31)

$$\bar{Q}_{j} (t + \Delta t) = \frac{{\bar{Q}_{j} (t) -({\Delta t} \mathord{\left/ {\vphantom {{\Delta t} {\Delta x_i}}} \right. \kern-\nulldelimiterspace} {\Delta x_j }) \, \left( {H_{j+1/2}^{(2)} (t) - H_{j-1/2}^{(2)} (t) - gI_{2} (x_{j -1/2} ,x_{j +1/2} ) + gB_{x} (x_{j -1/2} ,x_{j +1/2} )} \right)}}{{1 + \Delta t\,G_{j}^{(2)} (t)}}$$

(32)

where $H_{j + 1/2}^{(1)}$ and $H_{j + 1/2}^{(2)}$ are the first and second component of the numerical fluxes (8).

3.5 Node of Junction

For modeling of the hydraulic conditions at a junction, we will use one of two approaches: (1) assuming that only the mass conservation equation holds at the junction, and (2) based on mass and momentum conservation.

Consider a node s that connects two or more channels, where J_s,in is a set of channels which inflow and J_s,out that outflow from it (Fig. 5). Denote by $x_{s,i}^{in}$ and $x_{s,j}^{out}$ the coordinates of the ends of channels connected to the node s. Channel cross-sections and bottom levels are known at these points. In the first case, water surface elevation Y_s, or water surface elevation Y_s and water discharge Q_s for the second case at the junction are also given.

Fig. 5

Sketch of a node that is the junction of channels

Approximating that water surface is horizontal, the continuity equation over the control volume around a junction s can be discretized as follows

$$\begin{aligned} F_s^{A} \left( {Y_s (t + \Delta t)} \right) & = \sum\limits_{{i \in J_{s,in} }} {\left[ {V(x_{i,n_{i} +{1/2}},x_{s,i}^{in} ;Y_{s} (t + \Delta t),Y_{s} (t + \Delta t)) - V(x_{{i,n_{i} + {1 / 2}}} ,x_{s,i}^{in} ;Y_{s} (t),Y_{s} (t))} \right]} \\ & \quad + \sum\limits_{{j \in J_{s,out} }} {\left[ {V(x_{s,j}^{out} ,x_{{j,{1 / 2}}} ;Y_{s} (t + \Delta t),Y_{s} (t + \Delta t)) - V(x_{s,j}^{out} ,x_{{j,{1/2}}} ;Y_{s} (t),Y_{s} (t))} \right]} \\ & \quad - \Delta t\sum\limits_{{i \in J_{s,in} }} {\left\{ {\frac{{a_{{i,n_{i} + {1/2}}}^{ + } Q_{{i,n_{i} + {1/2}}}^{ - } - a_{{i,n_{i} + {1/ 2}}}^{ - } Q_{{i,n_{i} + {1 / 2}}}^{ + } }}{{a_{{i,n_{i} + {1 / 2}}}^{ + } - a_{{i,n_{i} + {1 / 2}}}^{ - } }} + \frac{{a_{{i,n_{i} + {1 / 2}}}^{ + } a_{{i,n_{i} + {1/ 2}}}^{ - } }}{{a_{{i,n_{i} + {1 / 2}}}^{ + } - a_{{i,n_{i} + {1/ 2}}}^{ - } }}\left[ {A_{{i,n_{i} + {1 / 2}}}^{ + } - A_{{i,n_{i} + {1 / 2}}}^{ - } } \right]} \right\}} \\ & \quad + \Delta t\sum\limits_{{j \in J_{s,out} }} {\left\{ {\frac{{a_{{j,{1 /2}}}^{ + } Q_{{j,{1 / 2}}}^{ - } - a_{{j,{1 / 2}}}^{ - } Q_{{j,{1 /2}}}^{ + } }}{{a_{{j,{1 /2}}}^{ + } - a_{{j,{1 / 2}}}^{ - } }} + \frac{{a_{{j,{1 / 2}}}^{ + } a_{{j,{1 / 2}}}^{ - } }}{{a_{{j,{1 / 2}}}^{ + } - a_{{j,{1 / 2}}}^{ - } }}\left[ {A_{{j,{1 /2}}}^{ + } - A_{{j,{1 / 2}}}^{ - } } \right]} \right\}} = 0 \\ \end{aligned}$$

(33)

Accordingly, the difference scheme for the equation of momentum conservation can be written in the form

$$ \begin{aligned} F_{s}^{Q} \left( {\bar{Q}_{s} (t + \Delta t)} \right) & = \left( {\bar{Q}_{s} (t + \Delta t) - \bar{Q}_{s} (t)} \right)\left[ {\sum\limits_{{i \in J_{s,in} }} {(x_{s,i}^{in} - x_{{i,n_{i} + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )} + \sum\limits_{{j \in J_{s,out} }} {(x_{{j,{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - x_{s,j}^{out} )} } \right] \\ & \quad - \Delta t\sum\limits_{{i \in J_{s,in} }} {\left\{ {H_{{i,n_{i} + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{(2)} - gI_{1} (x_{s,i}^{in} ,h_{s} ) + gI_{2} (x_{{{{i,n_{i} + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} }} ,x_{s,i}^{in} ) - gB_{x} (x_{{{{i,n_{i} + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} }} ,x_{s,i}^{in} )} \right\}} \\ & \quad + \Delta t\sum\limits_{{j \in J_{s,out} }} {\left\{ {H_{{j,{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{(2)} - gI_{1} (x_{s,j}^{out} ,h_{s} ) - gI_{2} (x_{s,j}^{out} ,x_{{j,{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} ) + gB_{x} (x_{s,j}^{out} ,x_{{j,{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )} \right\}} \\ & \quad - \Delta t\,\bar{Q}_{s} (t + \Delta t)\left[ {\sum\limits_{{i \in J_{s,in} }} {G_{s,i}^{(2)} (x_{s,i}^{in} - x_{{i,n_{i} + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} )} + \sum\limits_{{j \in J_{s,out} }} {G_{s,j}^{(2)} (x_{{j,{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - x_{s,j}^{out} )} } \right] = 0 \\ \end{aligned} $$

(34)

The interface values $A_{{j + {1 / 2}}}^{ \pm }$ are calculated from the water surface elevation, which is reconstructed in a similar way as described in Sect. 3.3. The interface values $Q_{{j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}^{ \pm }$ are computed: in the first case as $Q_{{i,n_{i} + {1 / 2}}}^{ + } = Q_{{i,n_{i} + {1 / 2}}}^{ - }$, $Q_{{j,{1 / 2}}}^{ - } = Q_{{j,{1 / 2}}}^{ + }$, and in the second case as $Q_{{i,n_{i} + {1 / 2}}}^{ + } = Q_{{j,{1 / 2}}}^{ - } = \bar{Q}_{s}$.

The functions $F_{s}^{A}$ and $F_{s}^{Q}$ may be thought of as the residuals of the continuity and momentum conservation Eqs. (33)–(34), respectively. The water surface elevation $Y_{s} (t + \Delta t)$ and the water discharge $Q_{s} (t + \Delta t)$ are found from the solution of Eqs. (33)–(34). The Newton–Raphson method is used to solve the nonlinear Eq. (33). Note that the derivative of $F_{s}^{A} \left( {Y_{s} (t + \Delta t)} \right)$ with respect to $Y_{s} (t + \Delta t)$ can be expressed as

$$\begin{aligned}\frac{{\partial F_{s}^{A} }}{{\partial Y_{s} }} & = \frac{1}{2}\sum\limits_{{i \in J_{s,in} }} {\left( {x_{s,i}^{in} - x_{{i,n_{i} + {1/2}}} } \right)\left[ {\sigma_{{i,n_{i} + {1/2}}} \left( {h_{{i,n_{i} + {1/2}}}^{ \ge 0} } \right) + \sigma_{s,i} \left( {h_{s,i}^{ \ge 0} } \right)} \right.} \\ & \quad - \left( {B_{s,i} - B_{{i,n_{i} + {1/2}}} } \right)^{ - 2} \left\{ {\left[ {\sigma_{{i,n_{i} + {1/2}}} \left( {h_{{i,n_{i} + {1/2}}}^{ \ge 0} } \right) + \sigma_{s,i} \left( {h_{s,i}^{ \ge 0} } \right)} \right]\left( {B_{s,i} - B_{{i,n_{i} + {1/2}}} } \right)^{2} } \right. \\ & \quad - 2\left. {\left. {\int_{{(Y_{s} - B_{{i,n_{i} + {1/2}}} )^{ \ge 0} }}^{{(Y_{s} - B_{s,i} )^{ \ge 0} }} {\left[ {(Y_{s} - B_{s,i} - y)\sigma_{{i,n_{i} + {1/2}}} (y) + (y - Y_{s} + B_{{i,n_{i} + {1/2}}} )\sigma_{s,i} (y)} \right]dy} } \right\}} \right] \\ & \quad + \frac{1}{2}\sum\limits_{{j \in J_{s,out} }} {\left( {x_{{j,{1/2}}} - x_{s,j}^{out} } \right)\left[ {\sigma_{s,j} \left( {h_{s,j}^{ \ge 0} } \right) + \sigma_{{j,{1/2}}} \left( {h_{{j,{1/2}}}^{ \ge 0} } \right)} \right.} \\ & \quad - \left( {B_{{j,{1/2}}} - B_{s,j} } \right)^{ - 2} \left\{ {\left[ {\sigma_{s,j} \left( {h_{s,j}^{ \ge 0} } \right) + \sigma_{{j,{1/2}}} \left( {h_{{j,{1/2}}}^{ \ge 0} } \right)} \right]} \right.\left( {B_{{j,{1/2}}} - B_{s,j} } \right)^{2} \\ & \quad - 2\left. {\left. {\int_{{(Y_{s} - B_{s,j} )^{ \ge 0} }}^{{(Y_{s} - B_{{j,{1 /2}}} )^{ \ge 0} }} {\left[ {(Y_{s} - B_{{j,{1 /2}}} - y)\sigma_{s,j} (y) + (y - Y_{s} + B_{s,j} )\sigma_{{j,{1 / 2}}} (y)} \right]} dy} \right\}} \right] \\ \end{aligned}$$

(35)

where $h^{ \ge 0} = \hbox{max} (0,h)$. Note that for $h_{{i,n_{i} + {1/2}}} = h_{s,i}$ and $h_{{j,{1/2}}} = h_{s,j}$ the expressions in the corresponding curly brackets on the right-hand side of (35) are zero.

3.6 Boundary Conditions

For the numerical solution of the ODEs (7), we should specify values of the central-upwind numerical fluxes $H_{{j + {1/2}}}$ at an upstream, and if required downstream, boundary of the computational area. We consider the most commonly used boundary conditions in the modeling of a subcritical fluid flow with a free surface in channel networks. In the following, we discuss the boundary treatment at the left boundary. The boundary treatment at the right boundary can be done in a similar way.

Discharge boundary conditions Let a water discharge Q(t) be provided on the left boundary. We will consider the following approximation of the boundary numerical fluxes $H_{{{1 /2}}}$

$$\begin{aligned} Q = \frac{a_{1/2}^+ Q_{1/2}^- - a_{1/2}^- Q_{1/2}^+ }{a_{1/2}^+ - a_{1/2}^-} + \frac{a_{1/2}^{ - } a_{1/2}^{ + } }{a_{1/2}^{ + } - a_{1/2}^{ - } } \left( {A_{1/2}^{ + } - A_{1/2}^{ - } } \right) \hfill \\ \frac{Q^{2} }{A_{1/2}^{ - } } + gI_{1} (h_{1/2}^{ - } ) = \frac{a_{1/2}^+ {(Q_{1 /2}^ - )^2 \mathord{\left/ {\vphantom {(Q_{1 /2}^- ) A_{1/2}^- }} \right. \kern-\nulldelimiterspace} A_{1/2}^- } - a_{1/2}^{ - } (Q_{1/2}^+ )^{2} \mathord{ \left/ {\vphantom {(Q_{1 /2}^- ) A_{1/2}^- }} \right. } A_{1/2}^+ } {a_{1/2}^+ - a_{1/2}^- } + \frac{a_{1/2}^- a_{1/2}^+ } {a_{1/2}^+ - a_{1/2}^- } \left( {Q_{1/2}^+ - Q_{1/2}^- } \right) + g\frac{{a_{1/2}^+ I_{1} (h_{1/2}^- ) - a_{1/2}^- I_{1} (h_{1/2}^+ )}}{{a_{1/2}^+ - a_{1/2}^- }} \hfill \\ \end{aligned}$$

(36)

The values of $Q_{1/2}^{ - }$ and $h_{1/2}^{ - }$ are found from a solution of the system of Eqs. (36).

Surface boundary condition Let a water surface elevation Y(t) be specified at the left boundary. Then we will use the following approximation of the boundary fluxes.

$$\begin{aligned} Q_{1/2}^- = \frac{a_{1/2}^+ Q_{1/2}^- - a_{1/2}^- Q_{1/2}^+ } {a_{1/2}^+ - a_{1/2}^- } + \frac{a_{1/2}^- a_{1/2}^+ }{a_{1/2}^+ - a_{1/2}^- } \left( {A_{1/2}^+ - A_{1/2}^- } \right) \hfill \\ \frac{{(Q_{1/2}^- )^2 }}{A(h)} + gI_{1} (h) = \frac{{a_{1/2}^+ {(Q_{1/2}^- )^2 } \mathord{ \left/ {\vphantom {(Q_{1 /2}^- ) A_{1/2}^- }} \right. } {A_{1/2}^- } - a_{1/2}^- {(Q_{1/2}^+ )^2 } \mathord{ \left/ {\vphantom {(Q_{1 /2}^+ ) A_{1/2}^+ }} \right. } {A_{1/2}^+ }}} {a_{1/2}^+ - a_{1/2}^- } + \frac{a_{1/2}^- a_{1/2}^+ }{a_{1/2}^+ - a_{1/2}^- } \left( {Q_{1/2}^+ - Q_{1/2}^- } \right) + g\frac{a_{1/2}^+ I_{1} (h_{1/2}^- ) - a_{1/2}^- I_{1} (h_{1/2}^+ )}{a_{1/2}^+ - a_{1/2}^- } \hfill \\ \end{aligned}$$

(37)

where $h = Y - B_{1/2}$; $Q_{1/2}^{ - }$ and $h_{1/2}^{ - }$ are found from a solution of the Eqs. (37).

We will also consider a simpler approximation of the surface boundary conditions. In this case, we assume that

$$\begin{aligned} Q_{1/2}^{-} & = Q_{1/2}^{+} \hfill \\ h_{1/2}^{-} & = h \hfill \\ \end{aligned}$$

(38)

Outflow boundary condition For outflow, we use the following extrapolations on the boundary

$$\begin{aligned} Q_{1/2}^{ - } = Q_{1/2}^{ + } \hfill \\ h_{1/2}^{ - } = h_{1/2}^{ + } \hfill \\ \end{aligned}$$

(39)

3.7 Positivity Preserving

In this section, we show that the resulting central-upwind scheme is positivity preserving. A sufficient condition for this is to ensure that at each time step no more water outflows from a cell than there is at the moment. For the positivity, our result is similar to the one obtained in (Balbás and Hernandez-Duenas 2014).

For any computational cell k of the link p, we define $\Delta t_{p,k}$ as

$$\Delta t_{p,k} = \hbox{min} \,\left\{ {\frac{\Delta x_{p,k} }{a_{p,k - 1/2}^{+} - a_{p,k + 1/2}^{-} },\frac{\Delta x_{p,k} \,\bar{A}_{p,k}^{n} }{{a_{p,k + 1/2}^+ \, A_{p,k + 1/2}^- - a_{p,k - 1/2}^- \, A_{p,k - 1/2}^+ }}} \right\}$$

(40)

For any node s, which has two or more links, we denote by $\Delta t_{s}$

$$\Delta t_s = \hbox{min} \left\{ {\mathop {\hbox{min} } \limits_{i} \,\left[ {\frac{\Delta x_{s,i}^{in} } {a^+_{i,n_i + 1/2} },\frac{ - \Delta x_{s,i}^{in} \,\bar{A}_{s,i}^{in,n} }{{a_{i,n_{i} + 1/2}^{-} \, A_{i,n_{i} + 1/2}^{+} }}} \right],\,\mathop {\hbox{min} }\limits_{j} \,\left[ {\frac{ - \Delta x_{s,j}^{out} }{a_{j,1/2}^- },\frac{\Delta x_{s,j}^{out} \,\bar{A}_{s,j}^{out,n} }{a_{j,1/2}^+ \, A_{j,1/2}^- }} \right]} \right\}$$

(41)

where $V(x_{s,j}^{out} ,x_{{j,{1/2}}} ;Y_{s} (t),Y_{s} (t)) = \Delta x_{s,j}^{out} \,\bar{A}_{s,j}^{out,n}$ and $\Delta x_{s,j}^{out} = x_{{j,{1/2}}} - x_{s,j}^{out}$. The values of $\,\bar{A}_{s,i}^{in,n}$ and $\Delta x_{s,i}^{in} \,$ are defined in a similar way.

Theorem 1

Consider the semi-discrete central-upwind scheme (7)–(10), (33) with the piecewise linear reconstruction described in Sect. 3.3and the discretization of the source terms (14)–(16). Assume that the system of ODEs (7) is solved by the forward Euler method with an implicit treatment of the friction slope and for all computational cells and nodes$h_{j}^{n} \ge 0$. Then all$h_{j}^{n + 1} \ge 0$, if

$$\Delta t \le \hbox{min} \left\{ {\mathop {\hbox{min} }\limits_{p,k} \,\Delta t_{p,k} ,\Delta t_{s} } \right\}$$

(42)

where $\Delta t_{p,k}$ and $\Delta t_{s}$ are calculated from (40) to (41).

Proof

The cross-sectional area A(h) is a nonnegative increasing function of water depth. Therefore, our task is to obtain conditions that will ensure the non-negativity of the cross-sectional area at the next time step.

Note, as follows from (9), that $a_{{j + {1/2}}}^{ + } \ge 0$, $a_{{j + {1/2}}}^{ - } \le 0$, $a_{j + {1/2}}^{ + } - u_{j + {1/2}}^{ + } \ge 0$, and $u_{j + {1/2}}^{ - } - a_{j + {1/2}}^{ - } \ge 0$ for all j. Moreover, the following inequalities $0 \le \frac{u_{j + 1/2}^- - a_{j + 1/2}^- }{a_{j + 1/2}^+ - a_{j + 1/2}^- } \le 1$, $0 \le \frac{a_{j - 1/2}^+ - u_{j - 1/2}^+ }{a_{j - 1/2}^+ - a_{j - 1/2}^- } \le 1$ are satisfied for all j. Then for any computational cell k of link p from (31), we have

$$\begin{aligned} \bar{A}_{p,k}^{n + 1} & = \bar{A}_{p,k}^{n} - \frac{\Delta t}{\Delta x_{p,k} } \, a_{p,k + 1/2}^{+} \left( {\frac{u_{p,k + 1/2}^{-} - a_{p,k + 1/2}^{-} }{a_{p,k + 1/2}^{+} - a_{p,k + 1/2}^{-} }} \right)\times{A_{p,k + 1/2}^{-}} - \frac{\Delta t}{\Delta x_{p,k} } \, a_{p,k + 1/2}^{-} \left( {\frac{a_{p,k + 1/2}^{+} - u_{p,k + 1/2}^{+} }{a_{p,k + 1/2}^+ - a_{p,k + 1/2}^- }} \right)A_{p,k + 1/2}^{+} \\ & \quad + \frac{\Delta t}{\Delta x_{p,k} } \, a_{p,k - 1/2}^{+} \left( {\frac{u_{p,k - 1/2}^{-} - a_{p,k - 1/2}^{-} }{a_{p,k - 1/2}^{+} - a_{p,k - 1/2}^{-} }} \right)A_{p,k - 1/2}^{-} + \frac{\Delta t}{\Delta x_{p,k} } \, a_{p,k - 1/2}^{-} \left( {\frac{a_{p,k - 1/2}^{+} - u_{p,k - 1/2}^{+} } {a_{p,k - 1/2}^{+} - a_{p,k - 1/2}^- }} \right)A_{p,k - 1/2}^{+} \\ \end{aligned}$$

(43)

The third and fourth terms on the right-hand side of (43) are nonnegative, so that for $\bar{A}_{p,k}^{n + 1} \ge 0$ it is sufficient if the following inequality is satisfied

$$\Delta x_{p,k} \,\bar{A}_{p,k}^{n} - \Delta t\left( {a_{p,k + 1/2}^{+} A_{p,k + 1/2}^{-} - a_{p,k - 1/2}^{-} A_{p,k - 1/2}^{+} } \right) \ge 0$$

(44)

For any node s, we rewrite (33) in the following way

$$\begin{aligned} \sum\limits_{{i \in J_{s,in} }} {\Delta x_{s,i}^{in} \, \bar{A}_{s,i}^{in,n + 1} } + \sum\limits_{{j \in J_{s,out} }} {\Delta x_{s,j}^{out} \, \bar{A}_{s,j}^{out,n + 1} } \, & = \sum\limits_{{i \in J_{s,in} }} {\Delta x_{s,i}^{in} \, \bar{A}_{s,i}^{in,n} } + \sum\limits_{{j \in J_{s,out} }} {\Delta x_{s,j}^{out} \, \bar{A}_{s,j}^{out,n} } \, \\ & \quad + \Delta t\sum\limits_{i \in J_{s,in} } {\left\{ {a_{i,n_{i} + 1/2}^{+} \left( {\frac{u_{i,n_{i} + 1/2}^{-} - a_{i,n_{i} +1/2}^{-} }{a_{i,n_{i} + 1/2}^{+} - a_{i,n_{i} + 1/2}^{-} }} \right)\times{A_{i,n_{i} + 1/2}^{-} }+ a_{i,n_{i} + 1/2}^{-} \left( {\frac{a_{i,n_{i} + 1/2}^{+} - u_{i,n_{i} + 1/2}^{+} }{a_{i,n_{i} + 1/2}^{+} - a_{i,n_{i} + 1/2}^{-} }} \right)A_{i,n_{i} +1/2}^{+} } \right\}} \\ & \quad - \Delta t\sum\limits_{j \in J_{s,out} } {\left\{ {a_{j,1/2}^{+} \left( {\frac{u_{j,1/2}^{-} - a_{j,1/2}^{-} }{{a_{j,1/2}^{+} - a_{j,1/2}^{-} }}} \right)A_{j,1/2}^{-} + a_{j,1/2}^{-} \left( {\frac{a_{j,1/2}^{+} - u_{j,1/2}^{+} }{a_{j,1/2}^{+} - a_{j,1/2}^{-} }} \right)A_{j,1/2}^{+} } \right\}} \\ \end{aligned}$$

(45)

It is clear that the left-hand side of (45) will be nonnegative if

$$\Delta x_{s,i}^{in} \,\bar{A}_{s,i}^{in,n} + \Delta t\,a_{i,n_{i} + 1/2}^{ - } A_{i,n_{i} + 1/2}^{ + } \ge 0$$

(46)

$$\Delta x_{s,j}^{out} \,\bar{A}_{s,j}^{out,n} - \Delta t\,a_{j,1/2}^{ + } \,A_{j,1/2}^{ - } \ge 0$$

(47)

The statement of the theorem follows from (44), (46)–(47) and the CFL restriction.

From (41), one can see that the positivity preserving step depends on the ratios ${{A_{{j \mp {1/2}}}^{ \pm } } \mathord{\left/ {\vphantom {{A_{{j \mp {1/2}}}^{ \pm } } {\bar{A}_{j} }}} \right. \kern-0pt} {\bar{A}_{j} }}$. Due to bottom topography, irregular channel geometry and the possibility of a channel drying, these ratios can be much more than unity. In this case, the positive preserving time step restrictions are more severe than the CFL restrictions. In order to overcome these restrictions and guarantee the positivity preserving of our scheme, we adopt a draining time step technique from Bernstein et al. (2016), Bollermann et al. (2013), Bollermann et al. (2011). The basic idea of this approach is to reduce the time step locally only for the cell faces at which the inequalities (44), (46)–(47) do not necessarily hold.

Following (Bollermann et al. 2011), we introduce the draining time step

$$\Delta t_{p,k}^{drain} = \frac{\Delta x_{p,k} \,\bar{A}_{p,k}^{n} }{\hbox{max} (0,H_{p,k + 1/2}^{(1)} ) + \hbox{max} (0, - H_{p,k - 1/2}^{(1)} )}$$

(48)

$$\Delta t_{s,i}^{in,drain} = \frac{\Delta x_{s,i}^{in} \,\bar{A}_{s,i}^{in,n} }{a_{i,n_{i} + 1/2}^{-} \,A_{i,n_{i} + 1/2}^{+} } \left( {\frac{a_{i,n_{i} + 1/2}^{+} - a_{i,n_{i} + 1/2}^{-} }{a_{i,n_{i} + 1/2}^{+} - u_{i,n_{i} + 1/2}^{+} }} \right)$$

(49)

$$\Delta t_{s,j}^{out,drain} = \frac{\Delta x_{s,j}^{out} \, \bar{A}_{s,j}^{out,n} }{a_{j,1/2}^{+} \, A_{j,1/2}^{-} } \left( {\frac{a_{j,1/2}^{+} - a_{j,1/2}^{-} }{u_{j,1/2}^{-} - a_{j,1/2}^{-} }} \right)$$

(50)

which describes the time when the water contained in cell in the beginning of the time step outflows from it. Now the evolution step in (43), (45) is replaced by

$$\bar{A}_{p,k}^{n + 1} = \bar{A}_{p,k}^{n} - \frac{\Delta t_{p,k + 1/2} \, H_{p,k + 1/2}^{(1)} - \Delta t_{p,k - 1/2} \, H_{p,k - 1/2}^{(1)} } {\Delta x_{p,k}}$$

(51)

$$\begin{aligned} \sum\limits_{{i \in J_{s,in} }} {\Delta x_{s,i}^{in} \, \bar{A}_{s,i}^{in,n + 1} } + \sum\limits_{{j \in J_{s,out} }} {\Delta x_{s,j}^{out} \, \bar{A}_{s,j}^{out,n + 1} } \, = \sum\limits_{{i \in J_{s,in} }} {\Delta x_{s,i}^{in} \, \bar{A}_{s,i}^{in,n} } + \sum\limits_{{j \in J_{s,out} }} {\Delta x_{s,j}^{out} \, \bar{A}_{s,j}^{out,n} } \, \hfill \\ - \sum\limits_{{j \in J_{s,out} }} {\Delta t_{j,1/2} \, H_{j,1/2}^{(1)} } + \sum\limits_{{i \in J_{s,in} }} {\Delta t_{i,n_{i} + 1/2} \, H_{i,n_{i} + 1/2}^{(1)} } \hfill \\ \end{aligned}$$

(52)

The outflow of a cell is the inflow of its neighbor. Thus, to keep the conservativity of the scheme, we define $\Delta t_{p,k + 1/2}$ as

$$\Delta t_{{p,k + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} = \hbox{min} (\Delta t,\Delta t_{p,m}^{drain} ), \qquad m = k + \frac{1}{2} - \frac{{\text{sgn} (H_{p,k + 1/2}^{(1)} )}}{2}$$

(53)

where $\Delta t$ satisfies CFL conditions.

For consistency, the momentum conservation Eq. (32) should be rewritten in the form

$$\begin{aligned} \bar{Q}_{p,k}^{n + 1} = \frac{\Delta x_{p,k} \, \bar{Q}_{p,k}^{n} - {\Delta t_{p,k + 1/2}} \, H_{p,k + 1/2}^{(2,a)} + {\Delta t_{p,k - 1/2} } \, H_{p,k - 1/2}^{(2,a)} } { \Delta x_{p,k} \left( {1 + \Delta t \, G_{p,k}^{(2)} } \right) } \hfill \\- \frac{\Delta t}{\Delta x_{p,k} } \frac{H_{p,k + 1/2}^{(2,g)} - H_{p,k - 1/2}^{(2,g)} - g I_{2} (x_{p,k - 1/2} ,x_{p,k +1/2} ) + g B_{x} (x_{p,k - 1/2} ,x_{p,k + 1/2} )} {1 + \Delta t\,G_{p,k}^{(2)} } \hfill \\ \end{aligned}$$

(54)

where $H_{{p,k \pm {1/2}}}^{(2,a)}$ and $H_{{p,k \pm {1/2}}}^{(2,g)}$ are the advective and gravity-driven parts of the $H_{{p,k \pm {1/2}}}^{(2)}$, which are defined as

$$\begin{aligned} H_{p,k + 1/2}^{(2,a)} & = \frac{a_{p,k + 1/2}^{+} }{a_{p,k + 1/2}^{+} - a_{p,k + 1/2}^{-} }\frac{(Q_{p,k + 1/2}^{-} )^{2} }{A_{p,k + 1/2}^{-} } - \frac{a_{p,k + 1/2}^{-} }{a_{p,k + 1/2}^{+} - a_{p,k + 1/2}^{-} } \frac{(Q_{p,k + 1/2}^{+} )^{2} }{A_{p,k + 1/2}^{+} } \\ & \quad + \frac{a_{p,k + 1/2}^{+} \; a_{p,k + 1/2}^{-} }{a_{p,k + 1/2}^{+} - a_{p,k + 1/2}^{-} } \left[ {Q_{p,k + 1/2}^{+} - Q_{p,k + 1/2}^{-} } \right] \\ \end{aligned}$$

(55)

$$H_{p,k + 1/2}^{(2,g)} = g\frac{a_{p,k + 1/2}^{+} \, I_{1} (x_{p,k + 1/2} ,h_{p,k + 1/2}^{-} ) - a_{p,k + 1/2}^{-} \, I_{1} (x_{p,k + 1/2} ,h_{p,k + 1/2}^{+} )}{a_{p,k + 1/2}^{+} - a_{p,k + 1/2}^{-} }$$

(56)

3.8 Well-balancing

The system (1), (2) admits smooth steady-state solutions, satisfying

$$Q = uA = Const, \quad E = \frac{1}{2}\frac{{Q^{2} }}{{A^{2} }} + g(h + B) = Const$$

as well as nonsmooth steady-state solutions. One of the most important steady-state solutions is a trivial stationary one (lake at rest)

$$Q = 0, \quad h + B = Const$$

A numerical scheme that exactly preserves steady flow at rest is called well-balanced.

Theorem 2

Consider the semi-discrete central-upwind scheme (51)–(52), (54) with the piecewise linear reconstruction described in Sect. 3.3and the discretization of the source terms (14)–(16). Assume that the numerical solution U(tⁿ) corresponds to the steady state at rest. Then U(tⁿ⁺¹) = U(tⁿ), that is, the scheme is well–balanced.

Proof

As a result of the reconstruction, we have $w_{j \pm 1/2}^{ \mp } = w_{j}$ and $Q_{j \pm 1/2}^{ \mp } = 0$. It is clear that the equations of the mass conservation (51), (52) are satisfied. Thus, we have to show that in the momentum conservation Eq. (54), there has to be a balance between the flux gradient and the source term.

For the second component $H_{j + 1/2}^{(2)}$ from (55) to (56), we obtain

$$H_{j + 1/2}^{(2,a)} + H_{j + 1/2}^{(2,g)} = g\,I_{1} (x_{j + 1/2} ,h_{j + 1/2} )$$

(57)

as $h_{j + 1/2}^{ + } = h_{j + 1/2}^{ - }$ and $A_{j + 1/2}^{ + } = A_{j + 1/2}^{ - }$. On the other hand, from (14) to (16) we have

$$\begin{aligned}I_{1} (x_{j + 1/2} ,h_{j + 1/2} ) - I_{1} (x_{j - 1/2} ,h_{j - 1/2} ) - I_{2} (x_{j - 1/2} ,x_{j + 1/2} ) + B_{x} (x_{j - 1/2} ,x_{j +1/2} ) \hfill \\ = \int_{{x_{j - 1/2} }}^{{x_{j + 1/2} }} {\frac{d}{dx}} \int_{0}^{h} {(w - B - y)\sigma (x,y)dy} \, dx - \int_{{x_{j - 1/2} }}^{{x_{j + 1/2} }} {\int_{0}^{h} {(w - B - y)\frac{d}{dx}\sigma (x,y)dy} } \, dx \hfill \\ + \int_{{x_{j - 1/2}}} ^{{x_{j + 1/2} }} {\frac{dB}{dx}\int_{0}^{h} {\sigma (x,y)dy} } \, dx = 0 \hfill \\ \end{aligned}$$

(58)

which ensures a balance of source terms with the flux terms in steady states at rest.

4 Numerical Results

In this section, we examine the scheme accuracy on the problem with smooth solutions and compare numerical results with analytical and measured data and also with the results obtained earlier using some other schemes.

In all tests, unless otherwise mentioned, we assume that water flow is frictionless, the CFL number is 0.5, the acceleration of gravity g = 9.81, and $\Delta x_{j} = 0.005$.

4.1 Accuracy of the Scheme in Smooth Regions

We expect that the scheme (51), (52), (54) has first-order accuracy in time and second-order accuracy in space. We consider the test proposed in (Balbás and Hernandez-Duenas 2014) to check second-order accuracy of the scheme in space for smooth flows by evaluation of L¹ error over the successively thinner grids.

Let a flat channel have a trapezoidal geometry

$$\sigma (x,y) = 1 + 0.3y$$

Initial water surface and velocity are described by

$$w_{0} (x) = 1.6 + 0.1 \, cos\left( {\frac{\pi (x - 0.4)}{0.2}} \right),\quad u_{0} (x) = 1, \quad 0 \le x \le 1$$

The outflow boundary conditions are set at the channel ends.

We compute L¹ error for $\Delta x = {1 \mathord{\left/ {\vphantom {1 N}} \right. \kern-0pt} N}$, N = 80,160,320,640,1280 and 2560 at the final time T = 0.05. The reference solution was computed with N = 5120. In order to suppress the scheme error in time, we provide all calculations with temporal step $\Delta t = 1 \times 10^{ - 8}$, which is less than $(\Delta x)^{2} = ({1 \mathord{\left/ {\vphantom {1 {5120)^{2} }}} \right. \kern-0pt} {5120)^{2} }} \approx 0.38 \times 10^{ - 7}$. The L¹ error is shown in Table 1. The results indicate the second-order accuracy in space of the scheme.

Table 1 Accuracy test for the scheme in smooth regions. The errors at T = 0.05

4.2 Large Perturbation of Rest

The following test is also taken from (Balbás and Hernandez-Duenas 2014), but with another channel geometry. In this test, propagation of perturbation of a steady-state solution at the left part of the simulated area is considered. When the perturbation propagates, part of it reflects back and part is transmitted, inundating the right side and leaving the area through the right boundary. According to (Balbás and Hernandez-Duenas 2014), the bottom topography is described by a piecewise spline. At the beginning, we construct the cubic spline of defect one with knots (0,0.3), (0.05,0.3), (0.1,0.2), (0.15,0.5), (0.3,0.4), (0.4,0.6), (0.75,0.6) and zero-values of the second derivatives at the boundaries. Then the spline is rescaled by a factor of 0.5 for $x \le 0.53$ and of 1.3 for $x \ge 0.53$. Finally, we assume that near the boundaries the bottom is flat and equals to 0.12 for $x \le 0.0693$ and 0.8 for $x \ge 0.7386$.

The geometry is given by (Fig. 6)

$$\sigma (x,y) = \left( {1 + \frac{3}{4}\cos (\pi x)} \right)\left( {1 - \frac{y}{2}} \right), \quad 0 \le x \le 1, \quad 0 \le y \le 2$$

(59)

Fig. 6

3D view of the channel geometry

We consider the left boundary as open (Antonopoulos and Dougalis 2017), (Shiue et al. 2011), and outflow condition is set up on the right boundary.

The initial water surface elevation at rest is $w = 0.8$, and a piecewise-constant perturbation from water surface of size $\varepsilon = 0.3$ is applied on the interval [0.1,0.15]. Propagation of the perturbation at different times is shown in Fig. 7. The wave partially reflects back and partially transmits through the bottom ridge, overtopping it and then propagating through the shore. These computed results are qualitatively similar to results obtained in (Balbás and Hernandez-Duenas 2014).

Fig. 7

Propagation of the perturbation at different times for the conditions of the test proposed in (Balbás and Hernandez-Duenas 2014)

4.3 Convergence to a Smooth Subcritical Steady State

In this example from (Hernandez-Duenas and Beljadid 2016) the topography is described by the cubic spline of defect one with knots (0.2,0), (0.3,0.6), (0.4,0.4), (0.5,0.5), (0.6,0.2), (0.7,0) and zero-value of the second derivatives at the boundaries. The channel geometry is given by (59).

The initial data is

$$w(x,0) = 0.8 \quad {\text{and}} \quad u(x,0) = 0$$

At the left (inflow) boundary, the discharge $Q_{in} = 0.3343$ and at the right (outflow) boundary the water surface elevation $w = 0.8$ are specified. Hernandez-Duenas and Beljadid (2016) calculated that the subcritical steady-flow with the given conditions has constant energy E = 10.0748. In our simulations, we obtained almost the same value of the energy $E = 10.0307$ with accuracy to four decimal places.

Comparison of computed water surface elevation, discharge, and energy for the smooth subcritical steady state with their exact values are shown in Figs. 8 and 9. In Fig. 9, the vertical axis extends by 0.15% of the exact value in each direction. We also, as in (Hernandez-Duenas and Beljadid 2016), calculated the maximal error between computed and exact solutions, and relative error in the L² norm defined as

$$err = \sqrt {\frac{1}{b - a}\int_{a}^{b} {\left( {\frac{{f(x) - f_{exact} (x)}}{{f_{exact} (x)}}} \right)^{2} dx} } , \quad x \in [a,b]$$

(60)

Fig. 8

Left: Comparison of computed and exact water surface elevation for the smooth subcritical steady state. Right: 3D view of the channel geometry

Fig. 9

Comparison of computed and exact discharges and energies for the smooth subcritical steady state

We obtained that for the discharge, the maximal error is $3.82 \times 10^{ - 4}$ and the relative error is $1.84 \times 10^{ - 4}$. For the energy, the maximal error is $5.67 \times 10^{ - 4}$ and the relative error is $2.15 \times 10^{ - 5}$.

4.4 Convergence to a Transcritical Steady State With Shock Wave

This test is also taken from (Hernandez-Duenas and Beljadid 2016). The topography is given by

$$B(x) = \;\left\{ {\begin{array}{*{20}l} {\frac{1}{2} + \frac{1}{2}\cos \left( {\pi \frac{x - 0.5}{0.4}} \right),} \hfill & {x \in [0.5,0.9]} \hfill \\ {0,} \hfill & {x \in [0,1]\backslash [0.5,0.9]} \hfill \\ \end{array} } \right.$$

The geometry is

$$\sigma (x,y) = \;\left\{ {\begin{array}{*{20}l} {\frac{1}{2} + \left[ {\frac{3}{8} - \frac{1}{8}\cos \left( {\pi \frac{x - 0.7}{0.2}} \right)} \right]\sqrt y ,} \hfill & {x \in [0.5,0.9]} \hfill \\ {\frac{1 + \sqrt y }{2},} \hfill & {x \in [0,1]\backslash [0.5,0.9]} \hfill \\ \end{array} } \right.$$

The boundary conditions are Q_in= 2.5561 and w_out= 1.9968.

Figure 10 shows that a good agreement is obtained between the analytical solution and the computed water surface elevation. A comparison of the computed discharge and energy with the theoretical results is shown in Fig. 11.

Fig. 10

Left: Water surface elevation for steady transcritical flow over a bump with a shock. Right: 3D view of the channel geometry

Fig. 11

Water discharge and energy for steady transcritical flow over a bump with a shock

4.5 Dam-Break Problem in a Triangular Channel

In this test (Lai and Khan 2012), a dam-break problem in a frictionless, horizontal, triangular channel is considered. The channel is 1000 m long with a side slope of 1H:1V. The dam is located in the middle of the channel. Both wet-bed and dry-bed conditions downstream of the dam are simulated. The upstream water depth was 1 m for both cases, and the downstream water depth was 0.1 m for the wet-bed case. For the wet-bed case, we compute the solution with 400 computational cells. For the dry-bed, we used two computational grids with 400 and 1000 cells.

The exact solutions of dam-break problems for a horizontal triangular channel can be found in Chen et al. (2011), Wu et al. (1993), Wu et al. (1999)). A comparison of the numerical and exact solutions for the wet-bed dam break at 80 s after dam removal is shown in Fig. 12. The dry-bed numerical results for water surface and flow rate at 45 s after dam removal are given in Figs. 13 and 14. In the dry-bed case, one can see that the greatest error is at the front of the moving water (Fig. 15), which decreases as the cell size is reduced. These results are similar to that obtained by Sanders (2001) and Lai and Khan (2012) using TVD schemes of second-order accuracy.

Fig. 12

Comparison of water surface and flow rate for wet-bed dam break in a triangular channel (N = 400, t = 80 s)

Fig. 13

Comparison of water surface and flow rate for dry-bed dam break in a triangular channel (N = 400, t = 45 s)

Fig. 14

Comparison of water surface and flow rate for dry-bed dam break in a triangular channel (N = 1000, t = 45 s)

Fig. 15

Water surface for dry-bed dam break in a triangular channel for different cell sizes (t = 45 s)

4.6 Drain on a Non-Flat Bottom

In this example, taken from Balbás and Karni (2009), Hernandez-Duenas and Karni (2011), a symmetric reservoir is being drained through a rectangular channel with a parabolic contraction. Due to the symmetry, the flow is computed on half the domain. The contraction is described by the quadratic interpolant through the points (0.25,1.0), (0.5,0.8), and (0.75,1.0), where the first number is the x-coordinate, and the second one is the width of the channel. The bottom topography consists of one hump (Fig. 16)

$$B(x) = \;\left\{ {\begin{array}{*{20}l} {0.25\left[ {1 + \cos (\pi (x - 0.5)/0.1)} \right],} \hfill & {{\text{if }}\left| {x - 0.5} \right| < 0.1} \hfill \\ {0,} \hfill & {\text{otherwise}} \hfill \\ \end{array} } \right.$$

Fig. 16

3D-view showing the topography at the bottom and lateral channel walls

Due to symmetry of the domain at x = 0, we specify wall boundary conditions on the left boundary. The right boundary condition is an outflow condition on a dry bed. This boundary condition on the right side of the domain allows the water that was at rest to flow freely through the right boundary into the initially dry region.

In Fig. 17, we show the solution obtained for the initial water surface elevation at 0.8.

Fig. 17

Water surface elevation under drain on a non-flat bottom at different times for the test (Balbás and Karni 2009, Hernandez-Duenas and Karni 2011) (N = 200)

After drainage begins, the solution converges to a steady-state solution in which water exists only on left-hand side of the hump. The obtained simulated results are in good qualitative agreement with numerical results in Balbás and Karni (2009) for all times except t = 1. The numerical results for t = 1 differ from Balbás and Karni (2009) by the presence of small waves. We carried out additional numerical computations with N = 100 and N = 400. The obtained results are illustrated in Fig. 18 and anyone can see the absence of the waves on the rough computational grid and their presence on the detailed grid.

Fig. 18

Water surface elevation under drain on a non-flat bottom at t = 1 for the test (Balbás and Karni 2009, Hernandez-Duenas and Karni 2011) for different number of cells (left N = 100 and right N = 400)

4.7 Subcritical Dam-Break Flow in a Rectangular Channel With a Junction

We consider a subcritical dam-break flow in a frictionless, horizontal, rectangular channel, 34 m long and 3 m wide. The channel is divided into two parts by a dam located 15 m from the left end. Initially, the water at rest to the left of the dam was at a level of 0.5 m and 0.1 m to the right. Vertical walls are placed at the ends of the channel, and their height is enough to prevent any spilling of water out of the channel (Fig. 19). The channel cross-sections for numerical solution comparisons G1 and G2 were respectively located 4.4 m and 5.9 m downstream from the dam.

Fig. 19

Outline of dam-break flows in a rectangular channel with a junction

After abrupt dam failure, the water flows downstream until it reaches the vertical wall at the right channel end and reflects from it. The left vertical wall reflects the flow too. The simulations are carried out for a uniform $\Delta x = 0.1$ m and $\Delta t = 0.01$ s. Numerical simulations are performed for a computational grid with and without channel junction that is located 5 m downstream from the dam. We apply two types of junction treatment for a control volume around the channel junction: (1) based on mass conservation (CUJ1) and (2) based on the Saint–Venant equations (CUJ2). We consider the numerical solution on the computational grid without junction (CU) as the reference solution.

Simulated water profiles for the two times, at 5 s and 18 s, are given in Fig. 20. For different junction treatments, comparisons of the simulated water depths at two channel cross-sections G1 and G2 are shown in Fig. 21. Note that the numerical solution with CUJ1 (red line in Figs. 20, 21) locally is slightly different from the reference solution CU (black line). Therefore, applying only the mass conservation equation for a channel junction can lead to small errors in the modeling of a subcritical open water flow.

Fig. 20

Simulated water surface profiles of the entire channel at t = 5 s and 18 s

Fig. 21

Comparisons of the simulated water depths and discharges with different junction treatments at selected cross-sections

4.8 Dam-Break Flow in a Rectangular Channel With a Loop

The purpose of this test is to compare different treatments of an open channel junction for subcritical and supercritical flows. Similar to the previous section, we consider a dam-break problem in a frictionless, horizontal, rectangular channel with a loop (Fig. 22). There are two vertical walls at the ends of the channel located at x = − 15 m and x = 19 m. The channel width before and after the loop is 3 m. The width of the lower loop channel (channel B in Fig. 22) is 2 m, and the upper (channel C) is 1 m. The channel is divided into two parts by a dam at x = 0 m. Initially, the upstream water depth is 1 m (supercritical flow), or 0.5 m (subcritical flow), and the downstream water depth is 0.1 m. The dam is then breached, and instantaneously water discharges from the higher to the lower level as a downstream-directed bore while a depression wave propagates upstream. The flow is reflected when the water front reaches the walls at the ends of the channel. Four channel cross-sections for comparison of numerical solutions G1, G2, G3, and G4 were located at different parts of the channel as shown in Fig. 22.

Fig. 22

Sketch of modeling area for dam-break problem in a rectangular channel with a loop

The simulations are performed with $\Delta x = 0.1$ m and $\Delta t = 0.01$ s. The 2D model COASTOX (Zheleznyak et al. 1992), (Kivva et al. 2018), (Zheleznyak et al. 2016) and 1D model are used to simulate dam-break problem. COASTOX is an unstructured finite volume model for solving the shallow water equations by using the Godunov scheme. COASTOX numerical solutions are the reference for comparison of the 1D numerical solutions. In 1D simulations, the junction treatment in the channel network is represented by two kinds of models: (1) mass conservation (CUJ1), and (2) mass and momentum conservation (CUJ2). These two junction models are compared at the four cross-sections with each other and with respect to the 2D numerical solution that averaged over the cross-section.

Comparison of the simulated water depths by using two junction models at the four cross-sections for the subcritical flow is presented in Fig. 23. Relative errors in L¹ norm and L² norm calculated by relation (60) between 1D and 2D solutions at the cross-sections for the subcritical flow are given in Table 2.

Fig. 23

Comparison of the simulated water depths at channel cross-sections G1, G2, G3, and G4 for the subcritical flow

Table 2 Relative errors in L¹ and L² norms between 1D and 2D solutions at the cross-sections for the subcritical flow

CUJ1 produces similar relative errors as CUJ2 does for the subcritical flow in this test. At the same time, the results of CUJ1 much better describe the dynamics of the first wave in the cross-section G2. For the supercritical flow in a channel network CUJ1 does not produce a physical numerical solution. Comparison between CUJ2 and COASTOX numerical solutions for the supercritical flow at the four gages is shown in Fig. 24. Agreement between 1D and 2D numerical results for subcritical and supercritical flows is generally satisfactory when only the CUJ2 approach is used for the supercritical flows.

Fig. 24

Comparison of the simulated water depths by CUJ2 and COASTOX at the channel cross-sections G1, G2, G3, and G4 for the supercritical flow

As one might suppose, the greatest error in the numerical solution of a 1D model in comparison with a 2D model is observed at the bore front under passing of the channel junction. Therefore, it is better to apply 2D or 3D models to simulate a bore propagation in a multiply-connected channel network.

4.9 Lower Danube River

As a test of the applicability of the above-presented numerical scheme for a large-scale natural open flow, we simulate the river network of the Lower Danube from 01/01/2000 to 31/12/2000 (Zheleznyak et al. 2016). The river network consists of 29 links and 26 nodes (Fig. 25, left). Blue and green triangles denote inflow nodes, yellow squares indicate outflow nodes and red circles mark channel junctions. Water gages D1–D10 are located at points indicated by red circles shown in Fig. 25, right. The total length of the main channels is more than 900 km. The total length of the main channels between water gages D1 and D10 is more than 500 km. For this part of the river network, we have only 10 channel cross-sections located at the water gages.

Fig. 25

River network of the Lower Danube River. Left: the inflow nodes of the river network are indicated by triangles (blue and green), and outflow nodes are indicated by yellow squares. Red circles denote the river channel junctions with their node numbers. Right: red circles denote the location of the water gage stations D1–D10. Yellow squares denote the junction nodes in which the water surface levels calculated by the CUJ1 and CHARIMA models are compared

The river network is discretized by a non-uniform grid with cell sizes varying from 1061.37 to 3948.45 m. The computational grid includes 321 cells and 338 cell interfaces. The CFL number is set to 0.8. To simulate water flow in the river network, we will use two models: (1) the described-in-this-paper central-upwind scheme with the junction treatment based on the continuity equation (CUJ1) and (2) an analog of the CHARIMA model (Holly et al. 1990). CHARIMA simulates water flow in a channel by applying the Preissmann implicit finite-difference scheme (Preissmann 1961). This scheme is a weighted four-point scheme. At a junction node, it is assumed that: (1) inflow discharges should be equal outflow discharges from all tributaries at the junction, and (2) the water levels at the ends of linked channels are equal at the junction.

The friction slope is evaluated from Manning’s Eq. (3). The hydraulic radius is calculated as A/P, where P is the wetted perimeter. Note that CUJ1 and CHARIMA have different difference approximations of the friction slope.

Water discharge is specified at the inflow nodes. At that, the water discharge at a node indicated by green in Fig. 25 (left) is simulated as a lateral inflow to the corresponding channel junction. Water surface elevation is set up at the outflow nodes (yellow squares).

The calibration of CUJ1 is performed by adjusting Manning’s coefficient n inside the river network by using the CMA Evolution Strategy (Suttorp et al. 2009). We assumed that $0.01 \le n \le 0.04$, and for its evaluation we adopted pCMALib code (Müller et al. 2009). The Nash–Sutcliffe model efficiency coefficient (NSE) is applied to assess the approximation of the measured water surface elevation at a water gage by the simulated water surface elevation. We took data of water surface level measurements at all water gages D1–D10 for the Manning’s coefficient evaluation. Thus, the objective functional is taken as

$$\Im = \sum\limits_{i = 1}^{10} {NSE_{i} = } \sum\limits_{i = 1}^{10} {\left[ {1 - \frac{{\sum\nolimits_{j = 1}^{m} {\left( {w_{i} (t_{j} ) - w_{i}^{0} (t_{j} )} \right)} }}{{\sum\nolimits_{j = 1}^{m} {\left( {w_{i}^{0} (t_{j} ) - \bar{w}_{i}^{0} } \right)} }}} \right]} \to \hbox{max}$$

where $\bar{w}_{i}^{0}$ is the mean of observed water surface elevations at an i-th water gage; $w_{i}^{0} (t_{j} )$ is the observed water surface elevation at an i-th water gage in time $t_{j}$; $w_{i} (t_{j} )$ is the modeled water surface elevation at an i-th water gage in time $t_{j}$; m is the number of times. When the functional $\Im$ reached a value of 9.5527, we interrupted the calibration. Values of the Nash–Sutcliffe coefficients at gages D1–D10 are given in Table 3.

Table 3 Values of the NSE at water gages D1–D10

Post-calibrated values of Manning’s n were used in simulations of the two models under the same computational grid, initial and boundary conditions. Comparisons of simulated and observed water surface elevations at the Borcea (D3), Izvoarele (D4), Gropeni (D8) and Braila (D10) water gages are shown in Fig. 26. Comparisons of the water surface elevations at junctions 3, 8, 17 and 22 (Fig. 25, right) calculated by the CUJ1 and CHARIMA are presented in Fig. 27.

Fig. 26

Comparison of simulated and observed water surface elevations at various gage stations

Fig. 27

Comparison of the water surface elevations at junctions calculated by CUJ1 and CHARIMA

A good agreement between the water surface elevations at the channel junctions calculated by the CUJ1 and CHARIMA models indicates that the continuity equation at a junction for a subcritical flow can be used instead a traditional internal boundary condition based on equality of the both discharges and water levels. That provides the possibility to apply the described central-upwind scheme CUJ1 for a multiply-connected channel network for subcritical flows through junctions.

4.10 Inundation of a Dry Channel Network

In this section, we simulate inundation of a dry channel network which consists of 9 inclined rectangular channels and four junction nodes (Fig. 28). We will denote each channel by the node names bounding this channel. In a channel name, the left letter will mark the upper end of the channel. Upper ends of channels AB, EF and HG are bounded by vertical walls. Characteristics of network channels are given in Table 4. Note that in the junction F channel ends connected with it have different elevations.

Fig. 28

Sketch of a channel network

Table 4 Characteristics of network channels

At the upper ends of the channels AB, EF and HG, there are water sources, rates of water inflow (m²/s) of which in the continuity equation are calculated by formulas

$$Q^{A} (t) = \;\left\{ {\begin{array}{*{20}l} {0.01\,\sin (\pi {{(t - 1)} \mathord{\left/ {\vphantom {{(t - 1)} {60}}} \right. \kern-0pt} {60}})} \hfill & {{\text{if }}1 \le t \le 61} \hfill \\ 0 \hfill & {\text{otherwise}} \hfill \\ \end{array} } \right.$$

$$Q^{E} (t) = \left\{ {\begin{array}{*{20}l} {0.0075\,\sin (\pi {{(t - 1)} \mathord{\left/ {\vphantom {{(t - 1)} {30}}} \right. \kern-0pt} {30}})} \hfill & {{\text{if }}1 \le t \le 31} \hfill \\ {0.05\,\sin (\pi {{(t - 150)} \mathord{\left/ {\vphantom {{(t - 150)} {60}}} \right. \kern-0pt} {60}})} \hfill & {{\text{if }}150 \le t \le 210} \hfill \\ 0 \hfill & {\text{otherwise}} \hfill \\ \end{array} } \right.$$

$$Q^{H} (t) = \left\{ {\begin{array}{*{20}l} {0.0075\,\sin (\pi {{(t - 50)} \mathord{\left/ {\vphantom {{(t - 50)} {80}}} \right. \kern-0pt} {80}})} \hfill & {{\text{if }}50 \le t \le 130} \hfill \\ 0 \hfill & {\text{otherwise}} \hfill \\ \end{array} } \right.$$

Such geometry of the channel network and the time distribution of the water source discharges were selected for this test to generate the consequence of the waves inundating and drying the network.

The channel network was discretized by a uniform grid with cell size of 0.2 m. Manning’s roughness coefficient was equal to 0.01 (s/m^1/3). The outflow boundary condition was specified at the node D. Numerical simulations were carried out for the CFL number equal to 0.9. Since flow along a dry bottom is supercritical we used the CUJ2 model to simulate open water flow in the channel network under inundation.

Open water surface elevations along various channels of the network at different times are shown in Figs. 29, 30, 31, 32 and 33. Time series of water surface elevations and the Froude number at the junction nodes B, C, F and G are given in Fig. 34, 35. In Figs. 29 and 30, anyone can see that during the passing of first wave along channel EF inflow discharge of water to the junction node F does not equal to outflow discharges from it, which are zero at this time. Thus, the condition that the sum of the discharges has to be zero at the junction, which is used in many hydrological models as the internal boundary condition for the junction treatment, is not valid in general case.

Fig. 29

Open water surface elevations at different times along channels EFB

Fig. 30

Open water surface elevations at different times along channels EFC

Fig. 31

Open water surface elevations at different times along channels HGB

Fig. 32

Open water surface elevations at different times along channels HGC

Fig. 33

Open water surface elevations at different times along channels ABCD

Fig. 34

Time series of open water surface elevations at the junction nodes

Fig. 35

Time series of the Froude number at the junction nodes

5 Conclusions

A central-upwind scheme has been applied to simulate shallow water flow in a multiply-connected channel network with arbitrary channel geometry and bottom topography. A new reconstruction algorithm of water surface elevation for partially flooded cells has been proposed. This reconstruction algorithm, that is a generalization of the idea from (Bollermann et al. 2013), with the exact integration of source terms of the shallow water equations, provides the well-balanced property and positivity preserving of the scheme. We considered two models of a channel junction treatment based on: (1) the continuity equation for a subcritical flow and (2) mass and momentum conservation equations for a supercritical flow. Using such channel junction treatment instead of the traditional internal boundary conditions allows a simpler way to apply other difference schemes instead of the Preissmann scheme and simplified forms of the Saint–Venant equations to simulate open water flow in a multiply-connected channel network. It is shown that for a subcritical flow the continuity equation at a channel junction is a generalization of the model of the equality of water surface elevations. Applying the local draining time approach from (Bollermann et al. 2011) to limit outflowing flux from the draining cell, we provide the positivity preserving without of a reduction of the CFL time step restriction. Implicit treatment of only a part of the friction slope according to Chertock et al. (2015) holds the scheme stability without additional time step restriction.

The set of the numerical tests demonstrates the scheme accuracy, positivity preserving and well-balancing, convergence of numerical results to steady-state solutions and their good agreement with exact solutions and experimental data, including measurements in the multiply-connected river channel network. We propose the new specialized test for the simulation of inundation and drying of the channel network by consequence of supercritical and subcritical flows. The test results demonstrate the stability of the scheme and the robustness of the new numerical algorithm for the simulation of the wetting/drying flow in the multiply-connected channel network.