Muskingum Method with Variable Parameter Estimation

Flood routing techniques are used for flood forecasting, reservoir design, watershed simulation, and comprehensive water resources projects. One category of flood routing is called the hydrology routing. The most well-known hydrology routing procedure is referred to as the Muskingum Method. McCarthy [1] first developed Muskingum flood routing for flood-control studies in the Muskingum River in Ohio, USA. Muskingum method has been continuously receiving attention from the researchers for several reasons. Firstly, it is a simple method, which can be used for flood routing without much complication as far as the procedural details are concerned. Its parameters can be calculated using the record of past historical floods. It does not require a knowledge of the river bed geometry as the phenomenon can be reproduced well enough on the basis of the calibration carried out using experimental data relative to the extreme sections of most significantly long reaches. Linear Muskingum models are (1) for continuity, and (2) for storage.

  $\frac{d S}{d t}=I-Q$ (1)

$S_{t}=K\left[x I_{t}+(I-x) Q_{t}\right]$  (2)

The variables S_t, I_t, and Q_t are simultaneous storage, inflow, and outflow, respectively, during the passage of a flood through the reach at time t; x and K are constants. Physically speaking, K is considered to represent the average reach travel time and will be equal to the time difference between the centroids of inflow and outflow hydrographs. The coefficient x is used to weigh the relative effects of inflow and outflow on reach storage and can take values from 0 to 0.5.

Traditionally, the parameters are estimated by plotting accumulated storage versus weighted flow of a given reach (Figure 1). Although some reaches may exhibit nonlinear loops, the value of x that gives a nearly collapsed loop, approximating a straight line as determined by visual judgment is considered the best value. The coefficient K is then the reciprocal of the slope of the straight line best representing the loop.

The trial and error graphical approach, which have been used for decades, is time consuming and prone to subjective interpretation. Furthermore, the visual judgement may not correctly identify the best among several nearly collapsed loops when all may appear acceptable.

1.png

Figure 1. Typical linear and nonlinear storage-weighted flow relationship

In case where the storage versus weighted flow relationship is not linear as implied in the original Muskingum model, using a linear form of Muskingum model may introduce considerable error. If the relationship for a reach is discovered to be nonlinear, as seen in Figure 1, the nonlinear models given by Eqns. (3) and (4) may be more appropriate [2-6].

$\mathrm{S}_{\mathrm{t}}=\mathrm{K}\left[\mathrm{x} \mathrm{I}_{\mathrm{t}}-(\mathrm{I}-\mathrm{x}) \mathrm{Q}_{\mathrm{t}}\right]^{\mathrm{n}}$    (3)

$\mathrm{S}_{\mathrm{t}}=\mathrm{K}\left[\mathrm{x} I_{t}^{n}-(\mathrm{I}-\mathrm{x}) Q_{t}^{n}\right]$       (4)

These models have an additional parameter n, which cannot be determined graphically from recorded inflow and outflow hydrographs, and alternative parameter estimation methods must be used.

In this study, we have studied the applicability of linear and nonlinear models on linear and nonlinear flood data. Subsequently, we have routed the flows based on the estimated parameters to compare the performance of these models.

3.1 Trial and error method

In the conventional graphical method, weighted flow given by Eq. (18) is plotted against the accumulated storage given by Eq. (19) for several selected values of x. The value of x which gives the narrowest loop is accepted as the Muskingum x. The narrowest loop is then represented by a straight line, and the reciprocal of the slope of this line is taken as storage constant K.

$S_{2}=S_{1}+(\overline{I}-\overline{Q}) \Delta t$    (19)

where,

$\overline{I}=\frac{I_{n}+I_{n+1}}{2}, \overline{Q}=\frac{Q_{n}+Q_{n+1}}{2}$ Weighted flow $=\mathrm{x} \mathrm{I}_{\mathrm{nt} 1}+(1-\mathrm{x}) \mathrm{Q}_{\mathrm{nt} 1}$    (20)

Instead of drawing the accumulated storage against weighted flow for each trial value of x, and then choosing the best fit line by visual judgement, I developed a computer program to compute the sum of the squares of the deviations from the best fit line. For each particular value of x, pairs of values for accumulated storage and weighted flow were computed using the given inflow and outflow hydrographs ordinates. The best fit line through these points was then determined by the method of least squares with the condition that this line should pass through the initial point when I₀ = Q₀. For each value of x, and the corresponding straight line giving the best fit line between weighted discharge and accumulated storage, the error sum of squares was recorded, and the value of x resulting in the least error sum of squares will be chosen as the Muskingum parameter x. The storage constant K will be the reciprocal of the slope of the best fit line. This procedure is summarized below.

Let  y = weighted discharge

b = slope of the straight line.

Then,

y = y₀ + b (S – S₀)    (21)

And

  $s u m=\sum\left[y-\left[y_{0}+b\left(S-S_{0}\right)\right]\right]^{2}=\sum\left[\mathrm{y}-\mathrm{y}_{0}-\mathrm{b} \mathrm{S}+\mathrm{b} \mathrm{S}_{0}\right.$    (22)

                      $\frac{d(s u m)}{d b}=0=\sum 2\left(y-y_{0}-b S+b S_{0}\right)\left(-S+S_{0}\right)-\sum y S+y_{0} \sum S+b \sum S^{2}-b S_{0} \sum S+S_{0} \sum y-n y_{0} S_{0}-$

$b S_{0} \sum S+n b S_{0}^{2}=0$    (23)

where n = number of data points except the initial point.

$b\left[\Sigma S^{2}-2 S_{0} \sum S+n S_{0}^{2}\right]=\sum y S-y_{0} \sum S-S_{0} \sum y+n y_{0} S_{0}$    (24)

$b=\frac{\sum y S-y_{0} \sum S-S_{0} \sum y+n y_{0} S_{0}}{\sum S^{2}-2 S_{0} \sum S+n S_{0}^{2}}$    (25)

If the initial storage is considered as zero,

sum = ⅀ (y – y₀ – b S)²    (26)

$\frac{d(s u m)}{d b}=0=\sum 2\left(y-y_{0}-b S\right)(-S)$    (27)

⅀ y S – y₀ ⅀S – b ⅀ S² = 0                       (28)

Then

$b=\frac{\sum y S-y_{0} \sum S}{\sum S^{2}}$    (29)

It should be noted that Eq. (29) is a special case of Eq. (25) for S_o= 0.

3.2 Least square 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:

$E=\sum_{j=1}^{N}\left[S_{0}(j)-S_{e}(j)\right]^{2} \rightarrow \min$    (30)

where S₀ (j) is the observed storage for the jth time interval; S_e (j) is the estimated storage for the jth time interval; and N is the number of observations. E is the error function to be minimized. Two cases are distinguished.

Case A: C ≠ 0 where C is the difference between absolute and relative storages. Eq. (21) can be written as (dropping j for brevity and assuming A = K x and B =K (1-x):

$\mathrm{E}=\sum_{j=1}^{N}\left[S_{0}-\mathrm{x} \mathrm{K} \mathrm{I}-\mathrm{K}(1-\mathrm{x}) \mathrm{Q}-\mathrm{C}\right]^{2}$    (31)

Thus we obtain:

$A=\frac{y_{1}}{y_{2}}-B\left(\frac{y_{3}}{y_{2}}\right)$, $B=\frac{y_{1} z_{2}-z_{1} y_{2}}{z_{2} y_{3}-y_{2} z_{1}}$ and $C=\frac{\sum S_{0}-A \sum I-B \sum Q}{N}$    (32)

where

$y_{1}=\sum S_{0} I-\frac{\left(\sum s_{0} \sum I\right)}{N}$, $y_{2}=\sum I^{2}-\frac{(\Sigma I)^{2}}{N}$ and $y_{3}=\frac{\sum Q I-\Sigma Q \Sigma I}{N}$    (33)

$z_{1}=\sum S_{0} Q-\frac{\left(\sum S_{0} \sum Q\right)}{N}$, $z_{2}=\sum I Q-\frac{(\Sigma I \Sigma Q)}{N}$ and $\Sigma Q^{2}-\frac{(\Sigma Q \Sigma Q)}{N}$    (34)

Therefore:

K = A+ B  and    X= A/(A+B)    (35)

Thus K, x and C can be determined objectively and conveniently.

Case B: C = 0. Solving for A and B as before:

$A=\frac{\left(\sum s_{0} I \sum Q^{2}-\sum s_{0} Q \sum I Q\right)}{D}$, $B=\frac{\left(\sum s_{0} I \sum I^{2}-\sum S_{0} Q \sum I Q\right)}{D}$ and $D=\sum I^{2} \sum Q^{2}-\left(\sum I Q\right)^{2}$    (36)  

And hence K and x can be determined from Eq. (34).

3.3 Direct Optimization (DO) method

The ditect optimization (DO) method determines directly the routing coefficients C₀, C₁, and C₂ without estimating K and x. The method is based on minimizing the difference between observed outflow hydrograph and computed hydrograph. Therefore, this method is also a least-squares optimization method and is in principle, equivalent to the least squares method defined previously.

There are, in fact, only two unknowns since the third is known from C₀ + C₁ + C₂ = 1. If we choose C₁, and C₂ to be the unknowns then:

$x=\frac{c_{1}+0.5 c_{2}-0.5}{c_{1}+c_{2}}$    (37)

$K=\frac{\Delta t\left(C_{1}+C_{2}\right)}{1-C_{2}}$    (38)

Further, rearranging

C₁ (I_n – I_n-1) + C₂ (I_n – Q_n-1) = I_n - Q_n    (39)

If we define:

R_n = I_n – Q_n; F_n = I_n – I_n-1; and G_n = I_n – Q_n-1    (40)

Then

R_n = C₁ F_n + C₂ G_n    (41)

The error function (dropping the subscript n for brevity) then follows:

$E=\Sigma_{1}^{N}\left(R_{0}-R_{e}\right)^{2} \rightarrow \min$    (42)

where subscripts 0 and e denote observed the estimated R values, respectively.

Following the usual procedure:

⅀R₀ F = C₁ ⅀F² + C₂ ⅀F G and ⅀R₀ G = C₁ ⅀F G + C₂ ⅀G²    (43)

Solving for C₁ and C₂:

$C_{1}=\frac{\left(\sum R_{0} F \sum G^{2}-\Sigma R_{0} G \sum F G\right)}{D E T}$ and $C_{2}=\frac{\left(\sum R_{0} G \sum F^{2}-\sum R_{0} F \sum F G\right)}{D E T}$    (44)

where

$D E T=\Sigma G^{2} \Sigma F^{2}-(\Sigma F G)^{2}$    (45)

Thus, once C₁ and C₂ are determined, K and x can be obtained.

3.4 Nonlinear method

Through an analysis of stream flow records at the upstream and downstream ends of a river reach, it is possible to construct a set of data comprising the inflow rates at the upstream end of the reach, the outflow rates at downstream end of the reach and the storage volumes in the reach at selected times over some historical period. These data may be denoted by:

S_t = actual or measured storage in the reach at time t

I_t = actual or measured inflow rate at the top of the reach at time t

Q_t = actual or measured outflow rate at the bottom of the reach at time t.

The standard, nonlinear Muskingum model for channel routing assumes a relationship between storage, inflow and outflow defined by Eq. (44). A standard selection procedure for parameters of the model is to minimize the sum of the squared deviations of model-predicted storage and actual storage. One way to define this procedure for the standard nonlinear model is with the mathematical model given by Eq. (47) and 48.

$S_{t}=K\left[x I_{t}^{n}+(1-x) Q_{t}^{n}\right]$    (46)

$N=\sum_{t}\left(S_{t}-P_{s t}\right)^{2}$    (47)

  $P_{s t}=K\left[x I_{t}^{n}+(1-x) Q_{t}^{n}\right]$    (48)

where P_st =model - predicted storage value;

N = error sum of squares.

One of the solution method for the mathematical program is: 1- select a group of values for n (e.g. n = 0.05, 0.10, ...); for each value of n, solve the resulting mathematical model using traditional calculus methods by setting all partial derivatives equal to zero and solving for the optimal values of K and x. 2- choose that value of n and the associated values of K and x that produced the minimum value of the objective function, N.

The standard nonlinear Muskingum model is a more general model and should perform better in minimizing the sum of the squares of deviations of actual storage from model-predicted storage values.

If we let

P = K x

then Eqns. (23), (24), and (25) yield

$N=\Sigma\left(S_{t}-P I_{t}^{n}-K Q_{t}^{n}+P Q_{t}^{n}\right)^{2}=\Sigma\left[S_{t}-P\left(I_{t}^{n}-Q_{t}^{n}\right)-K Q_{t}^{n}\right]^{2}$    (49)

$\frac{d(N)}{d P}=0$

$\Sigma 2\left[S_{t}-P\left(I_{t}^{n}-Q_{t}^{n}\right)-K Q_{t}^{n}\right]\left(Q_{t}^{n}-I_{t}^{n}\right)=\Sigma S_{t} Q_{t}^{n}-\Sigma S_{t} I_{t}^{n}+(k-2 P) \sum Q_{t}^{n} I_{t}^{n}+(P-K) \sum Q_{t}^{n}+P \sum I_{t}^{2 n}=0$    (50)

$\frac{d(N)}{d P}=0$

$\Sigma 2\left[S_{t}-P\left(I_{t}^{n}-Q_{t}^{n}\right)-K Q_{t}^{n}\right]\left(-Q_{t}^{n}\right)=0$    (51)

$\Sigma Q_{t}^{n} I_{t}^{n}-\Sigma S_{t} Q_{t}^{n}+(K-P) \Sigma Q_{t}^{2 n}=0$    (52)

Now let,

A₁ = ⅀ S_t$Q_{t}^{n}$, A₂ = ⅀ S_{t$I_{t}^{n}$}, A₃ = ⅀ $Q_{t}^{n}$$I_{t}^{n}$ , A₄ = ⅀ $Q_{t}^{n}$ and A₅ = ⅀ $I_{t}^{n}$    (53)

Then Eq. (28):

A₁ - A₂ + (K- 2P) A₃ + (P- K) A₄ + P A₅ = 0    (54)

Eq. (29):

PA₃ -A_I + (K- P) A₄ = 0    (55)

From Eq. (29) solve for P

$P=\frac{A_{1}-K A_{4}}{A_{3}-A_{4}}$    (56)

Substitution of Eq. (30) into (28) gives:

A₁ – A₂ + K (A₃ – A₄) + P (A₄ + A₃ -2A₃) = 0        (57)

Solving for K we get

$K=\frac{\left(A_{2}-A_{1}\right)-P\left(A_{4}+A_{5}-2 A_{3}\right)}{A_{3}-A_{4}}$    (58)

Substituting Eq. (24) into (25) yields

$K\left(A_{3}-A_{4}\right)+\left(A_{1}-A_{2}\right)+\frac{A_{1}\left(A_{4}+A_{5}-2 A_{1}\right)}{A_{3}-A_{4}}-\frac{K A_{4}\left(A_{4}+A_{5}-2 A_{3}\right)}{A_{3}-A_{4}}=0$    (59)

$K\left(\left(A_{3}-A_{4}\right)-\frac{A_{4}\left(A_{4}+A_{5}-2 A_{3}\right)}{A_{3}-A_{4}}\right)=\left(A_{2}-A_{1}\right)-A_{1} \times \frac{\left(A_{4}+A_{5}-2 A_{3}\right)}{A_{3}-A_{4}}$

Let

$B_{1}=\left(A_{3}-A_{4}\right)-\frac{A_{4}\left(A_{4}+A_{5}-2 A_{3}\right)}{A_{3}-A_{4}}$ and $B_{2}=\left(A_{2}-A_{1}\right)-\frac{A_{1}\left(A_{4}+A_{5}-2 A_{3}\right)}{A_{3}-A_{4}}$    (60)

Then:

$K=\frac{B_{2}}{B_{1}}$ and $x=\frac{p}{K}$    (61)

And

$N=S U M=\Sigma\left(S_{t}-P\left(I_{t}^{n}-Q_{t}^{n}\right)-K Q_{t}^{n}\right)^{2}$    (62)

The nonlinear model requires Q₁ to be solved by trial and error at every time step of flow routing. The routing equation for nonlinear model is given by the following equation.

$\left[\frac{I_{1}+I_{2}}{2}-\frac{Q_{1}+Q_{2}}{2}\right] \Delta t=S_{2}-S_{1}=K\left[\left(x I_{2}+(1-x) Q_{2}\right)^{n}-\left(x I_{1}+(1-x) Q_{1}\right)^{n}\right]$    (63)

I solved this equation by using Newton method as part of a Basic computer program.

Let

$F\left(Q_{2}\right)=\left[\frac{I_{1}+I_{2}}{2}-\frac{Q_{1}+Q_{2}}{2}\right] \Delta t-K\left[\left(x I_{2}+(1-x) Q_{2}\right)^{n}-\left(x I_{1}+(1-x) Q_{1}\right)^{n}\right]=0$    (64)

The procedure is as follows:

(1) At first choose an initial approximation for Q₂.

(2) Check $\frac{d f\left(Q_{2}\right)}{d Q_{2}} \neq 0$

If $\frac{d f\left(Q_{2}\right)}{d Q_{2}}=0$

then another approximation for Q₂ should be chosen.

(3) Calculate a new value of Q₂.

$Q_{2}(n e w)=Q_{2}-\left[\frac{F\left(Q_{2}\right)}{F \prime\left(Q_{2}\right)}\right]$    (65)

(4) If |Q₂ (new)-Q₂| < 0.0001, then

Q₂ =Q₂ (new)    (66)

However, if |Q₂ (new)- Q₂ | > 0.0001, then a new value of Q₂ should be calculated as described in step 3.

The data for the nonlinear model is given in table 3.

Table 3. Actual inflow and outflow data for nonlinear model*

*[19]

The parameters K and x, and various statistics are carried out as explained earlier and the results are shown in table 4.

Table 4. Parameter estimation and comparative statistics for the nonlinear method

6.png

Figure 6. Actual and estimated outflow relationship for Nonlinear method