( ƒj n + 1 + ƒj + 1n + 1 )                     ( ƒj n + ƒj + 1n ) ƒ (x,t )  ≈  θ ______________________   +  (1 - θ) ________________                                      2                                                2

(7)

∂                         ( ƒj + 1n + 1 - ƒj n + 1 )                       ( ƒj + 1n - ƒj n ) ____ ƒ (x,t )  ≈  θ ______________________   +  (1 - θ) ________________  ∂x                                   ∆ x                                               ∆ x

(8)

∂                        1                         _____ ƒ (x,t )  ≈  _____ [ ( ƒj n + 1 + ƒj + 1n + 1 ) - ( ƒj n + ƒj + 1n ) ]  ∂ t                      2∆t

(9)

in which ƒ = a dependent variable, and θ = weighting factor of the implicit scheme.

The application of Eqs. 7-9 to Eqs. 3 and 4 in the discrete x - t plane yields a set of simultaneous nonlinear algebraic equations. In order to solve for the values of the dependent variables at the advanced time level without resorting to an iterative process, a linearization scheme is devised. This scheme allows for the conversion of a system of nonlinear equations into a system of linear equations that can be solved directly, without iteration. The linearization scheme generates additional errors over and above those already imposed by the discrete solution. However, the size and growth of these errors can be effectively controlled by remaining within the limits of the computational method, i.e avoiding sharp changes in the values of the dependent variables. In terms of cost effectiveness, the linearized solution scheme compares favorably with respect to the nonlinear solution.

The linearization scheme is of the following form:

∆ƒj  ( ƒj n + 1 )δ ≈ ( ƒj n )δ + δ( ƒj n )δ-1∆ƒj       (valid for _____ < < 1 )                                                                                    ƒj n

(10)

in which δ = a real number.

The solution of the system of simultaneous linear equations is accomplished by the well-known double sweep algorithm (6).

4.  BOUNDARY CONDITIONS

At the upstream boundary, either the change in stage or bed elevation with time may be specified. Alternatively, the flow depth can be specified. This enables the modeling of a nonequilibrium sediment inflow hydrograph by relating the change in flow depth to the change in sediment inflow. For varying water and sediment discharge, the upstream boundary condition is formulated as follows (7):

qs1 n                 q1 n + 1  ∆d1T = ∆d1S + ∆d1Q = d1 [ (__________ )1/m(__________ ) - 1 ]                                                   qs1 n + 1                q1 n

(11)

in which ∆d_1T = change in flow depth at the upstream boundary; ∆d_1S = change in flow depth at the upstream boundary due to nonequilibrium sediment inflow; and ∆d_1Q = change in flow depth at the upstream boundary due to nonequilibrium water inflow. For time-invariant water discharge, Eq. 11 reduces to:

qs1 n                 ∆d1S = d1 [ (_________ )1/m - 1 ]                         qs1 n + 1

(12)

Equations 11 and 12 are based on the bed material transport equation (Eq. 6) and are strictly of a nonlinear type. A comparison between this formulation and another based on a linearization at the boundary [see Cunge and Perdreau (3)] shows the former to be inherently accurate, while the latter can result in gross errors at the boundary.

The downstream boundary condition is generally given as a stage-discharge relation. For the constant water discharge case, the downstream stage is usually kept constant.

5.  STABILITY AND CONVERGENCE

Every numerical model must satisfy certain stability and convergence requirements. Stability refers to the ability of the numerical scheme to inhibit the generation of error growth. Convergence is a measure of the smallness of the errors in the numerical solution due to improper discretization. In implicit schemes, stability is controlled by roundoff errors (9); convergence, by discretization errors. Stability is impaired if the discretization is made finer while the converse is true from the convergence viewpoint. Both properties are in conflict in numerical schemes: the more stable a scheme, the less convergent, and vice versa.

The first writer, et al. (10) have presented a comprehensive analysis of the numerical properties of implicit bed transient models. According to this study, the model behavior is controlled by three numerical and two physical parameters. The numerical parameters are: (1) The spatial resolution L /∆x ; (2) the temporal resolution T /∆t ; and (3) the weighting factor of the scheme θ. The spatial and temporal resolutions refer to the model discreteness in time and space, L and T being the characteristic length and duration of the transient under study, respectively. The ratio of the spatial to temporal resolution in a dimensionless value referred to as Courant number (1). For a certain L , T , and ∆x, the Courant number is a measure of the size of the time step ∆t.

The two physical parameters are: (1) the Froude number of the equilibrium flow F_o, defined as F_o = u_o / √ gd_o ; and (2) the dimensionless wave number σ_*, defined as σ_* = (2π / L)(d_o / S_o), in which u_o and d_o = equilibrium flow velocity and depth, respectively, and S_o = bed slope.

6.  NUMERICAL EXPERIMENTS

The numerical experiments are aimed at assessing the overall behavior of the model, measured by its ability to simulate hypothetical transients yielding physically realistic solutions. In addition, numerical experimentation is used to verify the theoretical analysis of stability and convergence.

Formation of Sand Waves. Sand waves are defined as transient sediment accumulations in an alluvial river bed. Various phenomena such as floods, landslides, and dam failure are known to cause the formation of sand waves. Once the transient has been formed, it slowly migrates downstream under subcritical flow conditions.

The mathematical model presented herein effectively simulates sand wave formation. A nonequilibrium bed material discharge hydrograph (in excess or defect of the equilibrium transport rate) creates a positive (aggrading) or negative (degrading) sand wave as it enters the channel reach. The process of sand wave formation is subject to a numerical constraint: the bed transient Courant number C must be equal to or greater than one. Otherwise, parasite oscillations are introduced in the model and eventually render the solution meaningless. Therefore, in modeling sand wave formation, it is of utmost importance to have prior knowledge of the celerity of the bed perturbation (4). With the celerity known at least approximately, the bed Courant number is checked to make sure that the numerical constraint is satisfied: C ≥ 1.

Migration of Sand Waves. Sand waves migrate downstream in subcritical flow, and are subject to dissipation. The amount of dissipation observed in the numerical model is usually the combination of physical and numerical dispersion. The amount of numerical dispersion is an indication of the convergence of the solution; the solution will be more convergent as the numerical dispersion is minimized. In general, this can be accomplished by using as high a spatial resolution L/∆x as practicable, a bed Courant number near 1, and a value of the weighting factor in the range 0.5 ≤ θ ≤ 0.75. However, stability may be impaired if the removal of all numerical dispersion is sought; therefore, in practice, a small amount of numerical dispersion may be necessary.

Erosional and Depositional Transients. The mathematical model may be used to study various erosional and depositional transients. The case of a sediment trapping dam is shown in Fig. 2. This transient phenomenon is simulated by increasing the downstream stage to reflect the presence of the dam. This causes a depth adjustment at upstream sections, flow velocities are reduced, and deposition occurs. A sand wave is deposited at the upstream end of the reach and it migrates downstream slowly, eventually filling the channel.

Fig. 2   Deposition upstream of sediment trapping dam.

Other examples of transient formation and propagation such as the process of degradation below dams (Fig. 3) and the behavior of dredge cuts and borrow holes (Fig. 4) can similarly be studied using the mathematical model.

Fig. 3   Degradation below dam.

Fig. 4   Migration of borrow hole.

Modeling Variable Water Discharge. When modeling a time-varying water discharge inflow at the upstream boundary, the hydrograph is discretized into a number of intervals. Each interval is made to correspond to one or more bed transient time steps. The simulation of the time-varying water discharge is accomplished by interfacing a backwater computation between bed transient time steps.

Sensitivity to Numerical Parameters. A series of test runs was designed with the purpose of assessing the sensitivity of the model to the numerical parameters. The problem is that of tracking the translation and attenuation of a sinusoidal sand wave as it travels downstream in a unit width channel. The physical parameters were set at F_o = 0.167, and σ_* = 132. A sinusoidal sand wave of 1 mile (1.6 km) length and an amplitude of 1 ft (0.305 m) was specified as initial condition. The numerical parameters were varied within the following ranges:

0.5 ≤ θ ≤ 1.0

(13)

L            10 ≤ _____  ≤ 100           ∆x

(14)

1 ≤ C ≤ 8

(15)

Figures 5(a) to 5(c) show the results of the numerical testing. Fig. 5(a) shows that the larger the value of θ, the more stable the model, i.e., the larger the numerical dispersion introduced in the model. A value of θ = 0.6 shows a weak pattern of instability. Fig. 5(b) shows that low values of L/∆x introduce numerical dispersion in the same way as do high values of θ. Fig. 5(c) shows that high values of C also introduce numerical dispersion in a manner similar to a high value of θ. In general, stability improves for increasing θ, increasing C, and decreasing L/∆x ; the converse is true from the convergence viewpoint. Therefore, a tradeoff between the conflicting demands of stability and convergence is necessary. An appropriate trade-off value for θ appears to be θ = 0.75 (Ref. 8). Values of L/∆x ≥ 40 are necessary in order to minimize the numerical dispersion introduced by taking θ = 0.75. The appropriate value of the Courant number C is a function of the spatial resolution L/∆x and the temporal resolution T/∆t. The second writer (7) has shown that values of T/∆t > 20 are necessary in order to remain within reasonable error bounds. Therefore, for a cenain T/∆t , the larger the L/∆x , the larger the Courant number that can be successfully used in the model.

Fig. 5   Effect of:(a) weighting factor θ on numerical solution;(b) spatial resolution L/∆x on numerical solution; (c) bed courant number C on numerical solution (sand bed profiles shown at intervals of 96 days).

Recommendations on Stability and Convergergence. Values of 0.5 ≤ θ ≤ 0.6 may result in a weakly stable or unstable solution. As a general rule, stability requires that θ > 0.6. A value of θ = 0.75 is recommended, provided it is used in conjunction with a sufficiently high spatial resolution (L/∆x ≥ 40) and as low a Courant number as possible. A value of C ≥ 1 is necessary to avoid oscillations due to ill-posing of the upstream boundary condition.

The second writer (7) has derived an approximate convergence criterion based on actual test runs. He defined a convergence parameter C_p as follows:

L                         ____                 ∆x   Cp  =  _______                       C θ

(16)

High values of C_p result in good convergence; conversely, low values of C_p are an indication of poor convergence. A lower limit of C_p may be dictated by accuracy considerations, while an upper limit may be necessary from the economic viewpoint. Eq. 16 can also be expressed as follows:

T                          ____                  ∆t   Cp  =  ________                          θ

(17)

Figure 6 shows the results of the convergence analysis. The degree of convergence is defined as the ratio of numerical to analytical celerity, for numerical solutions with varying T/∆t and θ. Figure 6 underscores the importance of the temporal resolution T/∆t as the relevant convergence parameter while depicting the secondary role of the weighting factor θ. A recommended range for the temporal resolution is as follows: 20 ≤ T/∆t ≤ 100. Midrange values will represent a tradeoff between accuracy and cost constraints.

Fig. 6   Degree of convergence as function of temporal resolution T/∆t and weighting factor θ.

7.  SUMMARY AND CONCLUSIONS

A mathematical model is developed to simulate bed transient formation and propagation in alluvial channels. The model is of the known-discharge type, i.e., within one time step, the water discharge is considered steady. This technique allows for the implicit numerical modeling of bed transient propagation using longer time steps than would otherwise be possible using the conventional sequential routing method.

Instability due to ill-posing of the sediment transport function is circumvented by using a power relation between bed material discharge and mean flow velocity. Errors in mass conservation due to a linearization at the boundary are eliminated by using a nonlinear formulation of the upstream boundary condition. The model can take account of time-varying water discharge by interfacing a backwater computation between one or a number of time steps. This feature allows for considerable flexibility in modeling long-term water surface and bed level transients.

Stability and convergence are analyzed, and criteria are established to ensure a successful model operation. In particular, convergence is shown to be primarily a function of the temporal resolution T/∆t, with the weighting factor θ playing a secondary but nonetheless important role.

APPENDIX I.   REFERENCES

1. Abbott, M. A., A. Damsgaard, and G. S. Rodenbuis. 1973. "System 21, 'Jupiter,' A Design System for Two-Dimensional Nearly-Horisontal Flows," Journal of Hydraulle Research, Vol. 11, No. 1, 1-28.

2. Chen, Y. H. 1973. "Mathematical Modeling of Water and Sediment Routing in Natural Channels," thesis presented to Colorado State University, at Fort Collins, Colo., in partial fulfillment of the requirements for the degree of Doctor of Philosophy.

3. Cunge, J. A., and N. Pedreau. 1973. "Mobile Bed Fluvial Mathematical Models," La Beadle Blanche, Grenoble, France, Vol. 28, No. 7, 561-580.

4. De Vries, M. 1973. "River Bed Variations-Aggradation and Degradation," Publication No. 107, Delft Hydraulics Laboratory, Delft, The Netherlands, Apr.

5. "HEC-6, Scour and Deposition in River and Reservoirs". 1973. Report 723-G2-L2470 , U.S. Army Corps of Engineers, Hydrologic Engineering Center, Davis, Calif., Mar.

6. Liggett, J. A., and J. A. Cunge. 1975. "Numerical Methods of Solution of the Unsteady Flow Equations," Unsteady Flow In Open Channels, K. Mahmood and V. Yevjevich, eds., Water Resources Publications, Fort Collins, Colo., 89-182.

7. Lopez-Garcia, J., "Mathematical Modeling of Alluvial Bed Transients," thesis presented to Colorado State University, at Fort Collins, Colo., in 1977, in partial fulfillment of the requirements for the degree of Master of Science.

8. Mahmood, K., and V. M. Ponce. 1976. "Mathematical Modeling of Sedimentation Transients in Sand-Bed Channels," Report CER75-76-KM-VMP28, Engineering Research Crater, Colorado State University, Fort collins. Colo., Apr.

9. O'Brien, G. G., M. A., Hyman, and S. Kaplan. 1951. "A Study of the Numerical Solution of Partial Differential Equations," Journal of Mathematics and Physics, Vol. 29, No. 4, 223-251.

10. Ponce, V. M., H. Indlekofer, and D. B. Simons. 1979. "The Convergence of Implicit Bed Transient Models," Journal of the Hydraulics Division, ASCE, Vol. 105.