Search This site The web












Back ] Up ] Next ]

Hydraulic Theory

Kinematic Flood Routing

From the MIDUSS Version 2 Reference Manual - Chapter 8
(c) Copyright Alan A. Smith Inc.

Flood routing methods can be classified as hydraulic - in which both continuity and dynamic equations are used - or hydrologic, which generally uses the continuity equation alone.  MIDUSS uses a method based on a kinematic wave equation and therefore falls into the second category.  For the type of conduits used in storm drainage systems, kinematic routing yields results of very acceptable accuracy.

The continuity equation is simply a statement that the difference between inflow and outflow must equal the rate of change of storage in the reach being considered.  Equation [8.20] is a general continuity equation in which lateral inflow is ignored.

If the flow in the channel can be assumed to be quasi-uniform then discharge is uniquely defined by the depth or stage, i.e.


 so that  

Now separating variables in [8.20] yields:

Equation [8.22] is seen to be a wave equation in which a function Q(x) is propagated with a celerity  c  which is given by  c  =  dQ/dA.  The variation of Q with respect to both time  t  and space  x  can be best described using a space-time coordinate system which defines changes in discharge Q in an elementary reach   Dx  and over a timestep   Dt. 

Figure 8.4 - An element of a space-time coordinate system.


Figure 8.4 shows such a situation.  For this element of the space-time system equation [8.22] can be expressed in terms of finite differences as follows.

This equation is applied around a 'nucleus' of the space-time element which is off-centre and defined by the weighting factors   a  and  b   which are applied to the  x  and  t  dimensions respectively as shown in Figure 8.4.

The finite difference quotients of equation [8.23] are expanded around the 'nucleus' point of Figure 8.4 in terms of the values of Q at points 5, 6, 7 & 8.  These in turn can be expressed as weighted averages of the values at points 1, 2, 3 & 4 and the weighting coefficients  a  and  b. Thus:







Substituting equations [8.24] in [8.23] yields:

The quantity (c.Dt/Dx) is a dimensionless time ratio which is equivalent to the Courant criterion for numerical stability.  Denoting this by  t   and substituting for Q5 etc. results in equation [8.26].

Now the process of flood routing usually involves a 'marching' solution in which the initial conditions are known at time  t  and it is required to predict conditions at time (t+Dt).  An upstream boundary condition is provided by the time-history of the inflow hydrograph at  x=0.  In the case of kinematic wave routing it is possible to advance the solution for all values of time t so that the solution advances over the whole time domain for each reach of channel.

In either case, the solution for the element of Figure 8.4 involves an estimate of Q4 in terms of the other three known values.  Collecting terms and casting Q4 as the dependant variable yields:


Equation [8.27] is a generalized form of the Muskingum flood routing method but for the special case of  b =0.5 this reduces to the more familiar form shown below.

Before equation [8.27] can be applied two additional pieces of information are required:

(i)            What values of  a  and  b  should be used to best represent the attenuation for a specific hydraulic condition?

(ii)           What conditions must apply for the computation to be numerically stable?

Evaluation of the Weighting Coefficients

Convergence is the condition in which the solution of a finite-difference equation for a finite grid size approximates the true solution of the partial differential equation which it represents.  It can be shown (Biesenthal (1975) and Smith (1980)) for the non-centred scheme of Figure 8.4 that as the coefficients  a  and  b  depart from a value of 0.5, truncation errors of the order of O(Dx) and O(Dt) increase respectively and independently.  This property is fundamental to the use and apparent success of kinematic routing methods in modelling attenuation.  It should be emphasized, however, that this attenuation results from truncation error and is a property of the numerical finite difference scheme and not of the physical system.  The trick is to find a way to make the numerical truncation error a close approximation to the attenuation which the flood wave will experience.

For a non-centred finite-difference scheme the error may be included in the partial differential equation as e in equation [8.29].

If only first order terms are included in the error term this becomes:


Equation [8.31] is in the form of a diffusion equation where  D  is the coefficient of diffusion.  In order to relate  D  to the physical characteristics of the channel an alternate diffusion equation can be developed using the continuity equation [8.20] with a simplified form of momentum equation in which the convective and temporal accelerative terms are assumed to be negligible, i.e.

where     h             =              water surface elevation

                Sf            =              friction gradient

This can be developed (Smith(1980) see References ) to yield the diffusion equation of [8.33].

where     K             =              channel conveyance defined by Sf = Q2/K2

Now from the definition of the conveyance K the total derivative is obtained as:

Comparison of the terms in equations [8.31] and [8.33] provides a means of evaluating the diffusion coefficient  D  and thus the weighting coefficients  a  and  b  in terms of the hydraulic characteristics of the channel.  Thus:

where the friction head loss over the reach being considered is given by


Equation [8.36] can be simplified by assuming an initial value of  b = 0.5 so that:-

If equation [8.37] yields a value for a  which is less than zero, MIDUSS sets  a = 0.0  and solves for  b   from [8.36].  This value will generally be in the range  0.5    b    1.0.  The generalized Muskingum coefficients of equation [8.27] can then be evaluated and a solution obtained for Q4.

In MIDUSS, evaluation of the diffusion coefficient differs depending on whether the conduit is a pipe or channel.  The Manning equation [8.12] can be differentiated to obtain an expression for dQ/dy as follows:

Substituting in the 2nd part of [8.35] yields an expression for D which can be evaluated in terms of the channel cross-section parameters and the channel gradient, thus:

For pipes, an alternative procedure is used in which a fitted polynomial represents the ratio Q/(dQ/dyr) as a function of the proportional discharge Qr which is the ratio of actual discharge to full pipe capacity.  With very acceptable accuracy this can be represented as follows:

from which the diffusion coefficient D can be found using equation [8.35].

Criteria for Numerical Stability in Flood Routing

Once values have been determined for the weighting coefficients a  and  b it is possible to carry out a check on the numerical stability of the process.  This involves the calculation of limiting values for the grid dimensions   Dand   Dt.

Biesenthal (1975) (see References ) obtained the following condition for numerical stability.

In general, this means that the nucleus of Figure 8.4 must lie above and to the right of a diagonal through the centre of the space-time element which has a slope of -(1/c). 

Figure 8.5 - Stability and numerical error characteristics of a space-time element.

Figure 8.5 shows a typical case in which the time step  dt is approximately half of the Courant time step given by   Dtc  =  Dx/c.  The heavy dashed lines parallel to the main diagonal form a family of lines each of which comprises the locus of points for which the nucleus will generate a numerical error  e  of a specific value.  It is significant that the same numerical attenuation can be produced by a set of ( a, b ) coordinates and MIDUSS makes use of this feature.

The shaded area of Figure 8.5 indicates locations of the nucleus which are numerically unstable.  Because the values of a  and  b are constrained to provide a required numerical error, the criterion for stability given by equation [8.38] must be satisfied by manipulating either the routing timestep  dt or the reach length  Dx.  The greater the coarseness of the space-time grid the greater the chance of numerical instability.  Thus we need to determine upper limits for both  Dx and the routing time-step  dt

Equation [8.35] shows a relationship between the diffusion coefficient D and the weighting coefficients  a and b.  If we assume initially that   b = 0.5 this reduces to:

Similarly if   b = 0.5 the inequality [8.41] becomes:

Taking the 1st and 2nd parts of [8.43] and substituting for  a   from [8.42] we obtain:



To get an upper bound for the routing time-step we use the 2nd and 3rd parts of [8.43] and obtain:



MIDUSS uses the limiting criteria of [8.44] and [8.45] to divide either the reach length L or the hydrograph time-step  ?t into sub-multiples which satisfy the conditions for stability. 

In the case of very long reaches the time step remains unchanged but routing is carried out over two or more subreaches.  The final outflow hydrograph is the only one presented to the user.

In the case of very short reaches only a single reach length is used but the routing time-step is set at   d t  =  Dt/n  (n = 2,3,4...).  Only flow values at intervals of  Dt are presented in the final outflow hydrograph.



(c) Copyright 1984-2023 Alan A. Smith Inc.
Terms and Conditions       Privacy Policy

Design makes
the difference
Latest Update:
v2.25 Rev473
On-line videos

There are videos available to help you use MIDUSS.  Click here.

Quick Story

How about a quick MIDUSS review in 180 seconds or less? We know your time is valuable. Click here for the fast story.

College or University faculty?

Please see our academic page for our low price offer.  Click here.

Theory anyone?

If you are interested in the hydrological or hydraulic theory built into MIDUSS click here.

You have a custom version of MIDUSS distributed through Winz Connections.
Find out more information about the custom Malaysia storms.  Click here.