Analytical verification
of Muskingum - Cunge routing

Victor M. Ponce,
A. K. Lohani, and C. Scheyhing

Online version 2015

[Original version 1996]


A verification of Muskingum-Cunge routing is accomplished by comparing theoretically calculated peak outflow and travel time with those generated using the constant-parameter Muskingum-Cunge method of flood rooting. The close agreement between analytical and numerical results underscores the utility of Muskingum-Cunge routing as a viable and accurate method for routine applications in flood hydrology.


The well-known Muskingum-Cunge method of flood routing (Natural Environment Research Council 1975; Ponce and Yevjevich, 1978) has for many years been successfully applied in flood modeling and prediction. The Muskingum routing parameters were theoretically derived by Cunge (1969) based on the approximation error obtained by a Taylor series expansion of the grid function. Although the parameters so derived have been shown to work reasonably well in a variety of practical applications (Hydrologic Engineering Center, US Army Corps of Engineers, 1990), their link to the analytical expressions for wave celerity (Seddon, 1900) and hydraulic diffusivity (Hayami, 1951) has yet to be demonstrated in an actual routing application.

This paper endeavors to test the Muskingum-Cunge model by comparing the analytically derived flood wave celerity and peak attenuation amounts with those obtained using an actual numerical application of the Muskingum-Cunge model. Agreement between analytical and numerical results is taken as an indication of the potential accuracy of the model, when used to route actual floods in a practical setting.


The Muskingum-Cunge model of flood routing (Cunge, 1969; Natural Environment Research Council, 1975) is an improvement of the classical Muskingum model (Chow, 1959), in which the routing parameters are based on readily measurable hydraulic data (stage-discharge rating, channel slope, channel width, and reach length), instead of being based on historic flood discharge data. This allows routing in ungaged streams with a reasonable expectation of accuracy, which would otherwise be possible only if the streams were gaged. As it is impractical to gage more than a few streams, a method that is both accurate and does not rely on stream gaging has to be acknowledged.

Ponce and Yevjevich (1978) improved the Muskingum-Cunge model by expressing the routing parameters as the dimensionless Courant and cell Reynolds numbers, where the former is the ratio of physical and numerical celerities (Seddon, 1900; O'Brien et al., 1950) and the latter is the ratio of physical and numerical diffusivifies (Hayami, 1951; Cunge, 1969).

       c Δt
C = _____


D = __________
         S0 c Δx


in which C is Courant number, D is cell Reynolds number, c is flood wave celerity, Δx is space step, Δt is time step, q0 is reference unit-width discharge, and S0 is bed slope.

With the Muskingum routing parameters defined in this way, the routing coefficients are (Ponce and Yevejevich, 1978):

        -1 + C + D
C0 = ___________
         1 + C + D


         1 + C - D
C1 = ___________
         1 + C + D


         1 - C + D
C2 = ___________
         1 + C + D


and the routing equation is:

Qj+1n+1 = C0 Qj n+1 + C1 Qj n + C2 Qj+1n


where j is a spatial index and n is a temporal index.

Cunge's formulation (Cunge, 1969) is based on diffusion wave theory, the applicability of which was confirmed by Ponce et al. (1978). Pioneering analytical work on the diffusion wave is due to Hayami (1951), who reasoned that the large-scale longitudinal mixing caused by channel irregularities could be effectively modeled as a diffusion process. Ponce and Simons (1977) used the theory of linear stability to derive celerity and attenuation functions for sinusoidal perturbations to the steady uniform flow.

In this paper we use the findings of Ponce and Simons (1977) to calculate theoretical (analytical) peak outflow and travel time for a series of sinusoidal hydrographs, and to compare these results with those obtained using the (numerical) Muskingum-Cunge model (Eqs. 1-6). The use and limitations of the Muskingum-Cunge model have been given by Ponce (1989).


The analytical diffusion wave celerity can be approximated as the Seddon speed (Ponce and Simons, 1977):

c =bu


in which c is flood wave celerity, u is mean velocity, and b is rating exponent, in

Q = αAb



q = α d b


in which Q is discharge, A is flow area, q is unit-width discharge, and d is flow depth.

The wave attenuation follows an exponential decay function of the following form:

qpo = q0 + (qpi - q0) exp(-βl*t *)


in which qpo is peak outflow, qpi is peak inflow, q0 is reference flow, βl* is dimensionless amplitude propagation factor, and t * dimensionless elapsed time.

The amplitude propagation factor is (Ponce and Simons, 1977)

β = __ σ * 2


in which σ * is the dimensionless wavenumber, defined as:

σ * = (2π / L) L0


in which L is sinusoidal wavelength, and L0 is reference channel length, defined as L0 = d0 / S0.

As L = cT, where T is wave period, taken as the time base Tb of a sinusoidal hydrograph, and t * = t (u0/L0), it follows that:

βl* t * = (______)2 νt


in which ν = q0 / (2S0) is the hydraulic diffusivity (Hayami, 1951).


We used an extension of Thomas' (1934) classical flood routing problem to test the agreement between analytically and numerically derived peak outflow and travel time. The Thomas problem routes a sinusoidal flood hydrograph of time base (i.e. wave period) Tb = 96 h, peak inflow qpi = 200 cfs-1 (cubic feet per second per foot), and basellow qb = 50 cfs ft-1 through a channel of length Lr = 200 miles, slope S0 = 1 ft mile-1, and rating curve

q = 0.688d 5/3


from which the Manning coefficent n = 0.0297 is calculated.

A series of test cases was developed based on Thomas' problem, using: (1) two wave periods Tb: 48 and 96 h; (2) three peak inflows qpi: 100, 200 and 500 cfs ft-1 , i.e. three peak inflow-to-base flow ratios qpi / qb: two, four and ten; (3) two channel lengths Lr: 200 and 500 miles. Twelve outflow hydrographs were generated using the constant-parameter Muskingum-Cunge model (Ponce and Yevjevich, 1978), where the reference flow is taken as:

                  qpi - q0
q0 = qb + __________


For Tb = 48 h, the chosen space and time intervals were Δx = 6 .25 miles, and Δt = 1.5 h, respectively; for Tb = 96h, the intervals were Δx = 12.5 miles, and Δt = 3 h. The chosen wave periods (48 and 96 h) satisfy diffusion wave applicability criteria (Ponce et al., 1978):

Tb S0 (g / d0)1/2 ≥ 30


where g is gravitational acceleration. In effect, for the less favorable case (Tb, = 48 h and qpi = 500 cfs ft-1), the criterion is:

Tb S0 (g / d0)1/2 = 31


which guarantees that the routed waves remain within the realm of diffusion waves.


Table 1 summarizes the comparison of analytical and numerical peak outflow qpo and travel time Tt, for the chosen series of sinusoidal hydrographs. Figure 1 shows inflow and oufflow hydrographs for a typical case (Tb = 96 h, Lr = 500 miles, and qpi /qb = 4). In this case, reference flow is q0 = 125 cfs ft-1, peak outflow qpo = 176.70 cfs ft-1, and travel time Tt = 79.8 h.

Table 1. Comparison between analytical and numerical peak outflow qpo (cfs ft-1) and travel time Tt (h)
for a series of sinusoidal hydrographs.

Fig 1. Inflow and outflow hydrographs for a typical case (Tb = 96 h, Lr = 500 miles, and qpi / qb = 4).

On the basis of the foregoing analysis, we conclude that Muskingum-Cunge routing accurately simulates flood wave propagation. In twelve (12) tests encompassing a wide range of conditions likely to be encountered in practice, the ratio of numerical to analytical peak outflows varied within the range 0.991-1.003, and the ratio of numerical to analytical travel times varied within the range 0.987-1.021. The agreement in indeed remarkable, considering that a constant reference flow was used in the computation of the analytical values.


The present study was performed in autumn 1994, while the second author was at San Diego State University, on leave from the Ganga Plains Regional Centre, National Institute of Hydrology, Patna (Bihar), India. His leave was funded by the United Nations Development Programme.


Chow, V.T., 1950. 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.

Hayami,S., 1951.On the Propagation of flood waves. Disaster Prey. Res. Inst, Kyoto Univ. Bull., 1: 1-16 Hydrologic Engineering Center, US Army Corps of Engineers, HEC-1, Flood hydrograph package. user's manual. Version 4.0, September, Davis, CA.

Natural Environment Research Council. 1975. Flood Studies Report, Vol. 5: Flood Routing Studies. NERC, London.

O'Brien,G.G., Hyman, M.A. and Kaplan, S., 1950. A study of the numerical solution of partial differential equations. J. Math. Phys., 29(4): 231-251.

Ponce, V.M., 1989. Engineering Hydrology, Principles and Practices. Prentice-Hall, Englewood Cliffs, NJ.

Ponce, V.M. and Simons,D.B, 1977. Shallow wave propagation in open channel flow. J. Hydraul. Div. ASCE. 103(.12): 1461-1476.

Ponce, V.M. and Yevjevich, V., 1978. Muskingum -Cunge method with variable parameter. J. Hydraul. Div. ASCE, 104(HY312): 1663-360.

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

Seddon, J.44,19110. River hydraulics. Trans. Am. Soc. Civ. Eng., 431 179-229.

Thomas, H.A., 1934. The hydraulics of flood movement in rivers. Engineering Bulletin. Carnegie Institute of Technology. Pittsburgh, PA.


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