A bivariable coupling model for river channel routing developed from the flow continuity equation and its application

In this study, a bivariable coupling model for river channel routing is presented. The proposed model is developed from the Priessmann 4-point implicit differential scheme with a weight coefficient of river f low continuity equation. It is based on the transformation of two different expression forms of river channel storage equation. Furthermore, we consider the impact of lateral inflow along the study river channel from another perspective. In this paper we deduct lateral inflow from the lower section instead of adding lateral inflow to the upper section. In order to be representative of geographical range, river channel characteristics, f lood magnitude, hydraulic characteristics and time, the proposed model is tested in 38 river channels of 6 river systems in China by using observed data during flood season. The rationality of model structure and the validity of model simulation are examined comprehensively. Comparison between the proposed model and Muskingum model shows that the proposed model can improve the simulation accuracy. The results show that the simulation accuracy and stability of the bivariable coupling model is much better than that of the Muskingum model.


INTRODUCTION
Flood forecasting based on hydrological models is an important non-engineering measure for flood control and disaster reduction, which has received increasing attention from public, government and academic communities (Al-Safi and Sarukkalige, 2017;Yu et al., 2014).Efficient reservoir operation, river management, flood control and warning depend on reliable and accurate real-time flood forecasts (Cloke and Pappenberger, 2009;Hartnett and Nash, 2017;Si et al., 2015), which can be achieved by using hydrological models (Hostache et al., 2018).However, both conceptual models and distributed models need an accurate river flow routing module to calculate flood processes (Barati et al., 2012;Xu et al., 2017).In addition, the river flow routing model can calculate water level of flood control section, because water level is an important index and reference basis for flood control and flood risk management (JING et al., 2014;Wu et al., 2014).
There are many methods for flood routing, which are mainly divided into two categories: one is the hydraulic method based on Saint-Venant equations (Molls and Molls, 1998); another is the hydrologic method based on the water balance equation and river channel storage equation (Kubo, 1985).The physical meaning of the hydraulic method is clear, but it needs detailed data of river section, roughness and water surface slope.The hydrologic method focuses on the relationship among hydrological elements, which can simulate the main characteristics of a flood in the river channel; the hydrologic method is simple, practical and operable (Bao and Zhao, 2011).Therefore, the hydrologic method is more widely used for flood routing than the hydraulic methods (Fread, 1981).
The hydrologic f lood routing method mainly includes the Muskingum method, linear regression method, runoff coefficient method, characteristic river length method and delay algorithm method.While the application of the Muskingum method is the most widely used due to the low requirement for river channel topography and roughness data (Singh and McCann, 1980;Tang et al., 1999), it has good performance in general river f lood calculation.In addition, the Muskingum model with variable parameters has been applied in ungauged basins, utilizing the Muskingum model with variable parameters of wave travel time KE and weight coefficient of discharge XE based on the physical characteristics of the river reach and f lood, including the reach slope, length, width, and f lood discharge (Song et al., 2011).
The Muskingum routing method has been improved by many researchers since it was proposed in 1939 (Koussis, 2009;McCarthy, 1939).Most of these theoretical studies on the Muskingum method include model validation (Afzali and Niazkar, 2015), analysis of model structure, model structure improvement (Haddad et al., 2015), and calibration of model parameters by considering the temporally varying factors.
Muskingum model structure is usually validated by the comparisons between the observed and simulated discharge in the downstream section.The simulated discharge in downstream section can be obtained using the observed discharge in upstream section via Muskingum routing method (Afzali and Niazkar, 2015).Then, the magnitude and characteristics of the simulation errors are examined to evaluate the rationality of the model structure and the accuracy of simulation (Baláž et al., 2010;Chatila, 2003;Nash, 1959;Singh and McCann, 1980).

200
Physical basis analysis of Muskingum model structure is conducted by analysing the relationship between the storage and release discharge of river channel, or by analysing the relationship between Muskingum river routing model and the differential model of Saint-Venant equations (Baymani-Nezhad and Han, 2013).Therefore, the physical property of the Muskingum model was indirectly proved by analysing the physical basis of model structure (Cunge, 1969;Diskin, 1967;Wang et al., 2006).
In order to consider the non-linear relationship between river channel storage and outflow, structural improvement of the Muskingum routing method is mainly about the nonlinearization of a linear storage-release relationship and the study of model structure simplification (Kumar et al., 2011;Ponce, 1979).
Research on model parameter calibration and time varying factors is conducted by establishing the relationships between model parameters and time varying hydraulic factors such as wave travel time KE and weight coefficient of discharge XE (Perumal and Raju, 2001;Ponce and Yevjevich, 1978;Saxena and Perumal, 2014).
Many achievements have resulted from the intensive studies on these original and improved Muskingum routing methods.The result reveal that most of these methods can only be applied in the river channel which has a stable stage-discharge curve.However, the traditional Muskingum model cannot provide accurate simulations when applied in river channels with unstable stage-discharge curve or in tidal channels without stage-discharge curve (Bao et al., 2009;Bao et al., 2010;Bao et al., 2007;Si-min et al., 2009).
There are great differences between hydrological methods and hydraulic methods in flood routing, especially in water level forecasting.For example, hydraulic methods need water level or discharge of upper and lower sections as boundary conditions.The hydraulic methods have no forecasting function, and even if they have forecasting function, the lead time will be lost.However, hydrological methods only require the cross-section survey data of flood control section, and can forecast the future water level process well (Bao et al., 2018;Zhang and Bao, 2012).
In order to solve the abovementioned problem relating to the Muskingum routing model, we propose a bivariable coupling model of universal applicability for river channel routing in this work.The method is based on the river channel storage equation and the differential scheme of river flow continuity equation with consideration of lateral inflow.As far as the structure and basic theory of the model are concerned, the model proposed in this paper belongs to the hydrological method.Therefore, the advantage of this proposed model is that it has the function of forecasting, and does not need highprecision river channel survey data as a prerequisite.

MODeL Differential flow continuity equation
If there is no lateral inflow, the one-dimensional flow continuity equation can be expressed by cross-sectional area A and discharge Q as: Using the Priessmann 4-point implicit differential scheme, the differential form of Eq. 1 can be expressed as follows: (2) where θ is the weight coefficient of Priessmann differential scheme; the subscript j and j + 1 represent the upper and lower section, respectively; ∆Q = Q i+1 -Q i ; ∆A = A i+1 − A i ; the superscript i and I + 1 represent the last and next time period, respectively.The ordinate t is time; the horizontal ordinate x is channel length.The meaning of superscripts and subscripts in Eq. 2 are shown in the schematic diagram in Fig. 1.
Lateral inflow comprises a large proportion of river channel storage during flood season.However, lateral inflow is not likely to be represented by a particular function; typically there is lack of data.In order to take into consideration lateral inflow in Eq. 2, a general treatment is to multiply the upper-stream section discharge and crosssectional area by an amplifier coefficient (Bao et al., 2010;Kumar et al., 2011).This means that the lateral inflow is considered by adding it into upper-stream section inflow.Actually, the lateral inflow is distributed along the river reach, and its impact on river storage is distributed along the river reach.Here, we propose to deduct the lateral inflow from the discharge of the lower cross-section, which is to multiply the lower-stream section discharge by a reducer coefficient.Then we can get a new differential equation of Eq. 1 as follows: (3) where β is the lateral inflow reducer coefficient, the other variables in Eq. 3 are as for Eq. 2.
The discharge and cross-sectional area of the lowerstream section include all of the information which can reflect the lateral inflow affected by river storage.Therefore, it is more reasonable to deduct the lateral inflow from the lower-stream section.

Relationship analysis
In order to get the solutions of Eq. 3, another equation is needed to describe the relationship between cross-sectional area and discharge.In hydrodynamics, the relationship between cross-sectional area and discharge is based on force balance.The simulation accuracy of the hydraulic method is not robust due to the large difference between the simulated and actual frictional resistance (Xiaoqin and Weimin, 2012;XiaoQin et al., 2009;Zhang and Bao, 2013).In hydrology, the storage discharge relation was proposed to develop the Muskingum routing method.However, this method is only appropriate for river channels with stable stage-discharge curve (BAO et al., 2009;Bao et al., 2010;Bao et al., 2007).In view of this, we tried to obtain the direct relationship between cross-sectional area and discharge by using the following two different expression forms of river channel storage equation; following this a bivariable coupling model for river channel routing can be established.
where W is the storage of river channel, L is the length of river channel, A is the average cross-sectional area of river channel, K is the flow propagation time, Q is the average discharge of river channel.Based on Eq. 4 and Eq. 5, an equation between average cross-sectional area A and average discharge Q can be obtained as below: Equation 3 can be solved once the relations between A and Q and the corresponding elements of upper and lower cross section are known.For this purpose, the following two assumed equations are introduced to solve this problem.
where α is the weight coefficient of average cross section, χ is weight coefficient of average discharge.The storage-discharge relationship of the river routing model can be obtained by taking Eq. 8 into Eq. 5.The rationality of the Muskingum storage-discharge relationship has been widely demonstrated (Baymani-Nezhad and Han, 2013;Kumar, Baliarsingh and Raju, 2011;Perumal and Raju, 2001).Then the relationship between cross section and discharge can be described by Eqs 6, 7 and 8.
The rationality of the above structure has been studied by Koussis (2009) and Kumar et al. (2011).

Model structure
The temporal differential form of Eq. 9 can be expressed as follows: (10) where ∆Q and ∆A are as defined for Eq. 1.
In both Eq. 3 and Eq. 10, all of the variables on the right-hand side of the equations are known, while the left hand values are unknown.The bivariable coupling model is established by using these two equations in this study.The coefficient determinant of these two equations can be expressed as: The above equation is tenable for any weighting coefficient value of α and χ within the range of (0, 1).The value of the determinant is > 0, which means that Eq. 9 is effective.The discharge of lower-stream section can be obtained by the simultaneous solution of Eqs 3 and 9.

Data description
In order to fully evaluate the rationality of model structure and the validity of model application, we selected river channels based on the following two principles in this study: • River channels should represent a wide range of hydraulic and cross-section characteristics.• In view of the simplified method for dealing with lateral inflow in this study, the lateral inflow area and the proportion of lateral inflow should not be too large.If the proportion of lateral inflow is too large, this simplified processing method will likely result in a large simulated error.According to the above principles, we selected river channels with different hydraulic and section characteristics for this study.These river channels are located in water systems all over the country (Yangtze River, Yellow River, Huaihe River, Liao River, Songhua River and an inland river channel of Tarim River).The river channel properties and data series used in this study are presented in Table 1, with river channels arranged in ascending order according to the river length.The locations of these study river channels are shown in Fig. 2. In order to avoid the influence of changes

202
in cross-section characteristics on data consistency, the selected data series should be as continuous as possible.
The selected data used for calibration and validation were all in flood season, with the duration of the flood event as long as possible (the shortest flood event lasted 10 days, while the longest is 20 days).The time interval of flood event data is 0.5 h.Since the time interval of the observed data we used was very short (0.5 h), it was difficult to collect data for many flood events (only 8 or 10 flood events in a river channel).

ResULTs AND ANALYsIs
The Muskingum river flow routing model can be expressed as follows: In order to test the performance of the bivariable coupling model, we used the same data series for both the Muskingum

203
model and the proposed model during the calibration period and validation period.Here, the Nash-Sutcliffe (NS) value was adopted to evaluate the model performance.
where Q -O is mean value of observed series, Q C is the computed series, Q O is the observed series, M is the number of data series.NS is used to evaluate the performance of the fit between calculated and observed flow.A larger value of NS indicates greater fitness of the model, the maximum value for which is 1.0.The collected data were divided into two parts, which were used for model parameter calibration and parameter validation, respectively.
The linearized parameter calibration method was adopted in this work.The linearized calibration method is a new optimization algorithm, which was developed to solve the theoretical problem of unrelated local optima produced in the nonlinear model parameter calibration by using the objective function based on error sum of squares.The calibrated parameter values are stable when using different initial parameter values via the linearized calibration method which can find the true parameter values without producing unrelated local optima.Furthermore, compared with the SCE-UA method and the simplex method, the linearized calibration method has higher calculation accuracy and convergence rate, and the parameter calibration results are also more stable, due to not being influenced by the different initial parameter values (Bao and Zhao, 2014).The linearized calibration method is an efficient, effective, and robust calibration method (Bao et al., 2013;Si et al., 2017).
In order to illustrate that there is no systematic deviation during the study by the calibrated parameters of the two models, we presented NS values of the traditional Muskingum model and the bivariable coupling model in the same scatter plot (Fig. 3).
In Fig. 3, NS-M and NS-B represent the average NS value of Muskingum model and bivariable coupling model, respectively.
From Fig. 3, it can be seen that the results of the calibration period and validation period are evenly distributed on both sides of the 45° line, which shows that the calibrated parameters of the two models are reasonable.
According to the analysis of the simulation results shown in Table 2 and Fig. 4, the results of the bivariable coupling model are better than those of the traditional Muskingum model, during both the calibration period and the validation period.In addition, Fig. 4 also shows that the stability of bivariable coupling model is better than that of the traditional Muskingum model.However, not all of the simulation results reach a higher NS value, the river channel length is treated as step length during practical application; a longer step length will cause a larger differential error in the continuous differential equation.In order to analyse this problem, we conducted a further statistical analysis and analysed the relationship between the mean NS value and the river channel length; the result is shown in Fig. 5.
Figure 5 shows that there is a significant correlation and trend between the simulation results and the river channel length.This shows that the interval inflow has a large influence on the simulation precision during flood routing.In general, the longer the river channel length, the larger the proportion of interval inflow, which is one of the important factors leading to the large difference error of the continuity equation.Although both the Muskingum model and the bivariable coupling model are influenced by the same difference error, simulation results of the proposed model are markedly better than the Muskingum method.
In order to show the application effect of the proposed model intuitively, the specific calculation results of a typical river channel (Wuzhou-Gaoyao) are presented in Table 3, and

204
the observed and simulated results by different models for the 10 flood events are shown in Fig. 6.The results in Table 3 and Fig. 6 show that the proposed model performs better than the Muskingum model.The overall performance indicates that the bivariable coupling model is more stable for river channels with different hydraulic characteristics, and the bivariable coupling model has a more reasonable structure, wider representativity and is more adaptive than the Muskingum model.Therefore, the structure rationality and simulation validity can be detected to a certain extent by using the same data series.

CONCLUsIONs
The bivariable coupling model for river channel routing is proposed in this paper.The proposed model is developed from the flow continuity equation based on the Pressimann 4-point implicit differential scheme with weight coefficient.In this work, the main innovation of the proposed model is the consideration of lateral inflow along the study river channel by deducting the lateral inflow from lower-stream section flow.

205
The proposed model was fully verified in 38 river channels of 6 river systems in China by using observed data during f lood season.The rationality of model structure and the validity of model simulation were estimated comprehensively.Comparisons between the proposed model and Muskingum model were also presented in this study.The results show that the simulation accuracy and stability of the bivariable coupling model is much better than that of the Muskingum model.The applicability of the bivariable coupling model is more general.In addition, the proposed model has the advantage of simple structure and stable performance.
Through this study, we also draw some conclusions on the applicability of the proposed model.We suggest that the river channel length should not be too long during the process of using this model.The river channel length is treated as step length during practical application; a longer step length will cause a larger differential error in the continuous differential equation.But even so, the proposed model is superior to the traditional Muskingum model.This conclusion also provides a theoretical basis and reference for the future application of the model.Further research will continue to improve the proposed model and demonstrate its applicability.

Figure 1
Figure 1The schematic diagram of superscripts and subscripts in Eq. 2

Figure 2
Figure 2Locations of study river channels

Figure 3
Figure 3 The relationship between mean NS value of calibration results and mean NS value of validation results of the two models

Figure 4
Figure 4 The average NS Value of the 38 river channels: (a) calibration result, (b) validation result

Figure 5
Figure 5 Changing trend of mean NS value with river channel length: (a) changing trend of Muskingum model in calibration period, (b) changing trend of Muskingum model in validation period, (c) changing trend of bivariable coupling model in calibration period, (d) changing trend of bivariable coupling model in validation period