Solute transport with equilibrium aqueous complexation and either sorption or ion exchange: Simulation methodology and applications

Journal of Hydrology, 90 (1987) 81-115 Elsevier Science Publishers B.V., Amsterdam - - Printed in The Netherlands

81

[2]

SOLUTE TRANSPORT WITH EQUILIBRIUM AQUEOUS COMPLEXATION AND EITHER SORPTION OR ION EXCHANGE: SIMULATION METHODOLOGY AND APPLICATIONS

FRANK M. LEWIS 1., CLIFFORD I. VOSS 1 and JACOB RUBIN 2

1 U.S. Geological Survey, National Center, Reston, VA 22092 (U.S.A.) 2 U.S. Geological Survey, Menlo Park, CA 94025 (U.S.A.) (Received M a r c h 19, 1986; revised and accepted J u n e 21, 1986)

ABSTRACT Lewis, F.M., Voss, C.I. and Rubin, J., 1987. Solute transport with equilibrium aqueous complexaeither sorption or ion exchange: Simulation methodology and applications. J. Hydrol., 90: 81 115. Methodologies that account for specific types of chemical reactions in the simulation of solute transport can be developed so they are compatible with solution algorithms employed in existing transport codes. This enables the simulation of reactive transport in complex multidimensional flow regimes, and provides a means for existing codes to account for some of the fundamental chemical processes t h a t occur among transported solutes. Two equilibrium-controlled reaction systems demonstrate a methodology for accommodating chemical interaction into models of solute transport. One system involves the sorption of a given chemical species, as well as two aqueous complexations in which the sorbing species is a participant. The other reaction set involves binary ion exchange coupled with aqueous complexation involving one of the exchanging species. The methodology accommodates these reaction systems through the addition of nonlinear terms to the transport equations for the sorbing species. Example simulation results show (1) the effect equilibrium chemical parameters have on the spatial distributions of concentration for complexing solutes; (2) t h a t an interrelationship exists between mechanical dispersion and the various reaction processes; (3) t h a t dispersive parameters of the porous media cannot be determined from reactive concentration distributions unless the reaction is accounted for or the influence of the reaction is negligible; (4) how the concentration of a chemical species may be significantly affected by its participation in an aqueous complex with a second species which also sorbs; and (5) t h a t these coupled chemical processes influencing reactive transport can be demonstrated in twodimensional flow regimes. INTRODUCTION

In recent years, research in subsurface contaminant transport has expanded to include the influence of various chemical reactions on solute migration (e.g., Rubin and James, 1973; Valocchi et al., 1981; Rubin, 1983). The principal objective in much of this research has been to develop numerical models that may be used to simulate the effect of specific equilibrium chemistries on * Present address: CH2M Hill, P.O. Box 4400 Reston, VA 22090 (U.S.A.)

0022-1694/87/$03.50

© 1987 Elsevier Science Publishers B.V.

82

reactive solutes during transport through porous media. To meet this objective successfully requires a strategy that efficiently couples two different processes, namely chemistry and mass transport. Each of these processes is characterized mathematically by fundamentally different equation forms. These forms are algebraic relationships that define equilibrium sorptive and chemical interactions, and differential equations that describe the physical movement of solute through porous media. In attempting to solve these mathematically mixed sets of equations, a variety of techniques have been developed. Most approaches may be classified following the general dichotomy in modeling reactive transport that is discussed by Jennings et al. (1982) and Cederberg et al. (1985). One such general approach involves incorporating the reaction mathematics directly into the partial differential equations (PDEs) that track solute-mass balances through advection and dispersion. Here, the governing equations all have the same form (differential), and thus only a single simultaneous solution algorithm is employed in the determination of both transport and chemical interaction. However, models based on this approach vary, particularly with respect to the chemical (mathematical) form of the transported variable. For example, the work of Kirkner et al. (1984) considers the transport of total dissolved concentrations and determines individual free ion concentrations from subsequent algebraic equations of chemical equilibrium which are independent from those for basic transport. On the other hand, models developed by Jennings et al. (1982) and Miller and Benson (1983) obtain the individual reactant concentrations directly. The alternative general approach to simulating reactive transport involves keeping the algebraic equations that define the chemistry separate from the PDEs that describe the mass balances. A solution is obtained by iterating between the two sets of equations: solving alternately for physical transport alone followed by chemical interaction until a new equilibrium is reached. This approach has been employed in studies by Grove and Wood (1979), Narasimhan et al. (1984), and Cederberg et al. (1985) among others. The approach makes it easier to treat chemical complexity because existing computer codes such as WATEQF (Plummer et al., 1976) and PHREEQE (Parkhurst et al., 1982) may be used to determine equilibrium calculations for large chemical systems. In the first approach the form of the mass-balance equations depends directly on the way in which the chemical mathematics are incorporated; consequently, numerical solution techniques applied to these equations are similarly dependent upon the nature of the chemistry and how it is incorporated. As a result, model development based on this approach has emphasized the specialized accommodation of particular chemistries within the transport equations as well as the subsequent solution of these equations. Previous work (e.g., Jennings et al., 1982; Miller and Benson, 1983; Kirkner et al., 1984) has resulted in new chemical transport models that solve specially derived transport equations. Each of these models employs a one-dimensional algorithm that has been developed specifically for the resulting equations. Accordingly, these models

83 may not be well suited to simulate complex multidimensional hydrodynamic systems that may be encountered in certain field applications. Methods based on the secohd approach may be more versatile in this regard, but they also have some limitations as pointed out by Jennings et al. (1982). A significant disadvantage may arise when chemical equilibria models such as those mentioned above are employed. These computer codes are large in order to be general, and may not be highly efficient inasmuch as the solution at each time step represents equilibrium in a physically static system. Their use in conjunction with transport models where the equilibrium code must be activated once for each node, on every iteration, and for all time steps, may lead to excessive computer time requirements, particularly if a fine grid is required to simulate the flow and transport system adequately. An alternative to the above approaches is to incorporate the chemistry so the resulting mass-balance equations remain in a form that is compatible with solution algorithms utilized by existing transport models. Such a formulation of the basic reactive-transport equations may be obtained by: (1) maximizing the number of linear mass-balance equations; (2) structuring the remaining mass-balance equations to be quasilinear in that the highest order derivative appears to the first power; and (3) formulating the mass-balance equations to have a single dependent variable whenever possible. In this way, the mathematical system describing the transport of reacting solutes consists of massbalance equations that are linear or quasilinear in a single dependent variable. Moreover, these equations may be solved independently, prior to the solution for the variables of the chemical algebra. The purpose of this paper is two-fold. First, to present a methodology utilizing the concept that facilitates the modification of existing, nonreacting, transport algorithms to simulate the effects of a few simple, but fundamental, types of chemical reactions; and second, to apply this methodology, and demonstrate the influence of various reactive processes on solute transport. Specifically, the reactions we consider are equilibrium-controlled aqueous complexation coupled with either sorption or binary ion exchange. Simulating generic forms of these reactions, individually or in various simultaneous combinations, provides some insight into the relationship between reactive effects and dispersive effects, and the role this relationship plays in determining the resulting distribution of solute concentrations. The methodology comprises a modular mathematical structure that essentially adds two new terms to the classic advection-dispersion equation, and should, therefore, be compatible with models that are based on this equation. Here, the methodology is applied to the model SATRA, a version of the U.S. Geological Survey flow and transport code SUTRA (Voss, 1984) for constantdensity fluid under fully saturated conditions. The result is a modified version of SATRA called SATRA-CHEM which is documented in Lewis et al. (1986). In this paper, however, SATRA-CHEM will be referred to simply as "the model". The methodology follows neither of the general approaches that were discussed above, but is actually a hybrid of the two. The method may be considered

84 to resemble the first approach because all nonlinear terms, which result from the particular simultaneous reaction combinations, are incorporated into coefficients of the mass-balance equation for the species that sorbs. All of the chemical equilibrium calculations, however, are performed externally to the basic solution of the PDEs. This is similar to the models based on the second method. The solution procedure that is used to solve the algebraic equations that define the aqueous chemistry is entirely separate from the equations that are used to solve for solute migration in porous media. Additionally, the coefficients that contain the nonlinear terms, which result from coupling sorptive and aqueous reactive processes, are determined explicitly and independently of the mass-balance equations. GENERAL PROBLEMDESCRIPTION, ASSUMPTIONS, AND NOTATION This paper concerns the simultaneous transport of multiple dissolved-inwater reaction participants. Specifically, these participants will be expressed in terms of tenads. In a given chemical system, a tenad is best defined as a chemical entity whose total mass is not influenced by the reaction process. In particular, the total mass of a tenad is always conserved within a given system, even if some of the substances containing the tenad are sorbed. If no substances of a particular tenad are sorbed, the definition implies that the total dissolved concentration of that tenad is transported conservatively. (A detailed explanation and specific examples of tenads in chemical systems may be found in Rubin, 1983.) In discussions that follow, tenads will be denoted as {Mi} where i = 1, 2, 3, 4, and serves to identify the particular tenad. The symbol M i refers to a dissolved chemical species composed of only one tenad {Mi} and MiMj represents a dissolved chemical species that contains two tenads {M~} and {M~} where i # j. ~ and M/Mj are the corresponding chemical species in the solid phase. The dissolved concentrations of each chemical species are designated by C~ (e.g., moles per unit volume of fluid) for Mi and C~j for MiMj.. Accordingly, adsorbed concentrations (e.g., moles per unit mass of solids) are referred to as ~ and Cij. The physical transport of each tenad is accomplished through the mechanisms of advection and hydrodynamic dispersion. These processes are represented in a general mass-balance relationship for an arbitrary tenad {M~} by:

c~dU~ + pbd5~t = L(U~) + Q(U* - U~)

(1)

where: L(U~) = V . ( ~ D . V ~ ) - ev. VU~

(2)

and e is porosity, pb is the bulk density, U~ is the total dissolved concentration of {M~} in the aquifer, ~ is the adsorbate concentration (defined above), D is

85 the hydrodynamic dispersion tensor, v is the fluid velocity, and Q is a fluid source term having a dissolved concentration of {M/} equal to Ui*. In eqn. (1) the only sorbed species involving {Mi } is/~i. The methodology that is developed in this paper will focus on solute transport in porous media involving two different chemical systems. In the first system, three tenads are considered in aqueous equilibrium-controlled reactions where one of the tenads is involved in sorption. The second system considers chemical species that involve four tenads that simultaneously interact through binary ion exchange and an aqueous equilibrium reaction. In both systems the aqueous reactions are considered to be Class I, i.e., homogeneous reactions, and the sorption/exchange reaction is considered to be a Class II, i.e., heterogeneous surface reaction (following the general classification in Rubin, 1983). As a result, these reaction systems may be termed Class I/Class II hybrid systems (Rubin, 1983). In order to facilitate the demonstration of our approach to simulating transport influenced by these reactions, the mathematics presented here represent the simplest physical conditions and relatively simple but significant chemistry. As a result, we will discuss only isothermal systems where the fluid density and viscosity are constant. These constraints, however, are not fundamental requirements. The porous medium is considered to be physically and chemically homogeneous. All of the reactions are reversible and sufficiently fast to assume that local chemical equilibrium is continually maintained. This implies that the rates of concentration changes due to reactions greatly exceed those due to the physical mass-transport process. Further assumptions involve stoichiometric and activity coefficients all of which, for simplicity, are set equal to one. Additionally, the chemical equilibrium constants (defined below) are independent of spacer and time. EQUILIBRIUM SORPTIONAND AQUEOUSREACTIONS The first reaction system we consider is defined by the following chemical equations: M1 + IVIx '

~ MIMx

(R1)

M1 + M2 ~

• M, M2

(R2)

M, + M4

~ M1M4

(R3)

'

where even subscripted participants may be viewed as anions and odd numbered ones as cations. Reaction (R1) describes the equilibrium sorption of M1 onto the surface of a reactive solid, /~x, and is governed by the functional relationship (isotherm) ~,

= f(C,)

(3)

Reactions (R2) and (R3) also involve M,, but represent the formation and dissociation of aqueous complexes of M1 with M2 and M~, respectively. Due to the local chemical equilibrium assumption, these reactions are controlled by

86 chemical equilibrium constants, which are defined through their respective chemical-relation equations: (4)

K'2 = C, C2

/(,4 -

C4

(5)

CC

Transport equation development

Simulating transport with reactions (R1)-(R3) implies that six unknowns are involved, namely the concentration of each reaction participant (C1, C,, C2, C4, C12, C14). The sorbed Ca, however, represents an immobile phase and thus is not directly transported. The value of ~1 can be computed from C, for any point in space and time by applying eqn. (3). Consequently, five equations, that account for mass transport as well as the interactions among all mobile reaction participants, are required to define the problem adequately. Due to the local equilibrium assumption, transport equations considered here have the same form as eqn. (1), i.e., the mass balances considered can be defined without kinetics-determined chemical source/sink terms. The number of transport equations is then dictated by the number of tenad. In this system the tenads are {M,}, {M2} and {M4}. The mass-balance equations, each one written for one of these tenads, are, respectively: ~C1

~C14

~C,2

~-~7 + ~ - ~ + ~-~- + ~

~CI

~t

= L(C, + C,2 + C t ) + Q(C* + C~ + C* - C - C,2 - C~4) ~C2 ~--~+~

~C2 Ot

= L(C2 + C~2) + Q(C* + C * -

0C4

~Ct

_

~-~- + ~ 0t

C2 - C~)

L(C4 + C4) + Q(C: + C* - C4 - C,4)

(6) (7) (8)

These transport equations and the two chemical-relation eqns. (4) and (5) provide the set of five equations necessary to determine the five unknowns of the system. However, even though this set contains the same number of equations as unknowns, it is difficultto obtain a solution with the equations in their present form. To facilitate the solution process, these equations may be combined to yield an alternative set of equations that can be solved sequentially. For the most part, this may be accomplished in one step by assigning a single variable to represent the mobile (dissolved) phase of each tenad. In other words, the respective Ui variables of eqn. (I) are defined for {MI}, {M2} and

{M4}: U, -- C, + C~2 + C,

(9)

87

us =

c~ + c ~

(10)

U4 = 64 + C~4

(11)

These definitions permit eqns. (6)-(8) to be simplified into: 8U~

8~i

8U2 8t

~u~

-

8t

_

L(UO + Q(U* -

U,)

(12)

L(U2) + Q(U~

Us)

(13)

L(U4) + Q(U* -

U4)

(14)

If the variable Ca could be expressed in terms of UI, our system would consist of the three mass-balance equations with three unknowns (U1, U2, U4), and five algebraic eqns. (4), (5), (9)-(11). The latter may be combined and solved simultaneously for the concentrations of the five reaction participants following the solution of the mass-balance equations for U,, Us, and U4. The transport eqns. (12)-(14) all have the same basic form as eqn. (1), but eqns. (13) and (14) represent transport for solutes that do not sorb. As a result, these equations are linear and their solutions may be obtained independently from all other equations. Equation (12), however, represents the transport of a tenad comprised of a substance that simultaneously participates in sorption and two aqueous complexations. This equation must therefore couple all three reactions together and account for the chemical processes which comprise this Class I/Class II hybrid reaction system. This coupling causes eqn. (12) to be nonlinear, and thus an iterative technique must be used to solve for U~. The complexity of the iterative approach is reduced if all of the partial derivatives in the equation can be expressed in terms of the same concentration variable. Moreover, if all of the partial derivatives in eqn. (12) are expressed with respect to the total dissolved concentration of {M~ } (i.e., U~), the same general solution approach can be applied to all three transport equations. The nonlinearities will remain, but they will be removed from all derivatives of U~ and placed into terms that can be determined explicitly at any given iteration. Transforming the adsorbate term in eqn. (12) involves two principal steps. The first step is to rewrite 8~1/(?t in terms of ~C~/(~t, and the second is to express 8C,/8t as a function of 8 U~/~t. We begin by redefining the sorbed concentration C1 in terms of the solute phase by applying the appropriate sorption model. From eqn. (3) we may obtain:

8t

-

8C, 8t

(15)

The derivative 8C~/8C~ mathematically represents the partitioning of/14, between the sorbed and solute phase, and is denoted by F. The function F is obtained from eqn. (3). The latter equation represents a variety of sorption

88 isotherms; however, a linear isotherm will be assumed in this study for simplicity. In linear sorption, F is constant and analogous to the equilibrium distribution coefficient, Kd, as defined in Freeze and Cherry (1979) and Reardon (1981) among others. The time derivative of the adsorbate concentration can now be expressed simply by:

~Cl 0t

~c~ -

F ~t

(16)

This allows eqn. (12) to be expressed entirely in terms of concentrations in the liquid phase as:

~-~U' + ~ F(~C'ot = L(U,) + Q(U* -

Vl)

(17)

We can now describe OC~/~)tas a function of chemical-relation eqns. (4) and (5) to read:

~U~/Otby

first rearranging the

C~2 = g,2 C~C2

(18)

C14 = g14 C1C4

(19)

Next, by combining eqn. (18) with (10), and (19) with (11), the following expressions for C2 and C4 are obtained:

U2

C2

(20)

(1 + K,2 C1)

g~

C4 -

(21)

(1 + Kl~C,)

The derivatives of these equations with respect to time yield:

c~C2 0t

( 1 ) =

c~U2

(1 + ~[~2C,)

~t

K12U2

c~C~

(22)

(1 + K~2C~)2 dt

and:

~C4 _ ( 1 ~t

) c~U4

(1 + K14C,)

~t

KI4U4

~C~

(23)

(1 + K,4C,)2 ~t

Similarly, introducing eqns. (18) and (19) into eqn. (9) yields: /-}1 = C~ + C~2 + C14 = C~(1 +

K12C2 + K14C4)

(24)

and:

~u~t

~c~t(1 +

oc=

K12C2 + K,4C4) + CIKa2 ~

Substituting the definitions of

~c,

+ C~K14 (~t

(~C2/Otand ~C4/~t, given

(25)

in eqns. (22) and (23),

89 into eqn. (25) results in: 0U~ _ ~ C ~ H + G ~t ~t

(26)

where: H

=

1 + K~2C2 + K,4C4 -

1 + KI2C~] U2C~ -

(

1 -~-K14C~] U~C~

(27)

and: a

=

1 + KI2C~]~

+

1 + K~4C1] ~t

Equation (26) may be rewritten explicitly for ~C1/~t as: 5C~ ~t -

1 ~U1 H ~t

V H

(29)

which, in turn, may be introduced into eqn. (17) such that the transport equation for U1 becomes: ~-+

pbF

~t

This equation may be rearranged to read: _ (E + pbk~) OV~ •t

L(U~) + Q ( U * -

[71) + pbk2

(31)

where: k~ = F / H

(32)

k2 = F G / H

(33)

Equation (31) represents the final form of the transport equation for {M1 }. All nonlinear components resulting from the coupling of the three reactions are contained in kl and k2. The variables G and H accommodate the aqueous chemical interactions, and sorption is incorporated through F. In this way the mathematical system describing solute interaction is arranged such that the mass-balance equations are linear or quasilinear in U1, and are solvable independently (by iteration, if necessary) before the variables of the chemical algebra are determined. The origins of the nonlinear components that comprise kl and k2 are clearly shown in eqns. (22)-(25). These equations also help identify mathematically how the chemical parameters F, K~2, and K~4, influence transport. The complex case occurs when the effects of both aqueous reactions (R2 and R3) are simulated. Without these reactions K12 = K~4 = 0, and eqns. (22)-(25) are linear.

90 Solution procedure Equation (31) along with eqns. (13) and (14) comprise the basic set of massbalance equations that describe the transport of total dissolved concentration of each of the three tenads in this system of sorption (R1) coupled with aqueous complexation (R2 and R3). These equations resemble the standard form of the advection-dispersion equation and hence are potentially compatible with solution algorithms used by many existing solute transport codes. These algorithms may be based on a variety of numerical methods. The model we use is based on finite elements. Values of the dependent variable are determined at discrete points (nodes) on a grid composed of several finite elements. Following solution of the mass-balance equations, the individual concentrations of the reaction participants can be determined by solving the algebraic eqns. (9)-(11), (4) and (5). The solution algorithm used to solve for the principal reaction participant concentrations, C1, C2, and C4, at specified points in space and time is summarized in Fig. 1. First, the linear transport eqns. (13) and (14) are solved at the end of a given time step (t k) for U2 and /-74, respectively. Equation (31), which is nonlinear due to the concentration-dependent variables G and H, is then solved for U1 at t k using a first-order (Picard) iteration scheme.

I

TRANSPORT U2

LINEAR EQUATIONS TRANSPORT U 4

CALCULATE GAND H

TRANSPORT U1.... 1

NONLINEAR EQUATIONS

I

SOLVEFORC1 I J

CALCULATE C 1 AND C 4

NO

NEXT TIME STEP

Fig. 1. Schematic diagram of the general solution algorithm for transport with sorption coupled with aqueous complexation.

91 The first step in solving for UI is to determine the values of G and H (eqns. (27) and (28) respectively) on a given iteration. These variables require U2 and U4 at tk, along with their derivatives with respect to time, and values of C,, C2, and C4. The derivatives are expressed simply as the linear change during the time step (e.g., ~U1/~t = (U~I - U~I-1)/At), and the values of C1, C2, and C4 are obtained from the previous iteration, or the previous time step if on the first iteration. With F known from the input data, all of the unknowns are accounted for and eqn. (31) can be solved for U,. Following the determination of U,, U2, and U4, the individual participant concentrations are obtained by subtracting eqns. (10) and (11) from (9), and then using the relationships (20) and (21) to arrive at a single equation with C~ as the only unknown: U1 -

U2 -

U4

:

C 1 -

C2 -

C4

=

C 1

U2

U4

(1 + K,2C,)

(1 + K14C1)

(34)

This equation can be converted into the following cubic expression: (KI2K~4)C~ + [K12K,4(U2 + U4 - U~) + (K,2K~4)]C~

+ [1 + KI2(U2 - U,) + K14(U, - U,)]C~ - U1 = 0

(35)

which may be solved for C1 analytically or with an iterative technique such as the method of Newton-Raphson. Values of C2 and C4 are subsequently obtained from eqns. (20) and (21). As a final step, the most recent value of U1 is checked against the value obtained in the previous iteration. If the difference between successive values of U1 is within a specified tolerance at every node, one proceeds to the next step. If not, another iteration is required and the solution procedure for U~ is repeated. When sorption and aqueous complexations do not occur simultaneously, the algorithm is simplified considerably. Specifically, without homogeneous aqueous reactions, eqns. (9)-(11) reduce to:

u, = C v2=

C~

v ~ = c4 and without sorption (either with or without homogeneous reactions), eqn. (32) becomes:

~v~ ~t

-- L(Vl) + Q(U* - UI)

(36)

Consequently, in both of these cases the resulting transport equations are linear and an iterative procedure is not required. However, in the latter case involving only complexation reactions (no sorption), the polynomial in eqn. (35) must be solved for each time step at which the concentrations of the reaction participants are desired.

92 From this outline of the general solution procedure it is obvious that by modifying only a few steps, our methodology would permit the addition of more complexation reactions involving C~. With each additional reaction in the form of R2 or R3, three adjustments must be made: (1) A linear transport equation in the form of eqn. (13) or (14) must be developed for the new tenad; (2) because U2 would be redefined, a new derivative corresponding to eqns. (25)-(28) must be determined; and (3) the polynominal in eqn. (35) must be expanded to the nth order, where n is the number of complexation reactions involving C1. EQUILIBRIUM ION EXCHANGE AND AQUEOUSREACTIONS The second reactive transport system we will examine is represented by the following simultaneous reactions in generic form:

M~ + M3Me Ma + M2

M, Me + M3 M,M.,

(R4) (R5)

Reaction (R4) describes binary exchange between M, and M3 in which a tenad in solution interchanges with another tenad in the sorbed phase. The common reactant for both tenads is the cation exchanger, ~re, which is also a tenad because its total mass does not depend on the reaction. Because interaction with a solid surface is involved (R4), like (R1), is a Class II surface heterogeneous reaction. On the other hand, the aqueous complexation (Rh), which is identical to (R2), is a Class I homogeneous reaction (Rubin, 1983). Inasmuch as these reactions are governed by the condition of local chemical equilibrium, the law of mass action yields the chemical-relation equations: g~3 K~

ClC C~3

C12

CC2

(37) (38)

for reactions (R4) and (Rh), respectively. Additionally, due to the assumption of chemically homogeneous porous media, the exchange capacity is assumed constant, and is defined as follows: CT = ~, + ~3

(39)

Transport equation development In this reaction system, the concentrations of six reaction participants are unknown (C~, C1, C2, C~2, C3, ~3), but only four of these are actually transported (C~, C2, C~2, C3). The system may also be described as containing four tenads {M,}, {/142}, {M3}, {/~e}, but only the first three are mobile. As a result, only three transport equations can be developed, even though the mass balance of

93 all four tenads must be considered. The mass-balance relationships for {M, }, {M.2}, {M3} and {/~e} are respectively: ~C1

~C,2

~ y + ~-~- +

Pb

(~t

0C2 ~Q~ L(C2 e-~-~-+~ ~t = ~C~ ~C3 ~ - - ~ - + pb Ot -

0C Pb-~-

-

L(C, + C,2 ) + Q ( C * + C * -

+ C,2) + Q ( C ~ + C* 2 -

L(C~) + Q ( C { -

C2 -

C~ -

C,2)

Ca)

C,2 )

(40) (41) (42)

~C3 ----- 0

-~- Pb ~t

(43)

These four expressions, together with relationships (37) and (38), comprise the basic set of six equations from which the six unknowns in this reaction system can be determined. As in the first reaction system, the solution to all of the unknowns, at the same point in time and space, is facilitated by transforming this set of equations to develop an alternative set that can be solved sequentially. It is desired to obtain a system of transport equations in U1 which can be solved independently. To achieve this, first, eqns. (42) and (40) are added together. This step takes advantage of the relationship in eqn. (43) and eliminates the adsorbate contribution of/t~3 from all of the mass-balance equations required for this system. The result is an expression containing three unknown solute concentrations:

~_~2C, ~--[[-+ ~--+

~

0C3 = L(C~ + C~2 + C3) + Q(C* + C~ ?t

+ C* - C, - C~2- C3)

(44)

To simplify the notation, we will express the combined dissolved concentrations associated with each of the three transport equations as a single variable. This allows eqns. (40), (41), and (44), to be written as follows:

~u~ + ph ~t

0v~

=

L(UI) + Q(U* -

U,)

(45)

= L(U~) + Q ( U ~ -

Us)

(46)

= L(U3) + Q ( U * -

U3)

(47)

~u~ where: U, =

C1 + C,2

(48)

U2 = C2 + C12

(49)

94

u~ = c , + c 2 + c ~

= v,+c~

(50)

Equations (45)-(47) now have the same form as eqns. (12)-(14) which describe the previously considered reactive transport system. Equations (46) and (47) are linear because they involve only the aqueous phase. Equation (45), on the other hand, is nonlinear due to the participation of/1//1 in both aspects of the Class I/Class II reaction hybrid. Although the transport equations have the same form as the previous case, two principal differences arise in accommodating the current set of reactions, namely the composition of the transported variables, and the conversion of the adsorbate term OCelOtin eqn. (45) into an expression that is a function of the transported concentrations/_71, U2, and/-73. In the previous case, with (R1), concentration C1 depends only on C~. Therefore, eqn. (16) can easily be applied to convert the adsorbate term into the liquid phase. But with ion exchange, concentration C1 depends on at least two variables, C~ and C3, as in eqns. (37) and (39). As a result, the steps required to eliminate ~C,/~t from eqn. (45) are slightly different than those for the case involving equilibrium sorption. The first step is to derive an explicit relationship for ~C~/~t in terms of the partial derivatives of the dissolved concentrations of the other reactants in (R4). The second step is to express these derivatives explicitly as functions of the transported variables U~, U2, and U3, taking into account the interrelationship of (R5) in the process. The details of these steps are in the Appendix. For convenience, only the final form of the transport equation for {/I//1} is presented here:

[~ + pbk~] c~U~ ~t _ L(U~) + Q(U~(- U~) + pbk2

(51)

where:

kl = (1/g)[£ + fl[(1/KI2)B + C~]]

(52)

k 2 = (l/g)

(53)

+ [~-~r-..

and: B = B(K12, UI, U2)

(54)

f, = f,(g13, CT, C, V,, V3)

(55)

f2 (K,3, ~T, C~, U,, U3)

(56)

f~=

g = g(g,3, C~, U,, U3)

(57)

The relationships (54)-(57) are fully defined in the Appendix. Equation (51) together with eqns. (46) and (47) comprise the basic set of transport equations for the system of binary ion exchange (R4) coupled with a homogeneous aqueous (complexation) reaction (R5). All nonlinear components resulting from coupling these reactions are contained in k, and k2. When ion exchange is

95

simulated without an accompanying aqueous reaction, the coefficients in eqn. (51) become: k, = (1/g)brl + f21

(58)

k2 = (llg)f2 ~t

(59)

The development of these expressions is explained in the Appendix. Note eqns. (51) and (31) are identical in form. The only difference is that kl and k2 have been redefined. In both chemical systems all solute interactions are condensed into these two terms. Solution procedure The procedure used to determine the concentration of all mobile reaction participants in this system is similar to that followed in the previous case involving equilibrium sorption. A schematic illustration depicting the solution process for C~, C2, and C3, is shown in Fig. 2. First, the linear transport eqns. (46) and (47) are solved at t~ for U2 and U4, respectively. Equation (51), which is nonlinear, is then solved for U, using a first-order (Picard) iterative process. TRANSPORT U 2 LINEAR EQUATIONS TRANSPORT U3

CALCULATE

CALCULATE k I AND k2

TRANSPORT

NONLINEAR EQUATIONS

U1

CALCULATE C2 AND C3

NO

NEXT TIME STEP

Fig. 2. Schematic diagram of the general solution algorithm for transport with ion exchange and aqueous complexation.

96 The solution to eqn. (51) requires the determination of all variables that comprise k, and k~, defined in eqns. (52), (53) and (58), (59). This may be accomplished by systematically solving the following explicit expressions in the order listed: 1.

A

=

U,-

Us-(1/K~)

2.

B

=

[A + 4(1/K~2)UI] '/2

3.

C~ = ~(A + B)

(60)

(61) (62)

4. g = K,3C~ + U3 - U1

(63)

5.

f2 = (g~3 C~ Ow)/g

(64)

6. f~ = K13(5T - f2)

(65)

In these calculations, values of U~ represent the previous iteration, or the previous time step if on the first iteration. When transport with only ion exchange is simulated, C1 = U1; and the solutions to the required explicit equations start with eqn. (63). With K,2, K~3, and CT known from the input data, eqns. (52) and (53) [or (58) and (59)] can be solved for kl and k2, and eqn. (51) is then ready to be solved for U1. If the difference between values of U, for successive iterations is greater than a specified tolerance at any node, new values of k~ and k2 are determined and the solution process for U, is repeated. If the difference is less than the tolerance, the model proceeds to the next time step. Once U, has converged, the concentrations C2 and C3 are calculated from eqns. (48) and (50). APPLICATIONS Simulating transport with Class I/Class II hybrid reaction combinations provides an opportunity to examine the combined influence of these processes on the spatial distribution of concentration for the reaction participants. For example, we can observe how local equilibrium-controlled sorption of a given chemical species influences the distribution of all other species involved in aqueous reactions at local equilibrium with that sorbing species. In addition, we are able to examine the results of varying the chemical parameters, e.g., F, K12, and K14, as well as observe the influence of mechanical dispersion on reactive transport. The results of several examples are presented in this section that demonstrate these concepts. For simplicity, all of the simulations presented involve steady-state flow fields and are representative of physical systems having dimensions that are typical of field rather than laboratory conditions. Although our model can simulate two-dimensional problems, important behaviors that occur during transient transport of these reactions are best illustrated in one dimension; therefore, we will start with examples presented in this perspective. A hypothetical two-dimensional flow system will be described later when multidimensional examples are discussed. For now, however, the region modeled may be

97 viewed as a r e p r e s e n t a t i v e s t r e a m tube w i t h i n a l a r g e r physical setting. The s t r e a m tube is 10,000m long and is discretized into 100 q u a d r i l a t e r a l finite elements, each 100 x 100m (Fig. 3). The system has a porosity (D of 0.20, and a solid-grain density (Qs) of 2.65 g cm -~. The h y d r o d y n a m i c dispersion t e n s o r (D) in eqn. (1) is r e p r e s e n t e d in o u r model by b o t h m o l e c u l a r diffusion (D M) and velocity-dependent dispersion components. The c o m p o n e n t s of D are defined in general (Bear, 1979) as:

Diy = a~lvlSij + (aL -- aT)ViVj/IVl -F D M ~ij

(66)

w h e r e aT is the t r a n s v e r s e dispersivity, a L is the l o n g i t u d i n a l dispersivity and 5ij is a K r o n e c k e r delta. In one dimension eqn. (66) reduces to: D

= aLv + DM

(67)

and in o u r simulations aL is 33.0 m. At this scale DM is negligible and is assumed to be zero. With respect to chemical parameters, all are in an a r b i t r a r y system of c o n s i s t e n t units, and are specified s e p a r a t e l y for each example. Initially, the porous medium is s a t u r a t e d with fluid t h a t c o n t a i n s zero c o n c e n t r a t i o n of all tenads, and h e n c e all r e a c t a n t species: U~(x, 0)

=

0.0 ~

C,(x, 0)

=

0.0

G ( x , O) =

o.o ~

G ( x , O) =

o.o

U,(x, o)

o.o ~

C,(x, O) =

o.o

=

The b o u n d a r i e s are defined by c o n s t a n t heads (Dirichlet), h(x, t), at both ends of the system (Fig. 3): h(0, t) =

H~

h(10000, t) = /-/2 The values of/-/1,/-/2, and the h y d r a u l i c c o n d u c t i v i t y , are chosen such t h a t the fluid v e l o c i t y (v) is c o n s t a n t at 3.865 x 10 6ms 1. The b o u n d a r y c o n d i t i o n for c o n c e n t r a t i o n is expressed as follows:

t x=xB =

~x +

x=xB

where XB is the l o c a t i o n of e i t h e r the u p s t r e a m or d o w n s t r e a m b o u n d a r y , and chemical equilibrium is assumed in the inflowing fluid. 10, 000 m

q

~.ILL~I,]Ii, !

NX

G

I~,TCT=

\ \ N\ \ \ \

I

XN

\

NN

T

Fig. 3. Finite-element mesh and boundary conditions for one-dimensional simulations.

98

The simulation commences as fluid-containing concentrations of the various reaction participants flow continuously into the system across the upstream constant head boundary. The assumption of chemical equilibrium also applies to the concentrations in the source fluid. Aqueous reactions such as R2 are assumed to occur not only in the simulated region but in the source fluid as well. The time step interval (At = 300 days) for all of the one-dimensional simulations is selected such that a particle of fluid will progress approximately one element (100m) per time step. Except where noted, all simulation results represent 1.2 x 104days (40 time steps) of elapsed time.

Equilibrium sorption and aqueous reactions Examples under this heading concern the chemical environment characterized by reactions R1-R3. The effects of these reactions may be simulated individually or in any collective combination; thus the influence of sorption and. complexation on transport may be observed separately or together.

Sorption Single species adsorption has been addressed extensively in the literature (e.g., Lai and Jurinak, 1972; Pickens and Lennox, 1976; Bear, 1979). The net effect of sorption on transport is to reduce the average pore velocity of the sorbing solute. The amount of this reduction is a function of the porosity, bulk density and the coefficient, kl, according to the relationship R = £ + Pbkl

(68)

where R is the retardation factor discussed by Lai and J u r i n a k (1972), and van Genuchten and Alves (1982) among others. In this case, without additional aqueous reactions, the calculation of R is straightforward because kl = F. The exact solution for both the sorbing and nonreactive fronts can be determined by noting that this problem may be described by rewriting eqn. (31) for a single arbitrary solute concentration, C, in one dimension to read: ~C

~2C

~C

R-zr•t = D-~--~vx.-v-~-xx

on

0 ~< x ~

10,000

subject to: C(x, o) = o ~c

- D - ~ x + vC]x~O = vCo ~C

~-~-(10,000, t) = finite

(Co = 1.0),

t > 0

(69)

99 The analytical solution to eqn. (69), as reported by van Genuchten and Alves (1982), is:

RxC(x, t) = ~erfc L-~Rt~7~jVvt]

f v2t ,~,12 ~_I (RX_4D_Rt vt)21j -

+ ~n--D-R) exp

1 (1 + vx + v2t

2_

~

~--~] exp

,x + (vx/D) erfc [L2-~D~t)~7~ j

(70)

where R = 1.0 for conservative transport. The effect of linear sorption on M1 without any accompanying aqueous reactions is shown in Fig. 4. The inflowing solution contains M1, and a nonreactive species, each at a concentration of 1.0. Figure 4 contains the analytical solutions for both fronts plotted as solid lines, and the corresponding simulation results represented by points. The results agree well with the exact solutions for both the nonreactive and sorbing distributions considering the relatively course temporal and spatial discretization.

Complexation The effect that equilibrium-controlled reactions have on the distribution of solutes and the influence that dispersion has on equilibrium-controlled reaction behavior are two interrelated processes that can be demonstrated by simulating transport with equilibrium-controlled aqueous complexation (R2 and/or R3). The results of simulations involving these reactions are shown in Figs. 5 and 6. The chemical system represented in Fig. 5 is defined only by a single reaction R2, with K12 arbitrarily set at 1.0. Inflow concentrations for the two reaction participants (C* and C~') are set at 1.0. For comparative purposes we assume that a nonreactive constituent (C~R) having a concentration of 1.0 is also introduced into the porous medium. The interrelationship between dispersion and transport with equilibriumcontrolled aqueous reactions is demonstrated in Fig. 5. With the same input concentrations, the concentration front for the reactive species M1 and M2 (C~ and C2) has higher values downstream from the source than that of the nonreactive constituent (CNR)- These reactive species have higher concentrations because at every point in the system they are in equilibrium with a compound that enters the medium with the input solution, and serves as an additional source of these species within the medium. Along the reactive solute

(M1M2) 1.0 [

* ~ k ~ \

F I O~ 0.0 I 0

Ml +Mx'~MIMx --Analytical

~

\ k ~

...... N~o_nreactive , ~

. . . . ivo, ,

2000 4000 6000 8000 10,000 DISTANCE FROM SOURCE, IN METERS

Fig. 4. Single species sorption [(R1), F = 0.1] and nonreactive transport after 1.2 × 104s.

100

2.0 Z 1.5 .............

MI +Me ~-MIM2KI2=I'O

1.0 i0.5 0.0 0

Nonreacti1 ~C~=C2 I

i

~

~

J

2000 4000 6000 8000 10,000 DISTANCEFROMSOURCE,INMETERS

Fig. 5. Transport with a single aqueous equilibrium reaction [(R2), K12 = 1.0] and a conservative solute (CNR)after 1.2 x 104s.

2.0,

~

CI

~) 1.5 ~

MI+M2 .~ MIM2 MI+M4~MIM4

Z 1.0 ~

Ktl4 5.0 :

0.5 O.C

i 2000 4000 6000 8000 10,000 DISTANCEFROMSOURCE,INMETERS

Fig. 6. Transport with two simultaneous aqueous equilibrium reactions [(R2), KI~ = 0.5, and (R3) K~4 = 5.0] after 1.2 × 104s. front, (i.e., w h e r e t h e c o n c e n t r a t i o n profile c h a n g e s from the i n p u t c o n c e n t r a tion to the a m b i e n t c o n c e n t r a t i o n ) r e d u c t i o n in the c o n c e n t r a t i o n of the r e a c t a n t s c a u s e s the r e a c t i o n (R2) to i n c r e a s e the a m o u n t of M1 and M2 t h r o u g h dissociation. T h e a c t u a l t r a n s p o r t e d q u a n t i t i e s , in t e r m s of the model, are the t o t a l dissolved t e n a d c o n c e n t r a t i o n s , /-71 = /-72, as well as CNR. T h e s e distributions, w h i c h a r e s h o w n in Fig. 5, a r e t r a n s p o r t e d c o n s e r v a t i v e l y b e c a u s e t h e y are i n d e p e n d e n t of a n y r e a c t i o n s . T h e individual c o n c e n t r a t i o n s C~ and C2, however, a r e a f u n c t i o n of r e a c t i o n (R2), and thus are not t r a n s p o r t e d conservatively. T h e solutes i n v o l v e d in the r e a c t i o n m u s t mix, i.e., m u s t be influenced by dispersion if the r e a c t i o n is to t a k e place. In a s y s t e m g o v e r n e d by p u r e a d v e c t i v e flow (no dispersion), the c o n c e n t r a t i o n d i s t r i b u t i o n s of b o t h the r e a c t i v e species a n d the n o n r e a c t i v e c o n s t i t u e n t would p r o g r e s s the s a m e d i s t a n c e into the p o r o u s medium, i.e., h a v e the s a m e s t r a i g h t v e r t i c a l front. In p u r e a d v e c t i v e s y s t e m s ( w i t h o u t decay) the c o n c e n t r a t i o n s of all t r a n s p o r t e d entities, at a n y p o i n t in s p a c e or time, are e i t h e r at the s o u r c e s t r e n g t h or zero by definition. Dispersion, therefore, e n a b l e s the c o n s e q u e n c e s of the r e a c t i o n to be seen d u r i n g t r a n s p o r t .

101

If the solute distributions are dispersion-affected, the equilibrium condition will, in turn, influence the nature of this distribution. This can be seen from the distribution of the reactive front (C~ = C2) which is different from the conservative CN~ distribution. Calculating the dispersion coefficient and the fluid velocity from the reactive distribution would yield higher values for both of these parameters; however, the dispersivity and velocity are the same for both the reactive and nonreactive distributions. Dispersivity is a scale-dependent parameter that is strictly a porous medium characteristic, and the velocity of water is similarly independent of reactions. Both parameters are assumed to be independent of concentration. If only the spatial distribution of the reactive species were used to obtain values of the dispersion coefficient, using conventional curve matching techniques, the result would be incorrect. Similarly, an incorrect value of v would also be obtained if the reactive distributions were used. The correct values of these parameters can only be determined from the distributions of conservatively transported solutes. In the system represented by Fig. 5, the conservatively transported solutes are U,, U2, and CNR.The variables U, and U2 move conservatively with the flow because the total dissolved concentrations of {M1 } and {M2} are independent of reaction R2. The individual concentrations C~ and C2, however, are a function of the reaction, and thus are not transported conservatively. These latter concentrations are determined from values of U1 and U2 through eqn. (35) which reduces to the following quadratic when K,4 = 0:

K12C~ + [1 + g12(U1 - U2)1C1 - Vl = 0

(71)

In this case, K12 = 1.0, and U~ = Us. As a result, eqn. (71) reduces to C~ + C~ - U1 = 0

(72)

Because the inflow concentrations of M~ and M2 are the same, C~ = C2; consequently, eqn. (72) also applies to C2 when substituted for C~. This nonlinear relationship yields values of C~ and C2 that are continually in equilibrium with C12 (shaded area in Fig. 5). Inasmuch as the spatial distributions of U1 and U2 as well as CNRare conservative, their values may be verified through the same analytical solution, eqn. (70), that was applied in the previous section. Additionally, from the boundary condition in this problem we know that: V1 :

U~

=

2CNR

This provides an analytical means of calculating C~ or C2. The only requirement is that U1 be known or obtained from either U2 or CNR. Once U~ is known, C~ (or C2) can be calculated using eqn. (72). Simulation results for a chemical transport system with both complexation reactions, (R2) and (R3), are shown in Fig. 6. Here, the inflow concentrations of M2 and M4 are set at 1.0, while C* is input at twice this level to maintain an electrostatic balance in the source fluid. (Note, however, the methodology does

102 not require electrical neutrality.) The equilibrium constants are arbitrarily set at K12 = 0.5 and K14 = 5.0 for (R2) and (R3), respectively. The concentration distributions in Fig. 6 indicate how the magnitude of the equilibrium constant affects the progress of the reactive fronts and the general nature of solute distribution. The participants in the reaction with the higher equilibrium constant have higher concentrations downstream. This result follows from the previous discussion of Fig. 5. A higher equilibrium constant implies the source fluid contains a higher concentration of the associated complex, given the same concentrations of reactants in both reactions. In this case, the equilibrium conditions imply that C,4 is everywhere greater than, or equal to, C,2, and that Ut is everywhere greater than, or equal to, Us. As a result, more C4 is available from dissociation (relative to C2) as both reactions (R2) and (R3) proceed to the left along their respective fronts. S o r p t i o n a n d one a q u e o u s reaction

The next few examples concern transport with linear sorption (R1) and a single aqueous complexation (R2) involving the sorbing species. The interrelationships between the two reactions and the equilibrium conditions imposed on the system cause the sorption of M1 to directly affect the resulting spatial distribution of C2. The nature of the effect depends on the value of the sorption coefficient (F), the equilibrium constant (K12) and the amount of elapsed time. In simulating transport with these reactions the inflowing species may be considered to be M,, M 6, M2, and Ms, wherein the inflow concentrations C*, C*, and C~R are all 1.0. M6 and Ms are nonreactive species present in the incoming solution accompanying M, and 3//2, respectively. As in the example shown in Fig. 5, the transport of CNR is included in these simulations for comparison purposes. The results for F -- 0.25 and K12 = 1.0 are shown in Fig. 7. The distribution of C2 is characterized by a rise and fall in concentration that peaks at almost twice the inflow level. These results are best explained by referring to the reactions (R1) and (R2). Because the system is at equilibrium, u2=c2+c12

2.0 i . . . . . . . . . . . . ~.

C2 Z

1.5

MI+Mx~ MIMx MI+M2"~-MIM2 F = 0.25

[...,

Kle=l ° Zm 1 . 0

ZOO 0.5

0.0 0

i f 2000 4000 6000 8000 10,000 DISTANCE FROM SOURCE, IN METERS

Fig. 7. Transport with equilibrium sorption [(R1), F = 0.25], one aqueous complexation [(R2), K12 = 1.0], and a nonreactive solute after 1.2 × 10~s.

103 as M~ is taken out of solution to form/~,, the reaction (R2) is forced to the left. This causes M1M2to dissociate producing both M, and/I//2. With the sorption of M 1, the system maintains equilibrium through the increased amount of M 2 in solution. Even when the inflowing Mi species represent ions, no electrical charge imbalance is produced when C2 increases and C~ decreases through sorption. Neutrality is preserved because ions M5 and Me present in the incoming solution balance the charges created. The magnitude and extent of the peak level of C2 are determined by the equilibrium condition that balances (R2), according to the degree of sorption in (R1), and the transport mechanisms of advection and dispersion. These factors combine to cause the level of C2 to peak at the point of maximum difference between the sorbing (C~) and the equivalent conservative concentration (CNR). The maximum amount of M 2 produced is limited by the total dissolved {M2 } concentration, U2, in the source fluid. For reference, the conservative distribution of U2 is included in Fig. 7. The concentrations C2 and U2 become the same at the point where C~ equals zero. Past this point, both reactions (R1) and (R2) do not occur because M1 is no longer available to react. If M, is highly sorbed, such that C1 reaches zero at a point where the conservative solute equals its inflow concentration, then C: will equal Us at the maximum possible level. If M 1 is not highly sorbed, however, then C2 will equal U2 at some point along the leading edge of the latter's concentration front (Fig. 8). The effect of progressively larger values of KlZ (for a constant value of F) on the distribution of C2 is shown in Fig. 9. A higher peak is associated with a higher K,2 because the level of C~z at equilibrium throughout the system is correspondingly higher. This causes U2, in turn, to be higher as well. Moreover, as U2 is increased, the potential limit to the peak in C2 is also increased. Additional C~2 implies that although C~ (not shown) is sorbed, it has higher values downstream. Consequently, the peak in C2 is also progressively shifted downstream with increased values of K12. In addition to being influenced by F and K,2, the peak level of C2 is time dependent. This is illustrated in Fig. 10 through successive distributions of C2 U2=C2+C12

2.o~-..... 7 - ~ F= 1.0 .'.

1.5

I F= O',O" ,o ",° ' °

o.ol 0

".:.2-,_ "-.

MI + M X ' ~ M I M x

M~+M2~-MIM2 .-. . .-

CI C

\

,"-..__,

2000 4000 6000 8000 10,000 DISTANCE FROM SOURCE, IN METERS

Fig. 8. The effect of sorption coefficient,F, on the spatial distributions of C~and C2 with K~2 = 1.0 after 1.95 × 104s.

104 2.0

Kl2 =1.25 MI+Mx~MIMx

Z 1.5 ©

MI+M2~-M~M2

F~

F=0,25 1.0

Z © 0.5

0.C

i

i

i

200 400 600 800 10,000 DISTANCE FROM SOURCE, IN METERS

F i g . 9. The effect of e q u i l i b r i u m c o n s t a n t , K,2, on the spatial d i s t r i b u t i o n of C2 w i t h F = 0.25 after 1.2 × 104s. 2.0

M I +~x_~ MIM x

,2 4/X/X \ \ \

/X/A/ \ \ \ \

.o 5

Z 0

O0

2000 4000 6000 8000 10,000 DISTANCE FROM SOURCE, IN METERS

F i g . 10. The effect of elapsed time on the d i s t r i b u t i o n of C= w i t h F = 0.25 a n d K12 = 1.0.

obtained using the same parameters employed in the example shown in Fig. 7. To conclude this reaction case, two footnotes should be considered. First, the steady-state solution for C2 in this example will show no increase in concentration above its inflow level. At sufficiently late time, C1 will reach its source concentration level at every point, and the sorption limit (as defined by the C1 boundary condition and eqn. (3)) will be attained as well. With no sorption, there is no production of M2 above its inflow concentration; thus the steady-state solution to this particular problem is the same as if no reactions or sorption occurred. Second, unlike the previous case in which only complexation is considered, the reactive processes involved in the hybrid case will be apparent in the absence of dispersion (pure advection). Sorption and two aqueous reactions Transport with all three reactions (R1)-(R3) is shown in Fig. 11. The results are for F = 0.25; K12 = 0.5; and K14 = 1.0. In addition, the source fluid contains M1, M2, and 2144 accompanied by non-reactive species M6, M3, and/145, respectively, with concentrations of C* = C* = 1.0; and C* = 2.0. In this case M1, while sorbing, also reacts with both M2 and M4. As a result, the formation of M, directly affects the distributions of C2 and C4. Both concentrations

105 2.5

MI+Mx~-MIMx /~C4

,.8

0.0

0

V~

2000

M I +M2 ~M1M 2

\t

4000

6000

8000

I0,000

D I S T A N C E F R O M S O U R C E , IN M E T E R S

Fig. 11. T r a n s p o r t w i t h e q u i l i b r i u m s o r p t i o n [(R1), F - 0.25] a n d two a q u e o u s e o m p l e x a t i o n s [(R2), K~2 - 0.5, a n d (R3), K~4 = 1.0] after 1.2 x 104 s.

increase above their inflow values, but differ as a function of their respective equilibrium constants. The distribution of C4 has the higher peak because it is associated with the larger equilibrium constant. As noted above, a higher value of K,4 implies th at there is a greater ratio of product to r e a c t a n t concentrations, and thus more M, M 4 is available to dissociate and form 341 and M 4. Sorption and aqueous reactions: some examples in two dimensions Inasmuch as the model to which we have applied our methodology can simulate flow and transport in two dimensions, the influence of (R1)-(R3) on solute migration can be demonstrated in more complex flow regimes. To demonstrate, we have selected a hypothetical two-dimensional flow field t hat we will examine through a cross-section (Fig. 12). The area simulated is 1000 m long by 20 m deep by 1 m thick, and has the following boundary conditions. Along the left edge a constant normal flux of 1.61 × 10-5 m 3s- 1representing regional flow enters the system. On the right side a constant head of 0.0m is imposed. Recharge is approximated through a constant normal flux of 1.61 x 10 -~ m 3s 1 RECHARGE (CONSTANT FLUX] OI]SOURCE

z lo~1~

8 NO FLOW 0

i 20O

400 600 DISTANCE, IN METERS

800

1000

Fig. 12. Flow field a n d g e n e r a l b o u n d a r y c o n d i t i o n s for t w o - d i m e n s i o n a l examples.

106 e n t e r i n g all across the top, r e p r e s e n t i n g an infiltration of 508 mm yr-1. Along the b o t t o m a no-flow c o n d i t i o n n o r m a l to the b o u n d a r y is prescribed. Initially, the system is s a t u r a t e d with zero c o n c e n t r a t i o n of all t r a n s p o r t e d solutes. Unless o t h e r w i s e specified, these solutes e n t e r the system t h r o u g h a 200m designated " c o n t a m i n a n t s o u r c e " a r e a in the top left corner. The b o u n d a r y c o n d i t i o n at this s o u r c e area is subject to the same prescribed flux as the r e c h a r g e on the u p p e r b o u n d a r y . The physical p a r a m e t e r s are the same as those used in the one-dimensional examples with the following exceptions. The dispersivities, a L and aT, are 25 m and 0.01 m, respectively; each element has the dimensions of 100 m h o r i z o n t a l l y and 0.5 m vertically; and h y d r a u l i c conductivity is isotropic with a value of 1.0 × 10 4ms -1. All of the two-dimensional results are shown a f t e r 1.6 x 103 days (40 time steps with At -- 40 days). The two-dimensional t r a n s p o r t of two r e a c t i v e species p a r t i c i p a t i n g in a c o m p l e x a t i o n r e a c t i o n (R2), are compared with the t r a n s p o r t of a n o n r e a c t i v e species in Fig. 13. The s o u r c e c o n t a i n s equal c o n c e n t r a t i o n s (1.0) of MI and M2. Both the r e c h a r g e across the top (excluding the source area) and the flow t h a t enters from the left side c o n t a i n zero c o n c e n t r a t i o n s for the r e a c t i v e and n o n r e a c t i v e solutes. The results in Fig. 13 are with K12 = 1.0. Because C~ = C2 both at the s o u r c e and initially in the simulated domain, t h e i r r e s u l t i n g spatial distributions are identical for all t > 0. As d e m o n s t r a t e d above (Fig. 5), equilibrium conditions cause the c o n c e n t r a t i o n s of the individual r e a c t i v e c o n s t i t u e n t s ( Q and C2) to gain from the dissociation of M1M2. Thus the distributions of C~ and C2 progress f a r t h e r t h a n those of n o n r e a c t i v e species whose s o u r c e c o n c e n t r a t i o n s are the same as C* and C* because of the p r e s e n c e of C*. However, the total dissolved t e n a d c o n c e n t r a t i o n s , U1 and U2 (not shown), t r a v e l c o n s e r v a t i v e l y and have the same distributions as a n o n r e a c t i v e solute. The n e x t example combines s o r p t i o n (R1) with a single c o m p l e x a t i o n (R2) in a similar m a n n e r to the problem, the results of which are shown in Fig. 7. A n o n r e a c t i v e solute is t r a n s p o r t e d to indicate e q u i v a l e n t c o n s e r v a t i v e distributions. The simulations are c o n t r o l l e d by the following conditions: the inflow c o n c e n t r a t i o n s of the r e a c t a n t s (M1 and M2), as well as the c o n s e r v a t i v e solute, are 1.0, F = 0.25, and K12 = 1.0. The results are shown s e p a r a t e l y for each

5

M~ + M e ~ M ~ M 2 KI2 =0 5 -

-

......

CI,C 2 Nonreactive

20 L

, 0

200

400

i 600

,

~ 800

i

J 1000

D I S T A N C E , IN M E T E R S

Fig. 13. Two-dimensional transport with one aqueous complexation [(R2), K12 = 1.0] and a nonreactive solute.

107

constituent in Fig. 14. All three diagrams in the figure are truncated at a depth of 10m for convenience. The bottom of each diagram does not, therefore, represent a prescribed boundary to flow. The distribution of the sorbing C~ is shown in the top diagram (Fig. 14A). The influence of F on the extent of C1 retardation can best be seen by comparing these results with the nonreactive distribution shown in Fig. 14C. Although it is not obvious from these figures, the distribution of C~ in Fig. 14A is ahead of where it would be, had the complexation reaction (R2) not been included. In the present simulation 3//1, while sorbing, is also being created by the dissociation of M~Ms. The simultaneous occurrence of both reactions (R1) and (R2) has the most noticeable influence on the transport of Ms. Figure 14B illustrates how the sorption of M1 in (R1) forces equilibrium to be met in (R2) through increased Cs. As in Fig. 7, the increase in C2 reaches levels higher than the input concentrations. The extent and magnitude of this increase is a function of the particular values of F and Kls chosen for the simulation. The effect of these parameters on the distribution of C2 is the same as that discussed previously (Fig. 7). The last example in a two-dimensional flow field involves transport with all three reactions (R1)-(R3). For this case, the concentrations of the three reactants are input somewhat differently than before. Here M~ and M4 enter the system all across the top through recharge with C* = C* = 1.0. Additionally, the source area in the upper left corner supplies/1//2 and additional M1 in equal concentrations (i.e., C* = C* = 1.0). The result is that C* has a value of 2.0 in

o:

CA}

C:

i0

(B)

E-

Z

E- I

(c)

1

i 0

i 200

L

i 400

DISTANCE,

i

L 600

,

i 800

i

i 1000

IN METERS

Fig. 14. T r a n s p o r t with (R1) and (R2); F = 0.25, Kl2 = 1.0. A. Two-dimensional distribution of C~. B. Two-dimensional distribution of C2. C. Two-dimensional distribution of n o n r e a c t i v e solute.

108 the source area and 1.0 elsewhere across the top; C* is 1.0 at the source area and zero elsewhere; and C* is 1.0 all across the top including the source area. The chemical parameters employed in this case are as follows: F = 0.25, K,2 = 0.5 and K,4 = 1.0. The results are shown in Figs. 15A-C. The top diagram (Fig. 15A) contains the sorbing C~ distribution. Lines of equal concentration form a bulge to the left because a higher value of C~ enters the system at the source than elsewhere across the top. Toward the right, away from the source, equal concentration lines gradually become parallel to the top boundary. The distribution of C2 is shown in Fig. 15B. The choice o f F and K12, coupled with the selection of the initial and boundary conditions for C2, in this example result in the same distributions as those shown in the previous case in Fig. 14B. Although C* is doubled at the source (from 1.0 to 2.0), K12 is halved (from 1.0 to 0.5). The value of C~2, therefore, remains the same (from eqn. 18), and C2 similarly does not change from one case to the other. The distribution of C4 is shown in Fig. 15C. The choice o f F and K,4, coupled with the selection of the initial and boundary conditions for C4, in this example result in a two-dimensional version of the type of distributions shown for the same chemical system in Fig. 1lB.

Equilibrium ion exchange and aqueous reactions The influence of the hybrid reaction system (R4)-(R5) on solute transport may also be demonstrated with a few examples. We will use the same physical parameters and boundary conditions as applied in the previous one-dimensional simulations. The only differences are the reaction participants, the nature of the reactions, and the initial conditions. Here, the porous medium is initially saturated with fluid containing the following r e a c t a n t concentrations: C~ = 0.0, C2 = 0.0, and C3 = 1.0. The electrostatic balance is assumed to be maintained by the presence of an additional nonreactive species. In the source fluid, concentrations are: C* = 1.0, C~ = 1.0 and C~ = 0.0.

Ion exchange Concentration distributions generated by the model for transport with binary ion exchange (R4) are shown in Fig. 16. The results are after 1.2 × 104 days with K,3 = 1.0 and ~-w = 0.02. The distribution of a nonreactive solute is included in the figure for comparison. Inasmuch as Ca is constant in the inflow solution and equal to the level of Ca initially present in the porous medium, the sum of these concentrations remains constant at every point and for all time (i.e., CT = Ca + C3). This relationship is valid because the total normality of the dissolved ionic species cannot be altered by the exchange process alone (Rubin and James, 1973). Fig. 16 shows that as C1 increases, C3 decreases equally so that their sum remains constant and equal to 1.0.

109

b-

INFLOW CONCENTRATIONS 2.0 1.0

(A)

10

(B)

0.0

0

Z

(C)

1.0

r~

i[

1

. . . . .

0

k

200

h

400 600 DISTANCE, IN METERS

i

800

i 1000

Fig. 15. Transport with (R1), (R2) and (R3); F = 0.25, Ku = 0.5, K14 = 1.0. A. Two-dimensional distribution of C1. B. Two-dimensional distribution of C2. C. Two-dimensional distribution of C4.

o.s

X",

k 8ooi

, 0

/ \ \ ,J, ~"~ac? 2000

4000

6000

<,=;oo e , 8000

10,000

D I S T A N C E F R O M S O U R C E , IN M E T E R S

Fig. 16. Transport with ion exchange [(R4); Kla = 1.0, 0 T = 0.02]. Ion exchange and one aqueous reaction Results of transport w h e n both reactions (R4) and (R5) occur s i m u l t a n e o u s l y are s h o w in Fig. 17. The distributions represent solute c o n c e n t r a t i o n s after 1.4 × 104 days w h e n Ow ---- 0.2 and K12 = /£13 = 1.0. As in the similar case of linear sorption, both C2 and Ca increase above their respective boundary conditions. Figure 17 illustrates the pattern of sorption and dissociation characteristic of Class I/Class II hybrid reaction systems. The sorption of M1 t h r o u g h the e x c h a n g e of Ma in (R4), leads to the production of both M1 and M2, the latter through the dissociation of M1 M2 in (R5). Addition of M1 leads to an increase in the proportion of Ow comprised of C~; c o n s e q u e n t l y , more Ma is in solution t h a n would be in the absence of (R5). N o t e the rise in Ca above its background c o n c e n t r a t i o n is due to the total dissolved a m o u n t of {M 1} in the system. In this

110 2.0[-

MI+MzMe-~MIMe+M;

]

MIM2-~MIM2

C T = 0.20 K22=1.00 K23= i. O0 J

i 0

2000

4000

6000

8000

10,000

DISTANCE FROM SOURCE. [N METERS Fig. 17. T r a n s p o r t with ion e x c h a n g e (R4) a n d one c o m p l e x a t i o n (R5) (K12 = 1.0, K,a = 1.0, (~, = 0.02).

case, the additional source of {/14,} is the M~M.2 complex; however, a distribution of Ca similar to that shown in Fig. 17 would also result if ion exchange (R4) occurred alone, and the inflow concentration of M, were higher than the initial background level of Ca (i.e., CT # constant). The important point is that with ~r constant, an increase in 3/1, regardless of origin, implies the ratio of C1 to O~r must increase. This, in turn, leads to an increase of M3 in solution. Therefore, the cause of the Ca peak is not the same as that for C2. The peak in the latter concentration distribution remains at the point of maximum difference between the exchanging (sorbing) front and an equivalent conservative distribution. However, for both 3/2 and M3, the peak concentrations are governed by the equilibrium constraints, as dictated by the values of K,2 and K,a, and the exchange capacity of the media. DISCUSSION AND CONCLUSIONS

The results of example simulations involving the various reactions that comprise the two hybrid chemical systems presented here, demonstrate some of the fundamental processes that occur during reactive solute transport. Understanding these processes, where only a few simple reactions are considered, is essential to comprehending the behavior of reactive solutes in more complex chemical environments. This understanding begins by identifying the interrelationship between the nature and consequence of the reaction and the hydrodynamics of the flow field. "Nature of the reaction" refers to two different kinds of reactions, sorptive (surface) and aqueous (homogeneous) reactions. "Consequence of the reaction" refers to the effect of the reactive process on the distribution of concentration during transport. The consequence of the reaction is, in turn, directly influenced by the magnitude of the chemical parameter (e.g., F, K12, K~a, K14) that controls it. If transport with only equilibrium aqueous complexation is considered, the

111 concentration front must be dispersed (i.e., mixing must occur) for the reactive process to affect the distribution of the front. With any deviation from pure advective flow, the resulting front depends on the equilibrium constant. The distribution of the reactive front then is not only a function of dispersion (i.e., dispersivity of the media and fluid velocity), but of the equilibrium chemistry as well. When a sorptive reaction involving one of the complexed species occurs in the system, the influence of the reaction process on the distribution of concentrations is apparent with or without dispersion. Values of dispersivity and true average fluid velocity determined from reactive concentration fronts without taking into account the influence of reactions may then be in error, the magnitude of which will depend on the governing chemical parameters (e.g., F and K12). In addition to providing some insight into the relationship between dispersion and equilibrium chemistry in reactive transport, analysis of example simulation results enables the following conclusions to be drawn. (1) The concentration of an individual (free) reactant chemical species (Mi) participating in a single aqueous, equilibrium-controlled, complexation reaction will be higher along the solute front (i.e., where the concentration profile changes from the input concentration to the ambient concentration) than a nonreactive species influenced by the same initial and boundary conditions. This occurs when the combined influence of boundary and initial conditions causes the complex to act as a source of ions due to the direction of the reaction. Conversely, the free ion concentrations along the front would be lower when the complex acts as a sink. Additionally, the difference in concentration between the two species is a function of the equilibrium constant for the particular reaction. The total dissolved tenad concentration (U~), however, travels conservatively and has the same kind of spatial distribution as a nonreactive solute. (2) Given two chemical species that participate in aqueous reactions under conditions of local chemical equilibrium, if one of these species is concurrently removed from solution through either sorption or ion exchange, the other will show a corresponding increase in concentration as dictated by the equilibrium condition. The concentration increase, which may temporarily exceed the boundary conditions specified for that species, is a function of the sorption coefficient (or the exchange properties of the medium), the equilibrium constant of the aqueous reaction, and the duration of the transport process. The maximum level of increase is governed by the total dissolved tenad concentration for the increasing reactant species. (3) If another aqueous equilibrium reaction containing the same sorbing species is added to the system described in (2) above, the resulting concentration increase will be the greatest for the non-sorbing species participating in the reaction with the highest equilibrium constant. The simulations on which these conclusions are based were made by applying our methodology for multisolute reactive transport to an existing groundwater flow/solute transport code. The code is modified to solve sequentially

112

three (with the potential for more) transport equations, followed by a set of algebraic expressions to obtain the concentration of all reaction participants at any given point in space and time. In this way, the basic solution algorithm for the transport equations is left intact and is minimally dependent on the chemistry. The compatibility of the methodology with existing transport codes enables the simulation of reactive transport in complex multidimensional flow regimes, and provides a means of easily adapting codes to account for a few simple, but fundamental, types of chemical reactions. APPENDIX The steps involved in c o n v e r t i n g t h e adsorbate term O~,/~t in eqn. (45) into a form defined explicitly in terms of the t r a n s p o r t e d variables are outlined as follows. (Note, equation n u m b e r s n o t prefaced by t h e letter A refer to equations stated previously in the text of t h e paper.) To express ~(~ ~Or in terms of the partial derivatives of the dissolved c o n c e n t r a t i o n s of the o t h e r r e a c t a n t s in (R4), we start with the chemical-relation eqn. (37) w h i c h may be differentiated as follows:

P~t

(K~3C'C3) = ~t (C,C~)

(A.1)

R e a r r a n g i n g this expression and applying the product rule of differentiation yields: K,3 C, - ~ - + K,3 C3 - - ~ - (~1 - ~

- C3 --Ot = 0

(A.2)

From eqn. (43), OC3/~t may be expressed in terms of 0C 1~or. This enables common coefficients to be combined such t h a t eqn. (A.2) becomes: ~C1 = K13(~3 0C1 (K, aC, + C3)--~ c~--~ - C, ~~C3

(A.3)

or: aq

~C37

1 f, GE_A_TEj

g

~t

(A.4)

where: g

= K,3C , + C3

(A.5)

fl = K~3C3

(A.6)

f2

(A.7)

C~

To be of use, however, t h e s e coefficients need to be expressed in terms of k n o w n or explicitly defined variables. Accordingly, eqn. (50) t o g e t h e r with (37) and (39) are used to replace C3, C3, and C1, such that: g = K13C1 + U~ - U1

(A.8)

f~ = (K,3C, Cr)/g

(A.9)

fl = K,a(~r - f2)

(A.10)

With the first step completed we now focus on arriving at explicit expressions to substitute for

113 t h e derivatives of the dissolved c o n c e n t r a t i o n in eqn. (A.4). The derivative of C3 is easily defined; from eqn. (50), it follows that: ~C3 9t

~U3 9t

8U~ ~t

(A.11)

Arriving at an explicit expression for PC 1/~t is slightly more involved. A r e l a t i o n s h i p for C 1in terms of the combined variables U1 and U2 must first be developed prior to d e t e r m i n i n g the differential. By combining eqns. (38) and (48), C 1 may be defined as:

u1

C~

(A.12)

(1 + K12C2)

Similarly, using eqns. (38) and (49):

u,

C~

(A.13)

(1 + K,2 C 1)

E q u a t i o n (A.12) can now be i n t r o d u c e d into eqn. (A.13) for C2, and t h e result expanded into the following quadratic equation:

K~2C~ + (1 + K,~U2 - K~2U~)Q - /_/1 = 0

(A.14)

w h i c h may be solved directly for C~. Applying t h e quadratic formula yields: C~ = 0.5[/./1 - U2 - ( l / g u ) ] + 0.5{[U~ - U2 - (lIKe2)] 2 + 4(1/K~2)U~]}~12 = 0.5A + 0.5B

= 0.5(A + B)

(A.15)

where: A

=

B =

[/1

[/2 - (1/K12)

(A.16)

[A 2 + 4(1/K~2)U1]'!2

(A.17)

The time derivative of C~ required for eqn. (A.4) is best obtained from t h e quadratic equation. Accordingly, differentiating (A.14) with respect to time results in: 2C, - - ~ + [(1/K12 ) + U2 - U,] - ~ -

~k 8t

8t /

(1/K~2)-~- = 0

(A.18)

2 C, ~U ~---~

(A.19)

or:

[2C, + (1/K~2) + U2 - U 1 ] ~C~ -~ =

[(1/K,2) + C ~ ] ~U~ -~-

This can be simplified in terms of B as: aQ B-~

=

~U1 aU2 [(1/K12) + C ~ ] ~ - C~-

(A.20)

which implies: aCI

8t

1

I~U1

B [(l/K'2) + C,] at

C1 ~ U 2

B ~t

(A.21)

By s u b s t i t u t i n g eqns. (A.11) and (A.21) for 8C3/~t and 9C~/St respectively, eqn. (A.4) may be r e w r i t t e n as: +

114

This latter expression may, in turn, be substituted into eqn. (42), and rearranged into: [~ + pbk~] (--~- = L(U1) +

Q(U*- U~) + pbk2

(A.23)

where:

k~ = (1/g) [f2 + f~[(1/K~2) + C~!] B

k= =

(l/g)

÷ r=

j

(A.24)

(A25)

Eqn. (A.23) together with eqns. (46) and (47) comprise the basic set of transport equations for the system of binary ion exchange together with a homogeneous aqueous (complexation) reaction. When binary ion exchange occurs without an accompanying aqueous reaction, the development must be adjusted to allow K12 to equal zero. Currently, k~ approaches infinity when K12 = 0 (eqn. A.24). In chemical systems defined only by (R4), U~ = C1 (eqn. 45); therefore, the quadratic eqn. (A.14) is not required, and the variables A and B, both of which contain the term (1/K~2), are not defined. Moreover, it follows that: ~t

~t

(A.26)

thus eqn. (A.21) is also unnecessary. As a result, when K~2 = 0, the coefficients (A.24) and (A.25) become: k~ = (1/g)[f~ + f2]

(A.27)

k2 = (1/g)f2 3-~

(A.28)

REFERENCES Bear, J., 1979. Hydraulics of Groundwater. McGraw-Hill, New York, N.Y., 567 pp. Cederberg, G.A., Street, R. and Leckie, J.O., 1985. TRANQL: A groundwater mass transport and equilibrium chemistry model for multi-component systems. Water Resour. Res., 21(8): 1095-1104. Freeze, A.R. and Cherry, J.A., 1979. Groundwater. Prentice-Hall, Englewood Cliffs, N.J., 604 pp. Grove, D.B, and Wood, W.W., 1979, Prediction and field verification of subsurface-water quality changes during artificial recharge, Lubbock, Texas. Groundwater, 17(3): 250-257. Jennings, A.A., Kirkner, D.K. and Theis, T.L., 1982. Multicomponent equilibrium chemistry in groundwater quality models. Water Resour. Res., 18(4): 1089-1096. Kirkner, D.J., Theis, T.L. and Jennings, A.A., 1984. Multicomponent solute transport with sorption and soluble complexation. Adv. Water Resour., 7: 120-125. Lai, S., and Jurinak, J.J., 1972. Cation adsorption in one-dimensional flow through soils: a numerical solution. Water Resour. Res., 8(1): 99-107. Lewis, F.M., Voss, C.I. and Rubin, J., 1986. Numerical simulation of advective~ispersive multisolute transport with sorption, ion exchange and equilibrium chemistry. U.S. Geol. Surv. Water Resour. Invest. Rep. No. 86.4022, 50 pp. Miller, C.W. and Benson, L.V., 1983. Simulation of solute transport in a chemically reactive heterogeneous system: model development and application. Water Resour. Res., 19(2): 381-391.

115 Narasimhan, T.N., White, A.F. and Tokunaga, T., 1984. Coupled hydrologic and multiple species chemical transport with precipitation and dissolution. EOS Trans., Am. Geophys. Union, 65(45): 889. Parkhurst, D.L., Thorstenson, D.C. and Plummer, L.N., 1982. PHREEQE - - A computer program for geochemical calculations. U.S. Geol. Surv. Water Resour. Invest. 80-96, 210 pp. Pickens, J.F. and Lennox, W., 1976. Numerical simulation of waste movement in steady groundwater flow systems. Water Resour. Res., 12(2): 171-180. Plummer, L.N., Jones, B.F. and Truesdell, A.H., 1976. WATEQF - - A FORTRAN IV version of WATEQ, a computer program for calculating chemical equilibrium of natural waters. U.S. Geol. Surv. Water Resour. Invest. 76-13, 63 pp. Reardon, E.J., 1981. Kd's - - can they be used to describe reversible ion sorption reactions in c o n t a m i n a n t migration? Groundwater, 19(3): 279-286. Rubin, J., 1983. 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. Rubin, J. and James, R.V., 1973. Dispersion-affected transport of reacting solutes in saturated porous media: Galerkin method applied to equilibrium-controlled exchange in unidirectional steady water flow. Water Resour. Res., 9(5): 1332-1356. Valocchi, A.J., Street, R.L. and Roberts, P.V., 1981. Transport of ion-exchanging solutes in groundwater: chromatographic theory and field simulation. Water Resour. Res., 17(5): 1517-1527. Van Genuchten, M. and Alves, W.J., 1982. Analytical solutions of the one-dimensional convect i v ~ i s p e r s i v e solute transport equation. USDA Tech. Bull., No. 1661. Voss, C.I., 1984. SUTRA: A finite-element simulation model for saturated-unsaturated, fluiddensity-dependent ground-water flow with energy transport or chemically-reactive singlespecies solute transport. U.S. Geol. Surv. Water Resour. Invest. Rep., 84-4369, 409 pp.