Variable-parameter Muskingum-Cunge method
revisited


Victor M. Ponce and P. V. Chaganti

Online version 2015

[Original version 1994]



ABSTRACT

This paper reexamines the variable-parameter Muskingum-Cunge method. Numerical experiments show that the loss of mass is small but perceptible throughout the wide range of peak inflows tested, from 200 to 1000 cfs ft-1 (18.58 - 92.9 m-2s-1). Modified three-point direct and four-point iterative variable-parameter methods are introduced. These methods have better mass conservation properties than comparable conventional methods. In practice, the small loss of mass experienced by the variable-parameter methods does not appear to constitute a significant drawback.


1.  INTRODUCTION

The Muskingum-Cunge method of flood routing is well established in the hydrologic engineering literature (Cunge, 1969; Ponce and Yevjevich, 1978; U.S. Army Corps of Engineers, 1990). It is a convenient method because the routing parameters are a function of channel properties and grid specification, which leads to results which are independent of grid size. The method has linear and nonlinear modes. In the linear mode, average flow values are used to calculate the routing parameters at the start of the computation, and these are kept constant throughout the computation in time. In the nonlinear mode, the routing parameters are recalculated for every computational cell as a function of local flow values.

Both modes of computation have their advantages and disadvantages. The linear mode conserves mass, but it requires a priori estimation of a reference flow on which to base the computation of the routing parameters. In addition, the routed flows lack the steepening exhibited by inbank flood waves. Conversely, the nonlinear mode does not require the specification of a reference flow and can simulate the wave steepening. Unfortunately, the nonlinear mode is saddled with a small but perceptible loss of mass.

Given the trade-offs involved, it seems certain that both routing modes will continue to be used in the future. The linear mode will be used where simplicity is desired, notwithstanding its lack of wave steepening. The nonlinear mode will be used where the simulation of wave steepening is judged to be important, albeit at the cost of a small loss of mass. A multilinear method, representing a compromise between linear and nonlinear methods, has been developed recently by Perumal (1992).

The Muskingum-Cunge is a viable alternative to the classical Muskingum method, particularly for the cases where hydrologic data (i.e., streamflow data) are not available, but where hydraulic data (cross-sectional data and channel slopes) can be readily ascertained. In many instances, the Muskingum-Cunge method is also an alternative to the more complex dynamic wave models, which lack robustness and have significant data requirements. The nonlinear features of the variable-parameter Muskingum-Gunge method makes it the method of choice in hydrologic flood routing.

This paper revisits the variable-parameter Muskingum-Cunge method of Ponce and Yevjevich (1978). The small but perceptible loss of mass is confirmed throughout a wide range of unit peak inflows, from 200 to 1000 cfs ft-1 (18.58-92.9 m-2s-1). Furthermore, a slight improvement in mass conservation is realized by an alternative way of calculating the variable parameters.


2.  BACKGROUND

The Muskingum-Cunge method is a variant of the Muskingum method (Chow, 1959) developed by Cunge (1969) and documented in the Flood Studies Report (Natural Environment Research Council, 1975). Ponce and Yevjevich (1978) expressed the routing parameters of the Muskingum-Cunge method in terms of the Courant and cell Reynolds numbers, two physically and numerically meaningful parameters. In addition, they developed a nonlinear mode of calculating the parameters, thereby enhancing the method's applicability to real-world routing problems. A three-point method and an iterative four-point method were suggested as a way to vary the parameters as a function of local flow values. The method has been adopted by the most recent version (Version 4.0) of the US Corps of Engineers' HEC-1 model (1990).

The routing equation of the Muskingum-Cunge method defined in the typical four-point grid configuration is:

   n + 1                 n + 1            n             n           
Q j + 1     =   C0Qj       +  C1Qj   +   C2Qj + 1

(1)

in which j is a spatial index, n is a temporal index, and:

C0 = (-1 + C + D) / (1 + C + D)

(2)

C1 = (1 + C - D) / (1 + C + D)

(3)

C2 = (1 - C + D) / (1 + C + D)

(4)

are the routing coefficients, with C the Courant number, defined as follows:

C = ct / Δx)

(5)

in which c is the wave cerelity, Δt is the time interval, Δx is the space interval, and D the cell Reynolds number, defined as follows:

D = q /(So c Δx)

(6)

in which q is the unit-width discharge and So is the bottom slope.

The wave celerity is defined as follows:

c = β (Q /A) = β (q / d )

(7)

in which β is the exponent of the rating, A is the flow area, and d is the flow depth.

The Muskingum-Cunge method matches the numerical diffusion of the discrete model with the physical diffusion of the analytical model. The advantage is that the routing results are independent of the grid specification, provided numerical dispersion is minimized. The numerical dispersion can be minimized by keeping the Courant number equal to or slightly greater than one (Cunge, 1969; Ponce and Theurer, 1982; Ferrick et. al., 1983).


3.  NUMERICAL EXPERIMENTS

In this paper, Thomas' (1934) classical problem was chosen to test the Muskingum-Cunge method under a wide range of flow conditions. The Thomas problem was chosen here to allow comparison with previous results (Ponce and Yevjevich, 1978). The original problem considered routing a sinusoidal flood wave with a 96-h period through a prismatic channel 500 miles long. A unit-width channel of baseflow qb = 50 cfs ft-1 (4.64 m2 s-1) and peak inflow qpi = 200 cfs ft-1 (18.58 m2 s-1) was specified. The bottom slope was set at 1 ft mile-1, and the discharge-depth rating was the following:

q = 0.688 d 5/3

(8)

This amounts to a Manning's n = 0.0297. It can be shown that the 96-h sinusoidal flood wave satisfies the diffusion wave applicability criterion (Ponce et. al., 1978).

The following methods are considered here:

  1. The constant-parameter method (CPMC), in which the routing parameters C and D are the same for all computational cells. They are based on an average unit-width discharge qa calculated as:

    qa = (qb + qpi) / 2

    (9)

    The average celerity ca is calculated with Eq. 7, using the discharge computed with Eq. 9.

  2. The three-point variable-parameter method (VPMC3) (Ponce and Yevjevich, 1978), in which the routing parameters C and D for each computational cell are based on the average unit-width discharge and celerity at the three known grid points:

                n         n                n + 1                      
    qa =  (qj  +   q j + 1   +    q j     ) / 3

    (10)

                 n         n              n + 1                      
    ca =  ( c j   +  c j + 1  +   c j      ) / 3

    (11)

  3. The iterative four-point variable-parameter method (VPMC4) (Ponce and Yevjevich, 1978), in which the routing parameters C and D for each computational cell are based on the average unit-width discharge and celerity at the four grid points:

                 n         n               n + 1        n + 1               
    qa =  (q j   +  q j + 1   +   q j      +   q j + 1    ) / 4

    (12)

                n         n               n + 1        n + 1               
    ca =  (c j   +  c j + 1   +   c j      +   c j + 1    ) / 4

    (13)

    To improve the convergence of the iterative procedure, the discharge qj + 1n + 1 obtained with the three-point method (VPMC3) is used as a first guess of the iteration. The celerity cj + 1n + 1 is calculated with Eq. 7, based on qj + 1n + 1. The iteration converges to a specified tolerance after a few trials.

  4. A modified three-point variable-parameter method (MVPMC3), in which the routing parameters C and D for each computational cell are based on the average unit-width discharge at the three known grid points:

                 n          n               n + 1                      
    qa =  (q j   +   q j   +  1  +  q j      ) / 3

    (14)

    Unlike the VPMC3, the average celerity c a is calculated with Eq. 7, based on the discharge computed with Eq. 14.

  5. A modified iterative four-point variable-parameter method (MVPMC4), in which the routing parameters C and D for each computational cell are based on the average unit-width discharge at the four grid points:

                  n        n                 n + 1      n + 1               
    qa =  ( q j   +  q j + 1   +    q j    +   q j + 1    ) / 4

    (15)

    Unlike the VPMC4, the average celerity c a is calculated with Eq. 7, based on the discharge qa computed with Eq. 15. To improve the convergence of the iterative procedure, the discharge qj + 1n + 1 obtained with the three-point method (MVPMC3) is used as a first guess of the iteration. The iteration converges to a specified tolerance after a few trials.

Two levels of spatial and temporal resolution were chosen:

  1. Resolution level I, the grid specification used by Thomas (1934): Δx = 25 miles (40.22 km), and Δt = 6 h; and

  2. Resolution level II, a finer grid specification used in this paper: Δx = 12.5 miles (20.11 km), and Δt = 3 h.

Three peak inflow to baseflow ratios qpi /qb were used:

  1. qpi /qb = 4, used by Thomas (1934), Dooge (1973), and Ponce and Yevjevich (1978);

  2. qpi /qb = 10; and

  3. qpi /qb = 20.

These peak inflow to baseflow ratios encompass a wide range of practical values. With qpi = 50 cfs ft-1 (4.64 m2s-1), the peak inflow range tested was 200 - 1000 cfs ft-1 (18.58 - 92.9 m2s-1).


4.  RESULTS

The mass conservation properties were assessed by integrating inflow and outflow hydrographs for the 500-mile reach. Table 1 shows peak outflow qpo, ratio of peak outflow to peak inflow qpo / qpi, time-to-peak tp, and mass conservation M for:

  1. Two resolution levels: I and II;

  2. Three peak inflow to baseflow qpi / qb ratios: 4, 10 and 20; and

  3. Five methods: (a) CPMC; (b) VPMC3; (c) VPMC4; (d) MVPMC3; and (e) MVPMC4.


TABLE 1.  Summary of routing results.
Resolution
level
(1)
qpi /qb
(2)
Method
(3)
qpo
(cfs ft-1)
(4)
qpo /qpi
(5)
tp (h)
(6)
M (%)
(7)
I 4 CPMC 176.231 0.881 126 100.000
VPMC3 175.097 0.875 120 99.213
VPMC4 175.860 0.879 120 99.464
MVPMC3 175.437 0.877 120 99.334
MVPMC4 176.084 0.880 120 99.559
I 10 CPMC 437.448 0.875 108 100.000
VPMC3 435.592 0.871 102 97.180
VPMC4 436.968 0.874 102 98.309
MVPMC3 435.784 0.872 102 97.904
MVPMC4 436.939 0.874 102 98.875
I 20 CPMC 876.675 0.877 96 100.000
VPMC3 876.628 0.877 90 95.370
VPMC4 878.236 0.878 90 97.444
MVPMC3 876.160 0.876 90 96.700
MVPMC4 877.786 0.878 90 98.441
II 4 CPMC 176.561 0.883 129 100.000
VPMC3 174.146 0.871 123 99.330
VPMC4 174.370 0.872 123 99.468
MVPMC3 174.189 0.871 123 99.368
MVPMC4 174.430 0.872 120 99.498
II 10 CPMC 438.666 0.877 105 100.000
VPMC3 428.704 0.875 102 97.582
VPMC4 429.937 0.860 99 98.224
MVPMC3 429.094 0.858 102 97.847
MVPMC4 430.580 0.861 99 98.429
II 20 CPMC 884.788 0.885 93 100.000
VPMC3 866.356 0.866 87 95.943
VPMC4 870.975 0.871 87 97.163
MVPMC3 868.557 0.869 87 96.477
MVPMC4 872.360 0.872 87 97.574
Note:  The peak outflows shown in Column 4 are those corresponding to the grid locations, and are not necessarily the actual peak outflows, which may lie in between.

The following conclusions are drawn from Table 1:

  1. All five methods give approximately the same peak outflow qpo and time-to-peak tp.

  2. The ratio of peak outflow to peak inflow (qpo /qpi) is approximately the same for all three peak-inflow-to-baseflow ratios. This suggests that wave diffusion is independent on the qpi /qb ratio.

  3. The constant-parameter method (CPMC) conserves mass exactly. In contrast, the variable-parameter methods are subject to a small but perceptible loss of mass. The loss of mass is cumulative, and can be assessed quantitatively by routing through a very long reach. In this case, the routing of a sinusoidal hydrograph through a 500-mile reach produced mass losses ranging from less than 2.5% for MVPMC4 to less than 5% for VPMC3.

  4. The loss of mass is greater for the three-point methods (VPMC3 and MVPMC3) than for the four-point methods (VPMC4 and MVPMC4).

  5. The loss of mass is greater for the conventional methods (VPMC3 and VPMC4) than for the modified methods (MVPMC3 and MVPMC4).

  6. The loss of mass increases with qpi /qb ratio, and does not appear to be strongly related to resolution level.

It should be noted that the mass conservation percentages shown in Column 7 of Table 1 are cumulative values for the 500-mile reach of the Thomas problem. A typical flood routing application would normally not consider such long reach without intervening lateral inflows which tend to mask the accuracy of the computation. Therefore, the percentages shown should be interpreted as an extreme or worst-case scenario.

Figure 1 shows a set of flood hydrographs for resolution level II and peak inflow to baseflow qpi /qb = 20. It is seen that the constant-parameter method (CPMC) lacks the nonlinear tendency to steepening of the rising limb. (It is worth noting that overbank flows have the opposite tendency, i.e., the flattening of the rising limb). On the other hand, all four variable-parameter methods show steepening of the rising limb, coupled with a corresponding flattening of the receding limb.


Fig. 1.   Outflow hydrographs for Resolution Level II and peak inflow to baseflow ratio qpi /qb = 20.


5.  CONCLUSIONS

The modified variable-parameter methods MVPMC3 and MVPMC4 developed in this paper result in a definite improvement in mass conservation, as compared with the conventional methods VPMC3 and VPMC4 (Ponce and Yevjevich, 1978). The improvement may be more marked when using a wide range of realistic discharges. Otherwise, the small loss of mass does not appear to impose a significant drawback. It is the price to be paid to model the nonlinear features of a flood wave within the framework of the Muskingum-Cunge method.


REFERENCES

Chow, V. T., 1959. Open-Channel Hydraulics. McGraw-Hill, New York.

Cunge, J. A., 1969. On the subject of a flood propagation computation method (Muskingum method). J. Hydraul. Res., 7(2) 205-230.

Dooge, J. C. I., 1973. Linear theory of hydrologic systems. U.S. Dep. Agric. Tech. Bull., 1468, 327 pp.

Ferrick, M. G., Bilmes, J. and Long, S. E., 1983. Modeling rapidly varied flow in tailwaters. Water Resour. Res., 20(2): 271-289.

Natural Environment Research Council, 1975. Flood Studies Report, Vol. III. Natural Environment Research Council, London.

Perumal, M., 1992. Multilinear Muskingum flood routing method. J. Hydrol., 133: 259-272.

Ponce, V. M. and Theurer, F. D., 1982. Accuracy criteria in diffusion routing. J. Hydraul. Div., ASCE., 108 (HY6): 747-757.

Ponce, V. M. and Yevjevich, V., 1978. Muskingum-Cunge method with variable parameters. J. Hydraul. Div. ASCE, 104(HY12): 1663-1667.

Ponce, V. M., Li, R-M. and Simons, D. B., 1978. Applicability of knematic and diffusion models. J. Hydraul. Div. ASCE, 104(HY3): 353-360.

Thomas, H. A., 1934. The hydraulics of flood movement in rivers. Eng. Bull., Carnegie Inst. of Technology, Pittsburgh, PA.

US Army Corps of Engineers, 1990. HEC-1, flood hydrograph package: User's manual, version 4, 1990. US Army Corps of Engineers, Hydrologic Engineering Center, Davis, CA.


150729 19:45

Documents in Portable Document Format (PDF) require Adobe Acrobat Reader 5.0 or higher to view; download Adobe Acrobat Reader.