Muskingum-Cunge Model

Basic Concepts and Equations

The Muskingum-Cunge routing method builds upon concepts within the Muskingum method which was previously presented. This method uses a combination of the continuity equation and a simplified form of the momentum equation. The Muskingum-Cunge method is sometimes referred to as a “variable coefficient” method since the routing parameters are recalculated every time step based upon channel properties and the flow depth.

This method begins by neglecting the convective and local acceleration terms within the momentum equation which yields:

1)	$\begin{array}{l}\displaystyle S_{f}=S_{o}-\frac{\partial y}{\partial x}\end{array}$

Equation 1 is then combined with the following equation:

2)	$\begin{array}{l}\displaystyle \frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=q_{L}\end{array}$

After combining Equation 2 and Equation 3, a linear approximation can be used to yield the convective diffusion equation (Miller & Cunge, 1975):

3)	$\begin{array}{l}\displaystyle \frac{\partial Q}{\partial t}+c \frac{\partial Q}{\partial x}=\mu \frac{\partial^{2} Q}{\partial x^{2}}+c q_{L}\end{array}$

where c = flood wave celerity, μ = hydraulic diffusivity, and q_L = lateral inflow. Flood wave celerity and hydraulic diffusivity can be expressed as:

4)	$\begin{array}{l}\displaystyle c=\frac{d Q}{d A}\end{array}$

5)	$\begin{array}{l}\displaystyle \mu=\frac{Q}{2 B S_{o}}\end{array}$

where B = top width of the water surface.

Recall that HEC-HMS solves the the following equation recursively to compute O_t given inflow (I_t and I_t_-1), an initial condition (O_t₌₀), K, and X when using the Muskingum method:

6)	$\begin{array}{l}\displaystyle O_{t}=\left(\frac{\Delta t-2 K X}{2 K(1-X)+\Delta t}\right) I_{t}+\left(\frac{\Delta t+2 K X}{2 K(1-X)+\Delta t}\right) I_{t-1}+\left(\frac{2 K(1-X)-\Delta t}{2 K(1-X)+\Delta t}\right) O_{t-1}\end{array}$

Using a finite difference approximation of the partial derivatives in Equation 3 and combination with Equation 6 yields:

7)	$\begin{array}{l}\displaystyle O_{t}=C_{1} I_{t-1}+C_{2} I_{t}+C_{3} O_{t-1}+C_{4}\left(q_{L} \Delta x\right)\end{array}$

The C₁, C₂, C₃, and C₄ coefficients are:

8)	$\begin{array}{l}\displaystyle C_{1}=\frac{\frac{\Delta t}{K}+2 X}{\frac{\Delta t}{K}+2(1-X)}\end{array}$

9)	$\begin{array}{l}\displaystyle C_{2}=\frac{\frac{\Delta t}{K}-2 X}{\frac{\Delta t}{K}+2(1-X)}\end{array}$

10)	$\begin{array}{l}\displaystyle C_{3}=\frac{2(1-X)-\frac{\Delta t}{K}}{\frac{\Delta t}{K}+2(1-X)}\end{array}$

11)	$\begin{array}{l}\displaystyle C_{4}=\frac{2\left(\frac{\Delta t}{K}\right)}{\frac{\Delta t}{K}+2(1-X)}\end{array}$

Within the previously mentioned Muskingum method, the X parameter is a dimensionless coefficient that lacks a strong physical meaning. Cunge (1969) evaluated the numerical diffusion that is produced through the use of Equation 6 and set this equal to the physical diffusion represented by Equation 3. This yielded the following representations for K and X (Ponce & Yevjevich, 1978):

12)	$\begin{array}{l}\displaystyle K=\frac{\Delta x}{c}\end{array}$

13)	$\begin{array}{l}\displaystyle X=\frac{1}{2}\left(1-\frac{Q}{B S_{o} c \Delta x}\right)\end{array}$

Since c, Q, and B can change during the passage of a flood wave, the coefficients C₁, C₂, C₃, and C₄ also change. As such, the C₁, C₂, C₃, and C₄ coefficients are recomputed each time and distance step (Δt and Δx) using the algorithm proposed by Ponce (1986).

Required Parameters

The parameters that are required to utilize this method within HEC-HMS are the initial condition, the reach length [ft or m], the friction slope [ft/ft or m/m], Manning’s n roughness coefficient, a space-time interval method and value, an index method and value, and a cross-section shape and parameters/dimensions. An optional invert can also be specified.

Two options for specifying the initial condition are included: outflow equals inflow and specified discharge [ft³/sec or m³/sec]. The first option assumes that the initial outflow is the same as the initial inflow to the reach from the upstream elements which is equivalent to the assumption of a steady-state initial condition. The second option is most appropriate when there is observed streamflow data at the end of the reach.

The reach length should be set as the total length of the reach element while the friction slope should be set as the average friction slope for the entire reach. If the friction slope varies significantly throughout the stream represented by the reach, it may be necessary to use multiple reaches with different slopes. If no information is available to estimate the friction slope, the bed slope can be used as an approximation. The Manning's n roughness coefficient should be set as the average value for the whole reach. This value can be estimated using “reference” streams with established roughness coefficients or through calibration.

The choices of space and time steps (Δx and Δt) are critical to ensure accuracy and stability. Three options for specifying a space-time method are provided within HEC-HMS: 1) Auto DX Auto DT, 2) Specified DX Auto DT, and 3) Specified DX Specified DT. When the Auto DX Auto DT method is selected, space and time intervals that attempt to maintain numerical stability will automatically be selected. Δt is selected as the minimum of either the user-specified time step or the travel time through the reach (rounded to the nearest multiple or divisor of the user-specified time step). Once Δt is computed, Δx is computed as:

14)	$\begin{array}{l}\displaystyle \Delta x = c\Delta t\end{array}$

When the Specified DX Auto DT method is selected, the specified number of subreaches will be used while automatically varying the time interval to take as long a time interval as possible while also maintaining numerical stability. When the Specified DX Specified DT method is selected, the specified number of subreaches and subintervals (rounded to the nearest multiple or divisor of the user-specified time step) will be used throughout the entire simulation.

Upon completion of a simulation, the minimum and maximum celerity of the routed hydrograph will be displayed as notes. Also, a reference space step, Δx_ref, will be computed using methodology presented in Engineer Manual 1110-2-1417 Flood-Runoff Analysis (U.S. Army Corps of Engineers, 1994):

15)	$\begin{array}{l}\displaystyle \Delta x=\frac{1}{2}\left(c \Delta t+\frac{Q_{o}}{B S_{o} c}\right)\end{array}$

where Q₀ = reference flow, which is computed from the inflow hydrograph as:

16)	$\begin{array}{l}\displaystyle Q_{o}=Q_{B}+\frac{1}{2}\left(Q_{p c a k}-Q_{B}\right)\end{array}$

where Q_B = baseflow and Q_peak = inflow peak. If Δx_ref is not equal to the Δx that was used during the routing computations, a note will be issued.

The index method is used in conjunction with the physical properties of the channel and the previously mentioned Δt and Δx interval selection to discretize the routing reach in both space and time. Appropriate reference flows and flood wave celerities are dependent upon the physical properties of the channel as well as the flood event(s) in question. Experience has shown that a reference flow (or celerity) based upon average values of the hydrograph in question (i.e. midway between the base flow and the peak flow) is, in general, the most suitable choice. Reference flows (or celerities) based on peak values tend to numerically accelerate the flood wave more than would occur in nature while the converse is true if a low reference flow (or celerity) is used (Ponce, 1983).

Six options are provided for specifying the cross-section shape: circle, eight-point, rectangle, tabular, trapezoid, and triangle. The circle shape is not meant to be used for pressure flow or pipe networks, but is suitable for representing a free surface inside a pipe. Depending upon the chosen shape, additional information will have to be entered to describe the size of the cross-section shape. This information may include a diameter (circle) [ft or m], bottom width (deep, rectangle, and/or trapezoid) [ft or m], or side slope (trapezoid and triangle) [ft/ft or m/m]. In all cases, cross-section shapes must be defined in such a way that all possible flow depths that will be simulated will be completely confined within the defined shape.

The tabular shape option allows for the use of user-defined elevation vs. discharge, elevation vs. area, and elevation vs. width relationships (which are all paired data objects). This option is typically used when relationships derived from hydraulic simulations are available. When the tabular shape is selected, no Manning's n roughness coefficients need to be entered.

When using the eight-point shape, a simplified cross-section (which is a paired data object) with eight station vs. elevation values must be selected. The cross-section is typically configured to represent the main channel plus left and right overbank areas. As such, separate Manning's n values are required for each overbank.

Many of the aforementioned parameters are typically estimated using GIS. However, field survey data may be necessary to accurately determine reach lengths, friction or bed slopes, and/or cross-section shape parameters.

A tutorial describing an example application of this loss method, including parameter estimation and calibration, can be found here: Applying the Muskingum-Cunge Routing Method.