The most successful and accepted procedure for solving the one-dimensional unsteady flow equations is the four-point implicit scheme, also known as the box scheme (see figure below). Under this scheme, space derivatives and function values are evaluated at an interior point, (n+\theta) \Delta t. Thus values at (n+1) \Delta t enter into all terms in the equations. For a reach of river, a system of simultaneous equations results. The simultaneous solution is an important aspect of this scheme because it allows information from the entire reach to influence the solution at any one point. Consequently, the time step can be significantly larger than with explicit numerical schemes. Von Neumann stability analyses performed by Fread (1974), and Liggett and Cunge (1975), show the implicit scheme to be unconditionally stable (theoretically) for 0.5 < \theta ≤ 1.0, conditionally stable for \theta = 0.5, and unstable for \theta < 0.5. In a convergence analysis performed by the same authors, it was shown that numerical damping increased as the ratio \lambda / \Delta x decreased, where \lambda is the length of a wave in the hydraulic system. For streamflow routing problems where the wavelengths are long with respect to spatial distances, convergence is not a serious problem.
In practice, other factors may also contribute to the non-stability of the solution scheme. These factors include dramatic changes in channel cross-sectional properties, abrupt changes in channel slope, characteristics of the flood wave itself, and complex hydraulic structures such as levees, bridges, culverts, weirs, and spillways. In fact, these other factors often overwhelm any stability considerations associated with \theta. Because of these factors, any model application should be accompanied by a sensitivity study, where the accuracy and the stability of the solution are tested with various time and distance intervals.

The following notation is defined:
1) |
\displaystyle f_j=f_j^n |
and:
2) |
\Delta f_j = f_j^{n+1} - f_j^n |
then:
3) |
f_j^{n+1} = f_j + \Delta f_j |
The general implicit finite difference forms are:
Time derivative
4) |
\displaystyle \frac{\partial f}{\partial t} \approx \frac{\Delta f}{\Delta t} = \frac{\Delta f_{j+1} + \Delta f_j}{2 \Delta t} |
Space derivative
5) |
\displaystyle \frac{\partial f}{\partial x} \approx \frac{\Delta f}{\Delta x} = \frac{\left( f_{j+1} - f_j \right) + \theta \left( \Delta f_{j+1} - \Delta f_j \right)}{\Delta x} |
Function value
6) |
\displaystyle f \approx \overline{f} = 0.5 (f_j + f_{j+1}) + 0.5 \theta (\Delta f_j + \Delta f_{j+1} ) |