A central-upwind scheme for open water flow in a wet/dry multiply-connected channel network

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 momentum conservation equation provide positivity preserving and well-balanced features of the scheme for various wet-dry states. We use two models based on 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 a subcritical and supercritical flow 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 initially dry channel network shows the merits of the new numerical algorithm to simulate the subcritical/supercritical open water flows in the network.


Introduction
The Saint Venant equations are 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 the hydraulic engineering (e.g. HEC-RAS [1], ISIS 1D [2], FEQ [3], CHARIMA [4], NETSTARS [5], FLDWAV [6] ) the Preissmann [7] 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 an efficient algorithm for the numerical solution of the Saint Venant equations in a multiply connected channel network. One of the disadvantages is that it 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 rivers flow. Steep bottom slopes can generate artificial "drying" and following crash of the calculations due to the negative values of depth in numerical solutions if Preissmann scheme or other schemes without the positivity preserving are used for modeling.
Many difference schemes for hyperbolic equations that preserve non-negativity of water flow depth were developed after the pioneer publication of Godunov [8]. We refer the reader to [9], [10], [11], [12], [13], [14], [15], [16], [17], [18] where different approaches are presented. Among them, central [19] and central-upwind [20] 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 [21], [22], [23]. Balbas and Karni [21] presented a central finite-volume scheme of second order accurate to simulate one-dimensional shallow water flows along channels with non-uniform rectangular cross-sections and bottom topography. In [22] Balbas and Hernandez-Duenas 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 [23] for shallow water flows in channels with arbitrary geometry and variable topography.
In this paper, we implemented a central-upwind scheme [24] for open water flow in multiplyconnected channel networks as a practical alternative to the Preissmann scheme. In the framework of the second order central-upwind schemes, we propose a novel reconstruction of water surface elevation that not requires additional correction in partially flooded cells. This reconstruction algorithm is a generalization of the reconstruction from [25].
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 momentum conservation equation includes terms that account for the forces due to changes in channel width and bed elevation. Assuming that the water and bed surface elevations and channel cross-section depend linearly on the water and bed elevations, and channel cross sections defined at cell faces, we can exactly calculate the source terms integrals on any internal cell interval. The exact values of the source terms integrals with our proposed reconstruction will guarantee us 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.
A treatment of channel junction in the modeling of open water flow in a multiply-connected channel network is one of the challenges. In many used hydraulic engineering codes, it is assumed that two groups of hydraulic conditions are held at a channel junction [26]. One is that the sum of the discharges has to be zero at the junction, and 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 [27], RIVICE [28] and ONE-D [29] use the equality of water surface elevations at the junction. Others as HEC-RAS and FEQ are energy based models. Usually, these hydraulic conditions are specified as interior boundary conditions. In general case, the assumption of equality of inflow and outflow discharges at a junction does not permit in principle to simulate open water flow in a channel network under wetting/drying conditions. Yoshioka et al. [30,31] use as the internal boundary conditions at the junction the continuity equation discretized on a dual mesh around it and the momentum flux at junction calculated as weighted linear combination of the momentum fluxes inflowing into it. Sanders et al. [32] and Bellamoli et al. [33] simulated open water flow in channel networks by nesting a 2D model at junctions. Contarino et al. [34] 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 equations for subcritical flow, and 2) based on mass and momentum conservation equations for supercritical flow. The first junction model is an analog of the equality of water surface elevations model.
Constraints on 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 [25], [35], we limit the outflow fluxes from 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 paper is organized as follows. Channel network definition is presented in Section 2. Open flow water equations in channels are given in Section 3. In Section 4.1 we shortly describe the centralupwind scheme from [24] implemented for the Saint Venant equations. In Section 4.2, we present relations for exact discretization of source terms in the Saint Venant equations. New reconstruction method we describe in Section 4.3. Slope friction treatment and channel junction models are given in Section 4.4 and Section 4.5, respectively. Boundary conditions, which usually applied in modeling a subcritical open flow in a channel, are presented in Section 4.6. Positivity preserving and wellbalancing properties of the scheme are proven in Sections 4.7 and 4.8, respectively. Section 5 contains different numerical tests illustrated the merits of the scheme. Section 6 provides some concluding remarks.

Channel network definition
The definition sketch of a channel network is shown in Fig.1 that is similar to CHARIMA [4]. Any channel network we will consider 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.  The link represents 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 nl at the lower end of the link. Remark that "lower" and "upper" generally, although not necessarily, mean "downstream" and "upstream", respectively. Each link should have at least one computational cell. Natural channel cross-section is replaced by its piecewise linear approximation (Fig.2) which assumed to be known for the each computational cell edge.

Channel water flow equations
The governing equations of one-dimensional shallow water flows for natural channel with irregular cross-section are given by the Saint Venant equations: 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; Sf is the friction slope due to the bed resistance; S0 is the bed slope; I1 is the hydrostatic pressure force; I2 is the wall pressure force.
The friction slope Sf is assumed to be given by the Manning equation where n is the Manning's roughness coefficient; R is the cross-section hydraulic radius.
The hydrostatic pressure force, wall pressure force, wetted cross-section area and bed slope are defined, respectively, as where h is the water depth; ( , ) xh  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 have its own unique number. The numbering of channels does not depend on the numbering of nodes.
In general, 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 in [24] is used for numerical discretization of the Saint Venant equations (1)- (2).
We divide the spatial domain into the grid cells 1 computed at time t Then, a semi-discretization of (1)-(2) can be written as the following system of ODEs where  is a small apriori chosen positive number. For consistency, the values of the discharge at cell interfaces are recomputed using are given by (10).
The interface values 12 j U   are obtained from the cell averages j U by a piecewise linear reconstruction which will be described below.

4.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 I1 and the wall pressure force I2.
We replace the bed function B and wetted area A on the 1 ( , ) Assuming that the open water surface w is linear in a computational cell   , LR xx , below we will use the following relations. Water volume V(x1,x2;wL,wR) between points x1 and x2 (Fig.3.) equals to 22 11 ( , ) 12 0 Changing the order of integration, we obtain 12 ( , ; , ) w the depth of water which 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 0 max(0, ) hh   .
into equation (13) and taking into account that  (18) which is a quadratic equation since We divide conditionally all computational cells into two groups: "wet" and "dry". A cell j will be called "wet" if Depending on combination of two neighboring cells, we use different formulas for approximating the forward and backward difference derivatives for w and Q. We approximate backward difference derivative   by the following way ( Fig.4): 1. j-1 and j cell are both "wet". Then 2. j-1 cell is "wet" and j cell is "dry" 3. j-1 cell is "dry" and j cell is "wet" (1 ) 4. j-1 and j cells are both "dry"  if (1 ) The forward derivatives   for w and Q are calculated in the same way.
The interface point-values 12 j w  are obtained by a piecewise linear reconstruction , 12 , , , According to [19] this reconstruction will be second-order accurate if the approximate values of the derivatives , xj w are at least first-order approximations of corresponding exact derivatives. To ensure the non-oscillary property in our numerical scheme we evaluate , xj w using the minmod limiter [13], [36], [37] , min mod( , ) The interface values

Friction slope treatment and time evolution
We calculate the second component of (2) 43 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 [38]. The explicit treatment of the friction term imposes a severe time step restriction, as the result of which time step should be in 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. In result, we obtained the following fully discrete central-upwind scheme Note that the scheme (31)-(32) is first-order accurate in time. A higher-order time discretization can be obtained using semi-implicit solvers proposed by Chertok et al. [39] which are based on the modification of explicit SSP Runge-Kutta methods [40]. The second-order semi-implicit method for a system ODE with a non-stiff term f(u,t) and stiff damping term g(u,t)u such that g(u,t) is a diagonal matrix with non-positive elements ( , ) ( , ) can be written as

Node of junction
For modeling of the hydraulic conditions at a junction, we will use two approaches: (1) assuming that only mass conservation equation is hold at the junction, and (2) based on mass and momentum conservations.
Consider a node s that connects two and more channels, Js,in is a set of channels which inflow and Js,out that outflow from it (Fig.5). Denote , in si x and , out sj x 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 Ys , or water surface elevation Ys and water discharge Qs for the second case at the junction are also given. Assuming that water surface is horizontal, the continuity equation over the control volume around a junction s can be written as follow  are found from the solution of equations (33)- (34). The Newton-Raphson method is used to solve the nonlinear equation (33). Note that the derivative of can be expressed as

Boundary conditions
For the numerical solution of the ODEs (7), we should specify values of the central-upwind numerical fluxes 12 j H  at an upstream and if required downstream boundaries 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. Boundary treatment at the right boundary can be done by similar way.
Discharge boundary conditions. Let a water discharge Q(t) has been provided on the left boundary. We will consider the following approximation of the boundary numerical fluxes 12 H   The values of 12 Q  и 12 h  are found from a solution of the system of equations (36).
Surface boundary condition. Let a water surface elevation Y(t) is specified at the left boundary. Then we will use the following approximation of the boundary fluxes   h  are found from a solution of the equations (37).
We will also consider a more simple approximation of the surface boundary conditions. In this case, we assume that 1

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 outflow from a cell than it is at the moment. For the positivity, our result is similar to the one obtained in [22].
For any computational cell k of the link p, we define , pk t  as , , , For any node s, which has two or more links, we denote by s t  , , ,  (7)- (10), (33) with the piecewise linear reconstruction described in section 4.3 and the discretization of the source terms (14)- (16). Assume that the system of ODEs (7) (40)- (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-section area at the next time step.
Note, as follows from (9) The third and fourth terms on the right-hand side of (43) are nonnegative so that For any node s, we rewrite (33) in the following way , , , , , , 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 ratios 12 jj AA  . Due to bottom topography, irregular channel geometry and 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 [41], [25], [35]. The basic idea of this approach is to reduce time step locally only for the cell faces at which the inequalities (44), (46)-(47) do not necessarily hold.
Follow to [35] we introduce the draining time step ,, , , , , , , where ( (54), there have to be a balance between the flux gradient and the source term.

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 as also with the results obtained earlier using some other schemes.
In all tests we assume that water flow is frictionless, the CFL number is 0.5, the acceleration of gravity g=9.81, and 0.005 j x  , unless otherwise mentioned.

Accuracy of the scheme in smooth regions
We expect that the scheme (51)-(52), (54) have the first-order accuracy in time and second order in space. We consider test proposed in [22] to check the second order of the scheme in space for smooth flows by evaluation of L 1 error over the successively thinner grids.
Let a flat channel has a trapezoidal geometry

Large perturbation of rest
The following test is also taken from [22], but with another channel geometry. In this test, propagation of perturbation of a steady state solution at the left part of 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 [22], the bottom topography is described by a piecewise spline. At the beginning we construct the cubic spline of defect one with knots (0,0. x  . Finally, we assume that near the boundaries the bottom is flat and equals to 0.12 for 0.0693 x  and 0.8 for 0.7386 x  . The geometry is given by (Fig.7) 3 ( , ) 1 cos( ) 1 42 We consider the left boundary as an open [42], [43] and outflow condition is set up on the right boundary.  Fig.8. The wave partially reflects back and partially transmits through the bottom ridge, overtopping it and then propagating through the shore. These computed results are similar to results obtained in [22].

Convergence to a smooth subcritical steady state
In this example from [23]  are specified. Hernandez-Duenas and Beljadid in [23] calculated that the subcritical steady-flow with given condition has constant energy E=10.0748, in our simulations we obtained almost the same value of the energy 10.0307 E  with accuracy to four decimal places.  Comparison of computed water surface elevation, discharge, and energy for the smooth subcritical state with their exact values are shown in Fig.9 and Fig.10. In Fig.10 the vertical axis extends on 0.15% of the exact value in each direction. We also as in [23] 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

Convergence to a transcritical steady state with shock wave
This test is also taken from [23]. The topography is given by 1 Fig.11 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.12.   Fig.11. Left: Water surface elevation for steady transcritical flow over a bump with a shock. Right: 3D view of the channel geometry.

Dam break problem in a triangular channel
In this test [44] 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 is simulated. The upstream water depth was 1 m for both cases, and the downstream water depth was 0.1m for the wetbed case. For the wet-bed case, we compute 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 [45], [46], [47]. Comparison of numerical and exact solution for wet-bed dam break at 80 s after dam removal is shown in Fig. 13. The dry-bed numerical results for water surface and flow rate at 45 s after dam removal are given in Fig.14 and 15. In the dry-bed case, one can see that the most error is at the front of the moving water (Fig.16) which decrease as the cell size is reduced. These results are similar to that obtained by Sanders [48] and Lai and Khan [44] using TVD schemes of second order accurate.

Drain on a non-flat bottom
In this example, taken from [21], [49], a symmetric reservoir is being drained through a parabolic contracting rectangular channel. Due to the symmetry, the flow is computed on half the domain.
The left boundary is assumed a line of symmetry of the domain and wall boundary conditions are applied. 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 freely flows through the right boundary into the initially dry region. In Fig.18 we show the solution obtained for the initial water surface elevation at 0.8.  [21,49].
After drainage beginning 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 agreement with numerical result in [21].

5.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 with 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 from 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 solutions comparison G1 and G2 were respectively located 4.4 m and 5.9 m downstream from the dam.
[32] 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). The numerical solution on the computational grid without junction (CU) we consider as the reference solution.  Simulated water profiles for 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. Applying only the mass conservation equation for a channel junction can lead to small errors in the modeling of a subcritical open water flows.

5.8.Dam-break flow in a rectangular channel with a loop
Purpose of this test is to compare different treatments of an open channel junction for subcritical and supercritical flows. Similarly, 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 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), 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. s. 2D model COASTOX [50], [51], [52] 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, junction treatment in channel network is represented by two kinds of model: 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 averaged over the cross-section numerical solution of the 2D model.
Comparison of the simulated water depths by using two junction models at the four cross-sections for the subcritical flow is present in Fig.23. Relative errors in L 1 and L 2 norms (60) between 1D and 2D solutions at the cross-sections for the subcritical flow are given in Table 2.   Table 2. Relative errors in L 1 and L 2 norms between 1D and 2D solutions at the cross-sections for the subcritical flow.
CUJ1 produce in similar but the smaller relative errors than CUJ2 for the subcritical flow in this test. The same time the results of CUJ1 much better describe the dynamics of the first wave in the section G2.
For the supercritical flow in a channel network CUJ1 produces not 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 CUJ2 approach is used for the supercritical flows. Fig.24. Comparison of the simulated water depths by CUJ2 and COASTOX at the channel crosssections G1, G2, G3, and G4 for the supercritical flow.
As one might suppose, 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.

5.9.Lower Danube River
As the test of the applicability of the presented above numerical scheme for the large-scale natural open flow we simulate the river network of the Lower Danube from 01/01/2000 to 31/12/2000 [52]. The river network consists of 29 links and 26 nodes (Fig.26, 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.26 (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 varied 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. We will use to simulate water flow in the river network two models: 1) the described in this paper central-upwind scheme with the junction treatment based on continuity equation (CUJ1) and 2) analog of the CHARIMA model [4]. CHARIMA simulates water flow in a channel by applying the Preissmann implicit finite-difference scheme [7]. 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 equation (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.26 (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 [53]. We assumed that 0.01 0.04 n  and for its evaluation we adopted pCMALib code [54]. The Nash-Sutcliffe model efficiency coefficient (NSE) is applied to assess 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      we interrupted the calibration. Values of the Nash-Sutcliffe coefficients at gages D1-D10 are given in Table 3.
Postcalibrated values of Manning's n were used in simulations of the two models under same computational grid, initial and boundary conditions. Comparison 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. Comparison 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 Table 3. Values of the NSE at water gages D1-D10.  A good agreement between the water surface elevations at the channel junctions calculated by CUJ1 and CHARIMA models indicates that 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 provide possibility to apply the described central-upwind scheme CUJ1 for a multiply-connected channel network for a subcritical flows through a junctions.

Inundation of a dry channel network
In this section we simulate inundation of a dry channel network which consists of eight inclined rectangular channels and four junction nodes (Fig.28). We will denote each channel by the node names bounded this channel. In 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.  Such geometry of the channel network and the time distribution of the water sources 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. The 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 equals to 0.9. We used 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 Fig.29-33. Time serious of water surface elevations at the junction nodes B, C, F and G are given in Fig.34. In Fig.29-30 anyone can see that during passing 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 in this time. Thus, the condition that the sum of the discharges has to be zero at the junction used in many hydrological models as internal boundary condition for junction treatment is not valid in general case.

Conclusions
A central-upwind scheme has been applied to simulate shallow water flow in a multiplyconnected channel network with arbitrary channel geometry and bottom topography. 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 [25], with the exact integration of source terms of the shallow water equations provide the well-balanced property and positivity preserving of the scheme. We considered two models of a channel junction treatment based on: 1) continuity equation for a subcritical flow and 2) mass and momentum conservation equations for a supercritical flow. It is shown that for a subcritical flow the continuity equation at a channel junction is generalization of the equality of water surface elevations model. Applying the local draining time approach from [35] to limit outflowing flux from 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 Chertok et al. [39] holds the scheme stability without additional time step restriction.