 

MuskingumCunge amplitude and phase portraits with online computation
Victor M. Ponce and Bavya Vuppalapati
30 April 2016


ABSTRACT
A comprehensive review of the amplitude and phase portraits of the MuskingumCunge method of flood routing is accomplished.
Expressions for the
amplitude and phase convergence ratios R_{1} and R_{2}, respectively, are developed as a function of:
(a) spatial resolution L/Δx; (b) Courant number C; and (c) weighting factor X.
It is concluded that the MuskingumCunge routing model is a good representation of the physical prototype, provided:
(1) the spatial resolution is sufficiently high,
(2) the Courant number is close to 1, and
(3) the weighting factor is high enough, better if in the range 0.3 ≤ X ≤ 0.5.
Two online calculators of the convergence ratios, one theoretical based on dimensionless parameters, and the other practical, based on actual hydraulic/hydrologic variables,
are developed and tested.

1. INTRODUCTION
The concept of amplitude and phase portraits in computational hydraulics
is due to Leendertse (1967). The amplitude portrait is a plot of the ratio
of numerical to analytical wave amplitudes, i.e.,
a convergence ratio with regard to wave amplitude,
as a function of relevant variables such as the spatial
resolution or temporal resolution, and Courant number. Likewise,
the phase portrait is a plot of the ratio
of numerical to analytical wave celerities (a convergence ratio with regard to wave phase)
as a function of relevant numerical parameters. Amplitude and phase portraits
for the MuskingumCunge method of flood routing were originally developed by Cunge (1969), and later
expanded for the generalized convection (advection) problem by Ponce et al. (1979a).
The objective of amplitude and phase portraits are to
assess the stability and convergence
properties of a given numerical scheme. This helps determine
the range of discrete parameters (spatial resolution, temporal resolution, and their corresponding ratio)
for which the scheme may be shown to be stable and convergent. Following O'Brien et al. (1950),
stability is related to roundoff errors, while convergence is related to
truncation (i.e., discretization) errors.
In this study, we revisit the amplitude and phase portraits of the MuskingumCunge method,
clarify their development, and further explain the details of the computation,
with aim to develop an online calculator.
2. MUSKINGUMCUNGE METHOD
In a seminal paper,
Cunge (1969) elucidated the kinematic wave basis of the Muskingum method of flood routing (Chow, 1959),
leading to the MuskingumCunge method
(Natural Environment Research Council, 1975; Koussis, 1978; Ponce and Yevjevich, 1978).
Cunge's procedure allowed for the calculation of the weighting factor X in terms of physical
and numerical parameters, instead of relying solely on the
discharge, as proposed by the original Muskingum method (McCarthy, 1938).
This revolutionary concept, based on the matching of the physical diffusivity of
Hayami (1951) with the numerical
diffusivity of Cunge (1969), enables the MuskingumCunge method to remain grid independent through a wide
range of relevant parameters (spatial resolution, Courant number, and weighting factor). It is worth noting that this fundamental property
of the
MuskingumCunge method is not shared by other flood routing methods
based on the kinematic wave (Ponce, 1986).
In his original paper (op. cit.), Cunge expressed the amplitude and phase portraits in terms of the temporal resolution
and a type of Courant number.
Ponce et al. (1979b) expanded Cunge's approach to the
entire set of numerical schemes that can be construed to solve the general convection (advection)
equation. Ponce et al. (op. cit.) chose to present their findings for the amplitude and phase
convergence ratios in terms of the spatial resolution (L/Δx) and the
Courant number C, defined as the ratio of physical celerity c (i.e., the kinematic wave celerity) to numerical celerity
(the grid ratio Δx/Δt).
In this study, we review the findings of Ponce et al. (op. cit.) and apply them
specifically to the MuskingumCunge method.
3. DERIVATION OF EQUATIONS
The kinematic wave equation is (Cunge, 1969; Ponce, 2014a):
∂Q ∂Q
^{____} + c ^{____} = 0
∂t ∂x
 (1) 
in which Q = discharge, c = kinematic wave celerity, and x and t are space and time variables, respectively.
The discretization of Eq. 1 on the xt plane (Fig. 1), assuming linearity (c = constant) leads to
(Ponce, 2014b):
X (Q_{ j }^{n+1}  Q_{ j}^{ n} ) + (1  X) (Q_{ j+1}^{n+1}  Q_{ j+1}^{ n} ) (Q_{ j+1 }^{n}  Q_{ j}^{ n} ) + (Q_{ j+1}^{n+1}  Q_{
j}^{ n+1} )
^{________________________________________________} + c ^{_______________________________________} = 0
Δt 2 Δx
 (2) 
Fig. 1 Spacetime discretization of the kinematic wave equation following the MuskingumCunge method.
Simplifying Eq. 2:
X (Q_{ j }^{n+1}  Q_{ j}^{ n} ) + (1  X) (Q_{ j+1}^{n+1}  Q_{ j+1}^{ n} ) + 0.5 C [ (Q_{ j+1 }^{n}  Q_{ j}^{ n} ) + (Q_{ j+1}^{n+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 wave length, and β = (2 π / T )
is the wave period (Ponce and Simons, 1977).
Replacing Eq. 5 into the various components of Eq. 3:
Q _{j}^{ n} = Q_{*} e ^{i}^{ [σj Δx  βn Δt ]}
 (6) 
Q _{j+1}^{n} = Q_{*} e ^{i}^{ [σ (j + 1) Δx  βn Δt ]}
 (7) 
Q _{j} ^{n+1} = Q_{*} e ^{i}^{ [σj Δx  β(n + 1) Δt ]}
 (8) 
Q _{j+1}^{n+1} = Q_{*} e ^{i}^{ [σj Δx  βn Δt ]}
 (9) 
Equations 6 to 9 can be respectively expressed as follows:
e ^{i} ^{σj Δx}
Q_{j} ^{n} = Q_{*} ^{___________}
e ^{i βn Δt}
 (10) 
e ^{i}^{ σj Δx} e ^{i}^{ σΔx}
Q_{j+1}^{n} = Q_{*} ^{_________________}
e ^{i βn Δt}
 (11) 
e ^{i}^{ σj Δx}
Q_{j} ^{n+1} = Q_{*} ^{_________________}
e ^{i βn Δt} e ^{i βΔt}
 (12) 
e ^{i}^{ σj Δx} e ^{i}^{ σ Δx}
Q_{j+1}^{n+1} = Q_{*} ^{___________________}
e ^{i βn Δt} e ^{i βΔt}
 (13) 
The first term of Eq. 3 is:
e ^{i}^{σj Δx} e ^{i}^{σj Δx}
X (Q_{ j }^{n+1}  Q_{ j}^{ n} ) = Q_{*} X [ ^{_________________}  ^{___________} ]
e ^{iβn Δt} e ^{iβΔt} e ^{iβnΔt}
 (14a) 
e ^{i}^{σj Δx}  e ^{i}^{σj Δx} e ^{iβΔt}
X (Q_{ j }^{n+1}  Q_{ j}^{ n} ) = Q_{*} X [ ^{__________________________} ]
e ^{iβn Δt} e ^{iβΔt}
 (14b) 
The second term of Eq. 3 is:
e ^{i}^{σj Δx} e ^{i}^{σΔx} e ^{i}^{σj Δx} e ^{i}^{σΔx}
(1X ) (Q_{ j+1 }^{n+1}  Q_{ j+1}^{ n} ) = Q_{*} (1  X) [ ^{__________________}  ^{_________________} ]
e ^{iβn Δt} e ^{iβΔt} e ^{iβnΔt}
 (15a) 
e ^{i}^{σj Δx} e ^{i}^{σΔx}  e ^{iβΔt} e ^{i}^{σj Δx} e ^{i}^{σΔx}
(1X ) (Q_{ j+1 }^{n+1}  Q_{ j+1}^{ n} ) = Q_{*} (1  X) [ ^{________________________________________} ]
e ^{iβn Δt} e ^{iβΔt}
 (15b) 
The third term of Eq. 3 is:
(C / 2 ) (Q_{ j+1}^{n+1}  Q_{ j}^{ n+1} + Q_{ j+1}^{n}  Q_{j}^{ n} ) =
 
C Q_{*} e ^{i}^{σj Δx} e ^{i}^{σΔx}  e ^{iσjΔx} + e ^{iβΔt} e ^{i}^{σj Δx} e ^{i}^{σΔx}  e ^{iβΔt} e^{ i}^{σj Δx}
^{________} [ ^{___________________________________________________________} ]
2 e ^{iβn Δt} e ^{iβΔt}
 (16) 
Substituting Eqs. 14b, 15b, and 16 into Eq. 3,
dividing by Q_{*} e ^{i}^{σj Δx}, and eliminating the common denominator, leads to:
2X (1 e ^{iβΔt} ) + 2(1 X )(e ^{i}^{σΔx}  e ^{i}^{σΔx} e ^{iβΔt} ) + C (e ^{i}^{σΔx}  1 + e ^{iβΔt} e ^{i}^{σΔx}  e ^{iβΔt} ) = 0
 (17a) 
e ^{iβΔt} [ 2X  2(1X ) e ^{i}^{σΔx} + C e ^{i}^{σΔx}  C ] + 2X + 2(1X ) e ^{i}^{σΔx} + C e ^{i}^{σΔx}  C = 0
 (17b) 
e ^{iβΔt} [ 2X + 2(1X )e ^{i}^{σΔx}  C e ^{i}^{σΔx} + C ] = 2X + 2(1X ) e ^{i}^{σΔx} + C e ^{i}^{σΔx}  C
 (17c) 
2[X + (1X ) e ^{i}^{σΔx}] + C [e ^{i}^{σΔx}  1]
e ^{iβΔt} = ^{____________________________________}
2[X + (1X ) e ^{i}^{σΔx}]  C [e ^{i}^{σΔx}  1]
 (17d) 
Rearranging Eq. 17d and taking the reciprocal value of the lefthand side:
2[(1X) e ^{i}^{σΔx} + X]  C [e ^{i}^{σΔx}  1]
e ^{iβΔt} = ^{____________________________________}
2[(1X) e ^{i}^{σΔx} + X] + C [e ^{i}^{σΔx}  1]
 (18) 
A trigonometric identity is (Spiegel, 1971):
e ^{i}^{σΔx} = p + iq
 (19) 
in which
p = cos(σΔx) and q = sin(σΔx).
The spatial resolution is L /Δx.
The quantity
σΔx = (2π / L) Δx = 2π / (L/Δx)
is the number of grid points
per wavelength. Substituting Eq. 19 into Eq. 18:
2[(1X)(p + iq) + X]  C [p + iq  1]
e ^{iβΔt} = ^{_____________________________________}
2[(1X)(p + iq) + X] + C [p + iq  1]
 (20a) 
{p [(1X)  C/2] + X + C/2} + iq [(1X)  C/2]
e ^{iβΔt} = ^{______________________________________________}
{p [(1X) + C/2] + X  C/2} + iq [(1X) + C/2]
 (20b) 
Defining for the sake of simplicity:
S = R  C = X  (C/2)
 (22) 
Substituting Eqs. 21 and 22 into Eq. 20b:
[p (1  R) + R ] + iq (1  R )
e ^{iβΔt} = ^{______________________________}
[p (1  S) + S ] + iq (1  S )
 (23) 
Multiplying the righthand side of Eq. 23 by the conjugate of the denominator:
[p (1  R) + R ] + iq (1  R ) [p (1  S) + S ]  iq (1  S )
e ^{iβΔt} = ^{______________________________} _{•} ^{______________________________}
[p (1  S) + S ] + iq (1  S ) [p (1  S) + S ]  iq (1  S )
 (24) 
{[p(1R) + R ] [p(1S) + S ] + q ^{2}(1R) (1S)} + iq {(1R) [p(1S) + S ]  (1S) [p(1R) + R ]}
e ^{iβΔt} = ^{____________________________________________________________________________________________}
[ p (1S )
+ S ] ^{2} + [ q (1S )] ^{2}
 (25) 
Now examine the components of Eq. 25. First consider the real part [RPN]
of the numerator of the righthand side of Eq. 25:
[RPN] = [p (1R) + R ] [p (1S) + S ] + q ^{2}(1R ) (1S )
 (26) 
Operating on Eq. 26, and after some algebraic manipulation:
[RPN] = 1  (1p) [R + S  2RS]
 (27) 
Given that S = R  C, it follows that
R = S + C. Thus, in Eq. 27:
[RPN] = 1  [ 2(1  p)(1  S)S ]  [ (1p)(1  2S)C ]
 (28) 
Consider the imaginary part [IPN] of the numerator of Eq. 25:
[IPN] = q {(1R ) [p (1S ) + S ]  (1S )[p (1R ) + R ]}
 (29) 
Operating on Eq. 29, and after some algebraic manipulation:
Replacing R = S + C from Eq. 22 results in:
Consider the denominator [DN] of Eq. 25:
[DN] = [p (1S ) + S ]^{2} + [q(1S )]^{2}
 (32) 
After some algebraic manipulation:
[DN] = 1  2 (1  p)(1  S) S
 (33) 
Defining for the sake of simplicity:
P = 1  2 (1  p)(1  S) S
 (34) 
Therefore, Eq. 25 can be written as follows:
P  M  iN
e ^{iβΔt} = ^{______________}
P
 (37a) 
which may be expressed as follows:
P  M N
e ^{iβΔt} = ^{_________}  i ^{______}
P P
 (37b) 
In the lefthand side of Eq. 37b, separating β into its real and imaginary parts:
e ^{iβΔt} = e ^{i (βR + i βI) Δt} = e ^{i βR Δt} e ^{βI Δt}
 (38a) 
e ^{iβΔt} = e ^{βI Δt} e ^{i βR Δt} = a e^{iφ}
 (38b) 
in which a = e ^{βI Δt}, and e^{iφ} = e ^{i (βR Δt)}
In the righthand side of Eq. 38a, taking the natural logarithm:
ln (e ^{i βR Δt} e ^{βI Δt}) = ln e ^{βI Δt}  i β_{R} Δt
 (39) 
Equation 37b can be expressed in the form z = x + iy, in which:
z = e ^{iβΔt} = a e^{iφ}
 (40) 
P  M
x = ^{________}
P
 (41) 
The solution of Eq. 40 is (Spiegel, 1971):
ln(a e^{iφ}) = ln (x^{2} + y^{2})^{1/2} + i tan^{1}(y / x)
 (43) 
Substituting Eqs. 38b, 41 and 42 into Eq. 43:
{(P  M)^{2} + N^{ 2}}^{1/2} N
ln e ^{βI Δt}  i β_{R} Δt = ln [ ^{____________________} ]  i tan^{1}( ^{ _________} )
P P  M
 (44) 
Separating Eq. 44 into its real and imaginary components:
[(P  M)^{2} + N^{ 2}]^{1/2}
e^{ βI Δt} = ^{____________________}
P
 (45) 
N
tan (β_{R} Δt) = ^{_________}
P  M
 (46) 
Since Eq. 1 does not allow for wave damping (e^{δ} = 1, where δ = logarithmic decrement
(Ponce and Simons, 1977), it is concluded
that Eq. 45 is a measure of the numerical damping ( sign)
or amplification (+ sign) of the numerical solution of Eq. 3 after a time increment Δt.
4. CONVERGENCE RATIOS
Following Ponce et al. (1979c), a convergence ratio with regard to wave amplitude is:
Therefore:
[(P  M)^{2} + N^{ 2}]^{1/2}
R_{1} = ^{____________________}
P
 (48) 
The celerity of the numerical solution is:
β_{R}
c_{n} = ^{______}
σ
 (49) 
Using Eqs. 46 and 49, a convergence ratio with regard to wave phase is defined as follows:
c_{n} 1 N
R_{2} = ^{_____} = ^{________} tan^{1} ( ^{_________} )
c σ c Δt P  M
 (50) 
in which c_{n} is the wave celerity of the numerical solution, c is the kinematic wave celerity (in Eq. 1), and σ is the wavenumber.
From Eq. 4: C Δx = c Δ t. Therefore:
1 N
R_{2} = ^{_________} tan^{1} ( ^{________} )
C σΔx P  M
 (51) 
In summary, given only: 
Weighting factor X,
 Spatial resolution L/Δx, from which: σΔx
= (2π /L) Δx = 2 π / (L /Δx),
and
 Courant number C,
Equations 48 and 51 enable the calculation of the MuskingumCunge amplitude and phase
convergence ratios, respectively.
The calculation of the convergence ratios is accomplished
using the calculator
Online MuskingumCunge Convergence Ratios.
5. AMPLITUDE AND PHASE PORTRAITS
The plotting of convergence ratios R_{1} and
R_{2} as a function of Courant number C and spatial resolution L/Δx
leads to a set of amplitude and phase portraits, a pair for each value of
weighting factor X.
Figures 2 to 8 show the amplitude and phase portraits of the MuskingumCunge method (Eq. 3),
for Courant numbers varying in the range
0.1 ≤ C ≤ 4.0, spatial resolution 20 ≤ L/Δx ≤ 100, at intervals of L/Δx = 20,
and weighting factor 0.0 ≤ X ≤ 0.6, at intervals of X = 0.1 (7 × 2 = 14 graphs).
Fig. 2 (a) Amplitude portrait for X = 0.
Fig. 2 (b) Phase portrait for X = 0.
Fig. 3 (a) Amplitude portrait for X = 0.1.
Fig. 3 (b) Phase portrait for X = 0.1.
Fig. 4 (a) Amplitude portrait for X = 0.2.
Fig. 4 (b) Phase portrait for X = 0.2.
Fig. 5 (a) Amplitude portrait for X = 0.3.
Fig. 5 (b) Phase portrait for X = 0.3.
Fig. 6 (a) Amplitude portrait for X = 0.4.
Fig. 6 (b) Phase portrait for X = 0.4.
Fig. 7 (a) Amplitude portrait for X = 0.5.
Fig. 7 (b) Phase portrait for X = 0.5.
Fig. 8 (a) Amplitude portrait for X = 0.6.
Fig. 8 (b) Phase portrait for X = 0.6.
6. ANALYSIS
The examination of Figs. 2 to 8 leads to the following conclusions:
Perfect convergence is achieved for X = 0.5 and Courant number C = 1.
Convergence in both amplitude and phase improves with an increase in spatial resolution (L/Δx).
Convergence in both amplitude and phase degrades with an increase in Courant number C, particularly for C > 1.
Convergence in amplitude improves with an increase in X in the range 0 ≤ X ≤ 0.5.
Convergence in phase improves slightly with an increase in X in the range 0 ≤ X ≤ 0.5.
At low spatial resolutions [refer to red curves in Figs. 2 (b) to 8 (b)], convergence in phase degrades significantly with Courant numbers in the range C < 1, leading to instability.
For X > 0.5, the amplitude convergence
ratio R_{1} becomes positive (amplification), and convergence and stability degrade accordingly.
We conclude that the MuskingumCunge model is a good representation of the physical prototype, provided:
The spatial resolution is sufficiently high.
The Courant number is close to 1.
The weighting factor is high enough, in the range 0.0 ≤ X ≤ 0.5.
Based on the analysis accomplished herein, the recommended parameter values for adequate convergence are: (a) spatial resolution L/Δx ≥ 100; (b) Courant number C ≅ 1; and (c) weighting factor 0.3 ≤ X ≤ 0.5.
7. PRACTICAL APPLICATIONS
In practical applications, the convergence of the MuskingumCunge method depends on the values of Courant and cell Reynolds numbers, which in turn,
depend on the hydraulics of the flow and the grid size (Δx and Δt) (Ponce, 2014c).
A set of amplitude and phase convergence ratios may be calculated
in terms of hydraulic and hydrologic variables.
To this end, define the timeofrise t_{r} of the flood wave as the time it takes for the flood wave to reach the peak flow.
Assuming a sinusoidal perturbation as a first approximation, the time base, or wave period T is:
The wave celerity c is:
in which u is the mean flow velocity.
In turbulent
openchannel 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, 2014d).
The wavelength L of the perturbation is:
The unitwidth discharge q is:
in which d is the mean flow depth. The spatial resolution is L/Δx.
The Courant number (Eq. 4), repeated here with a new equation number, is:
Δt
C = c ^{______}
Δx
 (56) 
The cell Reynolds number is defined as follows ((Ponce, 2014e):
q
D = ^{___________}
S_{o} c Δx
 (57) 
By definition, the MuskingumCunge weighting factor X is:
1
X = ^{____} (1  D)
2
 (58) 
For a given case, Eqs. 52 to 58, together with Eqs. 48 and 51, enable the calculation
of the practical convergence ratios of the MuskingumCunge method,
based on actual flow data, i.e., on the hydraulics of the flow.
The calculation of the MuskingumCunge practical convergence ratios may be accomplished
using the calculator
Online MuskingumCunge Convergence Ratios Practical.
For example, 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; rating exponent β = 1.6, (5) flood wave timeofrise t_{r} = 24 hr;
(6) reach length Δx = 4,800 m, and (7) time interval Δt = 0.5 hr. The application of
Online MuskingumCunge Convergence Ratios Practical leads to: R_{1} = 0.99953; R_{2} = 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 lie within the recommended ranges for adequate convergence.
8. CONCLUSIONS
A comprehensive review of the MuskingumCunge method's amplitude and phase portraits is accomplished.
Expressions for the
amplitude convergence ratio R_{1} and phase convergence ratio R_{2} are developed as a function of:
(a) spatial resolution L/Δx; (b) Courant number C; and (c) weighting factor X.
An online calculator
Online MuskingumCunge Convergence Ratios is developed and tested using a realistic data set.
In practical applications, the convergence ratios may be expressed alternatively in terms of
relevant hydraulic variables (mean velocity, flow depth, channel slope, rating exponent,
and flood wave timeofrise). For this case, an online calculator
Online MuskingumCunge Convergence Ratios Practical is also developed and tested.
REFERENCES
Chow, V. T. 1959. Openchannel hydraulics. McGrawHill, 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, 205230.
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, 109115.
Leendertse, J. J. 1967. Aspects of a computational model for longperiod water wave propagation. Report RM5294PR, 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, 223251.
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, 14611476.
Ponce, V. M., and V. Yevjevich. 1978.
MuskingumCunge method with variable parameters.
Journal of the Hydraulics Division, ASCE, Vol. 104, No. HY12, December, 16631667.
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, 10791086.
Ponce, V. M. 1986. Diffusion wave modeling of catchment dynamics.
Journal of Hydraulic Engineering, ASCE, Vol. 112, No. 8, August, 716727.
Ponce, V. M., and P. J. Porras. 1995. Effect of crosssectional shape on freesurface instability. Journal of Hydraulic Engineering, ASCE, Vol. 121, No. 4,
April, 376380.
Ponce, V. M. 2014. Fundamentals of openchannel hydraulics.
Online text.
Spiegel, M. R. 1971. Advanced mathematics for engineers and scientists. Schaum's Outline Series, McGrawHill Inc., New York.
