Improved flood forecasting using geomorphic unit hydrograph based on spatially distributed velocity field

This paper presents an energy model for determining the overland flow velocity in order to improve the low accuracy problem in flow concentration simulation. It furnishes a novel idea for studying flow concentration in ungauged basins. The model can be widely applied in analysis of spatial velocity field, extraction of instantaneous geomorphic unit hydrograph and development of distributed hydrological model. A distributed flood-forecasting model is constructed for Lianyuan Basin in Hunan Province of China. In the proposed method, gravitational potential energy is transformed into kinetic energy via an analysis of energy distribution of water particles in the basin. Based on the kinetic energy equation, the overland flow velocity simulating the geomorphic unit hydrograph is computed. Rainfall-runoff simulation is then performed by integrating with runoff yield and concentration model. Results indicate that the model based on energy conversion leads to more accurate results. The model has the following advantages: firstly, the spatial distribution of the velocity field is appropriate; secondly, the model has only one parameter, which is easily determined; and finally, flow velocity results can be used for the computation of river network flow concentration.


INTRODUCTION
Floods of medium and small catchments are among the most costly natural disasters all around the world, posing serious risk of loss of life, physical injury, damage to infrastructure and disruption of economic and social activities (Shi et al. ). Last decades showed a steady increase of damages due to flooding, which highlighted the need of developing effective measures to reduce flood event impacts (Barbetta et al. ). Conceptual rainfall-runoff (CRR) models have become a basic tool for flood forecasting and become increasingly important for catchment basin management (Xu et al. ). However, accurate flood prediction in mountainous watersheds remains a significant challenge for flood-forecasting models. In order to overcome this challenge, a good knowledge of the relationship between rainfall and runoff is entailed. A unit hydrograph is a common method in flood estimation, which is not only applied in peak flow estimation but also in the creation of complicated flood hydrographs (Khaleghi et al. ) and has been widely used in rainfall-runoff computation since it was proposed by Sherman (). Clark () combined the two concepts of the isochrone method and linear reservoir to establish the instantaneous unit hydrograph. Although the Clark unit hydrograph can be a very valuable technique in flood hydrology, the lack of appropriate techniques to estimate the storage coefficient-R for ungagged watersheds has diminished the application and utility of this technique (Sabol ).
Obviously, there is a close relationship between the hydrological response of a river basin and its geomorpholo- With the continuous development and improvement, the construction of a geomorphic unit hydrograph based on the spatial velocity field has become a practical technology suitable for the flow concentration computation in ungauged basins. However, it is difficult to get a reasonable velocity distribution if only a single factor is considered. When considering multiple factors to estimate the velocity, it suffers from the limitation that the influence weight and parameter value of each factor are determined. In a flat area, there will be many grids with a slope of 0, which need to be treated with a special method. In fact, the transformation from gravitational potential energy to kinetic energy is the root cause of water flow. In the process of transformation, energy will be lost due to terrain, soil, vegetation, river diversion, water conservancy engineering and other factors.
The objective of this paper is to develop a method of constructing GIUH with clear physical meaning and study a method of determining the value of parameters.
Based on this method, the distributed flood-forecasting model of Lianyuan Basin in China is established. In the proposed hydrological model, the spatial velocity field and the instantaneous unit hydrograph can be extracted based on the principle of energy conversion. It has many advantages as mentioned below. Firstly, only one parameter of the model needs to be determined and the process of parameter determination is very simple. Secondly, the error caused by the gradient of 0 in the flat area is avoided and the velocity distribution obtained from the analysis is reasonable. Finally, the velocity results can be used for the simulation of both slope and river network flow concentration, which can ensure the consistency of flow concentration computation to a certain extent and reduce the influence of watershed threshold on model parameters when subwatershed is divided.
The paper is organized as follows: Section 'Methodology' describes the methodology proposed for the area rainfall calculation and runoff generation model and presents the flow concentration model in this study. Section 'Case study' gives details on the selected case study and model evaluation criteria. Section 'Results' presents and discusses the results obtained through the proposed methodology. Finally, we provide the conclusions of this work in the last section.

METHODOLOGY Areal rainfall computation
The areal rainfall is computed by an inverse distance square (IDS) model. It can be approximately considered that the rainfall distribution in the basin is uniform because the area of sub-basin is very small. The result of the interpolation of the central point of the subwatershed is taken as the areal rainfall of the subwatershed.
IDS takes the distance between the interpolation point and a sample point as a weight. A sample point near the interpolation point has a greater weight, which is inversely proportional to the distance. It can be expressed as: where Z is the estimated value of the interpolation point, Z i is the measured sample value, n is the number of measured samples involved in the computation, D i is the distance between the interpolation point and the ith station, p is the power of the distance, which significantly affects the interpolation results, and its selection standard is the minimum average absolute error. When p is taken as 0, the weights of all sample points are equal, and the method degenerates to an arithmetic mean model. When p is infinite, the result is approximately equal to the value of the nearest sample point, and the method is equivalent to the Voronoi model. In practical application, p is usually taken as 2 under the IDS model, which is the method adopted in this paper.

Runoff generation model
The Xinanjiang model is adopted for runoff generation

Overland flow concentration
A unit hydrograph is used to compute the runoff concentration according to the convolution formula.  Exponential parameter with a single parabolic curve, which represents the nonuniformity of the spatial distribution of the soil moisture storage capacity over the catchment -7 I m Percentage of impervious and saturated areas in the catchment -Runoff routing parameter 8 S m Areal mean free water capacity of the surface soil layer, which represents the maximum possible deficit of free water storage mm 9 E x Exponent of the free water capacity curve influencing the development of the saturated area -10 K g Outflow coefficients of the free water storage to groundwater relationships -11 K i Outflow coefficients of the free water storage to interflow relationships -12 C i Recession constants of the lower interflow storage -

C g Recession constants of the groundwater storage -
where m is the number of net rain periods, I is the average net rain period and q is the discharge of geomorphic unit hydrograph period.
According to D8 algorithm (Fairfield & Leymarie ), a water particle in a grid is assumed to flow to the adjacent grid in the direction of the maximum gradient in DEM. According to the size of each grid and the flow velocity in the grid, the detention time of water particles in each grid can be computed by the following formula: where L is the length of the edge of the mesh converted to the actual distance and v is the flow velocity in the grid.
Along the flow concentration route, the flow concentration time from each grid to the outlet of the drainage basin can be computed by the following formula: where S is the slope of two adjacent grid points, A is the grid catchment area, h is the net rainfall of the unit grid, i is the net rainfall intensity and a, b, c, d are some parameters, where parameter a summarizes the combined effect of all other factors not described in the equations.
In fact, the fundamental cause of flood flow is the transformation of gravitational potential energy to kinetic energy.
For any grid in the basin, assuming that the gravitational potential energy of rainfall falling on the basin is completely converted into kinetic energy without energy loss, then the kinetic energy of water particles flowing through the grid point satisfies the following formula: where E ideal is the total kinetic energy of flow in the grid under the ideal conditions; n is the total number of grids in the basin upstream of the target grid (including the target grid); E p is the gravitational potential energy of each upstream grid relative to the target grid; m is the mass of unit net rain in the grid; g is the gravity acceleration and H Avg is the average elevation difference between all upstream grids and the target grid.
Due to the influence of topography, soil, vegetation, river diversion and hydraulic structures, the energy of water flow is gradually lost during the movement. The ratio of gravitational potential energy to kinetic energy is denoted by energy coefficient μ. According to the kinetic energy formula, the following formula is derived: where E k is the kinetic energy of the flow from the target grid. The following equation is obtained by combining Equations (9) and (10): For the grid on the watershed, the kinetic energy of the flow velocity is completely transformed from the gravitational potential energy, and the energy loss is mainly caused by the surface friction. Referring to the force of rigid body motion on the slope surface, the frictional force is directly proportional to the decomposition force of gravity in the vertical direction of the slope when the frictional coefficient is constant. When the slope is steeper, the supporting and frictional forces are smaller, and the speed is faster; when the slope is more gentle, the supporting and frictional forces are larger, and the speed is slower. Let μ ¼ μ 0 sin(θ=2), and rewrite Equation (11) as follows: where μ 0 is the energy residual coefficient, θ is the slope angle of the grid outflow direction and Δh is the elevation difference between the target grid and the outflow grid. For grids that are not on the watershed, the kinetic energy of water flow comes from two parts, one is the kinetic energy of the inflow grids, and the other is the gravitational potential energy difference between the target grid and the outflow grid. The energy equation is as follows: where N is the number of inflow grids of the target grid. The following direct equation of v can be obtained by combining Equations (10) and (13): The following equation is equivalent to Equation (14) v ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi

River network flow concentration
The Muskingum model is used for river network flow concentration, Muskingum equation is derived from the joint solution of water balance equation and Muskingum channel storage curve equation. The discharge at the river outlet can be computed according to the following equation: where I tÀ1 and Q tÀ1 are the inflow and outflow of the channel at the beginning of the simulation, respectively. I t and Q t are the inflow and outflow of the channel at the end of the simulation, respectively. C 0 , C 1 and C 2 are the parameters of the simulation and they can be determined by the following formula: where dt is the time step of the simulation. In order to determine the parameter K, the wave velocity must be computed first. The spatial velocity field is used to compute the average velocity v, and a flow equation (Todini ) is used to compute the wave velocity of flood C.
The following formula is used to compute the propagation time of flood wave in the reach: The computation of x is given by Song et al. () x ¼ 1 2 À n 0:6 Q 0:3 5:14λJ 1:3 L where L is the length of the reach, n is the roughness, Q is the average of the maximum and minimum flow values in the inflow process and J is the gradient of the reach.

CASE STUDY Study area and data
As shown in Figure 3, the study area is located in Hunan

)
: where Q 0 (i) and Q s (i) are the measured and simulated runoff or corresponding runoff, respectively, and N is the number of data points involved.  The NSCE is a well-known performance criteria (Nash & Sutcliffe ), which is defined as: where Q is the mean value of the measured processes. it is qualified if the error of peak discharge or total runoff is less than 20%, and the error of peak time is less than one computation step or 1 h. The qualified rate is computed by the following formula: where QR is the qualified rate; n represents the total number of floods that satisfy the acceptable criteria relative to the peak discharge, peak time and total runoff volume, respectively; and m is the total number of the calibrated floods or The scheme is excellent when QR reaches 85%. The scheme is good when QR is greater than 75% and less than 85%. Otherwise, the results of the performances of parameter calibration are unsatisfactory for online flood forecasting.

RESULTS
29 floods are selected for model parameter calibration from 1979 to 2006. Table 3 gives the results of parameter calibration. According to three statistical ratios of acceptable national criteria relative to the peak discharge, peak time and total runoff volume among the calibrated and validated historical flood events for flood forecasting in China, Table 4 presents the detailed results of statistics performance using calibrated parameters during the calibration for flood simulation. It indicates that the qualified quantity is 24 and the ratio of qualifying simulation is 82.76% relative to the error of peak discharge; the qualified quantity is 27 and the ratio of qualifying simulation is 93.1% relative to the error of peak time; and the qualified quantity is 25 and the ratio of qualifying simulation is 86.21% relative to the error of total runoff volume. Ten floods are used for validation from 2010 to 2017. Table 5 gives the detailed results of statistics performance using the calibrated parameters during the validation for flood forecasting. It shows that the qualified quantity is 9 and the ratio of qualifying simulation is 90.0% relative to the error of peak discharge; the qualified quantity is 10 and the ratio of qualifying simulation is 100% relative to the error of peak time; and the qualified quantity is 9 and the ratio of qualifying simulation is 90.0% relative to the error of total runoff volume.
In order to demonstrate the comparative performance of the proposed Distributed Model based on Spatially Distributed Velocity Field using Equation (14)    Notes: The total number of floods is 29, which are qualificatory relative to the error of peak discharge, is 24 and the ratio of qualifying simulation is 82.76%, which are qualificatory relative to the error of peak time, is 27 and the ratio of qualifying simulation is 93.1%, which are qualificatory relative to the error of total runoff volume, is 25 and the ratio of qualifying simulation is In order to analyze the effect of different velocity formulas, we use the same runoff generation model and parameters, and use Equations (7) and (14) to generate Notes: The total number of floods is 10, which are qualificatory relative to the error of peak discharge, is 9 and the ratio of qualifying simulation is 90.0%, which are qualificatory relative to the error of peak time, is 10 and the ratio of qualifying simulation is 100%, which are qualificatory relative to the error of total runoff volume, is 9 and the ratio of qualifying simulation is 90.0%.  (14) based on the energy conversion method are obviously better than the other models. When Equation (7) is used to compute the flow velocity in the grid, the velocity is often considered to be positively correlated with slope and rainfall intensity, that is, under the same conditions, the greater the slope and rainfall intensity, the faster the velocity will be. In this paper, parameter b of velocity formula is taken as 0.5 according to Maidment et al. (), parameter d is taken as 0.4 according to Li and Wen () and the parameter a is calibrated as 0.5 based on the hydrological data. Although this method can also get good results, but its velocity distribution results are very different from the energy conversion method.   The velocity distribution considering slope and net rainfall intensity (in mm/h) is adopted by Equation (7) in Figure 6(b).
It can be seen from Figure 6 that when the rain intensity is 1 mm/h, the velocity is generally 5-10 m/s on the steep slope, generally 2-5 m/s on the flat land and generally not more than 0.3 m/s in the river channel (many grid values are affected by a slope of 0  (7) can be obtained when parameter c tends to 0. When c is greater than 0, it can improve the situation of the velocity difference between slope land and river channel in Equation (7). The larger the catchment area is, the larger the velocity gain of the grid point is, but it will cause a very complicated calibration process.
It is easy to know from Equation (15) Table 7 presents the average velocity of each flood.
The result of μ 0 using Equation (24) is 0.00498, which is basically consistent with the calibration result μ 0 ¼ 0.005.
The average elevation is 268 m and the outlet elevation is 137 m in the studied area. Assuming that the gravity potential energy is all converted into kinetic energy under is about 51 m/s. The energy conversion μ efficiency can be computed by Equation (25): When μ 0 ¼ 0.005, the energy conversion efficiency is 0.45‰, and it can be seen that a large amount of energy is lost in the process of water flow movement, denoting an inefficient conversion rate.

CONCLUSIONS
Hydrologic models are widely used tools for flood forecasting in the world, as the demand for accurate and reliable flood forecasts has increased in the management of water resources. This paper presents a new method for watershed concentration computation, which can be used to construct the hydrological model in areas without enough data. In order to simplify the computation, the residual energy coefficient is assumed to be constant in this paper and its value can be determined by calibration or measurement of the velocity. From the application point of view, the discretization of residual energy coefficient in spatial distribution is realized by measuring the velocity of different control nodes.
As a matter of fact, it translates the flow concentration problem into a measurement problem (we can measure the flow velocity directly, or compute the flow velocity according to the hydraulic formula by using the measured channel section), and more precise flow concentration parameters can be obtained by more measurement data.
To demonstrate the applicability of the proposed method, a real-world catchment, Lianyuan River Basin is selected. For a 0.5 h time step, 29 historical floods in 28 years