An Efficient Numerical Solution of the Transient Storage Equations for Solute Transport in Small Streams

by
Robert L. Runkel and Steven C. Chapra

Center for Advanced Decision Support for Water and Environmental Systems,
University of Colorado, Boulder

Disclaimer: The online document displayed is based on the final draft provided to the editor. Minor discrepancies between the online document and the published version may therefore exist.

Introduction

Several investigators have proposed surface water solute transport models that incorporate the effects of transient storage [Thackston and Krenkel, 1967; Thackston and Schnelle, 1970; Valentine and Wood, 1977; Nordin and Troutman, 1980; Jackman et al., 1984; Bencala and Walters, 1983; Bencala et al., 1990; Kim et al., 1990]. These transient storage (or `dead zone') models consider a physical mechanism wherein solute mass is exchanged between the main channel and an immobile storage zone. This quasi-two-dimensional process is represented by adding a storage term to the conventional advection-dispersion equation.

Transient storage occurs in small streams when portions of the transported solute become isolated from the main channel in small pools and in the gravel underbed. As a result of this storage phenomenon, the magnitude of a typical solute tracer pulse is attenuated and its position delayed. To account for attenuation, two distinct zones are considered. The first zone represents the main channel that is usually considered when modeling advective surface-water systems. Processes influencing solute concentrations in this zone include advection, dispersion, lateral inflow, transient storage, and chemical reaction. A second area, the storage zone, represents recirculating pools, underflow channels and other areas that are immobile relative to flow in the main channel. Only the processes of storage and chemical reaction are considered in the storage zone. The zones are linked by an exchange term that acts to transfer mass between the two regimes.

The transient-storage equations presented herein are generally applicable to streams and rivers in which one-dimensional solute transport may be assumed. As described by Fischer et al. [1979], one-dimensional analysis is valid for systems in which solute mass is uniformly distributed over the stream's cross-sectional area. Although this uniformity rarely occurs in nature, it is a reasonable assumption for streams of small to moderate width and depth.

Although numerous applications of the transient storage equations may be found in the literature, little attention has been paid to the numerical aspects of the approach. Of particular interest is the coupled nature of the equations describing mass conservation for the stream channel and the storage zone. This paper details an implicit finite difference technique that decouples the governing equations. Additional numerical issues are also discussed.

Governing Differential Equations &
Finite Difference Approximations

Conservation of mass for the stream and storage-zone segments is given by (1) and (2), respectively [Bencala and Walters, 1983; Runkel and Broshears, 1991]:

(1)

(2)

where

A stream channel cross-sectional area [L2]
AS storage zone cross-sectional area [L2]
C in-stream solute concentration [M L-3]
CL solute concentration in lateral inflow [M L-3]
CS storage zone solute concentration [M L-3]
D dispersion coefficient [L2T-1]
Q volumetric flowrate [L3T-1]
qLIN lateral inflow rate [L3T-1L-1]
t time [T]
x distance [L]
alpha storage zone exchange coefficient [T-1]
lambda in-stream first-order decay coefficient [T-1]
lambda S storage zone first-order decay coefficient [T-1]

Boundary and initial conditions are given by:

(3)
(4)
(5)
(6)

where

l main channel length [L]
Cbc in-stream solute concentration at the upstream boundary [M L-3]
Cinit initial in-stream background concentration [M L-3]
CSinit initial storage zone background concentration [M L-3]

To solve (1) and (2) for the general case in which the parameters vary in time and space, numerical solution techniques are employed. In the development that follows, central divided-differences are used to approximate the spatial derivatives. It should be noted, however, that the following discussion is equally applicable when backward or arbitrary differences are used. In any event, (1) becomes:

(7)

(8)

where delta x denotes the length of a one-dimensional segment (i.e., control volume) and i-1, i and i+1 subscripts denote the upstream, central, and downstream segments, respectively. Note that CS in (7) (and the equations that follow) is the storage zone concentration in segment i; the i subscript is omitted as a notational convenience. Also note that Equation (8) is written for the special case of equal segment sizes (delta xi-1 = delta xi = delta xi+1). A similar equation for the general case of variable segment lengths is given by Runkel and Broshears [1991].

The Crank-Nicolson Method

Ordinary differential equations such as (2) and (7) may be solved using a variety of techniques. For the case of one-dimensional models, it is often expedient to use an implicit method such as Crank-Nicolson. This method has several noteworthy advantages. First, the Crank-Nicolson method is second-order accurate in both time and space. Second, the one-dimensional nature of the model results in the formation of a tridiagonal coefficient matrix that can be solved using an efficient formulation of the Thomas algorithm. Finally, Crank-Nicolson exhibits strong stability [Isaacson and Keller, 1966].

Application of the Crank-Nicolson approach proceeds as follows. To achieve second-order accuracy, difference equations are developed for the midpoint of a time step. Letting j denote an initial time and j+1 denote an advanced time, the main channel equation [Equation (7)] is formulated for time j+1/2. The right-hand-side of (7) at time j+1/2 is simply the average of the term at times j and j+1. In addition, the time derivative, dC/dt, is estimated using a centered difference approximation:

(9)

where delta t is the integration time step [T]. Equation (7) thus becomes:

(10)

where

(11)

Because (10) is dependent on the solute concentrations in the neighboring segments at the advanced time level (Ci-1, Ci+1 at time j+1), it is not possible to solve explicitly for Cij+1 (hence we have an `implicit' method). We can, however, rearrange (10) so that all of the known quantities (solute concentrations at time j) appear on the right-hand side and all of the unknown quantities (solute concentrations at time j+1) appear on the left. One exception to this rearrangement in that an unknown quantity, the storage zone concentration at the advanced time level (CSj+1), remains on the right-hand side. This exception is discussed in a subsequent section. Rearrangement yields:

(12)

This, in turn, may be simplified by collecting terms:

(13)

where

(14)
(15)
(16)
(17)

In developing (13) for each of the segments in the stream network, a set of linear algebraic equations is produced. These equations are solved simultaneously to obtain the in-stream solute concentration, Cj+1, in each of the stream segments. A hypothetical system of equations representing a five-segment network is shown below.

(18)

The concepts presented in the preceding paragraphs are also applied to the storage zone equation [Equation (2)]. This results in the following Crank-Nicolson equation:

(19)

In contrast to the stream equation, Equation (19) may be solved explicitly for the variable of interest, CSj+1. This yields:

(20)

where

(21)

Decoupling of the Stream and Storage Zone Equations

Equations (13) and (20) appear to be a set of coupled equations due to the presence of an unknown quantity, CSj+1, on the right-hand side of (13). (Recall that the right-hand side is to contain only known quantities.) This coupling suggests an iterative solution technique whereby (13) and (20) are solved in sequence until some desired level of convergence is obtained [Jackman et al., 1984]. This iterative predictor-corrector approach is inefficient in that (13) and (20) are solved more than once during each time step.

Fortunately, decoupling is easily achieved by noting a fundamental difference between the two equations. First, (13) is a Crank-Nicolson expression for a partial differential equation (PDE) containing spatial derivatives [Equation (1)]. As such, it is not possible to solve explicitly for the dependent variable, Cj+1. In contrast, (19) is a Crank-Nicolson approximation of an ordinary differential equation (ODE), wherein a spatial derivative is lacking [Equation (2)]. This allows one to solve explicitly for the dependent variable, CSj+1 [Equation (20)]. It is therefore possible to substitute (20) into (13), thereby eliminating CSj+1. This effectively decouples the Crank-Nicolson expressions. One may easily envision other sets of mixed equations (ODEs and PDEs) where this procedure is equally applicable.

Decoupling of the transient storage equations proceeds as follows. In Equation (20) the storage concentration is a function of two known quantities, Cj and CSj, and one unknown quantity, Cj+1. Substitution of (20) into (17) provides a new expression for R:

(22)

Although R' contains an unknown quantity, (22) is a much more convenient expression than (17). The unknown quantity is now Cj+1, a variable that already appears on the left-hand side of (13). We can therefore move the term involving Cj+1 to the left side of (13), creating new expressions for F and R:

(23)

(24)

Because R" involves only known quantities, (13) can be solved independently for the in-stream solute concentration, Cj+1. Having solved (13), the storage zone equation [Equation (20)] becomes a function of three known quantities, CSj, Cj, and Cj+1. We have thus decoupled the governing Crank-Nicolson expressions.

Thomas Algorithm

As given by (18), application of the Crank-Nicolson procedure results in the formation of a linear system of algebraic equations:

(25)

where

[M] tridiagonal coefficient matrix of dimension N by N
{C} vector of length N representing the unknown solute concentrations
{R} forcing function vector of length N, as given by Equation (24)
N
 number of segments in the stream network

Due to the tridiagonal nature of matrix [M], Equation (25) may be solved for {C} using the Thomas Algorithm [Thomas, 1949]. This approach exploits the tridiagonal structure of [M] by eliminating useless operations on zero elements. Although this leads to great computational savings, the algorithm is frequently presented in a format that confounds operations on the coefficient matrix [M] with those on the right-hand-side vector {R} [e.g. Cheney and Kincaid, 1985; Press et al., 1986). This leads to redundant computations when applied to cases involving an invariant coefficient matrix [M] and multiple right-hand-side vectors.

Increased efficiency can be obtained by formulating the Thomas Algorithm as an LU decomposition method [Chapra and Canale, 1988; Burden and Faires, 1988]. This method decomposes the coefficient matrix into `lower' and `upper' diagonal matrices. This `decomposition' is followed by `substitution' steps involving {R} that provide the solution for {C}.

The primary advantage of the LU Decomposition approach is the ability to efficiently evaluate multiple right-hand side vectors ({R}). Because the first step, decomposition, involves only the coefficient matrix, [M], it need not be repeated for each vector {R} . As a result, solutions for multiple right-hand side vectors may be obtained using a single decomposition step in conjunction with multiple substitution steps.

For the case at hand, there are multiple right-hand side vectors, as {R} is a function of the time-varying solute concentrations. Given a steady flow regime, the coefficient matrix remains constant throughout the simulation. A decomposition step is therefore not required for each time step. This considerably reduces the number of operations required to complete a given simulation.

For conservative substances, the coefficient matrix is a function of the model's physical parameters, i.e. [M] is not specific to a given solute. Due to this solute independence, the decomposition phase for the Thomas Algorithm is completed only once for each set of physical parameters (i.e. given steady flow conditions, only one decomposition step is required for the entire simulation). A substitution phase, meanwhile, is required for each solute at each time step, as {R} is a function of the solute concentrations at the current time level [see Equation (24)].

Unlike the conservative case, the coefficient matrix for nonconservative solutes is a function of both physical and chemical parameters. Matrix decomposition is therefore required for each solute being modeled, as the values of the chemical parameters vary between solutes. For steady flow regimes, the number of decomposition steps equals the number of solutes. As with conservative substances, the substitution step is required for each solute at each time step.

In summary, the solution of (25) is dependent on the types of solutes being modeled. Specifically, the coefficient matrix is solute-independent for conservative substances, while it is solute-specific for non-conservative solutes.

Benchmark Runs

In order to test the relative efficiencies introduced by the foregoing methods, three transient storage codes were developed. In the first version, the governing equations remained coupled and both steps of the Thomas Algorithm were completed at each time step. For the second version, the transport equations were uncoupled as described in the preceding text. A third and final version included the decoupled approach as well as the LU decomposition formulation of the Thomas algorithm.

Using the three computer codes, a series of benchmark runs were conducted. Simulation parameters for the benchmark runs are shown in Table 1. In each run, the transport of five conservative constituents was simulated under steady flow conditions (i.e. Q, A, qLIN and CL were time invariant). Solute concentrations were determined for a thirty-hour period using an integration time step of 0.01 hours. The five hundred meter stream was modeled as a linear reach comprised of one meter segments. All runs were completed on a SUN Sparc 2 workstation.

Figure 1 summarizes the results from the benchmark simulations. Two conclusions may be drawn from the figure. First, judicious use of the Thomas Algorithm leads to a small, but not insignificant, decrease in run-time. Second, decoupling of the transport equations leads to a significant decrease in computational expense. This decrease is directly attributable to the fact that decoupling eliminates the need for the iterative approach in which the Crank-Nicolson equations are solved more than once during each time step. Given the system configuration and the parameters chosen for the benchmark runs, the iterative approach converges after two iterations. As such, the decoupling results in a 50 percent decrease in run time. It is important to note that situations may arise in which more than two iterations would be required for the iterative scheme. In this instance the decoupling will lead to an even greater decrease in run time. Ongoing research is underway to evaluate this issue.

Conclusions

In the foregoing analysis we have described two operations that decrease the computational effort required for the solution of the transient storage solute transport equations. First, the governing equations have been decoupled so that the main channel and storage zone equations may be solved independently. This eliminates the need for the predictor-corrector approach, thereby reducing the required simulation time by approximately 50 percent. A second improvement is the efficient use of the Thomas Algorithm. By considering separate decomposition and substitution steps, run-time is reduced by an additional 14 percent. Both of these techniques may be easily incorporated into existing codes or new applications where simulation run-time is of concern.

As discussed previously, the decoupling procedure described here is relevant to other types of fate and transport models. Specifically, the method may be applied to the equations describing toxic substances in one-dimensional streams and estuaries underlain by stationary sediment beds. Thomann and Mueller [1987] have developed a similar decoupling approach for computing steady-state distributions with a matrix-based method. Our technique provides a means to efficiently simulate the dynamics of toxicants for such systems with stable implicit methods.

In addition, other applications are possible. For example, the approach is applicable to solute transport in porous media [e.g. Hornberger et al., 1992]. Research efforts are presently underway to generalize the decoupling procedure and to further define its performance and applicability.

Although computer run time may no longer seem to be the critical issue it was in the past, computations of the type described here still require efficiency for two reasons. First, the transient-storage solute-transport equations are presently being coupled with large chemical equilibrium models. For such applications, computational savings are critical for keeping total simulation times within reasonable bounds. Second, computer models are increasingly being integrated into decision-support systems. For such cases, the ability to implement interactive computations hinges on efficient model run times. The present paper has been directed toward making such applications feasible.

Acknowledgments

This work was carried out as part of a co-operative agreement with the United States Geological Survey (USGS) in Denver, Colorado. Funding for the work was provided by the USGS Toxic Substances Hydrology Program. The authors wish to thank Kenneth Bencala, Robert Broshears, Briant Kimball, Pete Loucks, Diane McKnight, and Pedro Restrepo for their review of this document and their input during software development. Finally, the authors wish to acknowledge two anonymous reviewers for their constructive comments on this manuscript.

References

Bencala, K.E. and R.A. Walters, Simulation of solute transport in a mountain pool-and-riffle stream: a transient storage model, Water Resour. Res., 19(3), 718-724, 1983.

Bencala, K.E., D.M. McKnight and G.W. Zellweger, Characterization of transport in an acidic and metal-rich mountain stream based on a lithium tracer injection and simulations of transient storage, Water Resour. Res., 26(5), 989-1000, 1990.

Burden, R.L. and J.D. Faires, Numerical Analysis, 4th Edition, pp. 373-376, PWS-Kent Publishing company, Boston, 1988.

Chapra, S.C. and R.P. Canale, Numerical Methods For Engineers, Second Edition, pp. 271-290, McGraw-Hill, New York, 1988.

Cheney, W. and D. Kincaid, Numerical Mathematics and Computing, pp. 232-234, Brooks/Cole Publishing Company, Monterey, CA, 1985.

Fischer, H.B., E.J. List, R.C.Y. Koh, J. Imberger and N.H. Brooks, Mixing in Inland and Coastal Waters, pp. 263-276, Academic Press, San Diego, 1979.

Hornberger, G.M., A.L. Mills, and J.S. Herman, Bacterial transport in porous media: Evaluation of a model using laboratory observations, Water Resour. Res., 28(3), 915-938, 1992.

Isaacson, E., and H.B. Keller, Analysis of Numerical Methods, pp. 508-512, John Wiley & Sons, New York, 1966.

Jackman, A.P., R.A. Walters and V.C. Kennedy, Transport and concentration controls for chloride, strontium, potassium and lead in Uvas Creek, a small cobble-bed stream in Santa Clara County, California, U.S.A., 2. Mathematical Modeling, J. Hydrol., 75, 111-141, 1984.

Kim, B.K., A.P. Jackman and F.J. Triska, Modeling Transient Storage and Nitrate Uptake Kinetics in a Flume Containing a Natural Periphyton Community, Water Resour. Res., 26(3), 505-515, 1990.

Nordin, C.F. and B.M. Troutman, Longitudinal dispersion in rivers: the persistence of skewness in observed data, Water Resour. Res., 16(1), 123-128, 1980.

Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, Numerical Recipes: The Art of Scientific Computing, pp. 40-41, Cambridge University Press, Cambridge, U.K., 1986.

Runkel, R.L. and R.E. Broshears, One dimensional transport with inflow and storage (OTIS): A solute transport model for small streams, CADSWES Technical Report 91-01, University of Colorado, Center for Advanced Decision Support for Water and Environmental Systems, Boulder, Colorado, 1991.

Thackston, E.L. and P.A. Krenkel, Longitudinal mixing in natural streams, J. San. Engr. Div., ASCE, 93(SA5), 67-90, 1967.

Thackston, E.L. and K.B. Schnelle, Predicting effects of dead zones on stream mixing, J. San. Engr. Div., ASCE, 96(SA2), 319-331, 1970.

Thomann, R.V., and J.A. Mueller, Principles of Surface Water Quality Modeling and Control, p. 581, Harper & Row, New York, 1987.

Thomas, L.H., Elliptic problems in linear difference equations over a network, Watson Scientific Computing Laboratory Report, Columbia University, New York, 1949.

Valentine, E.M. and I.R. Wood, Longitudinal dispersion with dead zones, J. Hydraul. Div., ASCE, 103(HY9), 975-990, 1977.

Figures

Table 1. Simulation Parameters for Benchmark Runs

Parameter Value
Number of Solutes
5

Type of Solute Conservative
Flow Regime
Steady

Integration Time Step [hours]
0.01

Simulation Period [hours]
30

Number of Segments
500

Segment Length [meters]
1

Fig. 1. Benchmark runs for versions 1-3.
Version 1: Coupled equations and complete Thomas Algorithm.
Version 2: Uncoupled equations and complete Thomas Algorithm.
Version 3: Uncoupled equations with the LU Decomposition form of the Thomas Algorithm.