Fig. 1.  Open-book catchment schematization (20).

A source of uncertainty in this type of modeling is the adequacy of the frictional representation. Flow over the planes is usually of such small depths that laminar flow may prevail. In contrast, flow in the channel is likely to be of a turbulent nature. Yet in certain cases, especially near the downstream end of the planes, the flow may be in the transitional, or mixed laminar-turbulent, regime. For routine applications, friction coefficients are generally estimated in such a way that they account for laminar flow in the overland flow planes (6).

The kinematic wave equations are a common choice in the simulation of flow over catchments. The equations are well established and generally lead to physically realistic solutions without the computational and other complexities that usually plague full dynamic wave formulations. To derive the kinematic wave equation, the usual statement of conservation of mass in a control volume is coupled with a simplified form of conservation of momentum which accounts only for frictional and body forces. The statement of conservation of mass is the following:

∂Q          ∂A            _____   + ______ = qL   ∂x           ∂t

(1)

in which Q = flow rate, in cfs; A = flow area, in sq ft; q_L = lateral inflow, in cfs/ft.

The simplified momentum statement is S_o = S_f , in which S_o = bed slope and S_f = friction slope. This in fact a statement of uniform slow, which leads to an alternate way of expressing the simplied momentum statement:

Q = αAβ

(2)

in which α and β are coefficient and exponent, respectively. The values of α and β contain information on frictional and cross-sectional characteristics [for instance, β = 3 for a wide channel with laminar friction; β = 5/3 for a wide channel governed by Manning friction; β = 3/2 for a wide channel with Chézy friction; and β = 4/3 for a triangular channel with Manning friction (11)].

The speed of propagation of kinematic waves is obtained from Eq. 2:

∂Q            Q            c = _____   = β ____ =  βv         ∂A             A

(3)

in which c = celerity of kinematic waves; and v = mean flow velocity.

The kinematic wave equation is obtained by multiplying Eqs. 1 and 3 and making use of the chain rule to obtain:

∂Q             ∂Q            _____  + c  ______ = cqL   ∂t               ∂x

(4)

Equation 4 contains convection terms on the left-hand side and a source term on the right-hand side. Significantly, Eq. 4 does not contain a diffusion term. This equation may be solved by analytical or numerical means. A typical class of numerical solutions is that in which the values of ∂Q/∂t, ∂Q/∂x, c and q_L are expressed in terms of some or all of four discrete adjacent values of Q, c and q_L in space and time (Fig. 2). Alternatively, the value of c could be kept constant, leading to the linear kinematic wave solution.

Fig. 2.  Spatial and temporal discretization of kinematic wave equation.

A numerical solution of Eq. 4 following the aforementioned procedure is not without its pitfalls. Unless the spatial and temporal derivatives are perfectly centered in the rectangular grid, a certain amount of diffusion will be generated in the numerical solution (Fig. 2). This numerical diffusion can be eliminated by centering the spatial and temporal derivatives, but a small amount of numerical dispersion will still remain (16). The latter is most often simply a nuisance, but in certain cases it may be of such magnitude as to invalidate the whole solution. Furthermore, given the fact that there usually is a small amount of diffusion in the physical problem, the complete elimination of numerical diffusion is generally ill-advised.

Fully off-centered schemes remain popular among practitioners (6, 7, 10). Ponce et. al. (16) have shown that of the three possible fully off-centered schemes (see Fig. 3), schemes I and II are conditionally stable and convergent to the analytical solution as the Courant number (defined as the ratio of physical celerity c to grid celerity Δs/Δt) approaches 1. Scheme I, stable for Courant numbers less than or equal to 1, amounts to "upwind differencing." Scheme II, stable for Courant numbers greater than or equal to 1, is exactly the opposite image of scheme I. Scheme III is unconditionally stable, but it is nonconvergent for any value of Courant number (16).

Fig. 3.  Three fully off-centered numerical schemes for solving that kinematic wave equation (16).

Schemes I and II are complementary, with versions of them in wide use, for instance in the kinematic wave routing option of HEC-1 and HEC-HMS (6). Scheme III (and related versions), although nonconvergent, has also been favored because of its features of unconditional stability (7, 10). In any given practical application, local Courant numbers are likely to vary widely and be different from 1. If such is the case, none of these three schemes would solve the kinematic wave equation, certainly not scheme III given its demonstrated nonconvergence. In theory, for Courant numbers close to 1, as the discretization is made exceedingly small, schemes I and II may asymptotically approach the analytical solution of the kinematic wave equation. With today's powerful computers, this is certainly a possibility. However, the question still remains whether the complete elimination of numerical diffusion is the desirable strategy for this type of problem.

In summary, the hydrologic modeler is faced with the following dilemma: If he or she seeks to eliminate numerical diffusion, at whatever cost, the solution may be purely kinematic and therefore lack the natural diffusion property inherent in physically realistic flows. If, on the other hand, the modeler uses a fully off-centered scheme such as I, II or III, the end result will be an uncontrolled amount of numerical diffusion only vaguely resembling the actual diffusion process.

The way out of this difficulty is to judiciously "match" physical diffusion with numerical diffusion, a procedure which has been successfully used in stream channel routing and other similar applications (9, 14, 18). A clear advantage of this approach is that it renders the numerical solution essentially independent of grid size. In other words, convergence is satisfied at a much coarser grid resolution than with either schemes I or II. Given the natural variation in Courant numbers likely to be present in practical applications, a method that converges for a wide range of grid sizes is certainly one to be reckoned with.

3.  DIFFUSION WAVE FORMULATION

The soundness of the method of "matched diffusivity" is illustrated here by comparing the results of two solutions for catchment dynamics. The first is an application of scheme III to the kinematic wave equation, Eq. 4, with the associated "uncontrolled" numerical diffusion effect. The second is an application of the matched diffusivity concept. For simplicity, both solutions are formulated in a linear mode, with celerity and diffusivity parameters kept constant. The routing of linear waves vis-à-vis nonlinear waves has indicated that for small catchments the nonlinear features are generally small and can be disregarded when comparing the numerical accuracy of alternative schemes (14).

The application of scheme III to Eq. 4 in a linear mode leads to the following routing equation:

n + 1                 n + 1            n               Q j + 1     =   C0Qj       +  C2Qj + 1  + C0QL

(5)

in which C₀ = C/(1 + C); C₂ = 1/(1 + C); and Q_L is the average lateral inflow rate. C is the Courant number of the computational cell, C = c (Δt /Δs), where Δs is either Δx or Δy.

For overland flow computations, Q_L is the product of rainfall access times the area of the overland flow cell. For channel routing, Q_L, is the time-averaged lateral contribution integrated along the length of a cell. With an estimate of average wave celerity for plane and channel, together with the proper specification of initial and boundary conditions, Eq. 5 enables the calculation of catchment response for a given rainfall excess excitation.

The derivation of the matched diffusivity method parallels that of the Muskingum-Cunge method with the addition of lateral inflow (3, 14, 15). The routing equations are the following:

n + 1                n + 1          n            n               Q j + 1     =   C0Qj     +  C1Qj  +  C2Qj + 1   + C3QL

(6)

in which C₀ = (-1 + C + D)/(1 + C + D); C₁ = (1 + C - D)/(1 + C + D); C₂ = (1 - C + D)/(1 + C + D); C₃ = (2C)/(1 + C + D); and D is the cell Reynolds number defined as the ratio of physical diffusivity q/(2S_o) to numerical diffusivity (cΔs)/2, as follows (14):

q                  D  =  __________             So c Δs

(7)

in which q = unit-width discharge (q = vd; d = flow depth); S_o = bottom slope; c = kinematic wave celerity, and Δs = spatial grid size.

Note that there are major differences between Eqs. 5 and 6 with regard to their physical data requirements. Equation 5 requires only values of c for both plane and channel (or alternatively: mean velocity v and a suitable β value, see Eq. 3). On the other hand, Eq. 6 requires values of c, q and S_o for both plane and channel (or alternatively: β, v, d, and S_o ).

Since the physical diffusivity is v = q /(2S_o ), it follows that unlike Eq. 6, Eq. 5 does not contain any physical information on diffusion. Not surprisingly, however, the application of Eq. 5 does show a diffusion-like behavior. This diffusion, unrelated to the physical problem, is directly traceable to the grid size and is a function of it.

4.  DIFFUSION-CUM-DYNAMIC WAVE FORMULATION

The preceding analysis has shown the advantages of using diffusion waves in lieu of kinematic waves in the modeling of catchment response. Aside from the assumption of linearity, the diffusion wave can be further improved to include the dynamics of the wave phenomena. Following Dooge (4), the diffusivity of the diffusion-cum-dynamic wave formulation for the case of turbulent Chézy friction in sufficiently wide channels is:

q                F2     v  =  _____ ( 1 - _____)           2So             4

(8)

in which F is the Froude number, defined as F = v/(gd)^1/2, and g is the gravitational acceleration. Equation 8 has the important property that the physical diffusivity v is zero for F = 2, a threshold of neutral stability at which both kinematic and dynamic waves propagate at the same speed (for turbulent Chézy friction in sufficiently wide channels) (11, 12).

Equation 8 can be generalized for any frictional formulation and cross-sectional shape by recalling that the neutral stability Froude number is equal to 1/(β - 1). Therefore, the physical diffusivity expressed in terms of β is:

q                     v  =  _____ [1 - (β - 1)2F2]           2So

(9)

Given Eq. 9, a diffusion-cum-dynamic wave model can be formulated based on Eq. 6, with Eq. 7 modified to account for the Froude-number and β dependence of the physical diffusivity, Eq. 9. In effect, for the diffusion-cum-dynamic wave, Eq. 7 is modified to:

q                     D  =  ___________  [1 - (β - 1)2F2]            So c Δs

(10)

With the cell Reynolds number defined as in Eq. 10, Eq. 6 ought to represent as complete a description of the wave dynamics as is possible within the context of a linear formulation. However, the usefulness of Eq. 10 as compared to Eq. 7 will depend on whether the wave phenomena being simulated can be properly construed as a dynamic wave (13). Experience indicates that for a large class of problems of practical interest in hydrology, the wave scale may well be in the diffusion (or kinematic) range.

5.  ILLUSTRATIVE EXAMPLE

The methodologies described in the previous section are compared here by applying them to a typical problem of catchment dynamics. Following widespread usage, the scheme of Eq. 5 is referred to as "kinematic wave," although its appreciable and uncontrolled numerical diffusion is recognized. The "matched diffusivity" method, Eqs. 6 and 7, is herein referred to as "diffusion wave." The diffusion-cum-dynamic wave, Eqs. 6 and 10, although properly within the family of diffusion models, is herein referred to as "dynamic wave" (see Tables 2 and 3).

The example consists of a catchment conceptualized as an open book with two planes adjacent to one channel (Fig. 1). The planes are assumed to be impermeable, with rainfall excess (net rainfall) equal to total rainfall. The inflow to the planes is the rainfall excess, and runoff is by overland laminar flow in a direction perpendicular to the channel alignment. Inflow to the channel is by lateral contribution from the planes, and the outflow from the channel is the catchment response.

The dimensions are: plane length (in the direction of plane flow) = 120 ft (36.58 m); plane width = 240 ft (73.15 m); channel length = 240 ft (73.15 m). The rainfall is 3 in/hr (76.2 mm/h), which, when multiplied times the catchment area, leads to a flow concentration of 4 cfs (0.113 m³/s) at the outlet (assuming translation only and neglecting diffusion). The bottom slopes are the following: slope of planes S_op (in the x-direction) = 0.01; slope of channel S_oc = 0.01.

Laminar flow in the planes and turbulent flow in the channel are calculated following established hydraulic practice (2). Average flow conditions in the planes are: discharge per unit width q_p= 0.004 cfs/ft (0.00037 m³/s/m); mean velocity v_p = 0.5 fps (0.15 m/s); mean flow depth d_p = 0.008 ft (0.0024 m); rating curve exponent (laminar flow in wide channels) β_p = 3; and mean wave celerity c_p = 1.5 fps (0.46 m/s). Average flow conditions in the channel are discharge per unit width q_c = 1 cfs/ft (0.092 m³/s/m); mean velocity v_c = 3 fps (0.91 m/s); mean flow depth d_c = 0.333 ft (0.1 m); rating curve exponent (turbulent flow in triangular channel) β_c = 4/3; and mean wave celerity c_c = 4 fps (1.22 m/s).

These calculations lead to a time-of-concentration t_c = (120 ft/1.5 fps) + (240 ft/4 fps) = 140 s. Therefore, a rainfall duration of 3 min is chosen to ensure the attainment of equilibrium at the catchment outlet. The total rainfall volume is 0.15 in. (3.81 mm), or 720 cu ft (20.39 m³).

5.  RESULTS

The three methods were tested with regard to their behavior under various levels of grid resolution. Five resolution levels were chosen, varying from very coarse to very fine, as shown in Table 1. The results of the simulations are shown in Tables 2 and 3 expressed in terms of: (1) peak flow; and (2) time-to-peak.

Figure 4 shows the calculated hydrographs for: (1) kinematic method; and (2) diffusion method. As expected, the results of the kinematic method varied substantially with the grid resolution. With a very coarse resolution (A), the peak flow was 2.3507 cfs (0.0666 m³/s), while for a very fine resolution (E), the peak flow was 3.9641 of (0.1122 m³/s). The results of the diffusion method showed much better convergence properties than the kinematic method. With a very coarse resolution (A), the peak flow was 3.9336 cfs (0.1113 m³/s), while for a very fine resolution (E), it was 3.9971 cfs (0.1132 m³/s). It is concluded that the diffusion method provides grid independence for a wide range of resolution levels, while the kinematic method does not. The demonstrated grid independence of the diffusion method is clearly a significant asset of this formulation.

TABLE 2. Calculated peak flow by three methods (in cubic feet per second).
Resolution level (1)	Kinematic (2)	Diffusion (3)	Dynamic (4)
A	2.3507	3.9336	3.9390
B	3.0721	3.9797	3.9834
C	3.5921	3.9932	3.9954
D	3.8612	3.9964	3.9979
E	3.9641	3.9971	3.9985
Note:  See text for method description. 1 cfs = 0.0283 m3/s.

TABLE 3. Calculated time to peak by three methods (in seconds).
Resolution level (1)	Kinematic (2)	Diffusion (3)	Dynamic (4)
A	240	180	180
B	210	180	180
C	195	180	180
D	187	180	180
E	184	180	180
Note:  See text for method description.

Fig. 4.  Calculated outflow hydrographs as function of grid resolution (see Table 1 for explanation of curve parameter): (a) kinematic method; (b) diffusion method.

The results of the dynamic method were not plotted in Fig. 4 because of their close resemblance to those of the diffusion method (see Tables 2 and 3). This resemblance underscores the lack of a significant dynamic contribution in the tested example. The dimensionless criterion of Ponce et. al. (13) may be used to prove whether the tested wave is indeed a diffusion wave. According to this criterion, the wave is a diffusion wave if the following inequality is satisfied:

g    1/2                           TSo (_____)    ≥ 30             d

(11)

in which T is the representative time scale of the wave phenomena. In this case a time scale can be assumed to be equal to twice the rainfall duration: T = 360 s. For the channel flow, with S_o = 0.01 and d = 0.333 ft (0.1 m), Eq. 11 leads to a left-hand side having a value of 35.4, satisfying the inequality. This explains why the dynamic wave solution is not very different from the diffusion wave solution for this particular problem.

Mass conservation was monitored by integrating the outflow hydrographs for all runs reported in Tables 2 and 3. All methods and runs were 100% mass conservative. In addition, the diffusion method was used with spatial resolution levels satisfying attendant accuracy criteria (17).

In closing, it should be pointed out that had the discretization been made exceedingly small, the peak flows calculated by the diffusion and dynamic methods (see Tables 2 and 3) would have converged to the peak flow value obtained by concentrating the flow at the outlet 4 cfs (0.1132 m³/s). This is attributed to the fact that translation (i.e., convection) is the predominant physical mechanism for the tested example, with diffusion playing a relatively minor role. Therefore, a diffusion method (and diffusion-cum-dynamic extension) in which the diffusion is small and controlled with physical parameters is certainly to be preferred over a kinematic method in which the diffusion could be excessive and bear no relation to the physics of the problem.

6.  SUMMARY AND CONCLUSIONS

A diffusion wave method for catchment dynamics is formulated. The method has better convergence properties than kinematic wave models which use off-centered derivatives in their numerical formulations. Off-centering the derivatives introduces numerical diffusion in uncontrolled amounts, with the results of the simulation being highly dependent on the chosen grid resolution.

Unlike the off-centered schemes, the diffusion wave scheme is formulated by matching physical and numerical diffusivity. This results in an effective control of numerical diffusion and leads to simulations which are essentially independent of grid size.

The diffusion wave method is further extended to the realm of dynamic waves by including the demonstrated Froude-number dependence of the physical diffusivity. The resulting formulation is believed to represent as complete a description of the wave dynamics as is possible within the framework of diffusion wave theory.

The three methods: (1) kinematic, with uncontrolled numerical diffusion; (2) diffusion, with matched physical and numerical diffusivities; and (3) dynamic, with matched and Froude number-dependent diffusivity, are tested for convergence and overall accuracy. The diffusion method is shown to have much better convergence properties than the kinematic method. In particular, the demonstrated grid independence of the diffusion method is an asset to be reckoned with. At least for the example problem tested here, the solution of the dynamic method is shown not to constitute a substantial improvement over that of the diffusion method.

APPENDIX I.  REFERENCES