Fast and high-precision simulation of hydrodynamic and water quality process in river networks

To realize a fast and high-precision simulation of unsteady flows and water quality process in river networks, a 1-D mathematical modelling system for free-surface flows and water quality process is developed. In the hydrodynamic model, a θ semi-implicit method is used to discretize the free-surface gradient term in the momentum equation while the Eulerian-Lagrangian method is employed to solve the advection term. To achieve good mass conservation, the finite-volume method is used to discretize continuity equation. Moreover, the resulting hydrodynamic model is unconditionally stable with respect to the gravity wave speed and flow advection. In solving the scalar transport equation, the forward difference is used for discretization in time while the centre-difference is in space. A sub-cycling is introduced to divide the computational time step and improve the stability of water quality model. Therefore, large time steps are allowed in both the hydrodynamic and the transport models. The coupling of the branches of the river networks is solved using a predictioncorrection method. Then, a Gaussian pulse test and a rectangular wave test are employed to demonstrate the accuracy and the efficiency of the system, which will also provide powerful supports for ecological operations of cascade reservoirs in drainage basins.


Introduction
In river-basin planning and management, fast decisions are required to deal with the emergency events, such as extreme floods, emergency scheduling of reservoirs and water quality operation, in which fast and high-precision simulations of unsteady flows and scalar transport in channel networks are often involved. With the development and construction of engineering in drainage basin, the ecological problems have become more obvious. To obtain the governing variables in river channel networks which flow out of upstream reservoir, a coupled one dimensional mathematical modelling system is needed. The foundation of real-time simulations of coupled model to solve the scalar transport problems in large channel networks is the hydrodynamic model, which significant efforts have been devoted to develop an fast and high-efficiency solver since 1960s [1][2][3]. When acquiring the value of governing variables, the transport equation can be solved with the initial conditions of concentration in a finite difference scheme or finite element scheme or other corresponding methods.
Rather than explicit discretization methods, implicit discretization methods use larger time steps which require less runtimes. In general, when implicit methods are used to discretize the governing equations, the junction points are always put in the whole computational matrix which will be rich of iteration and cost much runtimes. Therefore, many significant efforts have been done to simplify the structure of computational matrix of channel networks. The special node-numbering method [4,5], the hierarchical solution method [1,2,6], and the multistep iterative local method [3] are three representative methods.
The above methods or schemes simplify the coupling structures of computation nodes and junction points. However, it still exists some problems that not all kind of networks can be solved. And the model stability is not unconditional which will depend on the initial conditions and selected time steps. For some complicated channel networks, much computer memory and runtime are still required. Hence, further studies on the handling method of simplification are necessary.
This paper presents a 1-D hydrodynamic and water quality coupled model for rapid simulation of free water surface flows and scalar transport problems in river channel networks. A prediction-correction solver is used for the proposed model which can adopt a large calculated step, and a sub-cycling of water quality part is used for improving the model stability. A θ semi-implicit method is used to discretize the hydrodynamic equations. A forward difference in time and a centre-difference in space are used for discretization of scalar transport equation. The model is tested using a looped channel network for hydrodynamic part and two kind of waves for water quality part.

1-D hydrodynamic mathematical equations
Based on the Saint-Venant principle, the Governing equations of 1-D hydrodynamic equations can be described as combination of a continuity equation and a momentum equation: where A is area of a cross-section; B is width of a crosssection; Q is discharge; q is lateral inflow; t is time; x is distance along the river channel; η is water level; n is Manning's coefficient; g is acceleration of gravity; and R is hydraulic radius.

1-D water quality mathematical equations
The general governing equation of 1-D water quality is defined as: where C is concentration of water quality component; E is longitudinal comprehensive diffusion coefficient; S in is internal source term; S R is the biochemical reaction term; and S c is point source or surface source term.

Disposition of river channel network
The river channel network is composed of a quantity of branches and junction points, which is presented in Fig.1(a). For a junction point (JP), a polygonal unit is generated and defined as a JP unit, for which connects arbitrary branches in the meanwhile. For non-JPs, to divide one branch into several reaches which is to set the cross-section in the centre, and to define each reach as an unit called water level control body (WLCB). Correspondingly, an unit with two cross-sections is defined as discharge control body (DCB), as shown in Fig.1(b). One non-JP has two unit faces (two discharge sides or two water level sides). The first or last non-JP unit of each branch is defined as boundary unit which will be set to inflow or outflow boundary condition. To gain a geographical description of river channel networks, the coordinate of centre point of each cross-section in the horizontal plane is necessary to input the proposed model. In a JP unit, it usually exists more than two sides. In the computation, width of the side located in main stream or the largest side is considered as the JP unit width. The size of a JP unit is depended on the size of the channel width of the surrounding unit that belongs to the main stream side or the largest side. In the topological structure, a straggled grid is arranged by WLCB and DCB with corresponding variables(as shown in Figure.1). The variable η is defined at WLCB centre while the variable Q and the variable C are defined at DCB centre. The variable NB is defined as the number of branches of river channel networks, and the unit numbers of each branch is denoted by N i (i = 1,2,…, NB). The variable nwl and variable nd are denoted the quantities of WLCB (or units) and DCB (or sides). The variable i (i=0,…, nwl-1) is used to sign the indexes of WLCB while the variable i+1/2 (i=0,…, nwl-1) is for DCB. The index of an unit is unique in the global topology which is seriation from upstream to downstream and from main stream to tributaries.

Discretization of governing equations
In 1-D hydrodynamic model, the θ semi-implicit method [7] is adopted. For the momentum equation, operator splitting technique is used to discretized, and three steps are taken to solve the governing equations. First, to calculate all the explicit term which covers the advection term and (1-θ) non-implicit part of velocity-pressure term in the momentum equation to obtain intermediate state of discharges (Q ELM ). For the advection term, a Eulerian-Lagrangian method (ELM) is used to solve the problem. The state of a particle in the DCB centre at time level n+1 is acquired using backtracking from the state of this particle at time level n [8]. And the expression of ELM step can be written as: Based on the discharges of ELM, the bed resistance term and the (1-θ) non-implicit part of velocity-pressure term can be solved by Eq. (5) Second, the velocity-pressure implicit part is solved to get new water levels which can couple Eq.(6) and Eq. (7) to make up a tridiagonal matrix as shown in Eq.(8).
( ) where D i (i=1,2,3) is coefficients of leading diagonal and minor diagonal in the matrix; RHS is the right-hand term. Third, the discharges are obtained from an inverse computation of the water levels in second step by the Eq.(6).
For 1-D water quality model, the forward difference and the centre-difference are used for discretization in time and in space to solve the scalar transport problem, respectively [9], which can be described as Eq.

Solution with prediction-correction method
A prediction-correction method (PCM) is used for solving the implicit part of velocity-pressure term of the free water surface in river channel networks, in which concludes prediction step and correction step [10]. Each branch is made up a tridiagonal matrix, which is a subsystem of the entire networks, i.e., NB branches have NB tridiagonal matrices. To define the end-unit boundaries, which are obtained from input files, as the external boundaries and those around JPs as internal boundaries. The provisional internal boundaries are obtained at the prediction step and are used for the new internal boundaries at the prediction step. Using PCM can avoid to solve a large global matrix for the river channel networks but solve several subsystems locally which can use a large step Δt. For each time level, the PCM of river channel networks can be written as Eq. (10) At the prediction step, the matrix and the right-hand term of each branch subsystem are solved independently. The calculated downstream boundary condition uses the water levels of end units at time level n. After solving the prediction step, the temporary solutions of governing variables are obtained and deposited. At the correction step, the same subsystems are constructed with updated internal boundaries. The final solutions are obtained in this method.

Sub-cycling for water quality model
To improve the stability of water quality model, a subcycling is introduced to divide the computational time step. The transport calculated step Δt wq is given by Eq.(12): where N wqt is the number of substeps of the sub-cycling.

Boundary conditions
In the hydrodynamic model, a group of discharge sequenes are used for the inflow boundary conditions of the start units of branches, and a group of water level sequenes are used for the outflow boundary conditions of the end units (not JP-units) of channel networks.
In the water quality model, for each simulated transport material, the concentration of inflow and outflow units are necessary.

Time stepping
In the preceding description of governing equations, numerical discretization and solution method, the calculation process of one time stepping is summarized as follows: Step 1: Update the boundary conditions. To get inflow/water-level/lateral-inflow boundaries at each time stepping.
Step 2: Main computation. First, to calculate the explicit terms in the momentous equation. Second, to solve the velocity-pressure term to obtain the free water surface of river channels. Third, a backtracking calculation of final discharges. Forth, to solve the water quality problem. Step 3: Update the states of dry and wet units.

Tests and results
A hypothetical case of hydrodynamic model of [6], and two hypothetical waves of Gaussian pulse and a rectangular pulse [11], denoted by "HM-1", "GP-WQ1" and "RP-WQ2", respectively, are tested to demonstrate the coupling model in simulating simple scalar transport problem.

Test of HM-1
Test of HM-1 contains 14 branches with bed slope of 0.00016-0.00047 and 6 JPs, which is shown in Fig.2. Totally 158 cross-sections (CS) or reaches are uniformly distributed in branches with a distance of 100 m. The shape of each CS is triangle or rectangle of 10-30 m with the side slope of 1:1 or vertical. The manning's coefficient is of 0.022 or 0.025. The implicit factor θ is set to 0.6. And the time step Δt is set to 60s. Upstream boundaries are specified at nodes 1-7, and downstream boundary is specified at nodes 14. The JPs at nodes 11 and 12 are selected to compare discharges and water levels [6]. In Fig.3 and Fig.4, the water level process and discharge process at node 11 and 12 are presented, respectively. In Fig.3 and Fig.4, we can find that the solutions of water level and discharge are close to the measured data, which demonstrates the high-precision of the new model in river hydrological process simulation.

Test of GP-WQ1
This test considers the transport of a 1-D Gaussian pulse with an initial configuration centred at start CS owning a standard deviation σ. In the numerical simulations, the data are chosen as [11] and an inflow and outflow Dirichlet boundary condition is used. In the test, the time step is 0.05, and the sub-step is 0.005. The grid size is 1/60. The solutions shown in Fig.5 are closer to the analytical results within a little fluctuation which is caused by the difference method.

Test of RP-WQ2
To observe the performance of proposed method, another example considering the transport of 1-D rectangle wave is used here. The initial condition is given by [11]. The grid size is 0.01. The time step is 0.05, and the sub-step is 0.0001. In the resulting figure (Fig.6), the solutions present sinusoidal form at the high-value part of analytical

Discussions
One test of looped river channel networks and two tests of Guassian pulse and rectengle wave for scalar transport problems have been made. Rather than those methods with solving junction point equation, the proposed model is free of iteration at junction points which can obviously improve the efficiency. Using θ semi-implicit method can simplify the matrix system which considers speed and accuracy in the meanwhile.
Overall, the calculated precision relatively presents in a high level, especially in hydrodynamic part. In water quality part, the precision is relatively lower than hydrodynamic part because of the explicit scheme and sub-step selection. Therefore, the proposed model can be further improved in discretized scheme or using ELM to fastly acquire the concentration of scalar transport materials. Considering that the coupled model realizes the basic solutions of water quality problems, the further studies of different kind of containments can be joint and solved with some chemical reactions.

Conclusions
A 1-D mathematical modelling system for free-surface flows and water quality process is developed to realize a fast and high-precision simulation of unsteady flows and water quality process in river channel networks. In the proposed model, a prediction-correction solver is used for which can adopt a large calculated step instead of solving junction point equation independently, and a sub-cycling of water quality part is used for improving the model stability. To improve the calculated efficient, ELM is used to solve the advection term in the momentums equation. Besides, using θ semi-implicit method to discretize the hydrodynamic equations, and using forward difference in time and a centre-difference in space to discretize the scalar transport equation. Finally, the model is tested using a looped channel network for hydrodynamic part and two kind of waves for whole coupled model which proves the resulting model is free of large-scale iteration and very efficient. Based on the fast and high-precision simulations, the ecological operation of reservoir can improve the elaborated management techniques.