1. INTRODUCTION
The Muskingum method (1) and its improved MuskingumCunge version (2,4) are well established in the flood routing literature. The Muskingum method is based on the assumption of a linear relationship between the inflow, I, the outflow, 0, and the reach storage, V, of the following form:
V = K [X I + (1  X ) O ]
 (1) 
in which K and X are the parameters. In the conventional Muskingum method
these parameters are determined by calibration using measured inflow and outflow
hydrographs. On the other hand, in the MuskingumCunge version K and X
are calculated using the formulas derived by Cunge (2).
The experience with the MuskingumCunge version is reported herein. It is
shown that the way of calculating the parameters has a definite bearing on
the overall accuracy of the method. As shown by Dooge (3), the assumption
of constant parameters makes the solution dependent on the reference values
chosen to evaluate these parameters. A more physically realistic approach is
to consider the parameters K and X to vary in time and space according to
the flow variability (6). Koussis (5) has considered a dischargedependent K
but has assumed X constant on the grounds that the computation is relatively
insensitive to this parameter. In general, however, it is desirable to allow both
K and X to vary with the flow.
2. MUSKINGUM METHOD
The formula for the Muskingum Method is (Fig. 1) (1):
_{n +1} _{n} _{n +1} _{n}
Q_{ j +1} = C_{ 1} Q_{ j} + C_{ 2} Q_{ j} + C_{ 3} Q_{ j +1}
 (2) 
in which
( Δt / K ) + 2X
C_{1} = ^{_____________________}
2(1  X ) + ( Δt / K )
 (3) 
( Δt / K )  2X
C_{2} = ^{______________________}
2(1  X) + ( Δt / K )
 (4) 
2 (1  X )  ( Δt / K )
C_{3} = ^{_______________________}
2 (1  X ) + ( Δt / K )
 (5) 
Fig. 1 Spacetime discretization of Muskingum method.


where Δt is the routing period (time interval). In the MuskingumCunge version,
the parameters K and X are calculated by (2,3,4);
1 q
X = ^{___} ( 1  ^{__________}
)
2 S_{o} c Δx
 (7) 
in which Δx = reach length (space interval); c = flood wave celerity; q =
unit width discharge; and S_{o} = channel bed slope. Substituting Eqs. 6 and
7 into Eqs. 35:
1 + C  D
C_{1} = ^{______________}
1 + C + D
 (8) 
1 + C + D
C_{2} = ^{_______________}
1 + C + D
 (9) 
1  C + D
C_{3} = ^{______________}
1 + C + D
 (10) 
in which C = c Δt /Δx is the Courant number; and D = (q /S_{o})/cΔx is a type of cell Reynolds number. Both C and D have physical and numerical significance, C being a ratio of celerities and D a ratio of diffusivities.
3. VARIABLE PARAMETERS
Usually, Δt is fixed, and Δx and S_{o} are specified for each computational cell consisting of four grid points (Fig. 1). Therefore, it is necessary to determine
the flood wave celerity, c, and the unit width discharge, q, for each computational
cell. The values of c and q at grid point (j, n) are defined by:
in which Q = discharge; A = flow area; and B = top width. The following
ways of determining c and q for use in calculating C and D are investigated:
Directly, by using a twopoint average of the values at grid points (j, n)
and (j + 1, n);
Directly, by using a threepoint average of the values at
grid points (j, n), (j + 1, n), and (j, n + 1); and
By iteration, using a fourpoint average calculation. To improve convergence,
the values at (j + 1, n + 1) obtained by the threepoint average method are used as the first guess of the iteration.
4. EXAMPLE
The MuskingumCunge method with variable parameters is applied herein
to the problem posed by Thomas (8) in his classical paper on flood routing.
The problem consists of tracking the travel and subsidence of a flood wave
of sinusoidal shape in a unit width channel with a steadystate rating curve given by:
The inflow hydrograph is defined by
π t
I(t) = 125  75 cos ( ^{____} ), 0 ≤ t ≤ 96;
48
 (14a) 
in which t is in hours.
Thomas applied an approximate method to route the flood wave through
a channel 200 miles (322 km) long, using a time interval of Δt = 12 hr. Thomas'
approximate method neglects the inertia terms; therefore, his results are directly
comparable to those of the Muskingum method (both the Thomas and Muskingum
methods can be considered as numerical analogs of the diffusion wave equation).
Figure 2 and Table 1 summarize the results of computations using the MuskingumCunge
method with variable parameters. For comparison purposes, the calculations using constant parameters for three reference values of discharge are
also shown.
Fig. 2 Calculated hydrographs as described in Table 1.


TABLE 1. Summary of results.

Hydrograph
(1) 
Station, in miles
(2) 
Method^{ a}
(3) 
Δx, in miles
(4) 
Δt, in hours
(5) 
Peak q, in cubic feet per second
(6) 
Time to peak, in hours
(7) 
Mass conservation, as a percentage
(8) 
A 
0 
 
 
 
200 
48 
 
B 
500 
MC/200 
25 
6 
178.5 
114 
100 
C 
500 
MC/125 
25 
6 
177 
128 
100 
D 
500 
MC/50 
25 
6 
173.5 
162 
100 
E 
500 
VPMC2 
25 
6 
171 
124 
85 
F 
500 
VPMC3 
25 
6 
175.5 
121 
98 
G 
500 
VPMC4 
25 
6 
176.5 
121 
99 
H 
200 
VPMC4 
25 
12 
190 
77 
100 
I 
200 
Thomas (8) 
25 
12 
189 
79 
 
^{a}
MC/200= constant parameter MuskingumCunge Method, reference discharge 200
cfs (5.66 m^{3}/s); MC/125 = constant parameter MuskingumCunge Method, reference
discharge 125 cfs (3.54 m^{3}/s); MC/50 = constant parameter MuskinguaCunge Method,
reference discharge 50 cfs (1.41 m^{3}/s), VPMC2 = twopoint variable parameter MuskingumCunge; and VPMC3 = threepoint variable parameter MuskingumCunge, VPMC4
= fourpoint variable parameter MuskingumCunge.

The examination of Fig. 2 enables the following conclusions to be drawn:
The MuskingumCunge method with constant parameters (hydrographs B,
C, and D) shows results that are dependent on the value of reference discharge
chosen to evaluate the constant parameters. The higher the reference discharge,
the faster the rate of travel and the lesser the subsidence of the flood wave
(3). The calculated outflow hydrographs show negligible distortion from the
initially sinusoidal shape, implying that the constant parameter assumption is
tantamount to an assumption of linearity.
The MuskingumCunge method with variable parameters (hydrographs E,
F, and G) shows results that fall within the range encompassed by the constant
parameter calculations. The noticeable steepening of the rising limb of the
calculated outflow hydrographs indicates that the nonlinearity of the phenomenon
is being taken into account. The threepoint and fourpoint methods give similar
results; however, the twopoint method shows a smaller peak and a somewhat
slower rate of travel. Furthermore, the twopoint method results in a significant
loss of mass, as indicated in Col. 8 of Table 1.
The results of the Thomas and fourpoint variable parameter methods are
comparable (hydrographs I and H). The threepoint variable parameter hydrograph
(not shown) is also very close to Thomas' results.
5. SUMMARY AND CONCLUSIONS
The MuskingumCunge method in which the parameters K and X are allowed
to vary in time and space is investigated. A threepoint approach and an iterative
fourpoint approach to the calculation of the variable parameters are shown
to be sufficiently accurate in the simulation of flood flows. A twopoint approach
is shown to be inaccurate in the calculation of peak discharge and travel time.
Furthermore, the twopoint method results in a significant loss of mass.
APPENDIX. REFERENCES
Chow, V. T. 1974. Handbook of Applied Hydrology, McGrawHill Book Co., Inc., New
York, NY.
Cunge, J. A. 1969. "On the Subject of a Flood Propagation Computation Method (Muskingum
Method), Journal of Hydraulic Research, Vol. 7, No. 2, 205230.
Dooge, J. C. I. 1973. "Linear Theory of Hydrologic Systems," Agricultural Research Service
Technical Bulletin No. 1468, Oct.
"Flood Studies Report. 1975. Vol. III: Flood Routing Studies, Natural Environment Research
Council, London, England.
Koussis, A. 1978. "Theoretical Estimation of Food Routing Parameters," Journal of the
Hydraulics Division, ASCE, Vol. 104, No. HY1, Proc. Paper 13456. Jan.,
109115.
Miller, W. A.. and J. A. Cunge. 1975. "Simplified Equations of Unsteady Flow," Unsteady
Flow in Open Channels, K. Mahmood and V. Yevjevich, eds., Water Resources
Publications, Fort Collins, Colo.
Roache, P. 1972.Computational Fluid Dynamics, Hermosa Publishers, Albuquerque, NM.
Thomas H. A. 1934. "The Hydraulics of Flood Movement in Rivers." Engineering Bulletin,
Carnegie Institute of Technology, Pittsburgh, PA.
