1.  INTRODUCTION

The simplified methods of channel routing enjoy a position of preeminence within the gamut of available routing methods. These methods are simple to formulate and use, and may well account for most of the flood wave phenomena when practical applications are considered. The more refined versions simulate the diffusion wave model; therein lies their inherent strength.

A simplified method which is receiving increased attention is a Muskingum method in which the routing parameters are calculated based on channel and grid characteristics (1), herein referred to as Muskingum diffusion model. Its numerical behavior, however, is not completely understood. For instance, there is no generally accepted spatial resolution criteria to preserve accuracy. Evidence of the continuing controversy is the number of recent contributions (2,7,10), but the matter still remains to be clarified.

A direct consequence of insufficient spatial resolution is the notorious negative outflows which have been documented in certain cases (3,10). Until the reason for this anomaly is identified, the simplified methods will suffer from a lack of credibility. This can only hamper their wide acceptance for practical channel routing applications.

The objective of this paper is, thus, twofold: (1) To review and examine recent contributions regarding spatial and temporal resolution in Muskingum diffusion models; and (2) to develop a criterion for spatial resolution to preserve the acccuracy of the method. Accuracy here is understood to mean that which is attainable within the framework of diffusion wave routing.

2.  RECENT CONTRIBUTIONS ON ACCURACY CRITERIA

Koussis (2) was among the first to recognize the effect of insufficient spatial resolution on the hydrograph calculated by the Muskingum diffusion model, and to tie this effect to the grid size. He proposed, "in order to obtain reasonable results," the following restriction:

2 X Δx Δt ≥ ________              c

(1)

in which Δt = time step; Δx = space step; X = Muskingum parameter; and c = flood wave celerity for the reference flow. Koussis did not elaborate on how Eq. 1 was developed, and it is attributed here to his experience with the method.

In their recent review of approximate flood routing methods, Weinman and Laurenson (10) focus on the Muskingum diffusion model and its variations. They point out the possibility of unrealistic negative outflows, but suggest that for practical purposes these negative outflows are small enough and sufficiently short-lived to be ignored if the following inequality is satisfied:

Δx           Tr _____ ≤  _______    c             2X

(2)

in which T_r = period-of-rise of the inflow hydrograph.

The striking similarity between Koussis' and Weinman and Laurenson's criteria is apparent. A basis for comparison can be implemented if T_r is expressed in terms of Δt. Assuming that in a given case T_r / Δt = n is an integer greater than 5, Eq. 2 reduces to:

Δx          n Δt _____ ≤  ________    c              2X

(3)

which indicates that Koussis' criterion is more conservative in the estimation of Δx than Weinmann and Laurenson's.

The writers elsewhere (7) reported the results of preliminary findings that seemed to indicate that if accuracy is to be preserved (i.e., avoidance of negative outflows), the following inequality should hold:

C D ≥ ξ

(4)

in which C and D are Courant and cell Reynolds numbers, respectively, defined as follows:

Δt C = c  ______              Δx

(5)

qo D = _________         So c Δx

(6)

in which q_o = reference unit-width discharge; and S_o = equilibrium water surface slope. Ponce and Theurer further stated that a value of ξ = 0.25 provided satisfatory accuracy in most cases.

Kundzewicz (3) has applied the convolution method of analysis to the Muskingum diffusion model, and concluded that Eq. 2 does not hold for the general case. He further stated that a search for a general condition to supress negative outflows may prove to be fruitless.

3.  CHOICE OF SPACE AND TIME STEPS

Accuracy criteria such as Eqs. 1, 2 or 4 provide guidelines for the selection of space and time steps if negative outflows are to be eliminated. Following Ponce and Yevjevich (5), space and time steps are used together with channel and flood wave descriptors to calculate the routing parameters C and D, Eqs. 5 and 6. In turn, these parameters are used in the routing equation (Fig. 1):

n + 1                     n                   n + 1                   n Qj + 1 = C1 Qj + C2 Qj    + C3 Qj + 1

(7)

Fig. 1  Spatial and temporal discretization in Muskingum diffusion model.

in which C₁, C₂ and C₃ are defined as follows:

1 + C - D C1 = ____________           1 + C + D

(8)

- 1 + C + D C2 = ____________           1 + C + D

(9)

1 - C + D C3 = ____________           1 + C + D

(10)

Choice of space step. Given Eqs. 1, 2 and 4, the choice of space step is not altogether clear. An upper limit on Δx is certainly a necessity in any numerical procedure, and the present method is no exception. Equation 1 would seem to indicate that the upper limit should be:

c Δt Δx ≤ ______            2X

(11)

However, this condition suffers from a significant drawback; it is undefined for X = 0. In Muskingum diffusion routing, X = 0 implies that Δx is fixed, and Eq. 11 could no longer be applicable.

The same type of criticism holds for Eq. 2. Equation 4, on the other hand, is expressed in terms of D rather than X, and an expression for the upper limit on Δx based on this equation is:

qo Δt Δx ≤ (_______)1/2              So ξ

(12)

The question of a lower limit on Δx is also one that has led to some confusion, primarily due to the existence of a functional relationship between X and Δx:

qo Δx = 1/2 (1 - _________)                        So c Δx

(13)

Historically, the parameter X has been interpreted as a weighting factor, and this led Miller and Cunge (4), among others, to suggest a lower bound for X: X ≥ 0. This is echoed in Weinmann and Laurenson's review (10), which states (p. 1529): "Since X may not realistically decrease below zero, there is a limit to the allowable degree of subdivision, which implies a lower limit on Δx." Given Eq. 13, this lower bound can be calculated by setting X ≥ 0, from which:

qo Δx ≥ ______           S0 c

(14)

The right-hand side of Eq. 14 is immediately recognized as the "characteristic reach length" first identified by Kalinin and Milyukov (4) in connection with the routing method bearing their name.

Further analysis and experience with the method has shown, however, that in fact there is no theoretical lower limit on Δx. Ponce and Theurer (7), as well as others (3,9), have reported experience with Δx values violating Eq. 14, and yet with no physical unreality whatoever attached to the results of the computations. The reason for this behavior is apparent when the moment-matching characteristic of Eq. 13 is recognized. As long as this equation is satisfied, Δx may be arbitrarily chosen to violate Eq. 14, implying the possibility of negative values of X. Practical computations using such a condition help dispell the myth of X being a weighting factor to be restricted in the range 0.0 ≤ X ≤ 0.5.

Choice of Time Step. The criterion for selecting the time step has also led to some confusion. An upper limit on Δt is clearly a necessity if the solution is to remain within a reasonable numerical framework. Unlike spatial resolution, temporal resolution is a concept more readily grasped in hydrologic applications, since hydrographs are depicted in the temporal domain. Established practice is to set an upper limit on Δt such that:

Tr Δt ≤ ______            n

(15)

in which n = an integer having a minimum value of 5. In practice, hydrographs are skewed to the right; therefore, Eq. 15 guarantees a satisfactory level of temporal resolution.

In addition to the previously given requirement of temporal resolution, there is a question as to whether the Muskingum diffusion model is bound by the Courant numerical stability criterion usually applicable to hyperbolic problems. If this is the case, then an additional constraint on Δt should be imposed:

Δx Δt  ≤  ______               c

(16)

which follows from:

C ≤ 1

(17)

Practical computations, however, have shown that violating Eqs. 16 and 17 does not lead to numerical instability. This behavior is explained by recalling that the method is essentially a finite difference formulation of the kinematic wave equation. In theory (6), this feature provides unconditional stability with respect to the Courant condition, thereby invalidating Eq. 17.

A lower limit on Δt is mainly a question of computational resources. A decision on this matter is left to the individual modeler.

Accuracy Criteria. From the preceding analysis, it follows that there is no theoretical justification for a lower limit on either Δt or Δx. Furthermore, while there is an upper limit on Δt given by the temporal resolution requirement, a similar upper limit on Δx is not clearly discernible. Given Eqs. 5, 6 and 13, Koussis' criterion is modified here and expressed in terms of C and D:

C + D  ≥ 1

(18)

In this form, the overdeterminacy of Eq. 11 is removed. In effect, fo X = 0, D = 1 and Eq. 18 reduces to C ≥ 0. Since X = 0 implies that Δx is finite (Eq. 6), the condition C ≥ 0 reduced to the trivial one Δt ≥ 0.

Weinmann and Laurenson's criterion, Eq. 3, leads to

nC + D ≥ 1

(19)

and similar comments to that of Koussis' apply. Ponce and Theurer's remains:

C D ≥ 0.25

(20)

The following section describes a program of numerical experiments designed to throw additional light onto the question of accuracy criteria.

Consistency. Closely related to accuracy is the subject of consistency of a simplified routing method such as the one treated here. In the context of this study, accuracy refers to the avoidance of physically unrealistic negative outflows; consistency refers to the abiliity of the routing method to produce the same results regardless of grid size. For instance, if a certain hydrograph is obtained at the end of a given channel reach using ten (10) space steps, a consistent method would produce the same hydrograph if instead only two (2) space steps had been used.

Consistency, as previously defined, is to be used within the framework of a simplified routing method in which the objective is to optimize the numerical diffusion by making it match the physical diffusion. This concept should not be confused with the property of consistency of a conventional finite difference scheme, in which the objective is to minimize the amount of numerical error in order to properly solve the differential equation under consideration.

4.  NUMERICAL EXPERIMENTS

A program of numerical experiments was carried out in order to evaluate the accuracy (and therefore, consistency) of the Muskingum diffusion routing described by Eqs. 7-10. The following considerations were used in the design of the program:

1. To provide a wide variation in Courant and cell Reynolds numbers.

2. To provide a wide variation in grid size: space step Δx and time step Δt.

3. To provide a representative sample of channel propertied, specifically unit-width discharge q_o, bed slope S_o, and flood wave celerity, c.

Six hypothetical channels were developed and tested in the course of this study. For each channel, two flow conditions were specified: (1) A low base flow; and (2) a high base flow. Each channel and flow condition was tested by routing the same hydrograph through 12 different grid specifications, making a total of 144 runs. The 12 grid specifications represented combinations of three (3) time steps and four (4) space steps. To simulate natural hydrographs, a gamma distribution was used as the inflow hydrograph. This function is of the following form:

t                     tp - t Qt = Qb + (Qp - Qb) ( ___ )m exp ( ________ )                                     tp                   tg - tp

(21)

in which Q_t = inflow at time t; Q_b = base flow; Q_p = peak flow; t_p = time from beginning of direct runoff to peak of hydrograph; t_g = time from beginning of direct runoff to center of gravity of hydrograph; and m is given by the following equation:

tp m = _______          tg - tp

(22)

Table 1 shows channel and hydrograph characteristics for the twelve (12) testing series. A detailed summary of the results of the numerical experiments is given in Ref. 8.

5.  ANALYSIS OF RESULTS

The results of the numerical experiments can be interpreted in the light of the two routing parameters: Courant number, C, and cell Reynolds number, D. To start with, it is noted that Δx is inversely proportional to both C and D. Therefore, the larger the Δx, the smaller the values of C and D.

Figure 2 is a logarithmic plot of the values of C and D associated with all runs in the testing program. Also shown in this figure are the lines for C₁ = 0, C₂ = 0, and C₃ = 0. An analysis of the numerical results strongly indicates that accuracy is dependent on the value of C₂. Negative values of C₂ are invariably associated with dips in the rising portion of the outflow hydrograph, indicating the tendency of outflows to become "negative" (i.e., values less than the base flow). This tendency increases as the value of C₂ increases in the negative range.

Fig. 2  Courant and cell Reynolds Numbers of test series to determine accuracy criteria.

Values of C₂ in the range 0.0-0.5 do not show dips in the outflow hydrograph. However, depending on the problem, small inaccuracies in the calculated hydrographs may still remain. >

Values of C₂ in excess of 0.5 are associated with hydrographs which do not show any of the previous anomalies. In addition, once accuracy is achieved, say for C₂ ≥ 0.5, consistency follows. In other words, hydrographs calculated in the C₂ ≥ 0.5 range show excellent reproducibility of peak flow, time to peak and hydrograph shape, with varying Δx and Δt.

6.  SPATIAL RESOLUTION ACCURACY CRITERIA

The experimental evidence suggests that C₂ and only C₂ is determinant of accuracy. In fact, values of C₁ and C₃ were monitored throughout the testing series, and negative values for either were shown not to affect the accuracy of the calculated hydrograph. This finding points to the fallacy of limiting all three routing coefficients to nonnegative values in order to effect a "proper weighting," as argued by Miller and Cunge (4), among others.

The condition C₂ ≥ 0 is further examined here. This leads to, with Eq. 9, C + D ≥ 1, which is the modified Koussis criterion developed in the previous section.

A stronger condition is given by C₂ ≥ 0.5, which leads to

C + D ≥ 3

(23)

In general, a condition such as C₂ ≥ ζ  is applicable, in which ζ  is a positive real number. This leads to:

1 + ζ C + D ≥ ______ = κ               1 - ζ

(24)

It is instructive to examine the physical significance of the preceding criterion. Given Eqs. 5 and 6, Eq. 24 can be expressed as follows:

1 Δx ≤ ___ ( ΔxC + ΔxD )           κ

(25)

in which

ΔxC = c Δt

(26)

is the Courant length, defined as the length of channel traveled by the flood wave in one time step Δt; and

qo ΔxD = _____            So c

(27)

in the "characteristic reach length," defined in terms of the frictional and cross-sectional characteristics of the channel.

It is further noticed that:

qo = uo do

(28)

in which u_o = reference flow mean velocity; and d_o = reference flow depth; and

c = β uo

(29)

in which β = the exponent in the discharge flow area relation Q = α A^β, from which:

do ΔxD = _____            β So

(30)

Given Eqs. 29 and 30, Eq. 24 alternatively can be expressed as follows:

1                       do Δx ≤ ___ ( β uo Δt + _____ )           κ                     β So

(31)

in which β, u_o, d_o, and S_o = channel and flow properties; and κ = an accuracy parameter determined experimentally. A value of κ = 2 is sufficient to preserve accuracy for practical applications.

7.  SUMMARY AND CONCLUSIONS

An analysis of the Muskingum diffusion routing is carried out. Numerical experiments are performed to throw additional light onto the criterion setting an upper limit to the space step Δx if accuracy is to be preserved. Computational experience indicates that for very large values of Δx, there is a tendency for physically unrealistic negative outflows to occur.

The results of the numerical experiments indicate that the value of C₂ is determinant of accuracy. The condition C₂ ≥ ζ  is postulated to preserve the accuracy (and consistency) of the routing method, in which ζ is a positive real number. This leads to an upper bound for Δx given by the following equation:

1 Δx ≤ ___ ( ΔxC + ΔxD )           κ

(25)

A value of ζ = 0.33 is recommended, from which κ = 2. The quantities Δx_C and Δx_D are defined in terms of channel and grid characteristics.

 ACKNOWLEDGMENTS

The research on which this paper is based was funded in part by the USDA Soil Conservation Service, National Engineering Staff, Lanham, Maryland. The assistance of George H. Comer, Soil Conservationist, SCS, is gratefully acknowledged.

APENDIX I.  REFERENCES

1. Cunge, J. A. 1969. "On the subject of a flood propagation computation method (Muskingum method)," Journal of Hydraulic Research, Vol. 7, No. 2, 205-230.

2. Koussis, A. 1976. "An approximate dynamic flood routing method," Proceedings of the International Symposium on Unsteady Flow in Open Channels, Paper L1, Newcastle-Upon-Tyne, England.

3. Kundzewicz, Z. W. 1980. Discussion of "Approximate flood routing methods: A review," by P. Erwin Weinmann and Eric M. Laurenson, Journal of the Hydraulics Division, ASCE, Vol. 106, No. HY12, Proc. Paper 15866, Dec., 2072-2075.

4. 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., 183-758.

5. Ponce, V. M., and V. Yevjevich, V. 1978. "The Muskingum-Cunge method with variable parameters," Journal of the Hydraulics Division, ASCE, Vol. 104, No. HY12, Proc. Paper 14199, Dec., 1663-1667.

6. Ponce, V. M., Y. H. Chen, and D. B. Simons. 1979. "Unconditional stability in convection computations," Journal of the Hydraulics Division, ASCE, Vol. 105, No. HY9, Proc. Paper 14807, Sept., 1079-1086.

7. Ponce, V. M., and Theurer, F. D. 1980. Discussion of "Approximate flood routing methods: A review," by P. Erwin Weinmann and Eric M. Laurenson, Journal of Hydraulic Division, ASCE, Vol. 106, No. HY11, Proc. Paper 15784, Nov., 1945-1947.

8. Ponce, V. M. 1981. "Development of an algorithm for the linearized diffusion method of flood routing," Research Reposrt, San Diego State University Civil Engineering Series, No. 81144, San Diego, Calif., Sept.

9. Weinmann, P. E. 1977. "Comparison of flood routing methods for national rivers," Report No. 2/1977, Department of Civil Engineering, Monash university, Clayton, Victoria, Astralia.

10. Weinmann, P.E., and F. M. Laurenson. 1979. "Approximate flood routing methods: A review," Journal of the Hydraulics Division , ASCE, Vol. 105, No. HY12. Proc. Paper 15057, Dec., 1521-1536.

11. Weinmann, P. E., and F. M. Laurenson. 1981. Closure to "Approximate flood routing methods: A review," Journal of the Hydraulics Division, ASCE, Vol. 107, No. HY9, Proc. Paper 16479, Sept., 1112-1114.