Return to Applying the Muskingum Routing Method

Overview

In this tutorial you will apply the HEC-HMS Muskingum-Cunge channel routing method to a modeling application.  Initial parameter estimates will be estimated using GIS and observed information and the model will be calibrated through trial and error.

HEC-HMS version 4.12-beta.2 was used to create this tutorial.  You will need to use HEC-HMS version 4.12-beta.2, or newer, to open the project files.

Download the initial project files here:

MuskingumCunge_v4.12_start.zip

Background

The Muskingum-Cunge routing method builds upon concepts within the Muskingum method.  This method uses a combination of the continuity equation and a simplified form of the momentum equation.  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 the Muskingum routing equation and set this equal to the physical diffusion represented by the convective diffusion equation (Miller & Cunge, 1975).  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.

Parameters that are required to utilize this method within HEC-HMS include 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.

Estimate Initial Parameter Values

Reach Length

The reach length should be set as the total length of the reach element.  This parameter is typically estimated using GIS information and field survey data, when available.  The total length of the reach in question is 68861 ft and is shown in Figure 1.

Figure 1.  Mahoning Creek Reach of Interest

Friction Slope

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.  This parameter is typically estimated using GIS information and field survey data, when available.  In this tutorial, invert elevations were extracted from five locations along the reach of interest and are shown in Figure 2.  The profile data for this information is detailed within Table 1.

Figure 2. Mahoning Creek Stream Profile

Table 1. Profile Data for Mahoning Creek

LabelLocationDistance above Punx Gage (ft)Elevation (ft)
1Conf. with Stump Creek688611282.8
2Robertsville457381259.8
3Cloe354861240.2
4Riker277261220.5
5Punxsutawney01200.8

Use Figure 1 and Table 2 to answer the following questions.


Question 1: Estimate a single, representative slope (in ft/ft) for the reach in question.

Slope = \frac{Rise}{Run}

Slope = \frac{1282.8 ft - 1200.8 ft}{68861 ft - 0 ft}

Slope = \frac{82 ft}{68861 ft}

Slope = 0.0012 \frac{ft}{ft}


Question 2: Divide the reach into three sub-reaches and estimate a slope (in ft/ft) for each.

Subreach 1: Point 1 (Conf. with Stump Creek) to Point 2 (Robertsville): slope = 0.001 ft/ft

Subreach 2: Point 2 (Robertsville) to Point 4 (Riker): slope = 0.0022 ft/ft

Subreach 3: Point 4 (Riker) to Point 5 (Punxsutawney): slope = 0.0007 ft/ft


Manning's n Roughness Coefficient

The Manning's n roughness coefficient(s) should be set as the average value for the whole reach.  Depending upon the cross section shape that is chosen, zero , one, or three Manning's n roughness coefficients are required.  This parameter is typically estimated using “reference” streams with established roughness coefficients or through model calibration.

Space-Time Interval

The choice of space and time steps (Δx and Δt, respectively) are critical to ensure accuracy and stability.  Three options for specifying a space-time method are provided within HEC-HMS:

  • Auto DX Auto DT
  • Specified DX Auto DT
  • Specified DX Specified DT

When the Auto DX Auto DT method is selected, space and time intervals will automatically be selected that attempt to maintain numerical stability.  Δ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:

\Delta x=c \Delta t

where c = flood wave celerity.  Upon completion of a simulation, the minimum and maximum celerity of the routed hydrograph will be displayed as notes.  Also, a reference space step, Δxref, will be computed using methodology presented in Engineer Manual 1110-2-1417 Flood-Runoff Analysis (U.S. Army Corps of Engineers, 1994).  If Δxref is not equal to the Δx that was used during the routing computations, a note will be issued.  For most applications, the Auto DX Auto DT method is adequate.

Index Method

The index method is used in conjunction with the physical properties of the channel and the previously mentioned Δx and Δt interval selection to discretize the routing reach in both space and time.  Two options for specifying the index method are provided within HEC-HMS: 

  • Celerity
  • Flow

When the Celerity index method is used, a reference celerity [ft/s or m/s] must be specified.  When the Flow index method is used, a reference flow [cfs or cms] must be specified and cross-section properties are then used to infer a celerity.  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).

For most applications, using the Celerity index method and a celerity of 5 ft/s is adequate.

Cross Section

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.

Cross section shapes and properties are typically estimated using GIS and field survey data, when available.  In this tutorial, cross sections were extracted at three locations along the reach of interest, as shown in Figure 3.  The station-elevation data for these cross sections are detailed within Figure 3, Figure 4, and Figure 5 in addition to Table 2, Table 3, and Table 4.

Figure 3.  Cross Section Locations

Figure 3. Cross Section near Stump Creek Confluence

Table 2. Cross Section near Stump Creek Confluence

Station (ft)Elevation (ft)
035.1
7.526.2
14.119.7
2015.1
39.415.1
61.413.1
75.514.4
106.614.8
120.115.1
124.37.9
128.94.3
134.81
1440.7
160.81.3
1890.3
2100
213.61
219.26.6
223.113.1
224.115.1
232.916.1
242.815.1
255.913.8
305.115.1
319.915.4
321.919.7
324.826.2
339.935.1

Figure 4. Cross Section at Cloe

Table 3. Cross Section at Cloe

Station (ft)Elevation (ft)
035.1
7.526.2
14.119.7
2015.1
39.414.4
75.813.1
121.414.4
156.513.8
220.114.1
227.79.8
229.74.3
234.91
255.90.7
2690
285.81.3
297.61
3152
321.56.6
328.113.1
340.615.1
360.914.8
393.714.1
426.516.1
47914.4
499.316.7
50221.7
541.326.2
549.935.1

Figure 5. Cross Section at Punxsutawney

Table 4. Cross Section at Punxsutawney

Station (ft)Elevation (ft)
035.1
7.526.2
14.119.7
2015.4
46.315.4
69.613.1
88.314.1
106.616.1
122.715.1
124.37.9
128.94.3
127.61
134.80
158.11
178.81.6
201.42
2113.6
219.26.6
223.113.1
225.115.1
244.815.4
255.915.1
273.314.4
305.115.1
315.916.4
321.919.7
338.926.2
345.135.1


A tutorial containing steps that can be taken to enter these curves within HEC-HMS can be found here: Creating Cross Section Paired Data Curves.

Modify the Existing HEC-HMS Project

Now that you've estimated initial parameters, begin modifying the existing HEC-HMS project.

  1. Open the Punxsutawney project.
  2. Note that there are two basin models within this project: 1) Mahoning Creek and 2) Mahoning Creek Multi Reaches.
    1. The Mahoning Creek basin model contains a single routing reach while the Mahoning Creek Multi Reaches basin model contains three reaches.
    2. You will make use of both basin models within this tutorial.
  3. Expand the Paired Data node within the tree and take note of the three 8-point cross sections.
    1. These 8-point cross sections were simplified from the previously shown cross sections near Stump Creek, Cloe, and Punxsutawney

      • Open the Paired Data Manager from the Components menu.
      • Select the Cross Sections paired data type and click the New… button.
      • Enter a Name for the paired data curve (e.g. “upper”) and click the Create button.
      • Navigate to the new paired data curve in the Watershed Explorer.
      • The paired data curve is located in the Cross Sections folder within the Paired Data folder.
      • Click on the newly-created curve to open the Component Editor.
      • Within the Component Editor, set the Data Source to Manual Entry and the units to FT:FT.
      • Enter the eight previously-determined station-elevation pairs from into the table provided on the Table tab.

       

  4. Open the Mahoning Creek Basin Model.
  5. Select the Mahoning Creek reach element.
    1. Change the routing method to Muskingum-Cunge.
    2. Select the Routing tab to open the Muskingum-Cunge Component Editor.
    3. Make sure the Initial Type is set to Discharge = Inflow.
    4. Enter the previously stated reach length of 68861 ft.
    5. Enter the friction slope which was estimated within Question 1.
    6. Using one of the previously shown cross sections, enter a Manning’s n roughness coefficient for the main channel.
    7. Select the Auto DX Auto DT space-time interval method.
    8. Select the Celerity index method and enter an index celerity of 5 ft/s.
    9. For the same cross section, enter the left and right Manning’s n roughness coefficients.
    10. Finally, select the cross section within the Cross Section drop down menu.
    11. The Muskingum-Cunge Component Editor should resemble Figure 6.
      Figure 6.  Muskingum-Cunge Routing Component Editor
  6. Create a new simulation run by clicking Compute | Create Compute | Simulation Run.
    1. Name the new simulation run "April1994_MuskCunge".
    2. Select the Mahoning Creek basin model.
    3. Select the NoRain meteorologic model.
    4. Select the April 1994 control specifications.
  7. Select the April1994_MuskCunge simulation run from the Compute toolbar, .
  8. Press the Compute All Elements button,  , to run the simulation.
  9. View the result graph, , and summary table, , for the Mahoning Creek reach element, as shown in Figure 7.
    1. Notice the computed peak discharge, time to peak discharge, and discharge hydrograph shape match the observed data reasonably well when using a single routing reach and the initial parameter estimates from GIS and field surveys.  This is confirmed within the summary table.
      Figure 7.  Initial Results
    2. However, an improved estimate can be made by subdividing the routing reach into three separate reaches.  This has already been done for you within the Mahoning Creek Multi Reaches basin model.  Further modifications will be made using that basin model.

Calibrate the Model

Calibrating a model to afford better agreement between computed and observed runoff oftentimes requires simultaneous changes to more than just one process (e.g. calibrate both transform and routing reach parameters at the same time).  However, within this tutorial, only the Muskingum-Cunge routing method parameters will be modified.

  1. Open the Mahoning Creek Multi Reaches Basin Model.
  2. Select the Upper reach element.
    1. Change the routing method to Muskingum-Cunge.
    2. Select the Routing tab to open the Muskingum-Cunge Component Editor.
    3. Make sure the Initial Type is set to Discharge = Inflow.
    4. Use Figure 2 and Table 1 to estimate an appropriate reach length and friction slope.  Enter these values within the Length and Slope entry fields.
    5. Use the near Stump Creek cross section and enter the corresponding Manning’s n roughness coefficient for the main channel.
    6. Select the Auto DX Auto DT space-time interval method.
    7. Select the Celerity index method and enter an index celerity of 5 ft/s.
    8. For the near Stump Creek cross section, enter the left and right Manning’s n roughness coefficients.
    9. Finally, select the near Stump Creek cross section within the Cross Section drop down menu.
    10. The Muskingum-Cunge Component Editor for the Upper reach element should resemble Figure 8.
      Figure 8. Upper Reach Component Editor
  3. Select the Middle reach element.
    1. Change the routing method to Muskingum-Cunge.
    2. Select the Routing tab to open the Muskingum-Cunge Component Editor.
    3. Make sure the Initial Type is set to Discharge = Inflow.
    4. Use Figure 2 and Table 1 to estimate an appropriate reach length and friction slope.  Enter these values within the Length and Slope entry fields.
    5. Use the Cloe cross section and enter the corresponding Manning’s n roughness coefficient for the main channel.
    6. Select the Auto DX Auto DT space-time interval method.
    7. Select the Celerity index method and enter an index celerity of 5 ft/s.
    8. For the Cloe cross section, enter the left and right Manning’s n roughness coefficients.
    9. Finally, select the Cloe cross section within the Cross Section drop down menu.
    10. The Muskingum-Cunge Component Editor for the Middle reach element should resemble Figure 9.
      Figure 9. Middle Reach Component Editor
  4. Select the Lower reach element.
    1. Change the routing method to Muskingum-Cunge.
    2. Select the Routing tab to open the Muskingum-Cunge Component Editor.
    3. Make sure the Initial Type is set to Discharge = Inflow.
    4. Use Figure 2 and Table 1 to estimate an appropriate reach length and friction slope.  Enter these values within the Length and Slope entry fields.
    5. Use the Punxsutawney cross section and enter the corresponding Manning’s n roughness coefficient for the main channel.
    6. Select the Auto DX Auto DT space-time interval method.
    7. Select the Celerity index method and enter an index celerity of 5 ft/s.
    8. For the Punxsutawney cross section, enter the left and right Manning’s n roughness coefficients.
    9. Finally, select the Punxsutawney cross section within the Cross Section drop down menu.
    10. The Muskingum-Cunge Component Editor for the Lower reach element should resemble Figure 10.
      Figure 10. Lower Reach Component Editor
  5. Create a new simulation run by clicking Compute | Create Compute | Simulation Run.
    1. Name the new simulation run "April1994_MuskCunge_multi".
    2. Select the Mahoning Creek Multi Reaches basin model.
    3. Select the NoRain meteorologic model.
    4. Select the April 1994 control specifications.
  6. Select the April1994_MuskCunge_multi simulation run from the Compute toolbar, .
  7. Press the Compute All Elements button,  , to run the simulation.
  8. View the result graph and summary table for the Mahoning Creek Gage at Punx sink element, as shown in Figure 11.
    1. Notice the computed peak discharge, time to peak discharge, and discharge hydrograph shape match the observed data better than when using a single routing reach.  This is confirmed within the summary table.
      Figure 11. April1994_MuskCunge_multi Results
  9. Continue adjusting the Muskingum-Cunge routing method parameters in an attempt to simultaneously best match the peak discharge flow rate, time of peak discharge, hydrograph shape, and discharge volume.  Record the various differences between the computed and observed discharge hydrographs in addition to the RMSE Std Dev, Percent Bias, and Nash-Sutcliffe statistical metrics as you adjust model parameters.

Download the final project files here:

MuskingumCunge_v4.12_finish.zip

Summary Questions

Question 3. When using the Mahoning Creek Multi Reaches basin model, are the differences between the computed and observed hydrographs significant?  To what physical processes could any differences be attributed?  What could be changed to provide better agreement between the Muskingum-Cunge results and the results obtained from HEC-RAS?

The results from HEC-RAS and Muskingum-Cunge are very similar when the same channel data are used for both models.  This is to be expected.  If the results were not similar, we might expect that the abrupt change in slope at point 4 might be causing a backwater effect to occur.  It can be noted from the data that the slopes are large enough to minimize errors from simplifying the full unsteady flow equations (i.e. Slope > 0.0006).  Only the full unsteady flow equations are robust enough to predict depths and/or velocities when slopes are less than 2 ft/mi or if extensive backwater effects are present.  If channel obstructions such as bridges and culverts are present, backwater effects may exist which Muskingum-Cunge cannot explicitly recreate.


Question 4. Do you consider the differences in the routed hydrographs to be significant?  How would you quantify the significance from a study perspective?

There is an observable difference between the results obtained through the use of one reach when compared against to those obtained through the use of three reaches.  However, the significance of the differences are dependent upon the study objective and the study's sensitivity.  In some cases, such as one that analyzes possible levee configurations, the significance could be evaluated by examining the corresponding differences in peak stage given a fixed stage-flow rating curve.  In other cases, such as storage facility analysis, the differences in discharge volume are more critical to the study's success.


Question 5. If the answer to the previous question is yes, you might be inclined to subdivide the larger routing reach into smaller subreaches.  How would you know when to stop subdividing the reach in subreaches?

You would only need to subdivide until the simulated hydrograph stopped changing significantly.  Again, the sensitivity of the particular study objective will determine the point at which the changes in results cease to be important.