Go to ⇒ Presentation format 
9.1 MUSKINGUM METHOD
Stream Channel Routing Stream channel routing uses mathematical relations to calculate outflow from a stream channel once inflow, lateral contributions, and channel characteristics are known. Stream channel routing usually implies open channel flow conditions, although there are exceptions, such as storm sewer flow, for which mixed open channelclosed conduit flow conditions may prevail. In this chapter, stream channel routing refers to unsteady flow calculations in streams and rivers. Channel reach refers to a specific length of stream channel possessing certain translation and storage properties. The hydrograph at the upstream end of the reach is the inflow hydrograph; the hydrograph at the downstream end is the outflow hydrograph. Lateral contributions consist of point tributary inflows and/or distributed inflows; e.g., interflow and groundwater flow. The terms stream channel routing and flood routing are often used interchangeably. This is attributed to the fact that most stream channel routing applications are in flood flow analysis, flood control design, or flood forecasting (Fig. 91).
Two general approaches to stream channel routing are recognized: (1) hydrologic and (2) hydraulic. As in the case of reservoir routing (Chapter 8), hydrologic stream channel routing is based on the storage concept. Conversely, hydraulic channel routing is based on the principles of mass and momentum conservation. Hydraulic routing techniques are of three types: (1) kinematic wave, (2) diffusion wave, and (3) dynamic wave. The dynamic wave is the most complete model of unsteady open channel flow. Kinematic and diffusion waves are convenient and practical approximations to the dynamic wave. An alternate approach to hydrologic and hydraulic routing has emerged in recent years. This approach is similar in nature to the hydrologic routing methods yet contains sufficient physical information to compare favorably with the more complex hydraulic routing techniques. This hybrid approach is the basis of the MuskingumCunge method of flood routing. At the outset of the study of stream channel routing, it is necessary to introduce a few basic modeling concepts. A typical hydrologic model consists of: (1) input, (2) system, and (3) output (Fig. 92). In surface water hydrology, the system is usually a catchment, a reservoir, or a stream channel. In the case of a catchment, the input is a storm hyetograph. For reservoirs and stream channels, the input is an inflow hydrograph. For all three cases, catchments, reservoirs, and channels, the output is an outflow hydrograph.
In general, modeling problems are classified into three types: (1) prediction, (2) calibration, and (3) inversion. In the prediction problem, input and system are known and described by properties or parameters, and the task is to calculate the output based on the knowledge of system and input. For instance, with known inflow hydrograph, lateral contributions, and channel reach parameters, the outflow hydrograph from a stream channel can be computed using routing techniques (Example 91). In the calibration problem, input and output are known, and the objective is to determine the properties or parameters describing the system. In the case of a stream channel, with known upstream inflow, lateral contributions, and outflow hydrograph, the routing parameters are calculated by a calibration procedure (Example 92). The inversion problem is the third type of modeling problem. In this case, system and output are known, and the task is to calculate the inflow or inflows. This is accomplished by reversing the routing process in a technique known as inverse channel routing. For instance, with known upstream inflow, outflow, and channel reach parameters, the lateral contributions can be calculated by inverse routing. The prediction problem is the more common type of modeling application; however, a calibration is usually required in advance of the prediction. Model verification is the process of testing the model with actual data to establish its predictive accuracy. To calibrate and verify a model, it is usually necessary to assemble two different data sets. The first set is used in model calibration and the second set is used in model verification. A close agreement between calculated and measured data is an indication that the model has been verified. A detailed discussion of these subjects is given in Chapter 13. Muskingum Method The Muskingum method of flood routing was developed in the 1930s in connection with the design of flood protection schemes in the Muskingum River Basin, Ohio (Fig. 93) [11]. It is the most widely used method of hydrologic stream channel routing, with numerous applications in the United States and throughout the world.
The Muskingum method is based on the differential equation of storage, Eq. 84, reproduced here:
In an ideal channel, storage is a function of inflow and outflow. This is in constrast with an ideal reservoir, in which storage is solely a function of outflow (see Eqs. 85 to 87). In the Muskingum method, storage is a linear function of inflow and outflow:
in which S = storage volume; I = inflow; O = outflow; K = a time constant or storage coefficient; and X = a dimensionless weighting factor. With inflow and outflow in cubic meters per second, and K in hours, storage volume is in (cubic meters per second)hour. Alternatively, K could be expressed in seconds, in which case storage volume is in cubic meters.
Equation 91 was developed in 1938 and has been widely used since then [11].
It is esentially a generalization of the linear reservoir concept (Eq. 87).
In fact, for X = 0, Eq. 91 reduces to Eq. 87.
In other words, linear reservoir routing is a special case of Muskingum channel routing for which To derive the Muskingum routing equation, Eq. 84 is discretized on the xt plane (Fig. 82), to yield Eq. 813, repeated here:
Equation 91 is expressed at time levels 1 and 2:
Substituting Eqs. 92 to 93 into Eq. 813 and solving for O_{2} yields Eq. 815, repeated here:
in which C_{0}, C_{1} and C_{2} are routing coefficients defined in terms of Δt, K, and X as follows:
Since C_{0} + C_{1} + C_{2} = 1, the routing coefficients may be interpreted as weighting coefficients.
For Given an inflow hydrograph, an initial flow condition, a chosen time interval Δt, and routing parameters X and K, the routing coefficients can be calculated with Eqs. 94 to 96, and the outflow hydrograph with Eq. 815. The routing parameters K and K are related to flow and channel characteristics, K being interpreted as the travel time of the flood wave from upstream end to downstream end of the channel reach. Therefore, K accounts for the translation (or concentration) portion of the routing (Fig. 93). The parameter X accounts for the storage portion of the routing. For a given flood event, there is a value of X for which the storage in the calculated outflow hydrograph matches that of the measured outflow hydrograph. The effect of storage is to reduce the peak flow and spread the hydrograph in time (Fig. 94). Therefore, it is often used interchangeably with the terms diffusion and peak attenuation.
The routing parameter K is a function of channel reach length and flood wave speed; conversely, the parameter X is a function of the flow and channel characteristics that cause runoff diffusion. In the Muskingum method, X is interpreted as a weighting factor and restricted in the range 0.0 ≤ X ≤ 0.5. Values of X greater than 0.5 produce hydrograph amplification (i.e., negative diffusion), which does not correspond with reality (under the Froude numbers applicable to flood flows). With K = Δt and X = 0.5, flow conditions are such that the outflow hydrograph retains the same shape as the inflow hydrograph, but it is translated downstream a time equal to K. For X = 0, Muskingum routing reduces to linear reservoir routing (Section 8.2). In the Muskingum method, the parameters K and X are determined by calibration using streamflow records. Simultaneous inflowoutflow discharge measurements for a given channel reach are coupled with a trialanderror procedure, leading to the determination of K and X (see Example 92). The procedure is timeconsuming and lacks predictive capability. Values of K and X determined in this way are valid only for the given reach and flood event used in the calibration. Extrapolation to other reaches or to other flood events (of different magnitude) within the same reach is usually unwarranted. When sufficient data are available, a calibration can be performed for several flood events, each of different magnitude, to cover a wide range of flood levels. In this way, the variation of K and X as a function of flood level can be ascertained. In practice, K is more sensitive to flood level than X. A sketch of the variation of K with stage and discharge is shown in Fig. 95.
Unlike reservoir routing, stream channelrouting calculations exhibit a definite (time) lag between inflow and outflow. Furthermore, in the general case (X ≠ 0), maximum outflow does not occur at the time when inflow and outflow coincide. Example 91 has illustrated the predictive stage of the Muskingum method, in which the routing parameters are known in advance of the routing. If the parameters are not known, it is first necessary to perform a calibration. The trialanderror procedure to calibrate the routing parameters is illustrated by Example 92.
The estimation of routing parameters is crucial to the application of the Muskingum method. The parameters are not constant, tending to vary with flow rate. If the routing parameters can be related to flow and channel characteristics, the need for trialanderror calibration would be eliminated. Parameter K could be related to reach length and flood wave velocity, whereas X could be related to the diffusivity characteristics of flow and channel. These propositions are the basis of the MuskingumCunge method (Section 9.4).
9.2 KINEMATIC WAVES
Three types of unsteady openchannel flow waves are commonly used in engineering
hydrology: Kinematic Wave Equation The derivation of the kinematic wave equation is based on the principle of mass conservation within a control volume. This principle states that the difference between outflow and inflow within one time interval is balanced by a corresponding change in volume. In terms of finite intervals (i.e., finite differences) it is:
in which Q = flow; A = flow area; Δt = time interval; and Δx = space interval.
In differential form,
which is the equation of conservation of mass, or equation of continuity. The equation of conservation of momentum (Eq. 422) contains local inertia, convective inertia, pressure gradient (due to flow depth gradient), friction (friction slope), gravity (bed slope), and a momentum source term (Section 4.2). In deriving the kinematic wave equation, a statement of uniform flow is used in lieu of conservation of momentum. Since uniform flow is strictly a balance of friction and gravity, it follows that local and convective inertia, pressure gradient, and momentum source terms are excluded from the formulation of kinematic waves. In other words, a kinematic wave is a simplified wave that does not include these terms or processes. As shown later in this section, this simplification imposes limits to the applicability of kinematic waves. Uniform flow in open channels is described by the Manning or Chezy formulas (Section 2.4). The Manning equation is:
in which R is the hydraulic radius in meters, S_{f} is the friction slope in meters per meter, and n is the Manning friction coefficient. A pictorial on n is given by Barnes (1967). The Chezy equation is:
in which C = Chezy coefficient. Notice that in unsteady flow, friction slope is used in Eqs. 910 and 911 in lieu of channel slope. The hydraulic radius is R = A/P, in which P is the wetted perimeter. Substituting this into Eq. 910, leads to:
Assume for the sake of simplicity that n, S_{f}, and P are constant. This may be the case of a wide channel in which P can be assumed to be essentially independent of A. Equation 912 can then be written as:
in which α and β are parameters of the dischargearea rating (see rating curve, Section 2.4), defined as follows:
In Eq. 913, differentiating Q with respect to A leads to:
in which V is the mean flow velocity. Multiplying Eqs. 99 and 916 and applying the chain rule, the kinematic wave equation is obtained:
or, alternatively
Equation 917 (or 918) describes the movement of waves which are kinematic in nature . These are referred to as kinematic waves, i.e., waves for which inertia and pressure (flow depth) gradient have been neglected [10]. Equation 917 is a firstorder partial differential equation. Therefore, kinematic waves travel with wave celerity dQ/dA (or βV) and do not attenuate. Wave attenuation can only be described by a secondorder partial differential equation. The absence of wave attenuation can be further explained by resorting to a mathematical argument. Since dQ/dA is the celerity of the unsteady (i.e. , wavelike) Q, it can be replaced by dx/dt. Therefore, in Eq. 917:
which is equal to the total derivative dQ/dt. Since the right side of Eq. 919 is zero, it follows that Q remains constant in time for waves traveling with celerity dQ/dA. Discretization of Kinematic Wave Equation Equation 918 (or 917) is a nonlinear firstorder partial differential equation describing the change of discharge Q in time and space. It is nonlinear because the wave celerity βV (or dQ/dA) varies with discharge. The nonlinearity, however, is usually mild, and therefore, Eq. 918 can also be solved in a linear mode by considering the wave celerity to be constant. The solution of Eq. 918 can be obtained by analytical or numerical means. The simplest kinematic wave solution is a linear numerical solution. For this purpose, it is necessary to select a numerical scheme with which to discretize Eq. 918 on the xt plane (Fig. 98). A review of basic concepts of numerical analysis is necessary before discussing numerical schemes.
Order of Accuracy of Numerical Schemes. The order of accuracy of a numerical scheme measures the ability of the scheme to reproduce (i.e., recreate) the terms of the differential equation. In general, the higher the order of accuracy of a scheme, the better it is able to reproduce the terms of the differential equation. Forward and backward finite differences have firstorder accuracy, i.e., discretization errors of first order. Central differences have secondorder accuracy, with discretization errors of second order. When solving Eq. 918 by numerical methods, firstorder schemes create numerical diffusion and numerical dispersion, while secondorder schemes create only numerical dispersion. A thirdorder scheme creates neither numerical diffusion nor dispersion. Numerical diffusion and/or dispersion are caused by the finite grid size and are not necessarily related to the physical problem. Secondorderaccurate Numerical Scheme. The discretization of Eq. 918 following a linear secondorderaccurate scheme, i.e., using central differences in space and time, leads to (Fig. 98):
in which βV has been held constant (linear mode), leading to:
in which
and C is the Courant number, defined as follows:
Note that Courant number is the ratio of physical wave celerity βV to grid celerity Δx /Δt. The Courant number is a fundamental concept in the numerical solution of hyperbolic partial differential equations.
The three cases of Example 93 illustrate the properties of kinematic waves. The secondorderaccurate scheme has no numerical diffusion. In addition, for Courant number C = 1, i.e., the wave celerity βV equal to the grid celerity Δx/Δt, the scheme has no numerical dispersion, with the hydrograph being translated downstream without change in shape. In other words, the numerical solution by Eqs. 921 to 925 is exact only for Courant number C = 1. For other values of C, the numerical solution exhibits perceptible amounts of numerical dispersion. Firstorderaccurate Numerical Scheme. The numerical solution of Eq. 918 can also be attempted using a firstorderaccurate scheme, i.e., one featuring forward or backward finite differences. The discretization of Eq. 918 in a linear mode, using backward differences in both space and time yields (Fig. 98):
from which
in which
and C = Courant number, defined by Eq. 925.
Convex Method. The convex method of stream channel routing belongs to the family of linear kinematic wave methods. Through the 1970s, it was part of the SCS TR20 model for hydrologic simulation (Chapter 13). The routing equation for the convex method is obtained by discretizing Eq. 918 in a linear mode using a forwardintime, backwardinspace finite difference scheme, to yield (Fig. 98):
from which
in which
and C = Courant number (Eq. 925), restricted to C ≤ 1 for numerical stability reasons. In the convex method, C is regarded as an empirical routing coefficient. Example 95 illustrates the application of the convex method.
The convex method is relatively simple, but the solution is dependent on the routing parameter C.
The latter could be interpreted as a Courant number and related
to kinematic wave celerity and grid size, as in Eq. 925.
However, for values of C other than 1, the amount of diffusion introduced in the numerical problem is unrelated to the true diffusion, if any, of the physical problem. Therefore, the convex method, as well as all kinematic wave methods featuring uncontrolled amounts of numerical diffusion, are regarded as a somewhat crude approach to stream channel routing.
Kinematic Wave Celerity The kinematic wave celerity is dQ/dA, or βV. A value of β = 5/3 was derived for the case of a hydraulically wide channel governed by Manning friction. The kinematic wave celerity is also known as the KleitzSeddon, or Seddon, law [8, 19]. In 1900, Seddon [19] published a paper in which he studied the nature of unsteady flow movement in rivers and concluded that the celerity of long disturbances was equal to:
in which dQ/dy = slope of the dischargestage rating (Q versus y),
and T = top width.
The quantity c is the kinematic wave celerity.
Since dA = T dy, the kinematic wave celerity is equal to dQ/dA
From Eq. 934 it is concluded that the kinematic wave celerity is a function of the slope of the dischargestage rating. This slope is likely to vary with stage; therefore, the kinematic wave celerity is not constant but varies with stage and flow level. If c = βV is a function of Q, then Eq. 918 is a nonlinear equation requiring an iterative solution. Nonlinear kinematic wave solutions account for the variation of kinematic wave celerity with stage and flow level. The simpler linear solutions, as in Examples 93 and 94, assume a constant value of kinematic wave celerity βV. Notice that there is a striking similarity between the linear kinematic wave solutions and the Muskingum method. This subject is further examined in Section 9.4. Theoretical β values other than 5/3 can be obtained for other friction formulations and crosssectional shapes. For turbulent flow governed by Manning friction, β has an upper limit of 5/3, and it is usually greater than 1. For laminar flow in wide channels, β = 3; for mixed or transitional flowbetween laminar and turbulent Manning, it is in the range 5/3 < β < 3. For flow in a hydraulically wide channel described by the Chezy formula, β = 3/2 (Section 4.2). The calculation of β as a function of frictional type and crosssectional shape is illustrated by the following example.
Kinematic Waves with Lateral Inflow Practical applications of stream channel routing often require the specification of lateral inflows. The latter could be either concentrated, as in the case of tributary inflow at a point along the channel reach, or distributed along the channel, as with groundwater exfiltration (for effluent streams) or infiltration (for influent streams). As with Eq. 99, a mass balance leads to:
which, unlike Eq. 99, includes the source term q_{L}, the lateral flow per unit channel length. For Q given in cubic meters per second and x in meters, q_{L} is given in cubic meters per second per meter [L^{2} T ^{1}]. Multiplying Eq. 942 by ∂Q / ∂A (or βV), as with Eq. 917 (or Eq. 918), leads to:
which is the kinematic wave equation with lateral inflow (or outflow). For q_{L} positive, there is lateral inflow (e.g., tributary flow); for q_{L} negative, there is lateral outflow (e.g., channel transmission losses). Applicability of Kinematic Waves The kinematic wave celerity is a fundamental streamflow property. Flood waves which approximate kinematic waves travel with the kinematic wave celerity (c = βV) and are subject to very little or no attenuation. In practice, flood waves are kinematic if they are of long duration (Fig. 910) or travel on a channel of steep slope. Criteria for the applicability of kinematic waves to overland flow [20] (Section 4.2) and stream channel flow [14] have been developed. The stream channel criterion states that in order for a wave to be kinematic, it should satisfy the following dimensionless inequality:
in which t_{r} is the timeofrise of the inflow hydrograph, S_{o} is the bottom slope, V_{o} is the average velocity, and d_{o} is the average flow depth.
For 95% accuracy in one period of translation, a value of
9.3 DIFFUSION WAVES
In Section 9.1, the Muskingum method was used to calculate unsteady flows in a hydrologic sense. In Section 9.2, the principle of mass conservation was coupled with a uniform flow formula to derive the kinematic wave equation. Solutions to this equation have been widely used in hydrologic practice, particularly for overland flow and other routing applications involving steep slopes or slowrising hydrographs. The Muskingum method and linear kinematic wave solutions show striking similarities. Both methods have the same type of routing equation. The Muskingum method, however, can calculate hydrograph diffusion, whereas the kinematic wave can do so only by the introduction of numerical diffusion. The latter is dependent on the grid size and type of numerical scheme. Kinematic wave theory can be enhanced by allowing a small amount of physical diffusion in its formulation [10]. In this way, an improved type of kinematic wave can be formulated, a kinematicwithdiffusion wave, for short, a diffusion wave. A definite advantage of the diffusion wave is that it includes the diffusion which is present in most natural unsteady open channel flows. Diffusion Wave Equation In Section 9.2, the kinematic wave equation was derived by using a statement of steady uniform flow (i.e., friction slope is equal to bottom slope) in lieu of momentum conservation. In deriving the diffusion wave, a statement of steady nonuniform flow (i.e., friction slope is equal to water surface slope) is used instead (Fig. 911). This leads to:
in which the term [S_{o}  (d_{y}/d_{x})] is the water surface slope. The difference between kinematic and diffusion waves is in the term d_{y}/d_{x}. From a physical standpoint, the term d_{y}/d_{x} accounts for the natural diffusion processes present in unsteady open channel flow phenomena.
To derive the diffusion wave equation, Eq. 945 is expressed in a slightly different form:
in which m is the reciprocal of the square of the channel conveyance K, defined as:
With dA = T dy, in which T = top width, Eq. 946 changes to:
Equations 99 and 948 constitute a set of two partial differential equations describing diffusion waves. These equations can be combined into one equation with Q as dependent variable. However, it is first necessary to linearize the equations around reference flow values. For simplicity, a constant top width is assumed (i.e., a wide channel assumption). The linearization of Eqs. 99 and 948 is accomplished by small perturbation theory [4]. This procedure, while heuristic, has seemed to work well in a number of applications. The variables Q, A, and m can be expressed in terms of the sum of a reference value (with subscript o) and a small perturbation to the reference value (with superscript '): Q = Q_{o} + Q' ; A = A_{o} + A' ; m = m_{o} + m'. Substituting these into Eqs. 99 and 948, neglecting squared perturbations, and subtracting the reference flow leads to:
and
Differentiating Eq. 949 with respect to x and Eq. 950 with respect to t gives:
Using the chain rule and Eq. 949 yields:
Combining Eq. 952 with Eq. 953:
Combining Eqs. 951 and 954 and rearranging terms, yields:
Since by definition: mQ ^{2} = S_{f}, it follows that
and also
Substituting Eqs. 956 and 957 into Eq. 955, using the chain rule, and dropping the superscripts for simplicity, the following equation is obtained:
The left side of Eq. 958 is recognized as the kinematic wave equation, with ∂Q/∂A as the kinematic wave celerity.
The right side is a secondorder (partial differential) term that accounts for the physical diffusion effect.
The coefficient of the secondorder term has the units of diffusivity The hydraulic diffusivity is a characteristic of the flow and channel, defined as:
in which q_{o} = Q_{o}/T is the reference flow per unit of channel width. From Eq. 959, it is concluded that the hydraulic diffusivity is small for steep bottom slopes (e.g., those of small mountain streams), and large for mild bottom slopes (e.g., tidal rivers). Equation 958 describes the movement of flood waves in a better way than Eq. 917 or 918. It falls short from describing the full momentum effects, but it does physically account for peak flow attenuation. Equation 958 is a secondorder parabolic partial differential equation. It can be solved analytically, leading to Hayami's diffusion analogy solution for flood waves [7], or numerically with the aid of a numerical scheme for parabolic equations such as the CrankNicolson scheme [3]. An alternate approach is to match the hydraulic diffusivity with the numerical diffusion coefficient of the Muskingum scheme. This approach is the basis of the MuskingumCunge method [4, 12] (Section 9.4). Applicability of Diffusion Waves Most flood waves have a small amount of physical diffusion; therefore, they are better approximated by the diffusion wave rather than by the kinematic wave. For this reason, diffusion waves apply to a much wider range of practical problems than kinematic waves. Where the diffusion wave fails, only the dynamic wave can properly describe the translation and diffusion of flood waves. The dynamic wave, however, is very strongly diffusive, especially for flows well in the subcritical regime [14]. In practice, most flood flows are only mildly diffusive, and therefore, are subject to modeling with the diffusion wave. To determine if a wave is a diffusion wave, it should satisfy the following dimensionless inequality [14]:
in which t_{r} is the timeofrise of the inflow hydrograph, S_{o} is the bottom slope, d_{o} is the average flow depth, and g is the gravitational acceleration. The greater the left side of this inequality, the more likely it is that the wave is a diffusion wave. In practice, a value of N = 15 is recommended for general use.
9.4 MUSKINGUMCUNGE METHOD
The Muskingum method can calculate runoff diffusion, ostensibly by varying the parameter X. A numerical solution of the linear kinematic wave equation using a thirdorder accurate scheme (C = 1) leads to pure flood hydrograph translation (see Example 93, Part 1). Other numerical solutions to the linear kinematic wave equation invariably produce a certain amount of numerical diffusion and/or dispersion (See Example 93, Part 2). The Muskingum and linear kinematic wave routing equations are strikingly similar. Furthermore, unlike the kinematic wave equation, the diffusion wave equation does have the capability to describe physical diffusion. From these propositions, Cunge [4] concluded that the Muskingum method is a linear kinematic wave solution and that the flood wave attenuation shown by the calculation is due to the numerical diffusion of the scheme itself. To prove this assertion, the kinematic wave equation (Eq. 918) is discretized on the xt plane (Fig. 912) in a way that parallels the Muskingum method, centering the spatial derivative and offcentering the temporal derivative by means of a weighting factor X:
in which c = βV is the kinematic wave celerity.
Solving Eq. 961 for the unknown discharge leads to the following routing equation:
The routing coefficients are:
By defining the travel time
it is seen that the two sets of Eqs. 963 to 965 and Eqs. 94 to 96 are the same. Equation 966 confirms that K is in fact the flood wave travel time, i.e., the time it takes a given discharge to travel the reach length Δx with the kinematic wave celerity c. In a linear mode, c is constant and equal to a reference value; in a nonlinear mode, it varies with discharge.
It can be seen that for X = 0.5, Eqs. 963 to 965 reduce to the routing coefficients of the linear secondorderaccurate kinematic wave solution, Eqs. 922 to 924.
For X = 0.5 and C = 1 (C =
In practice, the numerical diffusion can be used to simulate the physical diffusion of the actual flood wave. By expanding the discrete function Q (jΔx, n Δt) in Taylor series about grid point (jΔx, n Δt), the numerical diffusion coefficient of the Muskingum scheme is derived (see Appendix B):
in which ν_{n} is the numerical diffusion coefficient of the Muskingum scheme. This equation reveals the following:
A predictive equation for X can be obtained by matching the hydraulic diffusivity ν_{h} (Eq. 959) with the numerical diffusion coefficient of the Muskingum scheme ν_{n} (Eq. 967). This leads to the following expression for X:
With X calculated by Eq. 968, the Muskingum method is referred to as MuskingumCunge method [12]. Using Eq. 968, the routing parameter X can be calculated as a function of the following numerical and physical properties:
It should be noted that Eq. 968 was derived by matching physical and numerical diffusion (i.e., secondorder processes), and does not account for dispersion (a thirdorder process). Therefore, in order to simulate wave diffusion properly with the MuskingumCunge method, it is necessary to optimize numerical diffusion (with Eq. 968) while minimizing numerical dispersion by keeping the value of C as close to 1 as practicable.
A unique feature of the MuskingumCunge method is the grid independence of the calculated outflow hydrograph, which sets it apart from other linear kinematic wave solutions featuring uncontrolled numerical diffusion and dispersion (e.g., the convex method).
If numerical dispersion is minimized, the calculated outflow at the downstream end of a channel reach will be essentially the same, regardless of how many subreaches are used in the computation.
This is because X is a function of Δx,
and the routing coefficients C_{0}, C_{1},
and An improved version of the MuskingumCunge method is due to Ponce and Yevjevich [15]. The C value is the Courant number, i.e. , the ratio of wave celerity c to grid celerity Δx/ Δt:
The grid diffusivity is defined as the numerical diffusivity for the case of X = 0. From Eq. 967, the grid diffusivity is:
The cell Reynolds number [18] is defined as the ratio of hydraulic diffusivity (Eq. 959) to grid diffusivity (Eq. 970). This leads to:
in which D = cell Reynolds number. Therefore:
Equations 971 and 972 imply that for very small values of Δx, D may be greater than 1, leading to negative values of X. In fact, for the characteristic reach length
the cell Reynolds number is D = 1, and X = 0. Therefore, in the MuskingumCunge method, reach lengths shorter than the characteristic reach length result in negative values of X. This should be contrasted with the classical Muskingum method (Section 9.1), in which X is restricted in the range 0.0 ≤ X ≤ 0.5. In the classical Muskingum, X is interpreted as a weighting factor. As shown by Eqs. 971 and 972, nonnegative values of X are associated with long reaches, typical of the manual computation used in the development and early application of the Muskingum method. In the MuskingumCunge method, however, X is interpreted in a momentmatching sense [2] or diffusionmatching factor. Therefore, negative values of X are entirely possible. This feature allows the use of shorter reaches than would otherwise be possible if X were restricted to nonnegative values. The substitution of Eqs. 969 and 972 into Eqs. 963 to 965 leads to routing coefficients expressed in terms of Courant and cell Reynolds numbers:
The calculation of routing parameters C and D, Eqs. 969 and 971, can be performed in several ways. The wave celerity can be calculated with either Eq. 916 or Eq. 934. With Eq. 916, c = βV; with Eq. 934, c = (1/T) dQ/dy. Theoretically, these two equations are the same. For practical applications, if a stagedischarge rating and crosssectional geometry are available (i.e., stagedischargetop width tables), Eq. 934 is preferred over Eq. 916 because it accounts directly for crosssectional shape. In the absence of a stagedischarge rating and crosssectional data, Eq. 916 can be used to estimate flood wave celerity. With the aid of Eqs. 969 and 971, the routing parameters may be based on flow characteristics. The calculations can proceed in a linear or nonlinear mode. In the linear mode, the routing parameters are based on reference flow values and kept constant throughout the computation in time. The choice of reference flow has a bearing on the calculated results [2, 15], although the overall effect is likely to be small. For practical applications, either an average or peak flow value can be used as reference flow. The peak flow value has the advantage that it can be readily ascertained, although a better approximation may be obtained by using an average value [15]. The linear mode of computation is referred to as the constantparameter MuskingumCunge method to distinguish it from the variableparameter MuskingumCunge method, in which the routing parameters are allowed to vary with the flow. The constant parameter method resembles the Muskingum method, with the difference that the routing parameters are based on measurable flow and channel characteristics instead of historical streamflow data.
Resolution Requirements When using the MuskingumCunge method, care should be taken to ensure that the values of Δx and Δt are sufficiently small to approximate closely the actual shape of the hydrograph. For smoothly rising hydrographs, a minimum value of t_{p} /Δt = 5 is recommended. This requirement usually results in the hydrograph time base being resolved into at least 15 to 25 discrete points, considered adequate for Muskingum routing. Unlike temporal resolution, there is no definite criteria for spatial resolution. A criterion borne out by experience is based on the fact that Courant and cell Reynolds numbers are inversely related to reach length Δx. Therefore, to keep Δx sufficiently small, Courant and cell Reynolds numbers should be kept sufficiently large. This leads to the practical criterion [16]:
which can be written as follows: 1 + C + D ≥ 0. This confirms the necessity of avoiding negative values of C_{0} in MuskingumCunge routing (see Eq. 974). Experience has shown that negative values of either C_{1} or C_{2} do not adversely affect the method's overall accuracy [16]. Notwithstanding Eq. 977, the MuskingumCunge method works best when the numerical dispersion is minimized, that is, when C ≅ 1. Values of C substantially less than 1 are likely to cause the notorious dips, or negative outflows, in portions of the calculated hydrograph. This computational anomaly is attributed to excessive numerical dispersion and should be avoided. Nonlinear MuskingumCunge Method The kinematic wave equation, Eq. 918, is nonlinear because the kinematic wave celerity varies with discharge. The nonlinearity is mild, among other things because the wave celerity variation is usually restricted within a narrow range. However, in certain cases it may be necessary to account for this nonlinearity. This can be done in two ways: (1) during the discretization, by allowing the wave celerity to vary, resulting in a nonlinear numerical scheme to be solved by iterative means; and (2) after the discretization, by varying the routing parameters, as in the variableparameter MuskingumCunge method [15]. The latter approach is particularly useful if the overall nonlinear effect is small, which is often the case. In the variable parameter method, the routing parameters are allowed to vary with the flow. The values of C and D are based on local q_{o} and c values instead of peak flow or other reference value as in the constantparameter method. To vary the routing parameters, the most expedient way is to obtain an average value of q_{o} and c for each computational cell. This can be achieved with a direct threepoint average of the values at the known grid points (see Fig. 911), or by an iterative fourpoint average, which includes the unknown grid point. To improve the convergence of the iterative fourpoint procedure, the threepoint average can be used as the first guess of the iteration. Once q_{o} and c have been determined for each computational cell, the Courant and cell Reynolds numbers are calculated by Eqs. 969 and 971. The value of bottom slope S_{o} remains unchanged within each computational cell. The variable parameter MuskingumCunge method represents a small yet sometimes perceptible improvement over the constant parameter method. The differences are likely to be more marked for very long reaches and/or wide variations in flow levels. Flood hydrographs calculated with variable parameters show a certain amount of distortion, either wave steepening in the case of flows contained inbank or wave attenuation (flattening) in the case of typical overbank flows. This is a physical manifestation of the nonlinear effect, i.e., different flow levels traveling with different celerities. On the other hand, flood hydrographs calculated using constant parameters do not show wave distortion. Assessment of MuskingumCunge Method The MuskingumCunge method is a physically based alternative to the Muskingum method. Unlike the Muskingum method where the parameters are calibrated using streamflow data, in the MuskingumCunge method the parameters are calculated based on flow and channel characteristics. This makes possible channel routing without the need for timeconsuming and cumbersome parameter calibration. More importantly, it makes possible extensive channel routing in ungaged streams with a reasonable expectation of accuracy. With the variableparameter feature, nonlinear properties of flood waves (which could otherwise only be obtained by more elaborate numerical procedures) can be described within the context of the Muskingum formulation. Like the Muskingum method, the MuskingumCunge method is limited to diffusion waves. Furthermore, the MuskingumCunge method is based on a singlevalued rating and does not take into account strong flow nonuniformity or unsteady flows exhibiting substantial loops in dischargestage rating (i.e., dynamic waves). Thus, the MuskingumCunge method is suited for channel routing in natural streams without significant backwater effects and for unsteady flows that classify under the diffusion wave criterion (Eq. 960). An important difference between the Muskingum and MuskingumCunge methods should be noted. The Muskingum method is based on the storage concept (Chapter 4) and, therefore, it is lumped, with the parameters K and X being reach averages. The MuskingumCunge method, however, is distributed in nature, with the parameters C and D being based on values evaluated at channel cross sections. Therefore, for the MuskingumCunge method to improve on the Muskingum method, it is necessary that the routing parameters evaluated at channel cross sections be representative of the channel reach under consideration. Historically, the Muskingum method has been calibrated using streamflow data. On the contrary, the MuskingumCunge method relies on physical characteristics such as rating curves, crosssectional data and channel slope. The different data requirements reflect the different theoretical bases of the methods, i.e., lumped storage concept in the Muskingum method, and distributed kinematic/diffusion wave theory in the MuskingumCunge method. 9.5 INTRODUCTION TO DYNAMIC WAVES
In Section 9.2, kinematic waves were formulated by simplifying the momentum conservation principle to a statement of steady uniform flow. In Section 9.3, diffusion waves were formulated by simplifying the momentum principle to a statement of steady nonuniform flow. These two waves, in particular the diffusion wave, have been extensively used in stream channel routing applications. The Muskingum and MuskingumCunge methods are examples of calculations using the concept of diffusion wave. A third type of open channel flow wave, the dynamic wave, is formulated by taking into account the complete momentum principle, including its inertial components. As such, the dynamic wave contains more physical information than either kinematic or diffusion waves. Dynamic wave solutions, however, are more complicated than either kinematic or diffusion wave solutions. In a dynamic wave solution, the equations of mass and momentum conservation are solved by a numerical procedure, either the method of finite differences, the method of characteristics, or the finite element method. In the method of finite differences, the partial differential equations are discretized following a chosen numerical scheme [9]. The method of characteristics is based on the conversion of the set of partial differential equations into a related set of ordinary differential equations, and the solution along a characteristic grid, i.e. a grid that follows characteristic directions (Fig. 416). The method of finite elements solves a set of integral equations over a chosen grid of finite elements. In the past four decades, the method of finite differences has come to be regarded as the most expedient way of obtaining a dynamic wave solution for practical applications [6, 9]. Among several numerical schemes that have been used in connection with the dynamic wave, the Preissmann scheme is perhaps the most popular. This is a fourpoint scheme, centered in the temporal derivatives and slightly offcentered in the spatial derivatives. The offcentering in the spatial derivatives introduces a small amount of numerical diffusion necessary to control the numerical stability of the nonlinear scheme. This produces a workable yet sufficiently accurate scheme. The stream channel is divided into several reaches for computational purposes (Fig. 912). The application of the Preissmann scheme to the governing equations for the various reaches results in a matrix solution requiring a double sweep algorithm, i.e., one that accounts only for the nonzero entries of the coefficient matrix, which are located within a narrow band surrounding the main diagonal. This technique leads to a considerable savings in storage and execution time. With the appropriate upstream and downstream boundary conditions (Fig. 913), the solution of the set of hyperbolic equations marches in time until a specified number of time intervals is completed.
In practice, a dynamic wave solution represents an orderofmagnitude increase in complexity and associated data requirements when compared to either kinematic or diffusion wave solutions. Its use is recommended in situations where neither kinematic nor diffusion wave solutions are likely to represent adequately the physical phenomena. In particular, dynamic wave solutions are applicable to flow over very flat slopes, flow into large reservoirs, strong backwater conditions and flow reversals. In general, the dynamic wave is recommended for cases warranting a precise determination of the unsteady variation of river stages. Relevance of Dynamic Waves to Engineering Hydrology Dynamic wave solutions are often referred to as hydraulic river routing. As such, they have the capability to calculate unsteady discharges and stages when presented with the appropriate geometric channel data and initial and boundary conditions. Their relevance to engineering hydrology is examined here by comparing them to kinematic and diffusion wave solutions. Kinematic waves calculate unsteady discharges; the corresponding stages are subsequently obtained from the appropriate rating curves. Usually, equilibrium (steady, uniform) rating curves are used for this purpose. Diffusion waves may or may not use equilibrium rating curves to calculate stages. Some methods, e.g., MuskingumCunge, use equilibrium ratings, but more elaborate diffusion wave solutions may not. Dynamic waves rely on the physics of the phenomena as built into the governing equations to generate their own unsteady rating. A looped rating curve is produced at every cross section, as shown in Fig. 914. For any given stage, the discharge is higher in the rising limb of the hydrograph and lower in the receding limb. This loop is due to hydrodynamic reasons and should not be confused with other loops, which may be due to erosion, sedimentation, or changes in bed configuration (Chapter 15).
The width of the loop is a measure of the flow unsteadiness, with wider loops corresponding to highly unsteady flow, i.e., dynamic wave flow. If the loop is narrow, it implies that the flow is mildly unsteady, perhaps a diffusion wave. If the loop is practically nonexistent, the flow can be approximated as kinematic flow. In fact, the basic assumption of kinematic flow is that momentum can be simulated as steady uniform flow, i.e., that the rating curve is singlevalued. The preceding observations lead to the conclusion that the relevance of dynamic waves in engineering hydrology is directly related to the flow unsteadiness and the associated loop in the rating curve. For highly unsteady flows such as dambreak flood waves, it may well be the only way to properly account for the looped rating. For other less unsteady flows, kinematic and diffusion waves are a viable alternative, provided their applicability can be clearly demonstrated (Eqs. 944 and 960). Diffusion Wave Solution with Dynamic Component A simplified approach to dynamic wave routing is that of the diffusion wave with dynamic component [2]. In this approach, the complete governing equations, including inertia terms, are linearized in a similar way as with diffusion waves. This leads to a diffusion equation similar to Eq. 958, but with a modified hydraulic diffusivity. The equation is [5]:
in which the hydraulic diffusivity (i.e., the coefficient of the secondorder term) is also a function of the rating curve parameter β and the Froude number, defined as:
with g = gravitational acceleration, and d_{o} = reference flow depth. Equation 978 can be expressed in terms of the Vedernikov number [17]:
With Eq. 980, Eq. 978 reduces to:
Equations 978 and 981 provide an enhanced predictive capability for the simulation of diffusion waves including a dynamic component. For instance, for β = 1.5 (i.e., Chezy friction in wide channels) and F_{o} = 2, the Vedernikov number V = 1 (Eq. 980), and the hydraulic diffusivity vanishes, which is in agreement with physical reality [10, 13]. On the other hand, the hydraulic diffusivity of the diffusion wave (Eq. 958) is independent of the Vedernikov number. Therefore, Eq. 981 is a better model than Eq. 958, especially for Froude numbers in the supercritical regime. Most natural flows, however, are in the range well below critical, with Eq. 958 remaining a practical model of unsteady open channel flow phenomena [7]. QUESTIONS
PROBLEMS
REFERENCES
SUGGESTED READINGS

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