Accuracy Criteria in Diffusion Routing


Victor M. Ponce and Fred D. Theurer



Online version 2018
[Original version 1982]



Abstract. An analysis of Muskingum diffusion model, in which routing paČrameters are calculated based on channel and grid chancteristics. is carried out. Numerical experiments are performed to throw additional light onto the criterion setting an upper limit to the space step if accuracy is to be preserved. Computational experience indicates that for very large values of the space step, there is a tendency fo physically unrealistic negative outflows to occur. The results of the numerical experiments point to the condition C2 ≥ ζ if accuracy is a real number. A value of ζ = 0.33 is recommended fro pratical applications. The spatial resolution criterion is experssed in terms of Courant and cell Reynolds numbers.


 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, and 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 on the subject 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).


 RECENT CONTIBUTIONS 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.

Weinman and Laurenson (10), in their recent review of approximate flood routing methods, 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 Tr = period-of-rise of inflow hydrograph.

The striking similarity between Koussis' and Weinman and Laarenson's criteria is apparent. A basis for comparison can be implemented if Tr is expressed in terms of Δt. Assuming that in a given case Tr / Δ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


CD ≥ ξ

(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 qo = reference unit-width discharge; and So = 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.


 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 udes 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)

mapa USA

Fig. 1  Spatial and Temporal Discretization in Muskingum Diffusion Model

in which C1, C2 and C3 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. Eq. 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 Eq. 4 , on the other hand, is expressed in terms of D rather than X, and 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, 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 can be arbitrarily chosen to violate Eq. 14, implying the possibility of negative values of X. Pratical 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)

Practicel 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


CD ≥ 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 10 space steps, a consistent method would produce the same hydrograph if instead only 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.


 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 numebers.

  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 qo, bed slope So, 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 3 time steps and 4 space steps. To simulate natural hydrographs, a gamma distribution was used as the inflow hydrograph. This function is of the following form:

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


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

(21)

in which Qt = inflow at time t; Qb = base flow; Qp = peak flow; tp = time from beginning of direct runoff to peak of hydrograph; tg = 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 12 testing series. A detailed summary of the results of the numerical experiments is given in Ref. 8.

Table. 1  Channel and Hydrograph Characteristics

mapa USA


 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 inverseley proportional to both C and D. Therefore, the larger the Δx, the smaller the values of C and D.

Fig. 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 C1 = 0, C2 = 0, and C3 = 0. An analysis of the numerical results strongly indicates that accuracy is dependent on the value of C2 . Negative values of C2 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 C2 increases in the negative range.

mapa USA

Fig. 2  Courant and Cell Reynolds Numbers of Test Series to Determine Accuracy Criteria

Values of C2 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 C2 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 C2 ≥ 0.5, consistency follows. In other words, hydrographs calculated in the C2 ≥ 0.5 range show excellent reproducibility of peak flow, time to peak and hydrograph shape, with varying Δx and Δt.


 SPATIAL RESOLUTION ACCURACY CRITERIA

The experimental evidence suggests that C2 and only C2 is determinant of accuracy. In fact, values of C1 and C3 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 C2 ≥ 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 C2 ≥ 0.5, which leads to


C +D ≥ 3

(23)

In general, a condition such as C2ζ  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 uo = reference flow mean velocity; and do = 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 β, uo, do, and So = channel and flow properties; and κ = an accuracy parameter determined experimentally. A value of κ = 2 is sufficient to preserve accuracy for pratical applications.


 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 C2 is determinant of accuracy. The condition C2ζ id 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 ΔxC and ΔxD are defined in terms of channel and grid characteristics.


 ACKNOWLEDGMENT

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., "On the Subject of a Flood Propagation Computation Method (Muskingum Method)," Journal of Hydraulic Research, Vol. 7, No. 2, 1969, pp. 205-230.

  2. Koussis, A., "An Approximate Dynamic Flood Routing Method," Proceedings of the International Symposium on Unsteady Flow in Open Channels, Paper L1, Newcastle-Upon-Tyne, England, 1976.

  3. Kundzewicz, Z. W., 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., 1980, pp. 2072-2075.

  4. Miller, W. A., and Cunge, J. A., "Simplified Equations of unsteady Flow," Unsteady Flow in Open channels, K. Mahmood and V. Yevjevich, eds., Water Resources Publications, Fort Collins, Colo., 1975, pp. 183-758.

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

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

  7. Ponce, V. M., and Theureer, F. D., 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., 1980, pp. 1945-1947.

  8. Ponce, V. M., "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., 1981.

  9. Weinmann, P. E., "Comparison of Flood Routing Methods for National Rivers," Report No. 2/1977, Department of Civil Engineering, Monash university, Clayton, Victoria, Astralia, 1977.

  10. Weinmann, P.E., and Laurenson, E. M., "Approximate Flood Routing Methods: A Review," Journal of the Hydraulics Division , ASCE, Vol. 105, No. HY12. Proc. Paper 15057, Dec., 1979, pp. 1521-1536.

  11. Weinmann, P. E., and Laurenson, E. M., Closure to "Approximate Flood Routing Methods: A Review," Journal of the Hydraulics Division, ASCE, Vol. 107, No. HY9, Proc. Paper 16479, Sept., 1981, pp. 1112-1114.


APENDIX II. - NOTATION

The following symbols are used in this paper:

    A = flow area;

    C = Courant number, Eq. 5;

    C1 = routing coefficient, Eqs. 7 and 8;

    C2 = routing coefficient, Eqs. 7 and 9;

    C3 = routing coefficient, Eqs. 7 and 10;

    c = reference flood wave celerity;

    D = cell Reynolds number, Eq. 6;

    do = reference flow depth;

    j = spatial index;

    m = exponent, defined by Eq. 22;

    n = temporal index, integer in Eq. 3;

    Q = discharge;

    Qb = base flow;

    Qp = peak flow;

    Qt = inflow at time t;

    qo = reference unit-width discharge;

    So = equilibrium water surface slope;

    Tr = period-of-rise of inflow hydrograph;

    t = time variable;

    tg = time from beggining of direct runoff to center of gravity of hydrograph;

    tp = time from beggining of direct runoff to peak of hydrograph;

    uo = reference flow mean velocity;

    X = Muskingum parameter;

    α = coeficient in Q = α A β;

    β = exponent in Q = α A β;

    Δt = time step;

    Δx = space step;

    ΔxC = Courant length

    ΔxD = characteristic reach length;

    ξ = accuracy parameter, Eq. 4;

    κ = accuracy parameter, Eqs. 24 and 25; and

    ζ = positive real number, in Eq. 24.


180829 10:30

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