Simplifying Eq. 2:

X (Q j n+1 - Q j n )  +  (1 - X) (Q j+1n+1 - Q j+1 n )   +   0.5 C [ (Q j+1 n - Q j n )  +  (Q j+1n+1 - Q j n+1 ) ] =  0

(3)

in which C = Courant number, defined as follows:

Δt       C =  c  ______                Δx

(4)

The solution of Eq. 3 is postulated in the following exponential form (O'Brien et al., 1950; Cunge, 1969):

Q = Q* e i (σx - βt )

(5)

in which σ = (2π /L) is the wavenumber, L is the wavelength of the disturbance, and β = (2 π /T ) is the wave period (Ponce and Simons, 1977).

4.  CONVERGENCE RATIOS

The derivation of the convergence ratios has been described in detail by Ponce and Vuppalapati (2016). The convergence ratio with regard to wave amplitude is:

[(P - M)2 + N 2]1/2       R1 =    ____________________                           P

(6)

in which:

P = 1 - 2 (1 - p)(1 - S) S

(6a)

M = (1-p)(1 - 2S)C

(6b)

N = - Cq

(6c)

wherein

p = cos(σΔx)

(6d)

q = sin(σΔx).

(6e)

S = X - (C/2)

(6f)

and C = Courant number (Eq. 4).

The convergence ratio with regard to wave phase is:

1                       N R2 =  _________ tan-1 ( ________ )            C σΔx                P - M

(7)

Therefore, given: (1) weighting factor X (in Eq. 2), (2) spatial resolution L/Δx [from which σΔx = (2π /L)Δx = 2π /(L/Δx)], and (3) Courant number C (Eq. 4), Eqs. 6 and 7 enable the calculation of the Muskingum-Cunge amplitude and phase convergence ratios, respectively. The calculation of the convergence ratios is accomplished using the calculator Online Muskingum-Cunge Convergence Ratios.

5.  AMPLITUDE AND PHASE PORTRAITS

The plotting of convergence ratios R₁ and R₂ as a function of Courant number C (in the abscissas) and spatial resolution L/Δx (as the curve parameter) leads to a set of amplitude and phase portraits R₁ and R₂, respectively, a pair for each value of weighting factor X in the range 0.0 ≤ X ≤ 0.5. We observe that this range of analysis for the weighting factor X is amply justified by kinematic wave theory (Ponce et al., 1979 a).

Figures 2 to 8 show the amplitude and phase portraits of the Muskingum-Cunge method (Eq. 3), for the following conditions: (a) Courant number 0.1 ≤ C ≤ 4.0, (b) spatial resolution 20 ≤ L/Δx ≤ 100, shown at intervals of L/Δx = 20, and (c) weighting factor 0.0 ≤ X ≤ 0.6, shown at intervals of X = 0.1 (7 × 2 = 14 graphs). Note that the value X = 0.6 (Fig. 8) is included here only for comparison purposes, i.e., to show the numerical effect of a value X > 0.5, which results in numerical amplification, R > 1; see Fig. 8 (a).

The selected range of Courant numbers shown (0.1 ≤ C ≤ 4.0) encompasses the central value C = 1, for which absolute accuracy is obtained for X = 0.5; see Figs. 7 (a) for amplitude portrait and 7 (b) for phase portrait. The range of spatial resolution shown in Figures 2 to 8 (20 ≤ L/Δx ≤ 100) depicts five (5) values, varying from low (L/Δx = 20; poor resolution) to high (L/Δx = 100; very good resolution). The spatial resolution is the number of discrete intervals Δx per wavelength L.

6.  ANALYSIS

The examination of Figs. 2 to 8 leads to the following conclusions:

1. Perfect convergence amounts to R₁ = 1 and R₂ = 1, i.e., implying that Eq. 2 is the exact solution of Eq. 1. It is achieved only for X = 0.5  [confirm with Fig. 7 (a)] and Courant number C = 1 [confirm with Fig. 7 (b)]. This behavior is due to the following reasons: (1) for X = 0.5, the numerical scheme (Eq. 2) is fully centered (Fig. 1) and, therefore, the second-order error vanishes; and (2) for C = 1, the third-order error also vanishes (Ponce et al., 1979 b).

2. For any X-C pair other than X = 0.5 and C = 1, the numerical solution of Eq. 2 is only an approximation of the analytical solution of Eq. 1. The accuracy of the approximation is a function of the following parameters: (a) spatial resolution L/Δx; (b) Courant number C; and weighting factor X of the numerical scheme (Fig. 1).

3. As expected, convergence in both amplitude and phase improves with an increase in spatial resolution. In Figures 2 to 8, the low resolution (L/Δx = 20) is depicted in red color, while the high resolution is shown in black. In particular, note the marked departure of the red-colored curves from the X-axis (R = 1), indicating a substantial lack of convergence for the low resolution case.

4. Convergence in both amplitude and phase degrades with an increase in Courant number as C increases above 1, particularly as C  ≫ 1; for example, for C = 4.

5. Convergence in amplitude improves substantially with an increase of X in the range 0 ≤ X < 0.5. At X = 0.5, amplitude convergence is perfect for all Courant numbers [see Fig. 7 (a)].

6. Convergence in phase improves slightly with an increase of X in the range 0 ≤ X < 0.5. At X = 0.5, phase convergence is perfect only for C = 1 [see Fig. 7 (b)].

7. At low spatial resolutions, refer to red curves in Figs. 2 (b) to 8 (b), convergence in phase (as depicted by R₂) degrades significantly as Courant numbers decrease below C = 1, leading to R₂ > 1, i.e., numerical instability [see Figs. 5 (b), 6 (b), y 7(b)].

8. For X > 0.5, as shown in Figs. 8(a) and 8(b), i.e., outside of the feasible range of X in Muskingum-Cunge routing, the amplitude convergence ratio R₁ becomes positive (amplification) and, accordingly, stability and convergence degrade in an explosive manner.

We conclude that the Muskingum-Cunge numerical model is a good representation of the analytical prototype, provided:

1. The spatial resolution L/Δx is sufficiently high, preferably of the order of L/Δx ≅ 100.

2. The Courant number C is sufficiently close to 1.

3. The weighting factor X is high enough in the range 0.0 ≤ X < 0.5; preferably between 0.3 and 0.5.

Based on the analysis presented herein, the recommended parameter values for reasonably adequate convergence are the following: (a) spatial resolution L/Δx ≅ 100; (b) Courant number C ≅ 1; and (c) weighting factor X between 0.3 and 0.5.

7.  APPLICATIONS

In practice, the convergence of the Muskingum-Cunge method depends on the values of Courant and cell Reynolds numbers, which in turn depend on the hydraulics of the flow and the chosen grid size (Δx and Δt) (Ponce, 2014 b). A set of amplitude and phase convergence ratios may be calculated in terms of hydraulic and hydrologic variables.

To this end, define the time-of-rise t_r of the flood wave as the time it takes it to reach the peak flow. Assuming a sinusoidal perturbation as a first approximation, the time base, or wave period T is:

T =  2 tr

(8)

The wave celerity c is:

c =  β u

(9)

in which u is the mean flow velocity.

In turbulent open-channel flow, the value of β varies from β = 5/3 for a hydraulically wide channel, down to β = 1 for an inherently stable channel (Ponce and Porras, 1995; Ponce, 2014 c). For a rectangular channel, β may be approximated as 5/3, for a triangular channel as 4/3, and for a hydraulically wide natural channel, a reasonably good approximation is β = 1.6.

The wavelength L of the perturbation is:

L =  c T

(10)

The unit-width discharge q is:

q =  u d

(11)

in which d is the mean flow depth. The spatial resolution is L/Δx.

The Courant number (Eq. 4), repeated here for convenience with a new equation number, is:

Δt       C =  c  ______                Δx

(12)

The cell Reynolds number is (Ponce, 2014e):

q       D =   ___________              So c Δx

(13)

By definition, the Muskingum-Cunge weighting factor X is:

1           X =  ____  (1 - D)           2

(14)

For a given case, Eqs. 6 and 7, with support from Eqs. 8 to 14, enable the calculation of the convergence ratios of the Muskingum-Cunge method. Two cases are explained here: (1) theoretical, and (2) practical.

The calculation of the theoretical Muskingum-Cunge convergence ratios may be accomplished using the calculator Online Muskingum-Cunge Convergence Ratios. Three examples are discussed here: (A) perfect convergence, (B) good convergence, and (C) poor convergence. For Case A, assume X = 0.5; L /Δx = 100; and C = 1. The results, shown in Fig. 9, are: R₁ = 1, and R₂ = 1, confirming the perfect convergence of the Muskingum-Cunge method for the case of X = 0.5, L /Δx = 100, and C = 1.

Fig. 9  Muskingum-Cunge convergence ratios for X = 0.5, L /Δx = 100, and C = 1:  R₁ = 1; R₂ = 1.

For Case B, assume X = 0.35; L /Δx = 60; and C = 1.1. The results, shown in Fig. 10, are: R₁ = 0.998195; R₂ = 0.999563, confirming the good convergence of the Muskingum-Cunge method for the case of X = 0.35, L /Δx = 60, and C = 1.1.

Fig. 10  Muskingum-Cunge convergence ratios for X = 0.35, L /Δx = 60, and C = 1.1:  R₁ = 0.998185; R₂ = 0.999563.

For Case C, assume X = 0.05; L /Δx = 20; and C = 2.5. The results, shown in Fig. 11, are: R₁ = 0.908286; R₂ = 0.944956, confirming the poor convergence (i.e., R ratios much smaller than 1) of the Muskingum-Cunge method for the case of X = 0.05, L /Δx = 20, and C = 2.5.

Fig. 11  Muskingum-Cunge convergence ratios for X = 0.05, L /Δx = 20, and C = 2.5:
 R₁ = 0.908286; R₂ = 0.944956.

The calculation of the practical Muskingum-Cunge convergence ratios may be accomplished using the calculator Online Muskingum-Cunge Convergence Ratios Practical. Two examples are discussed here: (D) small stream, and (E) large stream. For Example D, assume the following data set:

1. Mean velocity u = 1.75 m/s;

2. Mean flow depth d = 2.3 m;

3. Mean channel slope S = 0.0013;

4. Rating exponent β = 1.6;

5. Flood wave time-of-rise t_r = 24 hr;

6. Reach length Δx = 4800 m; and

7. Time interval Δt = 0.5 hr.

The application of Online Muskingum-Cunge Convergence Ratios Practical, shown in Fig. 12, leads to: R₁ = 0.99953; R₂ = 0.999915. In this case, the spatial resolution is L/Δx = 100.8; the Courant number is C = 1.05; and the weighting factor is X = 0.385. All three parameters are confirmed to lie within the recommended ranges for adequate convergence.

Fig. 12  Muskingum-Cunge convergence ratios for Example D:  R₁ = 0.99953; R₂ = 0.999915.

1. Mean velocity u = 0.5 m/s;

2. Mean flow depth d = 6.2 m;

3. Mean channel slope S = 0.018;

4. Rating exponent β = 1.65;

5. Flood wave time-of-rise t_r = 36 hr;

6. Reach length Δx = 2500 m; and

7. Time interval Δt = 1.5 hr.

The application of Online Muskingum-Cunge Convergence Ratios Practical, shown in Fig. 13, leads to: R₁ = 0.9996; R₂ = 0.999014. In this case, the spatial resolution is L/Δx = 85.536; the Courant number is C = 1.728; and the weighting factor is X = 0.458. All three parameters are confirmed to lie within the recommended ranges for adequate convergence.

Fig. 13  Muskingum-Cunge convergence ratios for Example E:  R₁ = 0.9996; R₂ = 0.999014.

8.  CONCLUDING REMARKS

Herein we endeavor to clarify some of the concepts elaborated in this article, with aim to encourage its wider use.

On the Courant number. Numerical convergence, which may be taken as numerical accuracy, depends on spatial resolution L/Δx and temporal resolution T/Δt, both of these being dimensionless expressions of grid size. In computacional hydraulics, the Courant number C is defined as the ratio of physical celerity to grid "celerity": C = c / (Δx /Δt) = c (Δt /Δx) = (L/T) (Δt /Δx) = (L/Δx) / (T/Δt). In other words, the Courant number is indeed the ratio of spatial to temporal resolution. Only two of the three variables: (1) spatial resolution, (2) temporal resolution, and (3) Courant number, are required to fully define the analysis. As stated in Section 2 of this article, we use spatial resolution and Courant number, as shown in Figs. 2 to 8.

On the limitations of the Muskingum method. The classical Muskingum method of flood routing (Chow, 1959) solves the kinematic wave equation numerically, thus, introducing a certain undefined amount ot numerical diffusion. In practice, the latter manifests itself as an appreciable amount of attenuation, or diffusion, of the flood wave being calculated. There is no relation of this introduced numerical diffusion with the true physical diffusion, if any, of the problem under consideration [Note: Recall that the true solution of a kinematic wave does not allow for any physical diffusion]. Therefore, the classical Muskingum method fails to solve the flood propagation problem in a meaningful way. The actual solution will be heavily dependent on the chosen grid size, which has not been explicitly related to the physical problem at hand.

On the strength of the Muskingum-Cunge method. Unlike the Muskingum method, the Muskingum-Cunge method (M-C method) is physically and numerically based. The weighting factor X is calculated by Eq. 14, which in turn is based on physical and numerical variables (Eq. 13). This remarkable feature of the M-C method is unique in the field of computational hydraulics, enabling it to become grid independent through a wide range of values of spatial resolution. Once the spatial grid interval Δx is selected or determined a priori, the weighting factor X is fixed by Eq. 14. In this way, the values of physical and numerical diffusion have been matched, and this procedure renders the method and its result grid independent!

On the meaning of the amplitude and phase portraits. Since the kinematic wave equation (Eq. 1) exhibits no physical diffusion (i.e., damping), Eq. 6 (amplitude portrait) is a measure of the numerical diffusion (i.e., a value less than 1) or, alternatively, numerical amplification (i.e., negative diffusion, a value greater than 1), after one time increment Δt (Ponce et al., 1979 d). Correspondingly, a similar conclusion for wave phase (either retardation or acceleration) follows from Eq. 7 (phase portrait). It should not escape our attention that the closeness to 1 of both portrait values shown in Figs. 9 to 13 reflects the fact that the values are applicable per time increment Δt.

9.  SUMMARY

A comprehensive review of the Muskingum-Cunge method's amplitude and phase portraits is accomplished. Expressions for the amplitude convergence ratio R₁ and phase convergence ratio R₂ are expressed as a function of the following basic numerical parameters: (a) spatial resolution L /Δx; (b) Courant number C; and (c) weighting factor X. An online calculator Online Muskingum-Cunge Convergence Ratios is used for convenience.

For practical applications, input may be expresssed in terms of relevant hydraulic variables, such as mean velocity, flow depth, channel slope, rating exponent, and flood wave time-of-rise. For this case, an online calculator Online Muskingum-Cunge Convergence Ratios Practical is developed.

REFERENCES

Chow, V. T. 1959. Open-channel hydraulics. McGraw-Hill, New York.

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.

Hayami, S. 1951. On the propagation of flood waves. Bulletin No. 1, Disaster Prevention Research Institute, Kyoto University, Kyoto, Japan, December.

Koussis, A. 1978. Theoretical estimation of flood routing parameters. Journal of the Hydraulics Division, ASCE, Vol. 104, No. HY1, Proc. Paper 13456, January, 109-115.

Leendertse, J. J. 1967. Aspects of a computational model for long-period water wave propagation. Report RM-5294-PR, The Rand Corporation, Santa Monica, California, May.

McCarthy, G. T. 1938. The unit hydrograph and flood routing. Unpublished manuscript, presented at a conference of the North Atlantic Division, U.S. Army Corps of Engineers, June 24, 1938.

Natural Environment Research Council. 1975. Flood Studies Report, Vol. III: Flood Routing Studies, London, England.

O'Brien, G. G., M. A. Hyman, and S. Kaplan. 1950. A study of the numerical solution of partial differential equations. Journal of Mathematics and Physics, Vol. 29, No. 4, 223-251.

Ponce, V. M. and D. B. Simons. 1977. Shallow wave propagation in open channel flow. Journal of the Hydraulics Division, ASCE, Vol. 103, No. HY12, December, 1461-1476.

Ponce, V. M. and V. Yevjevich. 1978. Muskingum-Cunge method with variable parameters. Journal of the Hydraulics Division, ASCE, Vol. 104, No. HY12, December, 1663-1667.

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, September, 1079-1086.

Ponce, V. M. 1986. Diffusion wave modeling of catchment dynamics. Journal of Hydraulic Engineering, ASCE, Vol. 112, No. 8, August, 716-727.

Ponce, V. M. and P. J. Porras. 1995. Effect of cross-sectional shape on free-surface instability. Journal of Hydraulic Engineering, ASCE, Vol. 121, No. 4, April, 376-380.

Ponce, V. M. 2014. Fundamentals of open-channel hydraulics. Online text.

Ponce, V. M. and B. Vuppalapati. 2016. Muskingum-Cunge amplitude and phase portraits with online computation. Online article.