The Muskingum flow routing method has been very well researched and established in the hydrological literature . Its modest data requirements make it attractive for practical use . The gives a general overview of the Flood routing concept and types , and then goes on to explain the Muskingum method in detail

Introduction to Flood Routing

to determine the flow hydrograph characteristics like shape and movement along a water course , and how these are affected by various factors like system storage and system dynamics on the shape and movement of flow hydrographs along a watercourse . In other words Flood routing can be described as a process of calculating outflow rates reservoir stage and storage volume from a stream channel once inflows and channel characteristics are known . The process of flood routing is used for the hydrologic analysis in flood forecasting , flood protection reservoir design and spillway design etc . The principle of routing is used here for predicting the temporal and spatial distribution of hydrograph , during the course of its travel through the various sections of a stream (Subramanya 2002

Basic Principles of Routing

All hydrologic routing methods use a common continuity equation as their common base . According to this equation , the difference between inflow and outflow rates is equal to the rate of change of storage Mathematically the equation can be written as below

(Gosh 1997

. 67

In the above equation , I is the rate of inflow , and at any time the corresponding outflow is O . dS is the storage that is accumulated during a very small duration of time dt . Figure below represents the pictorial relation between storage S and discharge Q

(Flood Routing

The above equation considers the losses due to seepage , evaporation and direct accretion to storage , as small enough to be ignored

Hydrological routing - These methods mainly use the continuity equation

Hydraulic routing - These methods combine the equation of continuity with the equation of motion for unsteady flow

(Subramanya 2002

. 271

Types of Flood Routing

In all the hydrologic analysis applications mentioned above , two categories of routing can be clearly recognised

Reservoir routing - In this type of routing , the effect of a flood wave entering a reservoir is studied . This is done by determining the volume-elevation characteristic of a reservoir in addition to the outflow-elevation characteristic of the spillway and also other outlet structures present in the reservoirs (Chadwick Morfett 1986 p318 . The results are used to predict the variation of reservoir elevation and outflow discharge with respect to time . This type of routing is necessary for

Designing the capacity of the spillway and other outlet structures

Determining the correct location and size of capacity of the reservoir pertaining to a particular requirement condition

Channel routing - In this type of routing , a study is made of the change in shape of a hydrograph as it travels down a channel . This done by considering a channel reach i .e . the specific length of the stream channel , and an input hydrograph at the upstream end of the stream . The results are used to predict the flood hydrograph at various sections of the reach (Chadwick Morfett 1986

. 322 . The output data obtained using this method is , the information on the flood-peak attenuation and the considered very important for

Flood-forecasting operations

Flood-protection related work

Hydrologic Channel Routing

In case of reservoir routing , the storage is a function of output discharge , whereas in case of channel routing , the storage is a function of both inflow and outflow discharges . This is the main reason why entirely different routing methods are needed for Channel routing

When a river is in flood , the flow can be characterized as gradually varied unsteady flow . In a particular channel reach the water surface as expected is not parallel to the channel both . Additionally it also varies with time

At the time of flood , the two categories Wedge storage - This term represents the wedge-like volume which is formed between the actual water surface pro and the prism storage surface In the downstream section of a river reach , the prism storage is observed to be constant , when the depth is fixed . However , the wedge storage changes from positive to negative depending on the type of flood . The wedge storage is positive at the time of advancing flood while it is negative in case of a receding flood (Subramanya 2002

br 282-283

(Flow Routing 2

.Muskingum Method

Introduction

Flood routing in open channels can be determined using a variety of modeling procedures . These methods follow a wide range of methodologies which can be categorized as

Simple like Muskingum-type approximations - Which have modest data requirements

Out of these the Muskingum and Muskingum-Cunge methods are well established in the hydrological literature , and the modest data requirements make these procedures attractive even though more rigorous hydraulic models are available for unsteady flow routing

(Birkhead James 2002

. 113

For Muskingum method which is a hydrologic method , the discharge measurements alone are sufficient for routing . This is because it is assumed that the parameters of the Muskingum model capture the combined flood-propagating characteristics of a river reach . When the water resources schemes to be built are in their initial planning stages , the river gauging system may remain either underdeveloped or insufficient to give precise and rigorous measurements of flow depths and discharges The Muskingum method is useful for predicting the preliminary outflow hydrographs required at the initial stages of planning spillway capacities . In addition , these outflow hydrographs can also aid the design of stream gauges for future use . Hence , the Muskingum model has a high significance for modern civil engineers (Das 2004

. 130

The Muskingum equation is frequently used for routing of floods in river channels

The Muskingum method for routing flood waves in rivers and channels has been widely used in applied hydrology , since its first use in connection with a flood control project in the Muskingum County of Ohio about fifty years ago . Since its development around 1934 by McCarthy , the Muskingum method has also been a subject of many investigations (Strupczewski Napiorkowski Dooge 2002

.235

The figure below shows the translational and storage processes in stream channel routing

(Gill 1979 22

Procedure for applying Muskingum Method

The Muskingum method as explained above is a widely used hydrologic method for routing flood waves in rivers and channels . The standard procedure for applying the Muskingum method involves two basic steps calibration and prediction

In the calibration step , inflow-outflow hydrograph data is taken from historical records to determine the parameter values for the Muskingum model of a river . This is used to solve the parameter-estimation problem . The prediction step is the solution of a routing problem in which the outflow hydrograph for a given inflow hydrograph is determined by using the routing equations

Figure below shows an example of routing using Muskingum method

(Elizabeth 1994

Disadvantages of Hydrologic Routing

The primary disadvantage is that the dynamic flow effect is not taken into account . In addition it is assumed that stage and storage is a single-valued function of discharge - implying flow is changing slowly with time . This may not be true always , and is a case of oversimplifying things

Muskingum Equations

Muskingum equation is a spatially lumped form of continuity equation and linear-storage and discharge relationship for a specified river reach The first equation is common to all conceptual models . It is a physical continuity equation i .e . law of conservation of mass , which is later integrated over the entire river reach

------ (1

Muskingum method differs from other conceptual models with the second equations , which relates storage in the river reach , inflow and outflow

------ (2

The equation (2 ) given above can be obtained by assuming that the storage in the river reach is a weighted sum of the storages corresponding to steady-state flow rates in terminating cross-sections and that a linear relation exists between the prismatic storage and the flow rate

As we have seen in the sections above the two types of storages prism and wedge , and the functions which they depend on . Revisiting the same we get The above equations can be rewritten as The above two storages In the equations above K is the travel time of flood wave in a reach in hours and x is the weighting factor , its values lies in the range 0 0 .5

The Muskingum equation based on equation (2 ) and the water balance equation can be written as

Where , i 1 , 2 . ------ (2a

Where

------ (2b

------ (2c

------ (2d

is the inflow at the beginning of the time interval (Guang-Te Singh 1992 57-58

Variable Parameters K and x - The parameter K , substituted in the equations above , is usually referred to as the storage constant and has the dimension of time . K is the first moment of the impulse response of the Muskingum model . This means that it measures the delay between the centre of gravity or centroid of the input wave and the centre of gravity of the model output wave . However , studies have shown that this parameter cannot be identified with the transport delay in the sense of theory of dynamic systems , since the Muskingum model never yields ideal translatory properties (Strupczewski Kundzewicz 1980

. 328 Physically speaking , the parameter K represents the average reach travel time , and is the function of length of the river reach and of the wave velocity . Its value and also the value of the parameter x can be analytically evaluated by matching the moments of impulse response of the Muskingum model and of the diffusion analogy or the linearized dynamic wave model (Overton 1966

. 186

While all the Hydrologists agree over the nature of the parameter K , but the views on the nature of the parameter x , a dimensionless quantity and the range of values it varies over are not consistent . According to Yevjevich and Barnes , the parameter x relates to the difference between flow rates in the cross-sections terminating the modeled river reach . It is thus the measure of shape of the wave passing through the reach and in principle , different for different types of waves . In the literature two typical hypotheses for the range of variability of x can be found The most popular opinion is that x as a weight coefficient should have a positive value from the range (0 , 1 , although some hydrologists like Cunge propose the choice of this parameter from the narrower range of values (0 , 0 .5 . One scientist Dooge even permits the possibility of negative values of x . Strupczewski Kundzewicz 1980

. 328-329

Theoretical basis of the variable parameters

(Guang-Te Singh 1992 59-60

In the year 1969 , Cunge derived the equation (2 ) by using a finite approximation of St . Venant 's equations , described in the sections below , but neglected the inertia term , and expressed the variable parameters K and x in terms of physical river characteristics as give below

------ (2e

------ (2f

is the unit width discharge , o is the discharge , B is the top width , S is the channel bed slope and ? is an exponent

The passage of a flood hydrograph through a reservoir or a channel reach is an unsteady-flow phenomenon . The basic differential equations used in the hydraulic routing are known as St . Venant equations which give a good of the unsteady flow

The equation of continuity is as below

------ (3

And the equation of momentum is

------ (4

In the above equations , A is the flow 's cross-sectional area , V is the flow 's average velocity , Q is the discharge (Q AV , Z is the level of water , g is the acceleration due to gravity , R is the hydraulics radius and C is the Chezy 's coefficient

For most of the flow situations , the inertia and acceleration terms in the above equation is much smaller than the friction slope and thus can be neglected . Equation (4 ) is then reduced as below

is the wetted perimeter of the flow cross section . It can be approximated for wide concave channels as

------ (8

Here , a is a constant and the exponent variable m varies between 1 .2 and 2 .0

------ (9

, we get

------ (10

Substituting equation (7 ) in (10 ) we get

------ (11

The equation can be further modified by applying the concepts pf hydraulics of open channel flow . According to this , and ? is an exponent . Its value is 0 .5 for Chezy formula and 2 /3 for Manning formula

Substituting these values in equation (11 ) and simplifying we get

------ (12

is an exponent

Equation (12 ) can be rewritten as

------ (13

becomes 1 .0 , the equation (13 ) becomes the same as equation (2e

X (Dimensionless Coefficient

(Guang-Te Singh 1992 61-62

The variable parameter x , for the Muskingum method can be represented as below

------ (14

Here , l is the characteristic length which can be written as below

------ (15

is the steady-flow discharge

For a steady state flow , the relation between the water level and the discharge is

------ (16

is an exponent . If partial-derivative is taken of the above equation we get

------ (17

------ (18

Substituting equation (18 ) in equation (14 ) we get

------ (19

are related as below , hence equation (19 ) becomes

------ (20

, we get

------ (21

Equation (21 ) is the same as equation (2f

Estimation of variable parameters

The estimation of variable parameters poses a tricky problem , which has attracted the attention of many researchers . The early method for the linear model is based on a trial-and-error graphical approach . A refinement of the trial-and-error-based graphical procedure is the use of the least squares method proposed by Gill in 1978 . For nonlinear models , the earlier researchers like Gill , Gavilan , Houc , Tung , Yoon Padmanabhan , Mohan - have proposed methodologies based on the concept of regression analysis . The regression is performed either between the observed and actual storages or between the observed and actual outflows (Das 2004

. 140

A software called MUPERS (Muskingum parameter estimation and flood routing system , was developed to select between linear and nonlinear forms of Muskingum models for the given data , and to route flows based on the estimated parameters . Eleven methods of estimating parameters are included in the MUPERS . Out of this eight are meant for linear form - TRERR , LSQREG , OTHREG , LEAST , FBOPT , FLTDAT , LP , and QP two for piecewise linearized form - CONSEG and DISSEG and one for nonlinear form of the model - NONLR , are included in MUPERS (Yoon Padmanabhan 1993

. 602

(Yoon Padmanabhan 1993

. 602

Some of the basic methods are as described below

becomes as close as possible to a straight line . The slope of the straight line fitted through the loop gives k . A sample is shown in the figure below (Flow Routing

The value of x that gives a straight line instead of a loop is considered the best x-value , and the storage constant K is estimated from the reciprocal of the slope . This trial and error method is subjective and has no objective-selection criteria for the proper x value . It is also very time consuming

A careful study of this method and its results shows that the least-squares method is a numerical expression of the graphical method hence , the two methods are equivalent . Therefore , the least-squares method should constitute a natural replacement for the graphical method

Mathematically the concept can be explained as below

, then according to covariance formula is the covariance of S and Y ?s is the standard deviation from S and ?y is the standard deviation from Y

It follows numerically that : is the mean of Y

The above expressions shows that for this method , the problem of determining K and x is maximizing r C as a function of x . In fact the same method is followed in the least-squares method , as described in the next section

(Singh McCann 1980

. 357-358

Least-squares method - This method is based on minimizing the sum of squares of deviations between observed storage and computed storage for a given inflow-outflow sequence . Mathematically it can be represented as

------ (22a

is the estimated storage for the jth time interval , N is the number of observations , E is the error function to be minimized

Here It can be seen that K , x and C can be determined objectively and conveniently , even with a calculator

Case B : C 0 , Solving for A and B , in the same way as above , we get

(Singh McCann 1980

. 356

Method of moments - This method is based on finding the moments of the IUH of the Muskingum reach . After finding the moments , they are expressed in terms of I and O

When equations (1 (2 ) are combined we get

------ (23

To solve this equation some initial conditions must be specified . One of these is If equation (23 ) is solved with the above condition we get Where , h (t ) is the instantaneous (IUH ) of the Muskingum reach and can be expressed as below

------ (24

function or impulse inflow

When moments of equation (24 ) are taken we get

------ (25

is the second moment of h (t ) about its centroid . The moments of IUH can be determined from the moments of I and O using Nash theorem . We get Here , U denotes the moment of the particular item . The superscript denotes the of U . Superscript denotes the of U . No superscript denoted the moment about the centroid . Thus , k and x can be obtained objectively and also conveniently by computing by computing the first two moments of inflow and outflow . The first moment is about the origin and other is about the centroid

(Singh McCann 1980

. 357-358

Method of Cumulants - From the above section it can be deduced that Here , K1 , and K2 are the first and second cumulants of IUH for equation (24

can be represented as represents the ith cumulant of quantity

Hence , an objective determination of k and h is possible using this method

As is seen from the above two methods , in this particular case the methods of cumulants and moments are equivalent (Singh McCann 1980

br 357

without estimating k and x , based on minimizing the difference between observed hydrograph and computed hydrograph . The difference can be expressed by an error function defined in a least-squares sense or differently . In this study we will employ the least-squares error function . Therefore , this method is also a least-squares optimization method and is , in principle , equivalent to the least-squares method defined previously

are chosen as the unknowns then then

------ (26

------ (27

is so small that it can be easily performed using a calculator

(Singh McCann 1980

. 359

Translatory solutions using Muskingum method

Scientists Kulandaiswamy , Cunge Gill etc . have taken part in researching translatory characteristics of the Muskingum method . Cunge in 1969 derived a partial differential equation for the Muskingum method , using the Muskingum assumption of a unique stage-discharge relationship in St . Venant equations . Later , he showed that the Muskingum method did not allow for wave damping , and that wave attenuation was due to errors arising from numerical approximation of the Muskingum equation . This is valid if it is assumed that the lumped system can be transformed into a distributed system , the Muskingum method has a dynamic basis and that this basis lies in St . Venant equations . This assumption is no more valid than the Muskingum hypothesis itself and has no particular basis other than to try to predicate a physical basis for a purely empirical assumption (Singh McCann 1980

. 347

The Muskingum method is not a translatory solution . Further , it can be shown that there is no unique dynamical analog of the Muskingum method even if the assumption of its dynamic basis is made

To write the standard flow equation in a finite-difference form using the ?- ?- ?- ?- br

. This is an approximate solution

For any other choice of k , x and ?t even the approximate translatory solution is not possible . Even for x 0 .5 , the approximate translatory solution is possible only if I and Q can be approximated linearly over At . This implies that At should be very small so that this assumption is not violated . Consequently , k will be very small and thence the channel reach under consideration will have to be small enough to have such a small approximate travel time (Gill 1979

.20

If we consider a continuous solution for the Muskingham method , then it is immediately clear that regardless of the choice of k and x this solution is not translatory

(Singh McCann 1980

. 348

Theoretical Analysis

Equation (23 ) derived in sections above is as below To solve this equation some initial conditions must be specified . One of these is If equation (23 ) is solved with the above condition we get

------ (29

This equation was found out by Kulandaiswamy in 1966 and also by Diskin in 1967

Now if we consider

And substitute in equation (29 ) we get for any

Kulandaiswamy additionally pointed out that the equation (29 ) will give translatory solution only when x 0 .5 and I (t ) is such that its third and higher derivates are very small

Also I can be written as below at any interval is unique . This hypothesis is unrealistic , however even if it does satisfy this hypothesis , there may not still be a translatory solution

As an illustration , considering , and substituting equation (29 ) we get This means In the above equation if x 0 .5 and k , we get Now when t k , we get are far from being zero . Obviously their third- and higher- derivatives do not vanish . Therefore , the assumptions are not valid

The above discussion shows that the Muskingum method is not a translatory solution

(Singh McCann 1980

. 347

Non-Linear Muskingum Model

The storage--flow relationship is not always expressible by the linear formula . The Muskingum equation is modified from its original form in this case . If the element of nonlinearity is appreciable , the following equation may be fitted to the data

------ (30

is as given in figure below

(Flow Routing

The constants , x ? and m cannot be obtained by Least Square Method because the resulting normal equations are implicit in character and require a tedious trial-and-error solution . Hence , the value of x should be determined by trial and error

is prepared and a smooth curve is drawn through the mean positions of the plotted data . The method of three points may be used for determining the constants ? and m Substituting these values in equation (30 ) we get

------ (31

------ (32

------ (33

(Gill 1978 356

Rearranging the equations and taking log on both sides we get

------ (34 C

o

p B

k C

D br

p W

D

p KH

D

W

br U

F

KH

jW

h

h

h

h

?i ?i ?y

?y

?y

ja

j

ha

hr

hr

hr

hr

j

p jn

jYA

j

j

j7

h

jeI

j K

h

joZ

j ?I

j ?I

jo /ooI

jCHY

jIY

j Y

jmI

j )K

jYY

joYA

h

jNYA

juYA

h

h

h

h

h .S

h .S

h .S

jNo

jCn

jyp

jw

j ?b

h

h

h

h

h

h

h

h

h

h

h

h

h

j (u

h

h

h

j ?t

h

h

h

h

h

h

j-f

h

j 'u

h

h

h

j`kh

jvkh

h

jQkh

h

h

h

h

j

j

h

h

jo /ooF

j :U

jGsch

j5sch

jzsh

jOsh

jNsch

j ?sch

jhsch

j 'w

j g

j d

------ (35

(Gill 1978 357

The coefficient ? can be determined by trial and error from equation (34 . The exponent m can then be determined from eq .35 . The coefficient ? can be determined using the calculated values of ? and m in any of the equations 31 , 32 and 33

(Gill 1978 357

As can be seen above , using the nonlinear relationship , the routing becomes more difficult . In 1985 a procedure was developed by Tung to perform channel routing using non linear Muskingum model References