of four-point implicit
water wave models

Victor M. Ponce,
Horst Indlekofer,
and Daryl B. Simons

Online version 2019

[Original version 1978]


The numerical solution of the equations of open channel flow (the Saint Venant equations) is an established subject in the hydraulic literature. Implicit numerical schemes are generally preferred over explicit schemes (9) because they allow larger time steps and thus reduce the computational effort. Among the implicit schemes, the four-point scheme, also referred to as the Preissmann scheme (2, 5), has received wide attention. However, the stability and convergence properties of this scheme are as yet not fully understood.

The Saint Venant equations constitute a set of quasilinear partial differential equations. At present, there exists no method for analyzing the stability and convergence of numerical solutions of equations of this type. On the other hand, there are various methods for analyzing the stabiity and convergence of linear partial differential equations (11), among them, notably, the heuristic von Neumann approach (7). Traditionally, a way out of the mathematical difficulty imposed by the quasilinear nature of the Saint Venant equations has been to linearize them, assuming that the behavior of the linear system approximates that of the more complex quasilinear system. In practice, the linear analysis serves, if not as a quantitative, at least as a qualitative description of the nonlinear phenomena.


The essential feature of a numerical model is the replacement of a derivative such as ∂ƒ/∂s by a ratio of finite differences such as Δƒ/Δs. In doing so, the numerical model must satisfy certain requirements of stability and convergence. Stability refers to the ability of the numerical scheme to march in time without generating unbounded error growth. Convergence refers to the ability of the scheme to reproduce the terms of the differential equation with sufficient accuracy. Stability [in the sense of O'Brien, et al. (7)] is governed by round-off errors; convergence is governed by discretization errors. In general, stability is impaired if the discretization is made finer, while the converse is true from the convergence viewpoint.

In practice, stability is a necessary condition for model operation, since an unstable model is of little or no use. In implicit models, stability is usually achieved at the expense of convergence. Therefore, convergence is the focal issue: Once the model is shown to be stable, then it becomes necessary to assess its convergence properties.

The numerical properties of four-point implicit schemes have been studied by Liggett and Cunge (5) and Fread (3). Liggett and Cunge used a simple system of linear partial differential equations to derive conclusions on the numerical properties of the four-point implicit scheme. While they fully realized the limitations of their analysis, they nevertheless pointed out that the Saint Venant equations would be expected to behave in a similar manner. Actual computations have shown, however, that the 'theoretically better" value of the weighting factor θ = 0.5 is impractical from the stability viewpoint, and that in practice a θ ≥ 0.6 may be required (2).

Fread (3) has analyzed the numerical properties of the four-point implicit scheme. He used a linearized version of the Saint Venant equations, but simplified the equation of motion by neglecting the convective inertia and bed slope terms, and the equation of continuity by neglecting the wedge storage term. The complete linear analysis has been made possible by the study of two of the writers (8), who formulated the analytical solution for the propagation of sinusoidal perturbations in open channel flow on the basis of a linearized version of the Saint Venant equations (6).

The present study aims to determine the convergence properties of the four-point implicit numerical models of unsteady open channel flow. The celerity and attenuation of the numerical solution will be compared with those of the analytical solution (8). The von Neumann technique (7,11) will be used in this regard. The conclusions relate to the assessment of convergence as a function of the relevant physical and numerical parameters.


One-dimensional unsteady open channel flow, per unit of channel width, can be expressed by a set of two partial differential equations (4):

Water Continuity Equation:

  ∂d            ∂u           ∂d
____  +  d ____  +  u ____  =  0
  ∂t            ∂x            ∂x

Equation of Motion:

  1       ∂u           u      ∂u            ∂d
____  _____  +   ____  _____  +   _____  +  (Sf  -  So )  =  0
  g       ∂t            g      ∂x            ∂x

in which u = mean velocity; d = flow depth; g = acceleration of gravity; Sƒ = friction slope; So = bottom slope; x = space variable; and t = time variable.

Equations 1 and 2 must satisfy the unperturbed flow (equilibrium or steady uniform flow) for which u = uo and d = do as well as the perturbed flow for which u = uo + u' and d = do + d'. The terms with a superscript represent small perturbations (unsteady flow components) to the steady uniform flow. Introducing the perturbed variables in Eqs. 1 and 2 yields, after linearization (neglecting second-order terms) (6):

  ∂d'             ∂u'              ∂d'
____  +  do _____  +  uo _____  =  0
  ∂t               ∂x               ∂x

  1       ∂u'          uo      ∂u'           ∂d'                    u'         d'
____  _____  +   ____  _____  +   _____  +  So ( 2 ____  -  ____ )  =  0
  g       ∂t            g       ∂x            ∂x                     uo        do


Two of the writers have shown (8) that the analytical solution of Eqs. 3 and 4 is characterized by the propagation celerity and attenuation of sinusoidal perturbations. For the celerity, they chose the dimensionless celerity c* defined as follows:

                     (____ )        
           c            T   
c* =  ____  =  _______
          uo           uo

in which L = wavelength of the sinusoidal perturbation; and T = wave period. For the attenuation, they chose the logarithmic decrement δ defined as follows:

ato+T  =  atoeδ (6)

in which ato and ato+T are the wave amplitudes at the beginning and at the end of a period T.

Both c* and δ are functions of the Froude number of the equilibrium flow Fo and the dimensionless wave number σ*, in which Fo = uo / ( gdo )1/2 and σ* = (2π / L) do /So (Appendix Ι).


The four-point implicit scheme is based on the following formulas (5), Fig. 1:

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

                  θ                                           (1 - θ)             
ƒ (x,t )  =  ___ ( ƒj+1 n + 1 + ƒj n + 1 )  +   _______ ( ƒj+1 n + ƒj n )
                  2                                               2


   ∂                        1                        
_____ ƒ (x,t )  =  _____ ( ƒj+1 n + 1 - ƒj + 1n + ƒj n+1 - ƒj n )
  ∂t                      2∆t                                          


   ∂                        θ                                       1 - θ
_____ ƒ (x,t )  =  _____ ( ƒj+1 n + 1 - ƒj n+1 ) + ______ ( ƒj+1 n - ƒj n )
  ∂x                      ∆x                                        Δx


in which ƒ(x,t) is any dependent variable; Δx and Δt are the space and time increments; and θ is the weighting factor of the implicit scheme. The weighting factor θ is introduced in the function and its space derivative to allow greater flexibility in the actual operation of the model. While the stable range is 0.5 ≤ θ ≤ 1.0, a strongly stable range is 0.6 ≤ θ ≤ 1.0.

The substitution of Eqs. 7-9 into 3 and 4 leads to discretized equations in u' and d'. The solution is postulated in the following exponential form:

____  =  u*  exp {i* x*  -  β* t* )}

____  =  d*  exp {i* x*  -  β* t* )}

in which x* = xSo / do; t* = tuoSo / do and β*R = 2πdo / TuoSo, and u* and d* are dimensionless amplitude functions.

The substitution of Eqs. 10 and 11 into the discretized equations in u' and d' yields:

          Δt*                                          Δt*                               ξ
u* [ i _____ (θξ + 1) tan α ] + d* [ i _____ (θξ + 1) tan α + ____ ]  =  0
          Δx*                                         Δx*                              2

                ξ                   Δt*                               
u* [ Fo2 _____ + i Fo2 ______ (θξ + 1) tan α + Δt* (θξ + 1) ]
                2                  Δx*

              Δt*                               Δt*
+ d* [ i ______ (θξ + 1) tan α - _____ (θξ + 1) ]  = 0
              Δx*                               2


in which i = √-1; ξ = [exp (-iβ* Δt*) - 1 ]; α = σ*Δx* / 2; Δt* = Δt uo So / do and Δx* = Δx So / do.

Equations 12 and 13 constitute a set of two homogeneous equations in u* and d*. For the solution to be nontrivial, the determinant of the coefficient matrix must vanish. Accordingly:

                                    C                                          C           α          ξ2 Fo2
[ (θξ + 1)2 (1 - Fo2) ( ____ ) 2 tan2α + ξ (θξ + 1) ( ____ ) ( ____ ) + _______ ]
                                    c*                                          c*          σ*             4

                                           C                               C             α
+ i tanα [ ξ (θξ + 1) Fo2 ( ____ ) + 3 (θξ + 1)2 ( ____ ) 2 ( ____ ) ] = 0
                                           c*                               c*             σ*


in which C / c* = Δt* / Δx*. The value C is referred to as the Courant number, and it is defined as the ratio of the physical celerity c to the "grid celerity" Δx / Δt. Equation 14 is a complex quadratic in ξ. It can be reduced to:

ξ2 (M2 + iN2) + ξ (M1 + iN1) + (Mo + iNo) = 0 (15)

in which:

Mo = (1 - Fo2)P2 (16)

No = 3PQ (17)

M1 = 2θMo + Q (18)

N1 = Fo2P + 2θNo (19)

M2  =  _____ + θQ + θ2Mo

N2 = θFo2P + θ2No (21)

P  =  ( ____ ) tanα


             C           α
Q  =  ( ____ ) ( ____ )
             c*          σ*


A = Q2 - P2 Fo2 (24)

B = -PQ Fo2 (25)

C = (A2 + B2)1/2 (26)

              A + C           
D  =  ( ________ ) 1/2

              A - C           
E  =  ( ________ ) 1/2

and solving the quadratic Eq. 15, the two roots are:

ξ1,2 = M + iN (29)

in which:

              -M2 ( M1 ± E ) + N2 ( -N1 ± D )
M  =  ( _______________________________ )
                         2( M22 + N22 )


              M2 ( -N1 ± D ) + N2 ( M1 ± E )
N  =  ( _______________________________ )
                         2( M22 + N22 )

Since ξ = [exp ( -i β* Δt* ) - 1], it follows that:

( β*R Δt* ) =  - ________
                         1 + M


exp ( β*I Δt* ) = [( 1+ M )2 + N 2 ]1/2 (33)

The dimensioness celerity of the numerical solution c*n is defined as follows (8):

             β*R               1                        -N
c*n  =  _______  =  ________ tan -1 ( _________ )
               σ*            σ* Δt*                  1 + M

or likewise:

                     1                            -N
c*n  =  _____________ tan -1 ( _________ )
                      C                         1 + M
            2α ( _____ )

The logarithmic decrement of the numerical solution δn is defined as follows (8):

                 β*I            π ln [( 1 + M )2 + N 2]
δn  =  2π ______  =  _______________________
                 β*R                              -N
                                   tan -1 ( _________ )
                                                  1 + M

Since c* is a function of Fo and σ*, and M and N are functions of Fo, σ*, α, C, and θ, it follows that c*n and δn are also functions of Fo, σ*, α, C and θ. Fo and σ* are the physical parameters, Fo being the Froude number and σ* a dimensionless wave number (σ* = 2πdo / LSo). The values α, C, and θ are the numerical parameters, in which α is the spatial resolution (α = π Δx / L); C is the Courant number (C = c Δt / Δx); and θ is the weighting factor of the implicit scheme.


The convergence analysis is based on the following ratios:

R1 = exp ( δn - δ ) (37)

R2  =  ______

The value R1 is the attenuation ratio and R2 is the translation ratio. For R1 > 1, the numerical damping is smaller than the physical damping; for R1 < 1, the numerical damping is greater than the physical damping. For R2 > 1, the numerical translation is greater than the physical translation; and for R2 < 1, the numerical translation is smaller than the physical translation. For R1 = R2 = 1, there is an exact coincidence between analytical and numerical solutions.

The sensitivity of the convergence ratios to the physical and numerical parameters has been studied by varying the parameters within a specified range. For the Froude number Fo, the following range has been studied: 0.01 ≤ Fo ≤ 0.99. This excludes the critical and supercritical regimes. For the wave number σ*, the following range has been studied: 0.01 < σ* < 1,000. This range includes both the kinematic and diffusion waves (σ* ≅ 0.01), and the inertia-pressure waves (σ* ≅ 1,000). The spatial resolution has been varied between π / 100 ≤ α ≤ π / 10. This range includes both a very poor resolution (α ≅ π / 10) and a very high resolution (α ≅ π / 100). The Courant number has been varied in the range 0.25 ≤ C ≤ 4.0. The weighting factor θ has been varied in the range 0 ≤ θ ≤ 1.

The analysis has been performed for the primary wave only (8), by taking the minus sign associated with D and E in Eqs. 30 and 31. Space limitations preclude a complete analysis of both primary and secondary components of the wave propagation.

Fig. 2  Attenuation ratio R1 as function of θ, σ*, Fo and C for: (a) α = π / 10; (b) α = π / 40; (c) α = π / 100.

Figures 2(a)-2(c) portray the attenuation ratio R1 as a function of Fo, σ*, α, C, and θ. The examination of theses figures leads to the following conclusions:

  1. For kinematic / diffusion waves (σ* ≅ 0.1) and inertia-pressure waves (σ* ≅ 1,000), the Froude number has a negligible effect on R1. For dynamic waves (σ* ≅ 10), however, the effect of Fo on R1 is very marked.

  2. For a given α, the attenuation ratios for kinematic/diffusion waves are not unlike those for inertia-pressure waves. In fact, Fig. 2 shows that for a given α, the R1, values for σ* = 0.1 and σ* = 1,000 are indistinguishable.

  3. The higher the spatial resolution (the lower the α value) and the lower the Courant number, the closer is R1 to 1 for all θ.

  4. For kinematic / diffusion and inertia-pressure waves, the convergence of the wave amplitude (R1 → 1) is obtained for values of θ equal to or slightly greater than 0.5. For dynamic waves, a value of θ much greater than 0.5 may be required to avoid an insufficient amount of numerical damping (R1 > 1) which is tied to numerical amplification (instability).

The following general conclusions are obtained from Fig. 2:

  1. The simulation is reasonably good for kinematic / diffusion waves and inertia-pressure waves, provided sufficient care is taken to minimize the amount of numerical damping. Large values of θ (θ approaching 1 from the left) can be used in conjunction with high spatial resolution. For poor spatial resolution, smaller θ values (θ approaching 0.5 from the right) may be required in order to avoid an excessive amount of numerical damping.

  2. For dynamic waves (waves with both resistance and inertia), the accuracy of the simulation is highly dependent on the value of θ. The optimum value of θ is a function of σ* and α, and may be difficult to determine in practice.

  3. Values of θ < 0.5 always cause numerical amplification; values of 0.5 ≤ θ < 1 may cause numerical amplification or attenuation; a value of θ = 1 always causes numerical attenuation (R1 < 1). The existence and amount of numerical attenuation in the range 0.5 ≤ θ < 1 is a function of Fo, σ*, α, and C.

The analysis of the translation ratio R2 shows that this ratio is nearly unity for a wide range of parameters Fo, σ*, α, and C. Predictably, for midrange dimensionless wavenumbers (σ* = 10) and poor spatial resolution, the ratios R2 depart from unity. As an illustration of this trend, Fig. 3 shows the values of R2 for varying Fo, σ*, C, and θ, and α = π / 10.

Fig. 3  Translation ratio R2 as function of θ, σ*, Fo and C for: α = π / 10.


A comprehensive theoretical treatment of the convergence of the four-point implicit numerical model of shallow water waves is presented. The propagation celerity and attenuation factor of the numerical analog are derived, and convergence is tested by establishing the ratio of attenuation and translation given by the numerical and analytical solutions. Convergence is shown to be a function of the complex interaction of five parameters: (1) The Froude number of the equilibrium flow Fo; (2) the dimensionless wave number σ** = (2π / L)do /So); (3) the spatial resolution α (α = π Δx / L); (4) the Courant number C (C = c Δt / Δx); and (5) the weighting factor θ of the implicit scheme.

The convergence analysis is carried out for Froude numbers in the subcritical regime. Furthermore, due to space limitations, only the primary wave convergence is analyzed.

The following conclusions are warranted by this study:

  1. The simulation is reasonably good for small σ* (kinematic / diffusion) waves and for large σ* (inertia-pressure) waves, provided sufficient care is exercised to minimize the amount of numerical damping.

  2. For intermediate waves σ* (dynamic waves, in which both friction and inertia predominate), the accuracy of the simulation is highly dependent on the correct value of the weighting factor θ. In practice, an optimum value of θ that will assure both stability and convergence may be difficult to determine.

  3. Values of θ < 0.5 always cause numerical amplification; values of 0.5 ≤ θ < 1 may cause numerical amplification or attenuation; a value of θ = 1 always causes numerical attenuation.


The findings of this study have been obtained by using a linearized version of the Saint Venant equations. Therefore, emphasis should be placed on the qualitative nature of the results. At present, the performance of a nonlinear model can only be assessed through carefully designed numerical experiments. The linear analysis presented herein provides an appropriate framework for the design of these experiments.


The research subject of this paper was supported in part by the Colorado State University Experiment Station. Horst Indlekofer was awarded a NATO fellowship to make possible his visit to Colorado State University.


The equations are: c* = 1 ± D; δ = -2π (BE)/(|1 ± D |), in which A = (1/Fo2) - B2; B = 1/(σ* Fo2); C = (A2 + B2)1/2; D = [(C + A)/2]1/2; E = [(C - A)/2]1/2; Fo = uo / (gdo)1/2; and σ* = (2π / L)(do / So).


  1. Chaudhry, Y. M., and D. N. Contractor. 1973. "Application of the Implicit Method to Surges in Open Channels," Water Resources Research, Vol. 9, No. 6, Dec., 1605-1612.

  2. Cunge, J. A., and M. Wegner. 1964. "Numerical Integration of the Saint Venant Equations by an Implicit Finite Difference Scheme," La Houille Blanche, Grenoble, France, Vol. 19, No. 1, Jan., 33-39.

  3. Fread, D. L. 1974. "Numerical Properties of Implicit Four-Point Finite Dfference Equations of Unsteady Flow," NOAA Technical Memorandum NWS HYDRO-18, National Oceanic and Atmospheric Administration, Mar.

  4. Liggett, J. A. 1975. "Basic Equations of Unsteady Flow," Unsteady Flow in Open Channels, K. Mahmood and V. Yevjevich, eds., VoL. I, Water Resources Publications, Fort Collins, Colo.

  5. 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., Vol. I, Water Resources Publications, Fort Collins, Colo.

  6. Lighthill, M. J., and G. B. Whitham. 1955. "On Kinematic Waves I, Flood Movements in Long Rivers," Proceedings, Royal Society of London, Vol. A 229, May, 281-316.

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

  8. 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, Proc. Paper 13392, Dec., 1461-1476.

  9. Price, R. K. 1974. "Comparison of Four Numerical Methods of Flood Routing," Journal of the Hydraulics Division, ASCE. Vol. 100, No. HY7, Proc. Paper 10659, July, 879-899.

  10. Quinn, F. H., and E. B. Wylie. 1972. "Transient Analysis of the Detroit River by the Implicit Method," Water Resources Research, Vol. 8, No. 6, Dec., 1461-1469.

  11. Roache, P. J. 1972. Computational Fluid Dynamics, Hermosa Publishers, Albuquerque, NM.


The following symbols are used in this paper:

A, B, C, D, E = parameters, Eqs. 24-28;

                   a = wave amplitude;

                  C = Courant number;

                   c = wave celerity;

                   d = flow depth;

                  Fo = Froude number;

                   ƒ = function;

                   g = acceleration of gravity;

                   i = √-1;

                   j = space discretization index;

                   L = wavelength;

M, Mo, M1, M2 = parameters;

N, No, N1, N2 = parameters;

                   n = time discretization index;

              P, Q = parameters;

                R1 = attenuation convergence ratio;

                R2 = translation convergence ratio;

                Sƒ = friction slope;

                So = bed slope;

                T = wave period;

                 t = time variable;

                to = initial time;

                 u = mean velocity;

                 x = space variable;

                 α = spatial resolution parameter;

                 β* = complex propagation factor;

                 β*i = imaginary part of β*;

                 β*R = real part of β*;

                 Δx = space step;

                 Δt = time step;

                  δ = logarithmic decrement;

                  θ = weighting factor; and

                σ* = dimensioinless wavenumber.


                  o = equilibrium flow; and

                  n = numerical.


                   ' = perturbed component; and

                  * = dimensionless parameter / variable.

220701 18:30

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