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

Solute transport models are concerned with the fate and transport of inorganic solutes including salts and trace metals [e.g., Zand et al., 1976; Bencala and Walters, 1983; Bencala, 1983; Kuwabara et al., 1984; McKnight and Bencala, 1989; Kim et al., 1992]. These models describe the physical processes of advection and dispersion and some specific chemical and biological reactions. A shortcoming of many solute transport models is the simplistic nature in which reactive chemistry is considered. Existing models rely primarily on the specification of kinetic rate constants and on simple partition coefficient representations of sorption phenomena. For the case of inorganic solutes, the database of kinetic rate constants is strikingly sparse, and many sorption reactions are thought to adhere to more mechanistic sorption models (e.g., surface complexation). While these transport models provide an accurate description of physical transport, they often do not include the degree of chemical sophistication needed to describe pH-dependent processes. Chemical equilibrium models, meanwhile, describe pH-dependent reactions in batch systems, but do not consider transport. Fortunately, many chemical reactions involving inorganic solutes are sufficiently fast so that local equilibrium may be reasonably assumed. It is therefore possible to develop a coupled model wherein a transport model is used to describe physical processes and a chemical equilibrium model is used to quantify pH-dependent reactions. These coupled models are developed by several authors for modeling solute transport in groundwater.

With respect to surface waters, few authors have taken the coupled approach. Chapman et al. [1982] studied a short stream reach using the model RIVEQL, a model that couples a semi-analytical solution of the advection-dispersion equation with the MINEQL equilibrium model. Later work by Chapman [1982] extended the RIVEQL model to include settling of precipitates and sorption processes. Although not explicitly for streams, a related model is the recent work by Wood and Baptista [1993], who developed a diagnostic model for estuaries that considers aqueous speciation and sorption.

In a recent review, Dzombak and Ali [1993] note that there is a need for general models that consider the suite of geochemical and physical processes affecting trace metals in streams. Three findings in their review are relevant to the work presented here. First, most existing models have been developed for specific applications. As such, there has been a tendency to focus on a relatively small subset of the potential chemical reactions affecting metals. Several models, for example, consider sorption but do not consider precipitation/dissolution [e.g., Bencala, 1983; Kuwabara et al., 1984; Wood and Baptista, 1993]. Second, existing models are predominately one-dimensional. Third, many existing models assume a steady flow regime. In this paper, we develop a general model for trace metal fate and transport that considers a variety of chemical processes. The model is pseudo-two-dimensional as it considers the physical process of transient storage. Finally, the model may be used to simulate solute transport under unsteady flow regimes [Runkel and Restrepo, 1993]. Our focus here is on the derivation of the governing differential equations and solution techniques that are unique to the problem of equilibrium-based transport in streams. Model applications are presented in Part II of this paper [Runkel et al., this issue].

Implicit in our approach is the "Local Equilibrium Assumption", wherein chemical reactions are considered sufficiently fast relative to hydrologic processes [Di Toro, 1976; Rubin, 1983]. An alternative to this approach is to model chemical reactions using kinetic rate equations. Although the kinetic approach is a useful paradigm for studying dynamic systems, it is not well suited for many practical geochemical problems. For our purposes, the equilibrium-based approach is preferred for two reasons. First, data obtained from pH modification experiments [McKnight and Bencala, 1989; Kimball et al., 1994] indicate that some pH-dependent reactions are rapid relative to physical transport, such that equilibrium may be assumed. Certain reactions (e.g. acid-base and complexation reactions) occur at very high rates, such that a study of kinetics is not relevant to the time scales of interest. Second, the kinetics of many reactions are poorly understood; this results in a lack of kinetic rate information for model development.

The conceptual and mathematical framework underlying MINTEQ and related models is well documented [Westall et al., 1976; Morel, 1983; Allison et al., 1991]. As a result, only the details essential to model development are given here. We begin by defining chemical "components" as the fundamental building-blocks from which all chemical "species" are derived. Chemical reactions involve two or more components that combine to form a chemical species. In general, components are selected such that: (1) the components combine linearly to form every possible species, and (2) no component may be formed as a combination of other components [Westall et al., 1976]. A species is simply a chemical entity that is formed by combining two or more chemical components. The chemical equilibrium problem entails solving for the unknown species concentrations at equilibrium. This is accomplished by developing mass-action equations to describe the species-producing reactions and mass balance equations for the chemical components.

The sequential iteration approach divides each time step into a "reaction" step and a "transport" step. During the reaction step, the equilibrium submodel is executed for each segment in the stream network. Each segment represents a batch reactor wherein chemical equilibrium is assumed. The equilibrium submodel thus determines the solute mass in dissolved, precipitated and sorbed forms. Based on this information, a transport step is taken in which the solute transport model physically transports the mobile phases of each solute. Because the transport and reaction steps neglect the coupling of the transport and chemistry, the procedure iterates until a specified level of convergence is achieved.

To consider the transport of solid-phase species, it is necessary to develop a set of governing equations and a solution algorithm that is specific to transport in streams. Several assumptions are embodied in the development of the equations that follow:

• Chemical Equilibria. Complexation, precipitation/dissolution and sorption reactions are in a state of local equilibrium.
  
• One-Dimensional Transport. Solute mass is uniformly distributed over the stream's cross-sectional area such that one-dimensional transport is applicable. Given one-dimensional transport, the physical processes of advection, dispersion, lateral inflow, transient storage and settling affect solute mass. All dissolved, precipitated and sorbed species resident in the water column travel at the same advective velocity.
  
• Mobile and Immobile Phases. Solute mass for each chemical component is distributed among five distinct phases. The first three phases represent dissolved, precipitated and sorbed mass that is present in the water column. These three phases are mobile, in that they are subject to transport. The final two phases represent precipitated and sorbed mass that resides on the streambed or is attached to stationary debris in the stream channel; these phases are immobile and are not transported.
  
• Precipitation. Dissolved mass in the water column may form precipitates if the solution becomes oversaturated with respect to the defined solid phases. Any precipitated mass initially resides in the water column and is subject to transport, until it settles to the streambed or redissolution occurs. Precipitation occurs in the water column exclusively, and precipitation directly to the immobile bed is excluded. Precipitated mass may accumulate on the bed, however, as transported precipitates are subject to the force of gravity and settle at a rate defined by a settling velocity. The settling rate of a particle is unaffected by other solids in solution.
  
• Dissolution. When the aqueous solution is undersaturated, dissolution occurs preferentially from the water column. All of the precipitate in the water column is allowed to dissolve before precipitate on the bed is considered for dissolution. This is a reasonable assumption, considering the intimate contact between precipitates in the water column and the flowing waters.
  
• Sorption/Desorption. Dissolved species may sorb to solid phases in the water column or to sorption sites on the streambed. Conversely, sorbed species may desorb from sites in the water column or on the streambed.
  

(1)

(2)

c	concentration of the uncomplexed component species (e.g. Fe3+) [Mass L-3];
xi	concentration of the ith complexed species [Mass L-3];
ai	stoichiometric coefficient of the component in the ith complexed species;
M	number of complexed species.

(3)

(4)

A summary of the processes considered for each phase is presented in Figure 1, where the system is represented as two compartments. The water column compartment contains the three mobile phases, C, P_w and S_w. Immobile substrate (i.e., the streambed or debris) constitutes the second compartment, containing the two immobile phases, P_b and S_b. The three mobile phases are subject to physical transport, as represented by the transport operator, L. The dissolved phase, C, takes part in precipitation/dissolution and sorption/desorption reactions that occur within the water column (interactions with P_w and S_w). The dissolved phase is also affected by dissolution of precipitate from the immobile substrate and by sorption/desorption from immobile sorbents (interactions with P_b and S_b). Finally, C may increase or decrease due to external sources and sinks, as denoted by s_ext. The precipitated and sorbed phases in the water column settle in accordance with settling velocities v_p1 and v_s1 [LT^-1], respectively.

A general mass balance equation for each component is developed by considering the mass associated with each of the five component phases within a control volume. An equation describing conservation of mass for each component is then developed by summing the equations for the individual phases. In the derivations that follow, the compartments depicted in Figure 1 are not treated as separate control volumes, but rather as a single control volume for which a macroscopic mass balance applies [Bird et al., 1960]. Note that this approach differs from the approach used in contemporary sediment-water models for toxic substances. These models are often developed for rivers and lakes in which significant volumes of sediment interact with the water column. In this instance, two or more control volumes are used to represent the sediments and the water column. For our purposes, we are concerned with streams where only a thin, immobile layer of precipitated and sorbed mass interacts with the overlying water column. As such, treatment of the system as a single control volume is an appropriate approach.

Mass balance equations for the five phases are developed below. To simplify the presentation, precipitated and sorbed phases for each component consist of a single species. The dissolved phase, meanwhile, is not limited by this assumption and may be composed of multiple aqueous species. The problem of multiple precipitated and sorbed species for a single component is revisited in a later section. Mass balances for the five phases are given by:

(5)

(6)

(7)

(8)

(9)

L	transport operator;
fw	source/sink term for precipitation/dissolution from the water column [Mass L-3T-1];
fb	source/sink term for dissolution from the immobile substrate [Mass L-3T-1];
gw	source/sink term for sorption/desorption from the water column [Mass L-3T-1];
gb	source/sink term for sorption/desorption from the immobile substrate [Mass L-3T-1];
sext	source/sink term representing external gains and losses [Mass L-3T-1];
vp1	settling velocity for the precipitated phase [LT-1];
vs1	settling velocity for the sorbed phase [LT-1];
d1	effective settling depth [L].

(10)

A	stream channel cross-sectional area [L2];
	instream concentration of an arbitrary component phase [Mass L-3];
	lateral inflow concentration of the arbitrary component phase [Mass L-3];
	storage zone concentration of an arbitrary component phase [Mass L-3];
D	dispersion coefficient [L2T-1];
Q	volumetric flowrate [L3T-1];
qLIN	lateral inflow rate [L3T-1L-1];
x	distance [L];
	storage zone exchange coefficient [T-1].

The mass balance equation for the total component concentration, T, is obtained by summing the mass balance equations for the five individual component phases. This yields:

(11)

As described by Runkel [1993], two solution techniques are available for the solution of the governing equations. The two techniques differ with respect to the primary variable used as input to the chemical equilibrium submodel. The first technique uses the "total water-borne component concentration" (= C + Pw + Sw) as the primary variable. Although this technique is useful for pedagogical purposes, it is unattractive as it often results in two calls to the equilibrium submodel for each computational segment during each iteration. Equilibrium computations make up a large portion of the computational expense required to solve the reactive solute transport problem for groundwater, where the equilibrium submodel is called once for each computational segment during each iteration. Implementation of a solution technique that can potentially double the number of equilibrium computations is therefore not an appealing prospect.

Here we present a second technique wherein the "total component concentration", T, is defined as the primary variable. A differential equation for T is presented as Equation (11). This equation is analogous to the explicit form of the groundwater equation [Yeh and Tripathi, 1989]. Here an implicit form is developed by combining Equation (1) with (11):

(12)

(13)

(14)

The equation set governing the problem consists of three PDEs for each component [Equations (12)-(14)] and the set of AEs representing chemical equilibria. This equation set is solved using a Crank-Nicolson approximation of the governing differential equations and the sequential iteration approach. Presentation of the solution technique requires additional nomenclature. Let n denote an initial time and n+1 denote an advanced time; time n is the previous time at which the state of the system is known, and time n+1 is the current time for which a solution is desired. In addition, k is a counter used to denote the iteration number. Finally, a caret is used to indicate that a given quantity is an estimate. For example, is an estimate of the dissolved concentration for the current iteration at the advanced time level.

The goal of sequential iteration is to solve the set of PDEs describing transport. In general, there is one equation in the form of (12) for each chemical component. Values of the state variables at the initial and advanced time levels are needed to solve for the total component concentrations at the advanced time level (Tⁿ⁺¹) using Crank-Nicolson. The state variables at time level n are available from the previous time step, while estimates of the state variables must be made for time level n+1. Specifically, estimates of P, P_b, S and S_b are needed, as well as the source/sink terms f_b, g_b and s_ext. As shown below, P and S are provided directly by the equilibrium submodel, P_b and S_b are provided via Equations (13) and (14), and the source/sink terms are developed algorithmically. Iteration is required because the values based on the equilibrium calculations are only estimates of the variables at the advanced time level. Solution of the reactive transport problem consists of four steps: initialization, equilibrium calculations, transport calculations and convergence testing (Figure 2).

Step 1: Initialization

Step 2: Equilibrium Calculations

The main task is to estimate the dissolution source/sink term, f_b, using information obtained from the equilibrium calculations. One way to do this is to consider the change in total precipitate, P, from one time step to the next. Re-examining the differential equations derived for each phase, the change in P with time is given by the sum of Equations (6) and (8):

(15)

An algorithm for determining f_b is presented in Figure 3. Using estimates of the total component concentrations, the equilibrium submodel is called to compute the total amount of precipitate present (). If the total amount of precipitate at time n+1 exceeds that for time n (> P ⁿ), the net amount of precipitate has increased, indicating that precipitation has occurred. For the case of precipitation, f_b equals zero. If, however, the net amount of precipitate had decreased ( < P ⁿ), dissolution has occurred and two situations are possible. First, if net dissolution (defined as P ⁿ - ) is less than the amount of precipitate present in the water column (P_wⁿ), all of the dissolved mass is taken from the mobile phase and no dissolution occurs from the bed. Here again f_b equals zero. Second, if net dissolution is not accounted for by mass residing in the mobile phase, mass has dissolved from the immobile phase. In this case, f_b is given by:

(16)

(17)

Step 3: Transport Calculations

Step 4: Convergence Test

(18)

(19)

Pbm	immobile precipitate concentration for precipitated species m [Mass L-3];
pm	total precipitate concentration for precipitated species m [Mass L-3];
vp1m	settling velocity for precipitated species m [LT-1];
fbm	source/sink term for dissolution of immobile precipitated species m [Mass L-3T-1].

(20)

This last point is an important assumption for the model presented here. Due to the empirical nature of the transient storage approach, the storage zone represents both open water (pool and eddies) and water within porous media (flow through the gravel bed and alluvium of the stream channel). This lumped system is not problematic for conservative solute transport models, as the solutes of interest are often dissolved-phase tracers that are transported through the porous media with the water. For the case of reactive solute transport, the modeled solutes are composed of both dissolved and mobile solid phases. If the storage zone is primarily open water, exchange of both dissolved and solid phases between the main channel and the storage zone is a plausible assumption. If the storage zone consists of water in porous media, exchange of solid phases between the two regimes is a questionable assumption, as solid mass entering the storage zone may not re-enter the main channel due to formation of bonds between the solid and the porous media.

The lumped nature of the storage zone therefore presents a problem for model development. If it is known a priori that the storage zone consists of open water, the correct model formulation involves a transient storage mechanism that affects both dissolved and solid phases. Conversely, if the storage zone is composed of porous media, the storage mechanism may be set up to affect the dissolved phase only. In reality, storage zones are usually some combination of open water and porous media. In the development that follows, the open-water scenario is assumed; i.e. dissolved and solid phases are subject to transient storage. This is a logical choice: it correctly models the open water portions of the storage zone, and is also applicable to the porous media component, if the solid particles are small relative to the pore space, such that solid-phase transport occurs.

To implement the transient storage approach, mass balance equations are developed to define the storage zone concentrations. The total component concentration in the storage zone, T_s, is given by:

(21)

(22)

(23)

(24)

(25)

(26)

fsb	source/sink for dissolution from the storage zone immobile substrate [Mass L-3T-1];
gsb	source/sink for sorption/desorption from the storage zone immobile substrate [Mass L-3T-1];
vp2	settling velocity for the storage zone precipitated phase [LT-1];
vs2	settling velocity for the storage zone sorbed phase [LT-1];
d2	effective storage zone settling depth [L].

The solution technique presented here was developed due to the need to consider solid-phase transport. As noted earlier, the groundwater models reviewed by Yeh and Tripathi [1989] assume that solid phases are not subject to transport. Although this is a reasonable assumption for many subsurface applications, it is not always valid as colloidal materials are known to be transported through porous media [McCarthy and Zachara, 1989]. The techniques developed herein may therefore be of use within groundwater models that incorporate solid-phase transport.

Although we have made significant progress, additional efforts are required before a truly general model for trace metals will be realized. First, our analysis has centered on the problem of precipitation/dissolution; future efforts should be directed towards the development of additional sorption algorithms. Second, a general model will need to consider both kinetic and equilibrium-based reactions. Flow velocities in many streams are such that the assumption of local equilibrium is not appropriate for certain reaction classes. Specific reactions that may require a kinetic approach include dissolution and sorption/desorption reactions between the water column and the streambed. Finally, field-scale model applications that consider multiple reaction types (e.g. sorption and precipitation/dissolution) and mixed (i.e., kinetic and equilibrium) reactions may require substantial computer resources. Innovative and efficient computational techniques such as parallel processing may be needed to lower execution times.

Bencala, K.E., Simulation of solute transport in a mountain pool-and-riffle stream with a kinetic mass transfer model for sorption, Water Resour. Res., 19(3), 732-738, 1983.

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.

Bird, R.B., W.E. Stewart, and E.N. Lightfoot, Transport Phenomena, John Wiley, New York, 1960.

Castro, N.M., and G.M. Hornberger, Surface-subsurface water interactions in an alluviated mountain stream channel, Water Resour. Res., 27(7), 1613-1621, 1991.

Cederberg, G.A., R.L. Street, and J.O. Leckie, A groundwater mass transport and equilibrium chemistry model for multicomponent systems, Water Resour. Res., 21(8), 1095-1104, 1985.

Chapman, B.M., Numerical simulation of the transport and speciation of nonconservative chemical reactants in rivers, Water Resour. Res., 18(1), 155-167, 1982.

Chapman, B.M., R.O. James, R.F. Jung, and H.G. Washington, Modelling the transport of reacting chemical contaminants in natural streams, Aust. J. Mar. Freshw. Res., 33, 617-628, 1982.

Chapman, B.M., D.R. Jones, and R.F. Jung, Processes controlling metal ion attenuation in acid mine drainage streams, Geochim. Cosmochim. Acta, 47, 1957-1973, 1983.

Di Toro, D.M., Combining chemical equilibrium and phytoplankton models - A general methodology, in Modeling Biochemical Processes in Aquatic Ecosystems, R.P. Canale, editor, Ann Arbor Science, Ann Arbor, Michigan, 224-243, 1994

Dzombak, D.A., and M.A. Ali, Hydrochemical modeling of metal fate and transport in freshwater environments, Water Poll. Res. J. Canada, 28(1), 7-50, 1993.

Engesgaard, P., and K.L. Kipp, A geochemical transport model for redox-controlled movement of mineral fronts in ground-water flow systems: A case of nitrate removal by oxidation of pyrite, Water Resour. Res., 28(10), 2829-2843, 1992.

Harvey, J.W., and K.E. Bencala, The effect of streambed topography on surface-subsurface water exchange in mountain catchments, Water Resour. Res., 29(1), 89-98, 1993.

Kim, B.K., A.P. Jackman, and F.J. Triska, Modeling biotic uptake by periphyton and transient hyporrheic storage of nitrate in a natural stream, Water Resour. Res., 28(10), 2743-2752, 1992.

Kimball, B.A., D.M. McKnight, G.A. Wetherbee, and R.A. Harnish, Mechanisms of iron photoreduction in a metal-rich, acidic stream (St. Kevin Gulch, Colorado, U.S.A.), Chem. Geol., 96, 227-239, 1992.

Kimball, B.A., R.E. Broshears, D.M. McKnight, and K.E. Bencala, Effects of instream pH modification on transport of sulfide-oxidation products, in Environmental Geochemistry of Sulfide Oxidation, C.N. Alpers and D.W. Blowes, editors, ACS National Meeting, August 1992, ACS Symposium Series 550, 224-243, 1994.

Konikow, L.F., and J.D. Bredehoeft, Ground-water models cannot be validated, Adv. Water Resour. Res., 15, 75-83, 1992.

Kuwabara, J.S., H.V. Leland, and K.E. Bencala, Copper transport along a Sierra Nevada stream, J. Environ. Engr., 110(3), 646-655, 1984.

Mangold, D.C., and C.F. Tsang, A summary of subsurface hydrological and hydrochemical models, Reviews of Geophysics, 29(1), 51-79, 1991.

McCarthy, J.F., and J.M. Zachara, Subsurface transport of contaminants, Environ. Sci. Technol., 23(5), 496-502, 1989.

McKnight, D.M., B.A. Kimball, and K.E. Bencala, Iron photoreduction and oxidation in an acidic mountain stream, Science, 240, 637-640, 1988.

McKnight, D.M., and K.E. Bencala, Reactive iron transport in an acidic mountain stream in Summit County, Colorado: A hydrologic perspective, Geochim. Cosmochim. Acta, 53, 2225-2234, 1989.

McKnight, D.M., K.E. Bencala, G.W. Zellweger, G.R. Aiken, G.L. Feder, and K.A. Thorn, Sorption of dissolved organic carbon by hydrous aluminum and iron oxides occurring at the confluence of Deer Creek with the Snake River, Summit County, Colorado, Environ. Sci. Technol., 26(7), 1388-1396, 1992.

Morel, F.M.M., Principles of Aquatic Chemistry, John Wiley, New York, 1983.

Rubin, J., Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions, Water Resour. Res., 19(5), 1231-1252, 1983.

Runkel, R.L., Development and application of an equilibrium-based simulation model for reactive solute transport in small streams, Ph.D. Dissertation, Univ. of Colo., Boulder, 1993.

Runkel, R.L., and R.E. Broshears, One dimensional transport with inflow and storage (OTIS): A solute transport model for small streams, Tech. Rep. 91-01, Center for Adv. Decision Support for Water and Environ. Syst., Univ. of Colo., Boulder, 1991.

Runkel, R.L., and S.C. Chapra, An efficient numerical solution of the transient storage equations for solute transport in small streams, Water Resour. Res., 29(1), 211-215, 1993.

Runkel, R.L., and P.J. Restrepo, Solute transport modeling under unsteady flow regimes: An application of the Modular Modeling System, in Water management in the `90s: a time for innovation, proceedings of the 20th anniversary conference of the Wat. Res. Planning and Mgmt. Div., ASCE, Seattle, Washington, May 1-5, 1993, K. Hon, editor, 1993.

Westall, J.C., J.L. Zachary, and F.M.M. Morel, MINEQL: A computer program for the calculation of chemical equilibrium composition in aqueous systems, Tech. Note 18, Water Qual. Lab., Dep. of Civ. Eng., Mass. Inst. of Technol., Cambridge, 1976.

Wood, T.M., and A.M. Baptista, A model for diagnostic analysis of estuarine geochemistry, Water Resour. Res., 29(1), 51-71, 1993.

Yeh, G.T., and V.S. Tripathi, A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components, Water Resour. Res., 25(1), 93-108, 1989.

Yeh, G.T., and V.S. Tripathi, A model for simulating transport of reactive multispecies components: Model development and demonstration, Water Resour. Res., 27(12), 3075-3094, 1991.

Zand, S.M., V.C. Kennedy, G.W. Zellweger, and R.J. Avanzino, Solute transport and modeling of water quality in a small stream, J. Res. U.S. Geol. Surv., 4(2), 233-240, 1976.

Figure 2 : The Sequential Iteration Approach for the equilibrium-based surface water model, where the total component concentration (T) is the primary variable.