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

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.


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 onedimensional 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 doublesweep 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 centralupwind (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 crosssections 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 Godunovlike 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 highresolution 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), 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 multiplyconnected river networks, we implemented a centralupwind 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 centralupwind 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 centralupwind 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 3422 S. Kivva et al. Pure Appl. Geophys. 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. (2015Yoshioka et al. ( , 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 longwave 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.

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 Vol. 177, (2020) Open Water Flow in a Wet/Dry Multiply-Connected Channel 3423 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. 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.
One-dimensional shallow water flow in a natural channel with irregular cross-sections is described by the Saint-Venant equations (Cunge et al. 1980): 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 0 is the bed slope; I 1 is the hydrostatic pressure force; I 2 is the wall pressure force.
The friction slope S f is assumed to be given by the Manning equation 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) where h is the water depth; rð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).

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.

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 x jÀ1=2 ; x jþ1=2 Â Ã of length Dx j , where x j is the center of a grid cell, and denote by U j ðtÞ the cell averages of the solution U ¼ A; Q ð Þ T of (1) and (2) computed at time t Then, a semi-discretization of (1) and (2) can be written as the following system of ODEs where H jþ1=2 are numerical fluxes at the cell interfaces x jAE1=2 , SðU; xÞ ¼ 0; gI 2 À gAS 0 ð Þ T and GðU; xÞ ¼ 0; ÀgAS f À Á T .
The central-upwind numerical fluxes H jþ1=2 are given by T are the right/leftsided values of the piecewise linear reconstruction U at x ¼ x jþ1=2 ; FðU; BÞ ¼ Q; Q 2 A þ gI 1 À Á T is the flux at the cell interface, and the one-sided local speeds are obtained using the eigenvalues of the Jacobian Here u ¼ Q=A and c ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi gA=r T p , and r 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) where e is a small a priori chosen positive number. For consistency, the values of the discharge at cell interfaces are recomputed using Q AE jAE1=2 ¼ A AE jAE1=2 u AE jAE1=2 , where u AE jAE1=2 are given by (10). The interface values U AE jþ1=2 are obtained from the cell averages U j by a piecewise linear reconstruction which will be described below.

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 1 and the wall pressure force I 2 .
We replace the bed function B and wetted area A on the interval x jÀ1=2 ; x jþ1=2 Â Ã by their piecewise linear approximation Vol. 177, (2020) Open Water Flow in a Wet/Dry Multiply-Connected Channel Assuming that the open water surface w is linear in a computational cell x L ; x R ½ , below we will use the following relations. Water volume V(x 1 , x 2 ; w L , w R ) between points x 1 and x 2 (Fig. 3) equals to Changing the order of integration, we obtain where 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 1 at a point x 1 can be calculated as Integration of the wall pressure force I 2 and AS 0 over an interval x 1 ; x 2 ½ yields the following expressions Remark 1. Note that for linear approximation of the water surface and bottom on an interval x jÀ1=2 ; x jþ1=2 Â Ã , if the values of w AE jÇ1=2 and B jAE1=2 are known at the cell faces x jÇ1=2 , then it is easy to calculate x Ã (Fig. 3) and integrals (13)-(16) for any internal interval x 1 ; x 2 ½ 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 x jÀ1=2 ; x jþ1=2 Â Ã .

Reconstruction of the Interface Values U AE jþ1=2
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 ¼ min mod ðD À U j ; D þ U j Þ where D À U j and D þ 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 ! 0 ¼ maxð0; hÞ. Substituting w 1 ¼ w 2 ¼ w st;j into Eq. (13) and taking into account that h jþ1=2 À h jÀ1=2 ¼ B jÀ1=2 À B jþ1=2 ¼ ÀDB, we obtain the following equation for finding w st;j 1 2 Equation (17) is a quartic equation for the piecewise linear approximation of the channel crosssection, 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.
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 D À by the following way ( Fig. 4): 1. j -1 and j cells are both ''wet''. Then 4. j -1 and j cells are both ''dry'' Vol. 177, (2020) Open Water Flow in a Wet/Dry Multiply-Connected Channel The forward derivatives D þ for w and Q are calculated in the similar way.
The interface point-values w AE jÇ1=2 are obtained by a piecewise linear reconstruction 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) where min mod ða; bÞ ¼ 0:5ðsgnðaÞ þ sgnðbÞÞ min ð a j j; b j jÞ.

Friction Slope Treatment and Time Evolution
We calculate the second component G x jÀ1=2 GðU; xÞdx using the midpoint rule and desingularization procedure where e 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 jþ1=2 are the first and second component of the numerical fluxes (8).

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 in s;i and x out s;j 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.
Approximating that water surface is horizontal, the continuity equation over the control volume around a junction s can be discretized as follows Accordingly, the difference scheme for the equation of momentum conservation can be written in the form Figure 5 Sketch of a node that is the junction of channels Vol. 177, (2020) Open Water Flow in a Wet/Dry Multiply-Connected Channel 3429 The interface values A AE jþ1=2 are calculated from the water surface elevation, which is reconstructed in a similar way as described in Sect. 3.3. The interface values Q AE jþ1=2 are computed: in the first case as s and F Q s may be thought of as the residuals of the continuity and momentum conservation Eqs. (33)-(34), respectively. The water surface elevation Y s ðt þ DtÞ and the water discharge Q s ðt þ DtÞ 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 A s Y s ðt þ DtÞ ð Þ with respect to Y s ðt þ DtÞ can be expressed as where h ! 0 ¼ 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.

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.
j;1=2 À gI 1 ðx out s;j ; h s Þ À gI 2 ðx out s;j ; x j;1=2 Þ þ gB x ðx out s;j ; s;j ðx j;1=2 À x out s;j Þ 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 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.
We will also consider a simpler approximation of the surface boundary conditions. In this case, we assume that Outflow boundary condition For outflow, we use the following extrapolations on the boundary

Positivity Preserving
In this section, we show that the resulting centralupwind 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 Dt p;k as Dt p;k ¼ min For any node s, which has two or more links, we denote by Dt s Theorem 1 Consider the semi-discrete central-upwind scheme (7)-(10), (33) with the piecewise linear reconstruction described in Sect. 3.3 and 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 n j ! 0.
where Dt p;k and Dt 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 Vol. 177, (2020) Open Water Flow in a Wet/Dry Multiply-Connected Channel 3431 non-negativity of the cross-sectional area at the next time step.
Note, as follows from (9), that a þ jþ1=2 ! 0, a À jþ1=2 0, a þ jþ1=2 À u þ jþ1=2 ! 0, and u À jþ1=2 À a À jþ1=2 ! 0 for all j. Moreover, the following inequalities 0 1 are satisfied for all j. Then for any computational cell k of link p from (31), we have A nþ1 The third and fourth terms on the right-hand side of (43) are nonnegative, so that for A nþ1 p;k ! 0 it is sufficient if the following inequality is satisfied For any node s, we rewrite (33) It is clear that the left-hand side of (45) will be nonnegative if 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 AE jÇ1=2 .
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 Following ( The outflow of a cell is the inflow of its neighbor. Thus, to keep the conservativity of the scheme, we define Dt p;kþ1=2 as Dt p;kþ1=2 ¼ minðDt; Dt drain p;m Þ; 3.8. Well-balancing The system (1), (2) admits smooth steady-state solutions, satisfying as well as nonsmooth steady-state solutions. One of the most important steady-state solutions is a trivial stationary one (lake at rest) 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.3 and the discretization of the source terms (14)-(16). Assume that the numerical solution U(t n ) corresponds to the steady state at rest. Then U(t n?1 ) = U(t n ), that is, the scheme is well-balanced.
Proof As a result of the reconstruction, we have w Ç jAE1=2 ¼ w j and Q Ç jAE1=2 ¼ 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 On the other hand, from (14) to (16) we have which ensures a balance of source terms with the flux terms in steady states at rest.

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 Dx j ¼ 0:005.

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 1 error over the successively thinner grids.
We compute L 1 error for Dx ¼ 1=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 Dt ¼ 1 Â 10 À8 , which is less than ðDxÞ 2 ¼ ð1 . 5120Þ 2 % 0:38 Â 10 À7 . The L 1 error is shown in Table 1. The results indicate the secondorder accuracy in space of the scheme.
The initial water surface elevation at rest is w ¼ 0:8, and a piecewise-constant perturbation from water surface of size e ¼ 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).

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 and 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 2 norm defined as We obtained that for the discharge, the maximal error is 3:82 Â 10 À4 and the relative error is 1:84 Â 10 À4 . For the energy, the maximal error is 5:67 Â 10 À4 and the relative error is 2:15 Â 10 À5 .

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 x 2 ½0; 1n½0:5; 0:9

(
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.

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.

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) 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. 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.

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.
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 Dx ¼ 0:1 m and Dt ¼ 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).

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),

Figure 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 1 and L 2 norms between 1D and 2D solutions at the cross-sections for the 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 downstreamdirected 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. The simulations are performed with Dx ¼ 0:1 m and Dt ¼ 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 dambreak 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 1 norm and L 2 norm calculated by relation (60) between 1D and 2D solutions at the cross-sections for the subcritical flow are given in Table 2. 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.
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.

Lower Danube River
As a test of the applicability of the abovepresented 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.
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 finitedifference 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 n 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 where w 0 i is the mean of observed water surface elevations at an i-th water gage; w 0 i ð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 = 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.
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.   Vol. 177, (2020) Open Water Flow in a Wet/Dry Multiply-Connected Channel 3451 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.

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.
At the upper ends of the channels AB, EF and HG, there are water sources, rates of water inflow (m 2 /s) of which in the continuity equation are calculated by formulas 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.

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 wellbalanced property and positivity preserving of the scheme. We considered two models of a channel junction treatment based on: (1) the continuity  (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 wellbalancing, 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.
JST through the Project ''Strengthening of the Environmental Radiation Control and Legislative Basis for the Environmental Remediation of Radioactively Contaminated Sites in Ukraine'' and by the Japan Society for the Promotion of Science, Grant (B) (18H03389) ''Long-term dynamics of radiocesium in aquatic ecosystems of Fukushima and Chernobyl contaminated areas''. Sergii Kivva and Oleksandr Pylypenko were supported within the part of the research under Grant of the National Academy of Sciences of Ukraine No 367/415.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.