Modeling Alluvial Channel
Bed Transients

Victor M. Ponce

Online version 2018

[Original version 1979]


Flow in an alluvial channel is of an unsteady nonuniform type. Accordingly, the water and sediment discharge vary in time and space. In an alluvial channel reach, equilibrium flows are defined as those in which the sediment inflow and outflow are balanced over a sufficiently long period of time. The resistance, bed material transport, and morphologic relations are usually based on equilibrium flow conditions. Deviations from the equilibrium condition cause the formation of bed surface perturbations that travel downstream in subcritical flow and are subject to dissipation. The regime associated with these perturbations, or the transitional regime between two equilibrium states, is referred to as a bed transient (8). Among the causes of bed transient formation are: (1) A variation in the upstream sediment supply to a channel reach; (2) a change in bed material transport capacity triggered by a change in hydraulic conditions; and (3) the effect of dredge cuts and borrow holes on hydraulic and bed material transport conditions.

The study of bed transient propagation can be approached from an analytical or numerical standpoint. The analytical approach (4) can be used to study certain types of transients, but the simplifying assumptions are too strenuous for the general case. On the other hand, numerical models can usually be built using a smaller number of assumptions and are more versatile than analytical solutions in the specification of boundary conditions. Therefore, a numerical model is regarded as a more effective tool in the study of bed transient propagation.

Relevant questions of stability, consistency, and convergence arise every time a numerical modeling job is attempted. Nevertheless, the experience with this type of modeling in the last two decades indicates that these difficulties can usually be brought under control.

A numerical model involves the solution of one or more governing partial differential equations by the use of finite difference techniques. In the bed transient problem, there are three equations: (1) water continuity; (2) sediment continuity; and (3) water motion. These equations constitute a set of quasilinear partial differential equations. Their numerical solution can be attempted in various ways. In the sequential routing method (2), the equations of water continuity and motion are solved simultaneously as a first stage of the calculations, and the changes in bed elevation are obtained by using the sediment continuity equation in a second stage. This type of formulation achieves good detail in the simulation of water surface transients but often results in an inefficient computation of the bed transients.

A practical alternative to the sequential routing method for the cases in which the bed transients are of primary interest is given by the known-discharge method. The rationale for this method of solution lies in the fact that the bed transients propagate much slower pace than the water surface transients. Therefore, within one time step, a steady water discharge can be assumed to hold. It should be noted that the assumption of steady water discharge in bed transient models is analogous to the assumption of rigid boundary in river flood routing models (4).

The steady water discharge assumption dictates that the water continuity equation be described by a simple algebraic relation. Therefore, only two differential equations remain: (1) the sediment continuity equation; and (2) the equation of motion. Furthermore, the neglect of the local acceleration term of the equation of motion follows directly from the steady water discharge assumption. As with the flood routing models based on a kinematic formulation, the unsteady nature of the problem is preserved through the rate-or-change term in the continuity equation.

The numerical solution of these two equations can be accomplished in two ways: (1) an uncoupled mode (5) may be used by first solving for the water surface profile and then adjusting the bed elevation using the sediment continuity equation; and (2) a coupled mode (3,8) may be used by simultaneously solving the sediment continuity equation and the equation of motion. Uncoupled models are based on explicit schemes, while coupled models lend themselves to computation by implicit schemes. Therefore, the coupled models are more flexible, since the maximum size of the time step is not governed by the strict stability criterion usually associated with explicit schemes, but rather, its size is solely based on accuracy considerations.

A known-discharge coupled bed transient model is presented in this paper. The model represents an improvement over earlier models of the same type in that instability due to ill-posing of the bed material transport function (6) is circumvented, and errors in mass conservation due to improper formulation of the boundary condition (3) are avoided by use of a strictly nonlinear formulation. Since the model simulates the slow bed transient phenomena in alluvial channels well in the subcritical range (Froude numbers less than 0.6), it can accept rather long time steps, and thus, its cost effectiveness is assured. In addition, the model can take account of time-averaged 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 transients. The price paid, of course, is the loss of detail in the simulation of the short-term water surface transients.

Three limitations of the mathematical model are apparent at the outset. The one-dimensional formulation does not permit the description of meandering. Also, since no provision is made for the routing of sediment by size fractions, no account can be made of armoring phenomena. In addition, the transport function does not effectively account for the condition of incipient motion; therefore, the model would work best when the flow conditions are such that the Froude number is greater than 0.1.

The remainder of the present paper treats the following subjects in detail: (1) governing equations; (2) finite difference scheme; (3) boundary conditions; (4) stability and convergence; and (5) numerical experiments.


The governing equations, expressed per unit of channel width, are continuity of sediment

  ∂qs         ∂                                       ∂z
_____  + ____ (Csd )  + (1 - p ) ρs g ____  =  qsL
  ∂x          ∂t                                       ∂t


and equation of motion

  1       ∂u            ∂         u 2             ∂d           ∂z 
____  _____  +   ____  (_____)  +   _____  +   _____  +  Sƒ    =  0
  g       ∂t            ∂x        2g             ∂x           ∂x


in which qs = bed material discharge, in pounds per second per foot; Cs = spatial bed material concentration, in pounds per cubic foot; d = flow depth, in feet; p = bed porosity; ρs = density of sediment particles; g = acceleration of gravity; z = bed elevation with reference to an arbitrary datum; qsL = lateral sediment inflow, in pounds per second per square foot; u = mean velocity; and Sƒ = friction slope.

According to Mahmood and the first writer (8), in river flow the temporal concentration variation term in Eq. 1 is very small when compared with the remaining terms; therefore, it is neglected here. In addition, the local acceleration term in Eq. 1 is neglected on the grounds that the steady discharge assumption does not warrant its inclusion (4). Making q = ud, in which q = steady water discharge, Eqs. 1 and 2 reduce to

  ∂qs                           ∂z
_____  + (1 - p ) ρs g ____  =  qsL
  ∂x                             ∂t



     q 2       ∂                    ∂d            ∂z 
   _____   ___ (d -2)  +   _____  +   _____  +  Sƒ    =  0
     2g      ∂x                   ∂x            ∂x                        


Two supplementary equations are necessary in order to solve the system of Eqs. 3 and 4: (1) a flow resistance formula; and (2) a bed material transport formula. No closed form equation for either is available, although empirical equations have been used in the past with some success. These equations have been derived for equilibrium flow conditions; their use for nonequilibrium flow conditions appears unwarranted at first. Nevertheless, the insight gained in the process of modeling far outweighs the inability of the model to properly account for all the physical details.

The resistance and transport functions used herein are:

            ƒr q 2           
Sƒ  =  _______
             gd 3     


                         Cg q m           
qs  = Cg u m =  ________
                            d m     

in which ƒr and Cg = empirical resistance and transport coefficients, respectively; and m = exponent of mean velocity. The dimensionless resistance coefficient ƒr is immediately recognized as ƒr = g / C 2, in which C = the Chézy coefficient.

For a wide channel, ƒr = f /8, in which f = the Darcy-Weisbach friction factor.


The four-point implicit scheme of Preissmann (6) is used here. The scheme is expressed by the following operators (refer to Fig. 1);

Dimensionless relative wave celerity vs dimensionless wavenumber

Fig. 1   Definition sketch for four-point implicit finite difference scheme.

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


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


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


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 n + 1 )δ ≈ ( ƒj n )δ + δ( ƒj n )δ-1∆ƒj (valid for _____ < < 1 )
                                                                             ƒj n


in which δ = a real number.

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


At the upstream boundary, either the change in stage or bed elevation with time can 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 (7)

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


in which ∆d1T = change in flow depth at the upstream boundary; ∆d1S = change in flow depth at the upstream boundary due to nonequilibrium sediment inflow; and ∆d1Q = 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           


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.


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 Fo, defined as Fo = uo /  gdo ; and (2) the dimensionless wave number σ*, defined as σ* = (2π / L)(do / So), in which uo and do = equilibrium flow velocity and depth, respectively, and So = bed slope.


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

Dimensionless relative wave celerity vs dimensionless wavenumber

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.

Dimensionless relative wave celerity vs dimensionless wavenumber

Fig. 3   Degradation below dam.
Dimensionless relative wave celerity vs dimensionless wavenumber

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 Fo = 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


10 ≤ _____  ≤ 100


1 ≤ C ≤ 8


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.

Dimensionless relative wave celerity vs dimensionless wavenumber

Fig. 5   Effect of:(a)weightins 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 Cp as

Cp  =  _______         
             C θ


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

Cp  =  ________         


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 θ. Fig. 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.

Dimensionless relative wave celerity vs dimensionless wavenumber

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


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.


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

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

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

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

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

  6. Liggett, J. A., and Cunge, J. A., "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., 1975, pp. 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 Ponce., V. M., "Mathematical Modeling of Sedimentation Transients in Sand-Bed Channels," Report CER75-76-KM-VMP28, Engineering Research Crater, Colorado State University, Fort collins. Colo., Apr., 1976.

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

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


The following symbols are used in this paper:

C = bed transient Courant number, Chézy coefficient;

Cg = bed material transport coefficient;

Cp = convergence parameter, Eqs. 16 and 17;

Cs = spatial bed material concentration;

d = flow depth;

do = equilibrium flow depth;

Fo = Fronde number of equilibrium flow;

ƒ = a dependent variable;

ƒr = friction factor;

f = Darcy-Weisbach friction factor;

g = acceleration of gravity;

j = space discretization index;

L = characteristic length of transient;

L / ∆x = spatial resolution;

m = exponent in bed material transport function;

n = time discretization index;

p = bed porosity;

q = steady water discharge;

qs = bed material discharge;

qsL = Lateral sediment inflow;

Sf = friction slope;

So = bed slope;

T = characteristic duration of transient;

T / ∆t = temporal resolution;

t = time variable;

u = mean velocity;

x = space variable;

z = bed elevation referred to arbitrary datum;

d1Q = change in flow depth at upstream boundary due to nonequilibritun water inflow;

d1S = change in flow depths at upstream due to nonequilibrium sediment inflow;

d1T = total change in flow depth at upstream boundary;

∆ ƒ = change in variable ƒ ;

t = time step;

x = space step;

δ = a real number;

θ = weighting factor of implicit scheme;

ρs = density of sediment particles; and

σ* = dimensionless wave number of transient.

181121 12:00

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