Application of a developed distributed hydrological model based on the mixed runoff generation model and 2D kinematic wave flow routing model for better flood forecasting

Due to the specific characteristics of semi‐arid regions, this paper aims to establish a Grid‐and‐Mixed‐runoff‐generation‐and‐two‐dimensional‐Kinematic‐wave‐based distributed hydrological model (GMKHM‐2D model) coupling a mixed runoff generation model and 2D overland flow routing model based on kinematic wave theory for flood simulation and forecasting with digital drainage networks. Taking into consideration the complex runoff generation in semi‐arid regions, a mixed runoff generation model combining saturation‐excess mechanism and infiltration‐excess mechanism is developed for grid‐based runoff generation and 2D implicit finite difference kinematic wave flow model is introduced for routing to solve depressions water storing for grid‐based overland flow routing in the GMKHM‐2D model. The GMKHM‐2D model, the GMKHM‐1D model with coupling the mixed runoff generation model and one‐dimension kinematic wave flow model, the Shanbei model and the Xin'anjiang model were employed to the upper Kongjiapo basin in Qin River, a tributary of the Yellow River, with an area of 1454 km2 for flood simulation. Results show that two grid‐based distributed hydrological models and the Shanbei model perform better in flood simulation and can be used for flood forecasting in semi‐arid basins. Comparing with GMKHM‐1D, the flood peak simulation accuracy of the GMKHM‐2D model is higher.


Introduction
Hydrological events always vary temporal and spatial distribution for nonlinear and interaction of hydrological process. Distributed hydrological model can describes hydrological mechanism and rainfall-runoff response with spatial precipitation input and the underlying surface information with GIS and RS technology. Coupled to the opening up of the technological feasibility of making such models possible has been the recognition and development of a demand for more distributed predictions as outlined by Freeze and Harlan (1969), Beven and O'Connell (1982), Abbott et al. (1986), Bathurst and O'Connell (1992). There are now a number of distributed hydrological models that are being used regularly for practical application, such as System Hydrological European (SHE) (Beven and O'Connell, 1982;Bathurst and O'Connell, 1992;Refsgaard and Storm, 1995), and the Institute of Hydrology Distributed Model (Beven, 1987;Calver and Wood, 1995). However, these models are not applied for flood operational forecast (Abbott et al., 1986). The Distributed Model Intercomparison Project (DMIP) was conducted to provide insights into the flood simulation and forecast capabilities of 12 distributed hydrological models and a widely used lumped hydrological model (Smith et al., 2004a). Result shows that calibrated distributed hydrological models can perform better than calibrated lumped hydrological models, and distributed hydrological models should be coupled with conceptual rainfall-runoff model and distributed flow routing model for flood forecast purpose (Smith et al., 2004b).
From the last part of the 20th century, a large number of distributed hydrological models or semi-distributed hydrologic models (Beven and Kirkby, 1979;Singh, 1995;Yu, 2000;Todini and Ciarapica, 2001;Yang et al., 2002;Xiong and Guo, 2004;Jia et al., 2005;Xu and Cheng, 2010;Ye et al., 2010;Xue et al., 2013) have been developed for flood simulation and flood forecast. Grid-based distributed hydrological models, which consider the spatial variability, have become one of the most important tools for the present hydrological investigations (Abbott and Refsgaard, 1996;Vieux, 2001). Large number of distributed hydrological models such as Distributed Hydrology-Soil-Vegetation Model (DHSVM) model (Wigmosta et al., 1994), TOPographical Kinematic Approximation and Integration (TOPKAPI) model (Ciarapica and Todini, 2002;Liu et al., 2005), The developed GMKHM-2D model for better flood forecasting 285 TOPgraphy based hydrological MODEL (TOP-MODEL, Beven and Kirkby, 1979) have been developed and applied in practical use. The runoff generation of the models is based on saturation-excess mechanism (Bao, 2006;Li et al., 2006;Wang et al., 2007;Yao et al., 2009Yao et al., , 2012Bao et al., 2010Bao et al., , 2011. These models perform well in humid regions and semi-humid regions, but do poorly in semi-arid basins and arid basins as the runoff generation mechanism is different in such regions (Bao et al., 2016(Bao et al., , 2017. The hydrological regime in semi-arid basins is extreme and highly variable, where flash floods from a single large storm can exceed the total runoff from a sequence of years (Wheater et al., 2008). In recent years, some distributed hydrological models based on infiltration-excess runoff mechanism such as the Grid and Holtan method based distributed hydrological (Grid-Holtan) model (Bao et al., 2016) have been developed for flood forecasting in semi-arid basins. Actually, there are saturation-excess runoff and infiltration-excess runoff in the semi-arid basin. The process of runoff generation on a hillslope with a complex topography should be calculated such that the process of flow concentration is adequately represented (Liu et al., 2004). Depressions storing water is usually ignored because of using D8 method in sinks filled module of distributed hydrological process modeling (Govindaraju et al., 1992;Singh, 1996Singh, , 2001Singh and Frevert, 2002;Mohammadian et al., 2004;Singh, 2006, 2007;Wang, 2010;Wang et al., 2010b;Ye et al, 2013).
In this paper, a conceptual mixed runoff generation model combining saturation-excess runoff and infiltration-excess runoff is introduced for runoff generation calculation, 2D implicit finite difference kinematic wave model is applied to solve catchment sinks storing water for grid-based concentration of overland flow routing and 1D implicit finite difference kinematic wave model is used for channel routing. A Grid-and-Mixed-runoff-generation-and-twodimensional-Kinematic-wave-based distributed hydrological model (GMKHM-2D model) coupling the mixed runoff generation model and 2D overland flow routing model based on kinematic wave theory was developed for flood simulation and forecast with digital drainage networks. In the GMKHM-1D model, 1D kinematic wave model is applied for flow routing and the modules of runoff generation and channel routing is the same to the GMKHM-2D model'. The GMKHM-2D model, the GMKHM-1D model, the Shanbei model and the Xin'anjiang (Zhao, 1983;1992;Zhao et al., 1995) model were applied in the upper reaches of the Qin River above Kongjiapo station, a tributary of the Yellow River, to investigate the reliability and improvement of the developed model.

Study basin
In order to investigate the applicability of the GMKHM-2D model, it was applied to the upper reaches of the Qin River above Kongjiapo station, a tributary of the Yellow River. The GMKHM-1D model, the Shanbei model and the Xin'anjiang model were also applied for comparison. The control area of Kongjiapo station is 1454 km 2 and is located between latitudes 36.18 ∘ N and 37.00 ∘ N and longitudes 111.95 ∘ E and 112.56 ∘ E (Figure 1). The basin average annual precipitation is 544 mm, 67.8% of which is within the period of the flood season. Three rainfall stations and one hydrological station are available in the Kongjiapo basin.

The GMKHM model
The GMKHM-2D model and the GMKHM-1D model are based on raster or grid data structures and use identical square DEM cells as primary computational elements for rainfall-runoff modeling, and adopt grids for the DEM as computational elements for distributed rainfall-runoff modeling with each element consisting of a runoff generation component and a grid-based flow routing component.

The mixed runoff generation model
Single point runoff generation mechanism is classified into, saturation-excess in humid regions and infiltration-excess in arid regions (Zhao, 1983). The mechanism of runoff is the mixture of both in semi-humid and semi-arid regions. The Xin'anjiang model, based on saturation-excess runoff generation, has established flood forecasting efficiency by long-term use in humid regions. However, in semi-arid region, the efficiency of the model is not always ideal. Therefore, a mixed runoff generation model with coupling the Green-Ampt infiltration-excess model (Green and Ampt, 1911) to the Xin'anjiang model runoff generation formation was developed in this paper.
According to the relation of soil moisture (W), water field capacity (W T ) and precipitation intensity (i) with infiltration rate (f ), we have the following governing equations. If If In which R S is a surface runoff, R G is groundwater runoff. The four equations are the complete runoff formula at a point of soil layer. The surface and groundwater runoff in saturated areas is calculated through the  Equations (2) and (3). The infiltration-excess runoff is calculated by Equation (1). The runoff formation on repletion of storage is opposite to the runoff formation in infiltration-excess, but the combination of both is the complete.
The Shanbei model (Zhao, 1983) is the simplest method simulating the infiltration-excess runoff. In the loess plateau, the runoff can be calculated by the infiltration curve: If During the process of a flooding event, the saturation area is varying. all precipitation is runoff ,which called the excess storage runoff in the saturation area. The runoff in unsaturated area is infiltration-excess overland flow. As flowing downhill through the saturation area, it can be as return flow, subsurface flow or groundwater flow. Therefore, the part of infiltration-excess overland runoff (the ratio is ) could be mixed with the runoff from saturation area. The other part (the ratio is 1− ) as overland runoff directly flows to channel network. The area ratio generation infiltration-excess overland flow is considered as parameter (i mf ). Therefore, the mixed runoff generation model can be summarized as follows; If PE = P − E > 0, then there is runoff, otherwise, not. Let R is the runoff amount calculated through the Xin'anjiang model. Then the runoff area is FR = R/PE and unsaturated is 1 − FR. A parabolic equation, which represents the area distribution of infiltration rate in the unsaturated area, defines as Where is the area in which infiltration rate is less than f and E f is a parameter. f mm = (1 + E f )f m . Here f m is the point maximum of f and f mm is the average maximum of f . The actual infiltration f can be calculated by the Green-Ampt infiltration model. i rs is runoff amount generation from the unsaturated area: If PE ≥ f mm then If PE < f mm then The amount of runoff (i rs ) flows through the saturated area and the amount i rs (1 − ) directly flow to channel networks.
If PE + AU + i rs FR < SMMthen The water storage in the free water reservoir is The other runoff generation and separation of runoff calculations are same with that of the Xin'anjiang model. In which, SM and SMM are respective the basin mean of free water storage capacity and free water storage capacity.

Routing model based on 2D kinematic wave
The influence of microtopography on the hydrological response of natural basins is significant (Singh, 2002). Some depressions in the overland storing water for some time during rainfall after ponding occurs. This water is retained on the surface and is ultimately evaporated or infiltrated (Cappelaere et al., 2003).
To solve this problem, this paper introduced 2D implicit finite difference kinematic wave model for overland flow routing. The errors contained in the DEM data are corrected with ANUDEM method (Hutchinson, 2004;Zhang et al., 2005). Once the retention storage in a grid cell is filled by excess rainfall, overland flow occurs. Figure 2(a) is the partly original DEM data of the upper reaches of the Qin River above Kongjiapo station. From the figure, the elevation of the central grid is the lowest in that of all grids. The elevation of the central one ought to be added with D8 method in sinks filled module of most distributed hydrological modeling, which results in grid flow can be out of the central grid with single direction (shown in Figures 2(b) and (c)). However, for the contribution for the hydrograph of the basin outlet, depressions water storing capability is neglected in this hydrological modeling process. In this paper, two-directional model (Singh, 2002) was introduced for taking into consideration sinks water storing. The sinks in DEMs were calculated as a reservoir, so the flow can routing out of the central grid while the water stage being above that of the other grid (shown in Figure 2(d)).
The Saint-Venant equations with the continuity equation and momentum equation were used to describe the governing equation of overland flow. Overland flow model based on kinematic wave theory can perform enough precision in flood simulation and forecasting (Singh, 2002;Rui et al., 2008;Rui and Jiang, 2010). In the GMKHM-2D model, 2D implicit finite difference kinematic wave model was applied overland flow routing. Expressed as the two partial differential equations in the form of Equations (13)-(15) ( Taky et al., 2009).
The governing equation for the overland flow is a 2D continuity equation: Where, h is the depth of overland flow, q x is the unit discharge in the x-direction, q y is the unit discharge in the y-direction, r e = (P e − f ),is excess rainfall rate , P e is the net rainfall rate and f is the infiltration rate.
Momentum equation: x-direction: y-direction: Where S o(x,y) is the land surface slopes in the xand-y-directions, S f (x,y) is the friction slopes in the x-and y-directions, u ,v are respectively the average velocity in the x-andy-directions, in units m/s, g is the acceleration of gravity, in units m/s 2 .
In accordance with the kinematic wave theory, ignoring the inertia term and the pressure term in the Saint-Venant equation, Equations (14) and (15) can be simplified to obtain kinematic wave equation, as following: Where S ox and S oy are the land surface slopes in the x-and y-direction, respectively, and are calculated from DEM data. S fx and S fy are the friction slopes in the xand y-direction, respectively.
The unit discharge in the x-direction q x is defined using a general resistance form: Similarly, in the y-direction q y : For the Manning equation is a constant (Maidment, 1993), and the terms x, y are defined as: 288 H. Bao et al. The GMKHM-2D model used a grid where each grid cell is assumed a homogenous unit with one representative value for any hydraulic or hydrological parameter. Figure 2(e) is the sketch of DEM grid in two-dimension kinematic wave model.
For a grid unit, its principle of mass conservation can be described that the variation of the grid cell water volume is the difference between inflow and outflow in the dt time.
Calculation of h in Equations (18) and (19) is as follows.
For x or y-direction, the continuity equation can be simplified based on kinematic wave theory.
where B is water width, L is channel length, P e (x,t) is the net rainfall rate.
Through to discretize Equation (12) with the four-point implicit finite difference scheme, the following equation can be shown.
( For j = 0, the following equations can be shown. Where N = L/Δx, is iteration number. According to Equations (24) and (25), D 1 and D 2 can be calculated. Substituting the upstream boundary condition into Equation (26), we get the following equation. ( Assuming the calculated runoff depth is z, so the following equation can be shown. h 1 1 can be calculated by bifurcated approach. In a similar way, according to Equation (26), runoff depth of every grid unit can be obtained and every grid discharge can be calculated by Manning formula in the studied basin.
Flow direction is based on friction slopes in the kinematic wave model. S 0x and S 0y can be calculated by DEM data.
The unit discharge in the x-direction between grid cell (j, k−1) and (j, k).

Routing model based on 1-D kinematic wave
One-dimensional implicit finite difference kinematic wave model is used for channel routing in the GMKHM-2D model and the GMKHM-1D model. Compared with routing model based on 2D kinematic wave, 1-D kinematic wave routing model is simplified as the single direction. Therefore, the governing equation and momentum equation for 1-D overland flow can be written as Calculation steps of numerical solution for 1D kinematic wave overland flow routing model are similar to that of 2D kinematic wave (from Equation (23) to Equation (31)), and same to that of 1D kinematic wave channel routing model with replacing h, q (Equation (32)) for channel.
According to the global DEM data and the land use data resolution of 30 ′′ × 30 ′′ provided by USGS (2015), the Kongjiapo basin was extracted. The types of land use are reclassified into four main types: forest, grassland, cropland, and water ,which respond to 0.1,0.17,0.035 and 0.0 of roughness coefficient value respectively (Yao et al., 2012).The interflow and groundwater runoff concentration calculations are same with that of the Xin'anjiang model.

Application and comparison with the GMKHM-1D model, the Shanbei model and the Xin'anjiang model
In order to verify the applicability of the GMKHM-2D model, it was applied to the upper reaches of the Qin River above Kongjiapo station for flood simulation.
The GMKHM-1D model, the Shanbei model and the Xin'anjiang model were used for comparison. Ten representative flood events from 1958 to 2001 are selected. Flood events from 1958 to 1988 were applied for models calibration and the other three events were used to verify the models in the studied basin. In this paper, the time step is 1 h due to the data constraints in Kongjiapo basin (Cundy and Tento, 1985;Wang et al., 2010a). Table 3 is simulated results of 10 representative flood events in Kongjiapo basin.
The statistic values of models results consist of the relative error of flood volume, the relative error of flood peak discharge, the peak time error and the Nash-Sutcliffe efficiency value (Schaefli and Gupta, 2007). The Shanbei model and the Xin'anjiang model parameters are calibrated with artificial optimization method based parameters physical characteristics and historical flood data (Zhao, 1983). For better comparison, the parameters of saturation-excess runoff generation are unified in the GMKHM-2D model, the GMKHM-1D model and the Xin'anjiang model, and the parameters of infiltration-excess runoff generation in the GMKHM-2D model agree with that in the GMKHM-1D model and the Shanbei model. The routing parameter is Muskingum coefficient, which is valued 0.2 in the Xin'anjiang model and the Shanbei model. The runoff generation parameters of GMKHM-2D model, the GMKHM-1D model, the Shanbei model and the Xin'anjiang model can be seen in the Tables 1 and 2. For the GMKHM-2D model and the GMKHM-1D model, most peak flows of simulated flood events were overestimated, probably because of the theory of runoff generation mechanism, since the two grid-based distributed hydrological models assumes the generation of overland flow as precipitation is over the infiltration rate. For the Xin'anjiang model, all of flood peak flows the simulated events were underestimated. The main reason is that the runoff generation calculation model is based on saturation-excess mechanism in the Xin'anjiang model, but infiltration-excess runoff generation is the main runoff in Kongjiapo basin. It caused that the Xin'anjiang model performed poor for peak occur time and flood volume simulation and Nash-Sutcliffe efficiency value. Therefore, the mean value of peak discharge simulated by the GMKHM-2D model and the GMKHM-1D model is better obviously than that of the Xin'anjiang model. From Table 3, for the relative error of flood volume, the peak time error and the Nash-Sutcliffe efficiency value, the two grid-based distributed hydrological models and the Shanbei model performed better than the Xin'anjiang model. This proves the advantage of mixed runoff generation mechanism in simulating flood hydrograph, especially peak flows simulation in semi-arid basin.

Comparison of the GMKHM-2D model and the GMKHM-1D model
For the advantages of GIS technology, DEM and excess infiltration mechanism, the GMKHM-2D model and the GMKHM-1D model performed well in simulating observed flood thin hydrograph in Kongjiapo basin. The two distributed hydrological models obtained the closed accuracy in the relative error of flood volume and the Nash-Sutcliffe efficiency value of all flood events, but the GMKHM-2D model performed significantly better in flood peak discharge and peak time simulation (Figure 3(a)). The reason is depressions storing water in the basin overland is neglected because of using single direction with D8 method in sinks filled module of hydrological process modeling, which results in flood peak discharge overestimating and peak time advances from zero to three hours. It is serious for flood simulation and forecasting in semi-arid basin. In the GMKHM-2D model, 2D implicit finite difference kinematic wave model was introduced to solve basin depressions storing water for grid-based concentration of overland flow routing, so the relative error of flood peak discharge and the peak time error is 7.3%, 0.5 h, respectively. Distributed hydrological models can represent the spatial map of the hydrological processes, which is one of the most representative feature in distributed hydrological modeling. From Figure 3(b), it can be clear that the peak (19th) time step simulated spatial map of the GMKHM-2D model runoff depth for FloodNo.19580802 is very reasonable and close to observation. Wherever, the relative error of flood peak discharge and the peak time error simulated by the GMKHM-1D model is 12.7%, 1 h, respectively. It proves the reliability and improvement of the developed GMKHM-2D model.

Conclusions
The GMKHM-2D model coupling a mixed conceptual runoff generation model and 2D overland flow routing model based on kinematic wave theory was developed for flood simulation and forecasting. The mixed runoff generation model with coupling the Green-Ampt infiltration-excess model to the Xin'anjiang model runoff generation formation on repletion of storage was applied for runoff generation calculation in semi-arid regions. Two-dimensional implicit finite difference kinematic wave model was introduced to solve depressions storing water of basins overland for grid-based concentration of overland flow routing. In order to investigate the reliability and improvement of the GMKHM-2D model, the GMKHM-1D model, the Shanbei model and the Xin'anjiang model were also employed for comparison in the Kongjiapo basin.
Depressions storing water of overland flow routing was taken into consider, so the GMKHM-2D model perform higher accuracy, especially in flood peak discharge and peak time simulation in semi-arid basin. Compared with the Xin'anjiang model, advantages of the two grid-based distributed hydrological models show obviously flood simulation and forecasting. However, compared with the Xin'anjiang model, five parameters are added to the GMKHM-2D model. The further research should focus on the parameter physical meanings and on how to estimate the parameters according to the physical aspect such as soil and land cover information.