6.9 Geochemical TransportJIÌÍ ŠIMæNEK, USDA-ARS, George E. Brown, Jr. Salinity Laboratory,
Riverside, California
ALBERT J. VALOCCHI, University of Illinois, Urbana, Illinois

6.9.1 Introduction
Chemicals migrating from industrial and municipal waste disposal sites, or intentionally applied to agricultural lands, are often subjected to a large number of simultaneous physical, chemical, and biological processes, including sorption–desorption, volatilization, precipitation–dissolution, complexation, and biodegradation.
Sections 6.2 through 6.5 considered the transport of only one chemical species andassumed that the behavior of this solute is independent of that of other species that
may be present in the soil solution. In reality, the soil solution is always a mixtureof many ions that may interact, create complexed species, precipitate, dissolve,
and/or compete with each other for sorption sites (van Genuchten & Šimçnek, 1996).
Many important environmental problems require analysis of the coupled transportand reaction of multiple chemical species; a few examples include acid minedrainage (Walter et al., 1994; Lichtner, 1996), radionuclide transport (Viswanathan
et al., 1998), fate of metal-organic mixed waste (Rittmann & VanBriesen, 1996;
VanBriesen, 1998), analysis of redox zone formation in organic-contaminatedaquifers (Abrams et al., 1998; Essaid et al., 1995), and reactive permeable barriers
for aquifer remediation (Fryar & Schwartz, 1994). There are also important soil science applications that require consideration of multiple chemical species, for example, evaluation of water management practices and irrigation techniques underarid and semiarid conditions, evaluation of water suitability for irrigation, andreclamation of sodic soils (Rhoades, 1997, Šimçnek & Suarez, 1997). This section
considers such complex problems that lead to multicomponent transport.
While many models have been developed to quantify solute transport insoils, they usually consider only one solute and thus severely simplify the different chemical interactions. For example, the complex processes of adsorption andcation exchange are often accounted for by invoking simplified linear (Huyakorn
et al., 1991) or nonlinear Freundlich isotherms (Yeh & Huff, 1985; Šimçnek & van
Genuchten, 1994; Liu & Narasimhan, 1989a, b) that lump all reactions between thesolid and liquid phases into a single distribution coefficient KD, either with or without a Freundlich type exponent. Other processes, such as precipitation, biodegradation, or volatilization, are frequently simulated by means of simple first- or zeroorder rate equations. While several other models have been developed to simulate
1511

1512

CHAPTER 6

multiple solutes, they focus primarily upon sequential first-order decay reactions(Gureghian, 1981; Wagenet & Hutson, 1987; Šimçnek & van Genuchten, 1994).
Aqueous equilibrium chemistry models, such as WATEQ (Truesdell & Jones,
1974), MINEQL (Westall et al., 1976), GEOCHEM (Sposito & Mattigod, 1980),
FITEQL (Westall, 1982), PHREEQE (Parkhurst et al., 1980; Parkhurst & Appelo,1999), MINTEQ (Felmy et al., 1984), GMIN (Felmy, 1990), MINTEQA2 (Allison et al., 1991), EQ3NR (Wolery, 1992), and The Geochemist’s Workbench(Bethke, 1994, 1996), which solve for solution equilibrium without considerationof transport, were originally developed independently of models for water flow andsolute transport. The most widely used equilibrium models have been recently reviewed by Mattigod and Zachara (1996). Only during the last two decades has there
been a significant effort to couple these models with transport models.
Recent reviews of the development of hydrogeochemical transport models involving reactive multiple components are given by Kirkner and Reeves (1988), Yeh
and Tripathi (1989), Mangold and Chin-Fu Tsang (1991), Lichtner (1996), andSteefel and MacQuarrie (1996). Kirkner and Reeves (1988) and Rubin (1983) presented an analysis of several methods for approximate solution of multicomponent
transport with homogeneous and heterogeneous chemical reactions, and discussed
how the nature of the chemistry may affect the choice of the numerical formulation and solution algorithm. Yeh and Tripathi (1989) provided a critical review ofmany computational methods that have been presented in the hydrologic literaturefor solving multicomponent, equilibrium-controlled transport.
Most modeling efforts involving multicomponent transport have thus far focused on the saturated zone, where changes in the water flow velocity, temperature,and pH are often much more gradual, and hence less important, than in the unsaturated zone. Consequently, most multicomponent transport models assume one- or
two-dimensional steady-state saturated water flow with a fixed value of the flow
velocity, temperature, and pH (Valocchi et al., 1981; Walsh et al., 1984; Cederberg
et al., 1985; Bryant et al., 1986; among many others). Only recently have several
multicomponent transport models been published which also consider variably saturated flow; these include DYNAMIX (Liu & Narasimhan, 1989a, b), HYDROGEOCHEM (Yeh & Tripathi, 1990), UNSATCHEM-2D (Šimçnek & Suarez,
1994), FEHM (Zyvoloski et al., 1997), MULTIFLO (Lichtner & Seth, 1996), HYDROBIOGEOCHEM (Yeh et al., 1998), and FLOTRAN (Lichtner, 2000). For example, the model DYNAMIX of Narasimhan et al. (1986) and Liu and Narasimhan(1989a) couples the chemical equilibrium submodel PHREEQE with an integrated
finite difference model for fluid flow in variably saturated porous media. Yeh andTripathi (1991) developed and demonstrated a two-dimensional finite element hydrogeochemical transport model, HYDROGEOCHEM, for simulating transport ofreactive multispecies chemicals. Yeh et al. (1998) have recently developed the moregeneral model HYDROBIOGEOCHEM that includes microbiological reactions, inaddition to geochemical reactions.
Geochemical models can be divided into two major groups: those with specific chemistry and general models. Models with specific chemistry are limited in
the number of species they can handle, and their application is quite restricted toproblems having a prescribed chemical system. They are, however, much easier to
use and can be more computationally efficient than general models. Typical ex-

MISCIBLE SOLUTE TRANSPORT

1513

amples of models with specified chemistry are those simulating the transport ofmajor ions, such as LEACHM (Wagenet & Hutson, 1987) and UNSATCHEM(Šimçnek & Suarez, 1994; Šimçnek et al., 1996), or reclamation models (Dutt etal., 1972; Tanji et al., 1972) (see Section 6.9.5). Models with generalized chemistry(DYNAMIX, HYDROGEOCHEM, MULTIFLO, FLOTRAN, and OS3D/GIMRT[Steefel & Yabusaki, 1996]) provide users with much more freedom in designing
their particular chemical systems; possible applications of these models are alsomuch wider.
In the following sections we will give an overview of the mathematical equations representing the major chemical reactions and governing transport processes,
briefly discuss how these equations can be implemented in reaction multispecies
transport models, and describe several possible applications relevant to soil scientists.

6.9.2 Geochemical Reaction Equations
Reversible chemical reactions may be described conceptually using either akinetic or equilibrium approach. The kinetic approach entails describing the temporal progress of the reaction and is usually based upon some model of the elementary reaction mechanisms. Mathematically this leads to a system of ordinary
differential equations. On the other hand, the equilibrium approach uses thermodynamic principles to describe the equilibrium state of the reaction; this approach
does not require a model for the reaction mechanism and mathematically leads toa system of nonlinear algebraic equations. Since groundwater and soil water are open
systems that are subject to fluxes of mass and energy, it is more appropriate to adopt
the concepts of “partial” and “local” equilibrium where thermodynamic equilibrium
applies to certain subsets of reactions (Knapp, 1989; Bethke, 1996; Lichtner, 1996).
Several investigators have studied the question of when local equilibrium conditions
apply in a dynamic transport system. For example, Valocchi (1985, 1988) and Jennings (1987) studied adsorption processes. General reactions, including mineral precipitation and dissolution, were studied by Knapp (1989), Bahr (1990). and Lichtner (1996). Dissolution of nonaqueous phase liquids was studied by Miller et al.
(1990) and Powers et al. (1991). In general, these studies show that local equilibrium conditions apply when the chemical reaction rates are much faster than therates of advection and dispersion processes. Lichtner (1993) demonstrated that local
equilibrium applies asymptotically, that is, after the system has evolved for a long(infinite) time, regardless of the values of the kinetic rate constants.
The term component is generally used when dealing with chemical equilibrium systems, while the term primary species is used for systems in partial equilibrium. Following Westall et al. (1976), components may be defined as a set of linearly independent chemical entities such that every species in the system can be
uniquely represented as the product of a reaction involving only these components.
No component can be represented as the product of a reaction involving only theother components. As a typical example, the chemical species CaCO30
Ca2+ + CO32− º CaCO30[6.9–1]

1514

CHAPTER 6

consists of the two components Ca2+ and CO32−. However, from the definition of acomponent it follows that any other combination of two species selected from theabove system as basis components, such as CaCO30 and Ca2+, or CaCO30 and CO32−,
may be used also, since the third species can always be represented by a linear combination of the two selected components.
Reversible chemical reaction processes are generally represented mathematically using mass action laws that relate thermodynamic equilibrium constants
to activities (the thermodynamic effective concentration) of the reactants and products (Mangold & Tsang, 1991; Bethke, 1996; Appelo & Postma, 1993). For example, the reaction
bB + cC º dD + eE

[6.9–2]where b and c are the number of moles of substances B and C that react to yield dand e moles of products D and E, is represented at equilibrium by the law of massactionK = (adDaEe )/(abBaCc )[6.9–3]where K is a temperature-dependent thermodynamic equilibrium constant, and ai
is the ion activity, being defined as the product of the activity coefficient (γi) andthe ion molality (mi), that is, ai = γimi. Single-ion activity coefficients may be calculated using either the Davies equation, an extended version of the Debye–Hückel
equation (Truesdell & Jones, 1974), or by means of Pitzer expressions (Pitzer, 1979)(see Section 6.9.2.6). Equation [6.9–3] can be used to describe all of the major chemical processes, such as aqueous complexation, sorption, precipitation–dissolution,
and acid–base and redox reactions, provided that the local chemical equilibrium assumption is valid. When Eq. [6.9–3] is used to describe heterogeneous reactions
involving solids and/or gases, then appropriate definitions for activities of solid or
gas phases are used rather than ion activities. However, mineral reactions involving precipitation–dissolution are extremely slow, as discussed below, and usuallyrequire a kinetic formulation. Acid–base reactions are usually formulated by using
the hydrogen ion H+ as a component (Yeh & Tripathi, 1989; Mangold & Tsang,1991). The essential difference between the proton and other components lies in
the meaning of the total aqueous concentration, which is defined with respect to auser-defined reference level. Different choices of the reference level can lead to anelectroneutrality condition (Yeh & Tripathi, 1989; Šimçnek & Suarez, 1994) or a
proton condition (Yeh & Tripathi, 1989; Rittmann & VanBriesen, 1996; Steefle &MacQuarrie, 1996; Tebes-Stevens et al., 1998). Redox reactions are often slow, as
exhibited by redox disequilibria in natural systems. Nevertheless, several models
treat redox equilibria by including the electron e− as a component. Since the electron doesn’t exist as a separate entity, these models invoke the notion of the “operational electron” as a conservative property that defines the redox state of the system (Yeh & Tripathi, 1989; Engesgaard & Kipp, 1992; Liu & Narasimhan, 1989a;
Keating & Bahr, 1998). A more detailed analysis of the role of the electron is given
by Lichtner (1998).

MISCIBLE SOLUTE TRANSPORT

1515

6.9.2.1 Complexation
Equations for the aqueous complexation reaction [6.9–1] can be obtainedusing the law of mass action as follows:or

aCaCO30
KCaCO30 = _____________
aCa2+ aCO32−
KCaCO30
x CaCO30 = ______
c 2+γ 2+c 2−γ 2−
γCaCO30 Ca Ca CO3 CO3[6.9–4]Equations for general aqueous complexation (i.e., speciation) reactions are then obtained similarly (e.g., Lichtner, 1985; Yeh & Tripathi, 1989):
Naxi = (Kix/γix) Π(γkack)aik
k=1x

i = 1,2,3,...,Mx[6.9–5]where xi is the concentration of the ith complexed species, Kix is the thermodynamic
equilibrium constant of the ith complexed species, γix is the activity coefficient of
the ith complexed species, Na is the number of aqueous components (or primary
species for systems in partial equilibrium), ck is the concentration of the kth aqueous component, γka is the activity coefficient of the kth aqueous component species,
axik is the stoichiometric coefficient of the kth aqueous component in the ith complexed species, Mx is the number of complexed species, and the subscripts and superscripts “x” and “a” refer to complexed species and aqueous components, respectively. The chemical equilibrium assumption often applies to most aqueous
complexation reactions since they do not involve any interphase mass transfer andthus are generally considered “fast”. Exceptions are oxidation–reduction reactions
that involve the transfer of one or more electrons for which disequilibrium is therule rather than the exception. These reactions must be treated kinetically (Lichtner, 1996). Aqueous speciation reactions can exert an important indirect effect upon
species mobility and fate, since sorption and degradation processes can depend
strongly upon the particular complexed form of a chemical.
6.9.2.2 Cation Exchange Reactions
Solute transport through soil is significantly affected by cation exchange reactions. These reactions govern the distribution of dissolved nutrients and otherchemical species between the soil solution and solid exchange phase, and thus influence the rate of movement of chemicals in soils. These reactions are particularly
important in practical problems of sodic soils reclamation (see example in Section
6.9.6).
Partitioning between the solid exchange phase and the solution phase can be
described by the general exchange equation (White & Zelazny, 1986)
zjAzi+Xzi + ziBzj+ º ziBzj+Xzj + ziAzj+
or by the widely used Gapon equation[6.9–6]

1516

CHAPTER 6z+izj+ º B zj+X + (1/z )Azi+
A1/z
1/zj
i
iX + (1/zj)B[6.9–7]

where A and B are chemical formulas for particular cations, X refers to an “exchanger” site on the soil, and zi and zj are the valences of species A and B, respectively. The mass action equation resulting from the exchange reaction (Eq. [6.9–6])
is
c–jzj+ zi _____
γiacizi+ zj
Kij = _____az+j— γj cj  — c–izi+ [6.9–8]
while the mass action expression resulting from the Gapon exchange equation (Eq.
[6.9–7]) is
c–jzj+ ________(γiacizi+)1/zi
Kij = ____
‰ c–izi+  (γjacjzj+)1/zj[6.9–9]Kij is the selectivity coefficient, and –ci is the exchanger-phase concentration of theith component (expressed in moles and equivalents per mass of solid in Eq. [6.9–8]
and [6.9–9], respectively). In the above equations all concentration quantities are
molar except for the adsorbed concentrations in Gapon equations, which are expressed as equivalents. As an example, for the case of exchange of the cations Ca2+and Na+,2+ º c2+ + 2c+–Ca
2 c–+Na + cCaNa[6.9–10]

we obtain the following Gapon equation
KCa/Na =2+
c–Ca
γNac+Na___________
2+)1/2‰ c–+Na  (γCacCa[6.9–11]

White and Zelazny (1986) provide a review of other general forms for cation selectivity coefficients. It is often assumed that the cation exchange capacity c–T is constant and independent of pH. When exchanger-phase concentration is expressed inequivalents per mass of solid, then
Mz

Σ c–i
c–T = i=1

[6.9–12]

where Mz is the number of ion exchange species.
Cation exchange reactions, recognized early by soil scientists and hydrologists to significantly regulate soil solution composition, have been implemented inmany models with both generalized and specific chemistry (e.g., Dutt et al., 1972;
Jury et al., 1978; Robins et al., 1980; Valocchi et al., 1981; Wagenet & Hutson, 1987;Šimçnek et al., 1996). Most of these models assume that equilibrium conditions
apply to the cation exchange reactions. Although the actual ion exchange reactions
written above are generally fast, in complex porous media (e.g., unsaturated soils)
the reaction rates may be limited by mass transfer from the flowing fluid to the exchange sites. A few transport models presented in the recent hydrology literature

MISCIBLE SOLUTE TRANSPORT

1517

include such kinetic effects for chemical species that undergo ion exchange reactions (e.g., Andre et al., 1998).
6.9.2.3 Adsorption Reactions
In addition to the ion exchange processes discussed above, dissolved chemical species may adsorb onto solid surfaces as a result of a variety of physical andchemical processes (Sposito, 1984; Brusseau & Rao, 1989; van der Zee & Destouni,
1992; van der Zee & Riemsdijk, 1994). Although the presentation here will focus
on equilibrium models, in many practical situations the time scales required for transfer of the dissolved chemical to the active sorption sites are large enough to require
use of nonequilibrium models. This is especially true for unsaturated soils when
there are dynamic changes in water and solute flux. These models are often based
upon diffusion and mass transfer concepts. Further elaboration of nonequilibrium
sorption modeling is widely available in the literature (e.g., see Selim et al., 1976;
Cameron & Klute, 1977; Sparks, 1989; Sardin et al., 1991; Harmon et al., 1989;
van Genuchten & Wagenet, 1989; and Section 6.3).
There are two approaches to describe equilibrium sorption reactions: isotherm
equations (leading to a KD approach for a linear case), and surface complexation
models. These are summarized briefly in the next two sections.
6.9.2.3.a Sorption Isotherms
Isotherms relate the adsorbed concentration of a particular species to its solution concentration. In the past, a linear equilibrium isotherm, s = KDc, where s isthe concentration in the adsorbed phase and KD is the distribution or partition coefficient, has been often used to describe solute interactions between the liquid andsolid phases of the soil. This leads to a constant retardation factor Rf in the governing
solute transport equation. While the use of a linear isotherm can greatly simplify
the mathematics of a transport problem, sorption and exchange reactions are generally nonlinear and usually depend also on the presence of competing species inthe soil solution. The solute retardation factor for nonlinear adsorption is not constant, as is the case for linear adsorption, but changes as a function of the slope of
the adsorption isotherm:
Rf = 1 + (ρb/θ)[ds(c)/dc][6.9–13]

where ρb is the bulk density and θ is the volumetric water content. It should also
be noted that the isotherm approach can be applied to equilibrium ion exchange reactions, in which case the multicomponent isotherm may be a nonlinear functionof all aqueous species that participate in the reactions (e.g., Valocchi et al., 1981;
Rubin & James, 1973).
A variety of models have been used in the past to describe the nonlinear nature of adsorption–desorption process. Van Genuchten and Šimçnek (1996) listed
some of the most commonly used sorption models relating the adsorbed concentration, s, to the total concentration, c, of an ion in solution. Although several of theseequations (e.g., the Langmuir and Freundlich equations) can be derived rigorously,
such as for the adsorption of gases onto solids, the expressions are generally used

1518

CHAPTER 6

only in a purely empirical fashion. Of the equations listed by van Genuchten andŠimçnek (1996), the most popular sorption models are the Langmuir, Freundlich,
and Temkin equations.
A general classification of adsorption as reflected by different features of theadsorption isotherm, such as the initial slope, the presence or absence of a plateau,or the presence of a maximum, was proposed by Giles et al. (1960). They divided
possible adsorption processes into four main classes—S, L (Langmuir), H (highaffinity), and C (constant partitioning) isotherms—and discussed mechanisms explaining the different types of isotherms. Increasing solution concentrations resultin increasing or decreasing adsorption rates for the S and L isotherms, respectively.
An H isotherm is characterized by extremely high affinities of the exchanger for
exchangeable ions, whereas a C isotherm results from constant partitioning of thesolute between the solution and adsorbed phases. The shape of an isotherm can havesignificant effects on transport predictions. For example, S and L isotherms lead toso-called unfavorable and favorable exchange situations, respectively (e.g., Valocchi et al., 1981). The latter (favorable) condition (e.g., for a Freundlich isothermwith exponent n1) produces a sharp constant-shape concentration front duringtransport in a soil profile, while the former (unfavorable) conditions produces a front
that continues to spread during transport (Reiniger & Bolt, 1972; van Genuchten& Cleary, 1979; Valocchi et al., 1981; van der Zee & van Riemsdijk, 1994).
Although many adsorption processes are more accurately described with themore sophisticated surface complexation models described below, isotherm models have been successfully applied to the environmentally significant class of neutral, relatively nonpolar organic compounds, such as chlorinated hydrocarbons andpesticides. In fact, for low solution concentration, the simple linear isotherm often
applies. In soils with significant fractions of organic carbon these compounds adsorb primarily to solid-phase organic matter as a result of hydrophobic interactions,and the KD of these compounds is often found to correlate directly with the fractional organic C content of the soil. Further details about hydrophobic sorption bysoils can be found in a number of references, including Karickhoff et al. (1979),
Hassett and Banwart (1989), Pignatello (1989), and Weber et al. (1991).
6.9.2.3.b Surface Complexation Models
Models based on the sorption isotherms are not sufficiently general to accountfor variations in sorption with pH, multiple oxidation states, electrostatic forces, andother factors. Surface complexation models, such as the constant capacitance, diffuse double layer, and triple layer models (Mattigod & Zachara, 1996), must be used
for these more complex conditions. These models were developed for adsorption
onto oxide surfaces, which are modeled as functional groups that undergo acid–base,metal binding, ligand exchange, and other reactions (Dzombak & Morel, 1990;
Stumm and Morgan, 1996). This enables a quantitative description of the impactof pH upon adsorption and the development of surface charge. The various surface
complexation models differ in their depiction of the interfacial region surrounding
an adsorbent, that is, the number of considered planes and the charge-potential relationships. The constant capacitance model includes only one plane, the diffuse
double layer model considers two planes, and the Stern and triple layer models con-

MISCIBLE SOLUTE TRANSPORT

1519

sider three planes. Detailed descriptions of different surface complexation models
are beyond the scope of this discussion, but are provided by several authors, including Papelis et al. (1988), Dzombak and Morel (1990), and Stumm and Morgan(1996).
6.9.2.4 Precipitation–Dissolution
Equations describing precipitation–dissolution reactions are also obtainedusing the law of mass action, but contrary to the other processes are represented byinequalities rather than equalities, as follows:Nap

Kip ≥ Qip = Π(γkack)aik
k=1

i = 1,2,...,Mp[6.9–14]

where Mp is the number of precipitated species, Kip is the thermodynamic equilibrium constant of the ith precipitated species (i.e., the solubility product equilibrium
constant), Qip is the ion activity product of the ith precipitated species, and apik is
the stoichiometric coefficient of the kth aqueous component in the ith precipitatedspecies. The inequality in Eq. [6.9–14] means that a particular precipitate is formed
only when the solution is supersaturated with respect to its aqueous components;
if the solution is undersaturated then the precipitated species (if it exists) will dissolve in order to reach equilibrium conditions. Equation [6.9–14] assumes that theactivity of the precipitated species is equal to unity.
Precipitation–dissolution reactions are often orders of magnitude slower thanother chemical reactions, while rates of dissolution of different minerals can also
differ by orders of magnitude. Therefore, precipitation–dissolution reactions usually have to be considered as kinetic, rather than equilibrium, reactions. Dissolution of dolomite or silicate minerals, such as feldspar, kaolinite, gibbsite, anorthite,and others, are typical examples of kinetic reactions. The precipitation–dissolution
reaction equations are often written with the rate depending on the difference between Qip and Kip, that is, on the disequilibrium driving force (see, e.g., Inskeep &Bloom, 1985, for calcite precipitation). Different mechanistic (Plummer et al.,1978) and empirical models can be used as well. For a more elaborate discussionof kinetics and related issues (the surface area, a moving boundary problem, aboundary layer, quasistationary states, and so on) see the excellent paper by Lichtner (1996).
6.9.2.5 Reactions with Organic Matter and Effects of Bacteria
There is increasing evidence that microorganisms, present throughout the subsurface, can catalyze many chemical reactions and thus affect the fate of organicand inorganic contaminants (Zysset et al., 1994b; Rittmann & VanBriesen, 1996).
Transport and biodegradation of anthropogenically introduced chemical contaminants such as hydrocarbons are mediated by subsurface aerobic or anaerobic bacteria. Bacteria catalyze redox reactions in which organic compounds act as the electron donor and inorganic substances (oxygen, nitrate, sulfate, or metal oxides)
serve as the electron acceptor. Bacteria, by catalyzing such reactions, gain energy
and organic C to build up new biomass (Zysset et al., 1994b). The transport of bacteria is more complex than that of solutes, since it is affected by various additional

1520

CHAPTER 6

mechanisms, such as filtration, straining, sedimentation, adsorption, and desorption, growth, and decay and death of bacteria (see Section 6.7).
The effect of bacterial growth and microbiological reactions on geochemical transport in the subsurface has been recently incorporated into the biogeochemical transport models HYDROBIOGEOCHEM (Yeh et al., 1998), FEREACT(Tebes-Stevens et al., 1998), BIORXNTRN (Hunter et al., 1998), and FLOTRAN(Lichtner, 2000). As an interesting illustration of the coupling between biochemical and geochemical processes, an example from Tebes-Stevens and Valocchi(1998) is presented. The example considers the fate of nitrilotriacetate (C6H6NO63−
or simply NTA), an organic chelating agent that has been detected along with metals and radionuclides in some contaminated groundwater. Although NTA isbiodegradable, there is some uncertainty regarding the specifics of the degradation
process (VanBriesen, 1998). For illustration purposes, it is assumed that only thesingle protonated form is biodegradable, which leads to the following overall stoichiometry under aerobic conditions (Tebes-Stevens & Valocchi, 1998):
HNTA2− + 1.62O2 +1.272H2O + 2.424H+→ 0.576C5H7O2N + 3.12H2CO3* + 0.424NH4+[6.9–15]

where C5H7O2N is the simple elemental formula for biomass (Rittmann & VanBriesen, 1996) and H2CO3* represents both aqueous CO2 and H2CO3. In the above
equation, electrons are transferred from the donor substrate, HNTA2−, to the acceptor,
O2, in order to produce the energy required for synthesis of new biomass. The overall rate of the above biological reaction is often described using a Monod-type kinetic model (Yeh et al., 1998; Rittmann & VanBriesen, 1996; Lichtner, 1998):
r = rr[S/(KS + S)][A/(KA + A)

[6.9–16]

where r and rr are the actual and maximum reaction rates, S and A are the concentrations of the electron donor substrate and the electron acceptor, respectively, andKS and KA are the half saturation constants for the substrate and the electron acceptor, respectively, that is, the concentrations at which the rate is reduced by 50%.
In systems with dynamically growing biomass, the maximum reaction rate is normally assumed directly proportional to the amount of biomass. The overall stoichiometry above clearly shows that biodegradation of NTA is coupled with otherpotentially important chemical changes; specifically, the increase in the pH and production of carbonic acid may exert significant influence over other reactions (e.g.,precipitation–dissolution). In the example considered by Tebes-Stevens and Valocchi (1998), which also included metal species that complex NTA, the change inpH also affected the biodegradation rate by changing the fraction of NTA that wasin the form of HNTA2−.
The production of carbon dioxide and consumption of oxygen represent important effects of microbiological activity in the soil environment in the unsaturated
zone. Microbiological respiration, combined with the plant root respiration and highwater contents, can change CO2 concentrations several orders of magnitude fromvalues that are at equilibrium with the CO2 content in the atmosphere (0.035%), up

MISCIBLE SOLUTE TRANSPORT

1521

to extreme values of about 20%. The solubility of many solid phases, such as carbonates and oxihydroxides, significantly changes within this range of CO2, primarily
because changes in soil CO2 produce changes in soil pH for all but acid soils(Šimçnek & Suarez, 1994). The UNSATCHEM models (Šimçnek & Suarez, 1993a,b; Šimçnek et al., 1996) consider the effects of CO2 producing microbiological activity and CO2 transport in the soil environment on geochemical transport. Suarez
and Šimçnek (1993) gave several examples of the seasonal simulations of the CO2
concentrations in the soil and the CO2 fluxes to the atmosphere for different crops
and compared their results with published field data.
6.9.2.6 Activity Coefficients and Thermodynamic Equilibrium Constants
6.9.2.6.a Thermodynamic Equilibrium Constants
Large databases of thermodynamic equilibrium constants for hundreds of reactions, that is, equilibrium constants for aqueous complexed species, solubility
products for solid phases, and the redox potentials for various oxidation–reduction
reactions, are usually available as part of general geochemical equilibrium models(Mattigod & Zachara, 1996). These constants are normally valid at reference temperature (25°C) and for zero ionic strength. Since soil solutions involve different
temperatures, and have finite ionic strength, the thermodynamic constants must be
corrected to reflect these factors.
6.9.2.6.b Activity Coefficients
The solute activity coefficients are formally defined as (Stokes, 1979)
ai = γi(mi/m0)[6.9–17]

where ai is the activity, mi is the molality (mol kg−1), m0 is the standard state unit
molality (i.e., 1 mol kg−1), and γi is the activity coefficient of the ith ion. Single-ion
activity coefficients may be calculated either using the Davies equation and an extended version of the Debye–Hückel equation (Truesdell & Jones, 1974) for diluteto moderately saline solutions, or by use of the Pitzer expressions (Pitzer, 1979),
which are considered more accurate for solutions having very high ionic strength
down to solutions having infinite dilution. The Davies equation has the following
form:%I& − 0.3I
lnγ = − Az2 ______— 1 + %I&


[6.9–18]

The extended version of the Debye–Hückel equation (Truesdell & Jones, 1974) is
given by
Az %I& + bI
lnγ = − ________
1 + Ba%I&2[6.9–19]

where A and B are constants depending only upon the dielectric constant, density,and temperature; z is the ionic charge (valence); a and b are two adjustable pa-

1522

CHAPTER 6

rameters accounting for the hydrated ionic radius of a species and for ion pair interactions, respectively (Lichtner, 1996), and I is the ionic strengthM

I = 0.5 Σ
z2mi
i=1 i

[6.9–20]

where M is the number of species in the solution mixture. The activities of neutral
species are calculated as (Sposito, 1981)
lnγ = a′I

[6.9–21]

where a′ is an empirical parameter.
At high ionic strength (solutions at high concentrations, such as in brines orin soils in arid regions), activity coefficients are no longer universal functions ofionic strength, but depend on the relative concentration of the various ions present
in solution (Felmy & Weare, 1986). The activity coefficients can then be expressedin terms of a virial-type expansion of the form (Pitzer, 1979)
lnγi = lnγiDH + Σj Bij(I)mj + Σj Σk Cijkmjmk + ...

[6.9–22]

where γiDH is a modified Debye–Hückel activity coefficient, which is a universal
function of ionic strength, and Bij and Cijk are specific coefficients for each ion interaction. The disadvantage of using Pitzer’s equations is their extensive demandon computer power.
6.9.2.6.c Temperature Dependence
Most thermodynamic equilibrium constants depend on the temperature andpressure of the system. The temperature dependence of the thermodynamic equilibrium constants is often expressed as a power function of the absolute temperature (Truesdell & Jones, 1974; Plummer & Busenberg, 1982)
logK = a1 + a2/T + a3T + a4logT + a5/T2[6.9–23]

where T is absolute temperature, and a1 through a5 are empirical constants. The pressure dependence is usually neglected for near soil-surface conditions. The temperature dependence of the equilibrium constants for which the constants of Eq.
[6.9–23] are not known is estimated with the enthalphy of reaction and the Van’t
Hoff expression (Truesdell & Jones, 1974)
lnK = lnKr − (∆Hr/R)[(1/T) − (1/Tr)][6.9–24]

where Tr is the reference temperature (298.15 K), R is the universal gas constant,
Kr is the equilibrium constant at the reference temperature, and ∆Hr is the enthalpy
change of the reaction.
6.9.3 Mass-Balance Transport Equations
Mass balance equations for Na aqueous components are written in terms ofthe total analytical dissolved concentration Cj, the total precipitated concentration

MISCIBLE SOLUTE TRANSPORT

1523

Pj, the total sorbed (adsorbed and ion-exchanged) concentration Sj, and the total analytical concentration Tj, of the jth aqueous component (Yeh & Tripathi, 1989):
Mx

Cj = cj + i=1
Σ aijx xi

j = 1,2,...,Na[6.9–25]

where the variables are as defined above when Eq. [6.9–5] was introduced, that is,
xi is the concentration of the ith complexed species, axij is the stoichiometric coefficient of the jth aqueous component in the ith complexed species (this is contrary
to ai which is the activity of the ith species), and cj is the concentration of the jthaqueous component.
My

Mz

Sj = i=1
Σ aijs si + i=1
Σ aijz c–i

j = 1,2,...,Na[6.9–26]

Mp

Pj = i=1
Σ apijpi
Tj = Cj + Sj + Pj

j = 1,2,...,Naj = 1,2,...,Na[6.9–27][6.9–28]

where pi is the concentration of the ith precipitated species, si is the concentrationof the ith sorbed species, and aijz and asij are the stochiometric coefficient of the jthaqueous component in the ith ion-exchange and sorbed species, respectively. Allconcentrations in the above equations are expresses in units of moles per unit volume of water in the system. When neglecting source terms, decay, and biodegradation reactions, we obtain the governing mass balance equations:
∂θT
∂
∂Cj
∂viCj
j + Sj + Pj)___j = ∂θ(C____________= __ θDik ___ − ____
∂t
∂t
∂xi ‰
∂xk 
∂xi

j = 1,2,...,Na[6.9–29]

where θ is the water content, vi are components of Darcy’s velocity, Dik are coefficients of the dispersion tensor D, t is time, and xi are spatial coordinates. Since thewater content can change rapidly in the soil environment, the total precipitated concentration, Pj, and the total sorbed concentration, Sj, are often expressed not on a
volume, but on a mass basis (P*,j S*)
j (Šimçnek & Suarez, 1994):∂θC∂(Sj* + P*)
∂
∂Cj
∂viCj
j
___ j + ρ _________= __ θDik ___ − ____
∂t
∂t
∂xi ‰
∂xk 
∂xi

j = 1,2,...,Na[6.9–30]

The second term on the left side of Eq. [6.9–30] is zero for components that
do not undergo sorption, ion exchange, or precipitation–dissolution. Transport Eq.
[6.9–29] or [6.9–30] have to be coupled with a numerical solution of the Richards
equation (or even better with the two-phase flow of gas and water equations) for
spatially and temporally variable water contents and fluxes. If the reactions are assumed to be equilibrium controlled, then the transport Eq. [6.9–29] or [6.9–30] arealso coupled with a geochemical speciation model that essentially solves Eq.
[6.9–25] through [6.9–28] (along with the mass action equations) for the individual aqueous, sorbed, and precipitated species (i.e., ci, si, pi); this step is needed tocompute the mass balance for the sorbed, ion exchange, and precipitated species(Yeh & Tripathi, 1991; Engesgaard & Kipp, 1992).

1524

CHAPTER 6

The governing equations given above assume that all reactions are equilibrium-controlled. The more general case with mixed kinetic and equilibrium reactions is presented by many investigators, including Chilakapati (1995), Lichtner(1996, 1998), Tebes-Stevens et al. (1998), Yeh et al. (1998), and Chilakapati et al.
(1998). This more general formulation includes kinetic source–sink terms in thetransport Eq. [6.9–30], plus separate mass balance equations for immobile species(e.g., sorbed, precipitated, and biomass species) governed by kinetically controlled
reactions.
6.9.4 Numerical Implementation
The governing equations consist of partial differential equations for transport(i.e., Eq. [6.9–29] or [6.9–30]) plus a set of nonlinear algebraic and ordinary differential equations describing the equilibrium and kinetic reactions. Since the reaction of one species depends upon the concentration of many other species, theseequation sets are coupled. For complex geochemical systems, consisting of many
components and multidimensional transport, numerical solution of these coupledequations is challenging.
The chemical equilibrium reactions can be formulated and solved mathematically by means of two different but thermodynamically equivalent approaches,that is, (i) by minimizing the Gibbs free energy of the system subject to the constraints of mass balance, or (ii) by using equilibrium constants. For the approachbased on the Gibbs free energy, the same mass is maintained throughout the iteration process, while equilibrium is iteratively reached by transferring mass from one
species to another (Felmy, 1990). The solution procedure for the second approach
begins with an initial guess for a set of components from which the minimum Gibbs
free energy composition can be readily calculated using equilibrium constants, after
which the mass balance equations are solved by iteration (Westall et al., 1976). Although the first approach tends to be more robust (Steefel & MacQuarrie, 1996),
the second approach is a more popular method when dealing with chemical equilibrium systems (Truesdell & Jones, 1974; Westall et al., 1976). A Newton–Raphson method is usually employed to solve the final set of nonlinear algebraic equations. Steefel and MacQuarrie (1996) list several packages that are available forsolving sets of ordinary differential and differential algebraic equations that can beused to solve systems of both kinetic and equilibrium reactions.
In a review, Yeh and Tripathi (1989) identified three different approaches formathematically solving multicomponent transport problems: (i) a mixed differential and algebraic approach, (ii) a direct substitution approach, and (iii) a sequential iteration approach. In the first approach, the sets of differential and algebraic
equations describing the transport processes and chemical reactions, respectively,
are treated simultaneously (Miller & Benson, 1983; Lichtner, 1985). In the second
approach, the algebraic reactions representing the nonlinear chemical reactions aresubstituted directly into the differential mass balance transport equations (Rubin &James, 1973; Valocchi et al., 1981; Jennings et al., 1982). The third approach (theiterative two-step method) divides the numerical solution into two independent steps
that are solved separately. The first step provides solutions of the partial-differential equations governing convective–dispersive transport, whereas the second step

MISCIBLE SOLUTE TRANSPORT

1525

involves algebraic equations describing the different chemical equilibrium reactions,as well as ordinary differential equations describing kinetic reactions. The two steps
are solved sequentially and iteratively (Walsh et al., 1984; Cederberg et al., 1985;
Bryant et al., 1986; Yeh & Tripathi, 1991; Šimçnek & Suarez, 1994; Tebes-Stevens
et al., 1998). Based on a study of computer resource requirements, Yeh and Tripathi(1989) suggested that only the third method (sequential iteration) is feasible for realistic multidimensional applications. In order to speed up the calculations, the iteration part between the chemical and transport steps is often avoided, which leads
to a sequential noniterative approach. This is also called the operator-splitting
method, and it has many advantages in terms of computational speed and programming flexibility (Steefel & Yabusaki, 1996; Steefel & MacQuarrie, 1996).
There are, however, numerical errors associated with this noniterative approach, particularly for fast reactions. These errors can be reduced somewhat by selecting
higher-order time splitting approaches (Valocchi and Malmstead, 1992; Zysset etal., 1994a; Steefel & MacQuarrie, 1996).
6.9.5 Effects of Solution Composition on Hydraulic Properties
and Reclamation Models
High concentration of soluble salts accumulated in the soil can significantly
decrease agricultural productivity by imposing osmotic stress on growing crops, resulting in decreased crop yields and possibly even crop failure (Šimçnek & Suarez,
1997). Moreover, accumulation of monovalent cations, particularly Na, may leadto clay dispersion, swelling, and overall poor soil physical properties. These
processes have an adverse effect on the hydraulic properties of a soil, including hydraulic conductivity, infiltration rates, and soil water retention (Dane & Klute, 1977;
Shainberg & Levy, 1992). These negative effects are usually explained on the basisof the diffuse double layer theory. Monovalent cations satisfy the negative chargeof clays at a greater distance away from the surface than divalent cations, thus leading to a more extensive diffuse double layer for monovalent cations. The same istrue for a dilute as compared with a high-electrolyte soil solution. A consequenceof the greater extension of the diffuse double layer associated with monovalent ionsis the greater repulsion force or swelling pressure between neighboring clay
platelets, a situation that results in swelling and also dispersion.
The hydraulic conductivity can be reduced several orders of magnitude bysoil swelling, while the effect upon water content is much smaller (Russo, 1988).
Models that describe the effects of sodicity and salinity on the hydraulic conductivity can be divided into two broad groups, (i) theoretical models based on double
layer theory (Lagerwerff et al., 1969; Russo & Bresler, 1977; Russo, 1988) and (ii)semiempirical models (Yaron & Thomas, 1968; McNeal, 1968, 1974; Dane, 1977).
McNeal’s work is based on a simple clay-swelling model, where clay swelling isrelated to decreases in soil hydraulic conductivity.
Initial attempts to describe and quantify the reclamation processes (i.e., removal of excessive salts from the soil profile and exchanging Na for Ca in the sorbed
complex) were based mostly on chromatographic models (Dutt et al., 1972; Tanjiet al., 1972). Models of this type assume that the soil profile is stratified chemicallyand physically into a series of horizontal layers, with water moving from one layer

1526

CHAPTER 6

to another by piston flow. Chemical processes between the solution, exchanger, andsoil mineral phase within a particular layer are assumed to be instantaneous equilibrium processes. The models by Dutt et al. (1972) and Tanji et al. (1972) both consider major ion species, several complexation reactions, cation exchange, and dissolution of gypsum and calcite. Since these two models did not consider explicitly
the soil hydraulic conductivity and water flow velocities, they cannot predictwhether an amendment is concentrated enough to significantly improve the infiltration rate. Hence, they cannot accurately predict the time required for completing the reclamation process. More advanced models, such as LEACHM (Wagenet& Hutson, 1987) or UNSATCHEM (Šimçnek & Suarez, 1993a; Šimçnek et al.,1996), couple major ion chemistry with variably saturated water flow and solutetransport. These models are powerful tools for quantitative predictions of theamount of water and amendment required in order to reclaim soil profiles to desired levels of salinity and ESP (exchangeable Na percentage). Of the reclamation
models, only the UNSATCHEM model considers dynamic changes in CO2 concentrations that can significantly alter pH and the solubility of CaCO3, while also
accounting for kinetic dissolution and precipitation processes for calcite and dissolution of dolomite (Šimçnek et al., 1996; Šimçnek & Suarez, 1997). The UNSATCHEM model also takes into account effects of solution composition on hydraulic properties. Example results using this model for sodic soil reclamation will
be presented in the next section.

6.9.6 Applications
There are a large number of possible applications of reactive multispecies
solute transport models in the soil environmental sciences. For example, Yeh andTripathi (1991) used the HYDROGEOCHEM model to simulate the release of trace
metals and acidity from an acidic uranium mill tailings pile. Water leaching fromuranium mills tailing piles often contains, in addition to U, many trace metals including Mo, Se, As, and Cr (Brookins, 1984). Using the DYNAMIX model,
Narasimhan et al. (1986) also studied groundwater contamination from an inactive
uranium mill tailings pile. Similar studies were carried out by Walter et al. (1994)
and Gerke et al. (1998).
Abandoned mines and mine tailing piles represent significant environmental problems. Water, rich with O2, infiltrates through these structures, which often
contain pyrite. Pyrite oxidation leads to the precipitation of iron oxide and the formation of acid that, in turn, mobilizes heavy metals with consequent pollution ofenvironment (Lichtner, 1998). Multicomponent transport models are crucial tools
in analyzing the potential environmental hazard of these sites.
Alkaline floods, that is, the process of oil recovery by the subsurface injection of high pH solutions that involve both precipitation–dissolution and cation exchange, were simulated by Bryant et al. (1986) using a multicomponent model. The
injected alkaline solution neutralizes organic acids in the oil, thereby mobilizing
the oil by a variety of mechanisms.
Safe disposal and long-term isolation of radioactive waste, especially highlevel nuclear waste, represents another significant challenge. Since prediction has

MISCIBLE SOLUTE TRANSPORT

1527

to be done over geological times, numerical models of reactive transport are indispensable tools to handle such complex problems. For example, Viswanathan etal. (1998) used a multicomponent model to study the migration of Np between therepository and the groundwater table after a hypothetical repository breach at thepotential high-level nuclear waste repository at Yucca Mountain, Nevada.
Sodic soil reclamation (i.e., removal of excessive salts from the soil profileand replacement of exchangeable Na+ with Ca2+) represents a typical application
for geochemical transport models in agriculture and soil science. This complex problem involves variably saturated water flow, solute transport, dissolution of solid
phases (i.e., reclamation amendments, such as gypsum or calcite), cation exchange(replacement of Na+ with Ca2+), complexation reactions, effects of solution composition on hydraulic conductivities and infiltration rates, and effects of elevated
CO2 on pH and consequently calcite dissolution. We will briefly summarize here
results of the study of sodic soil reclamation by Šimçnek and Suarez (1997).
Šimçnek and Suarez (1997) considered six different reclamation strategies.
In the first two cases, high-quality (i.e., low ion concentrations) and gypsum-saturated waters were applied to the soil without additional amendments. In the third
example, high-quality water was applied to soil where gypsum was incorporated
into the top 20 cm of the soil profile. In the next three cases, high-quality water andwaters at two different acid concentrations were applied to the surface of a calcareous
soil. The solution composition of the water initially present in the soil profile was
assigned to be that of the following highly sodic soil water system: CaT = 0.2, MgT= 0.2, NaT = 4.8, ClT = 4.8, alkalinity = 0.4 mmolcL−1 (alkalinity is defined as 2CO3T+ HCO3T + OH–H; for high H concentrations, alkalinity can be negative). The cationexchange capacity was set at 200 mmolckg−1 and the ESP (exchangeable Na percentage) at 60 (Ca = 40.0, Mg = 40.0, Na = 120.0 mmolckg−1). Four irrigation waters were used for infiltration: high-quality dilute water (CaT = 1.5 mmolcL−1, MgT= 0.5, NaT = 2.0, ClT = 1.0, SO4T = 2.5, Alkalinity = 0.5), water high in Ca2+ andSO42− (almost gypsum saturated, CaT = 32.0, MgT = 0.5, NaT = 2.0, ClT = 1.0, SO4T
= 33.0, Alkalinity = 0.5), and two acidified waters, one at pH = 2.05 (CaT = 1.5,
MgT = 0.5, NaT = 2.0, ClT = 3.0, SO4T = 11.0, Alkalinity = −10.0) and the other at
pH = 1.09 (CaT = 1.5, MgT = 0.5, NaT = 2.0, ClT=93.0, SO4T = 11.0, Alkalinity =−100.0). Model-predicted CO2 concentrations were used for the last three cases,
where CO2 has an important effect on calcite dissolution. The hydraulic conductivity reduction due to solution composition (see Section 6.9.5) was considered for
all calculations. The examples assume that the soil surface was flooded with water,
with the water kept at a constant level of 1 cm above the soil surface. Šimçnek andSuarez (1997) give additional information about all relevant input parameters andthe different reclamation strategies. Figure 6.9–1 shows Na adsorption ratio, SAR[= Na/{(Ca+Mg)/2}0.5], profiles for all six reclamation strategies. Reclamation
strategies were evaluated based on the considerations of the time required forreclamation to be completed and the quantity of water needed. Times to profile saturation for six strategies for selected conditions were calculated to be 358, 15.3, 2.0,
10, 6.4, and 3.2 d, respectively; times to profile reclamation were 86,870, 195.6,
116.5, 91, 64, and 15.5 d, respectively; and amounts of infiltrated water needed forreclamation were 7230, 559, 470, 600, 485, and 197 cm, respectively (Šimçnek &Suarez, 1997).

1528
CHAPTER 6

Fig. 6.9–1. Sodium adsorption profiles for six reclamation strategies: (A) irrigation with high quality water and no amendments, (B) irrigation with gypsum-saturated water,(C) irrigation with high quality water and gypsum incorporated in the top 20 cm, (D) irrigation with high quality water and calcite throughout the soil profile, (E) irrigation with acid water at pH = 2.05 and calcite throughout the soil profile, and (F) irrigation with acid water at pH = 1.09 and calcite throughout the soil profile.

MISCIBLE SOLUTE TRANSPORT

1529

Reclamation simulations using the UNSATCHEM model demonstrated thecritical importance of including the effects of solution composition on soil hydraulicparameters in the model predictions (Šimçnek & Suarez, 1997). Only models withthis feature can successfully estimate the time required for the reclamation process.
This information is critical to decisions about choice of strategy to be employed,
and its economic viability. Šimçnek and Suarez (1997) also showed using the multicomponent transport model that parameters, such as quantities of water, as wellas the amounts of reclamation materials (and depth of incorporation), are dependent upon the CO2 concentration, thus requiring a model capable of predicting CO2
concentrations under dynamic field conditions. Flooding, combined with CO2 production, elevates CO2 concentrations in the root zone and results in a very effective reclamation practice. Application of acid or dilute water to calcareous soil withelevated CO2 concentrations, incorporating gypsum into the top of the soil profile,and application of gypsum-saturated water were shown to be decreasingly less effective reclamation strategies. All these conclusions were consistent with previously
reported field experiments.
Since sodic soil reclamation is a very complex process, only a model thatquantitatively addresses and integrates all the involved processes can help in exploring issues that were impossible to be studied earlier and in drawing conclusions
that could not be obtained using simpler models.
6.9.7 Use of Geochemical Transport Models
In recent decades, much of the research in the fields of soil science and groundwater hydrology has focused on understanding the fundamentals of pollutant fate
and transport processes. The ultimate practical goal is to utilize such fundamental
knowledge to develop cost-effective, technically sound strategies for resource management and pollution remediation and prevention. As society continues its technological development, and as our scientific understanding matures, the types ofpollution problems that pose the more significant environmental threats become increasingly complex. Problems such as radioactive waste disposal, acid mine pollution, heavy metal transport, bioremediation of organic contaminants, and reclamation of saline soils all require understanding the coupling of multiplehydrogeological, geochemical, and microbiological processes. These are the problems that have motivated development of the class of sophisticated multicomponent reactive transport models presented here. It is the continually increasing speed
and power of modern computers that has enabled such models to become a convenient tool for analysis of complex geochemical systems.
It is an open question whether we can expect to achieve the practical goal ofusing multicomponent reactive transport models as true predictive tools that can beused with confidence to make the quantitative forecasts necessary for remediation
design. At present, there is even controversy regarding the predictive capability ofsimpler models of saturated water flow and single component transport (see, e.g.,
Konikow & Bredehoeft, 1992; de Marsily et al., 1992; Oreskes et al., 1994). The
precision of the obtained predictions depends to a large extent upon the accuracyof available model input parameters and upon a proper conceptualization of soil heterogeneity and other system complexities. As discussed in Chapter 3, it is techni-

1530

CHAPTER 6

cally difficult to accurately determine the soil hydraulic function parameters characterizing soil water retention and permeability properties, which are the most important input variables for models based on the numerical solution of the variably
saturated flow (Richards) equation. Adding multicomponent geochemical and microbiological processes only confounds the problem further. There may exist uncertainty about some of the equilibrium constants contained in thermodynamic data
sets, and also about the composition of mixed solid phases that are present in natural soils. Many kinetic rate laws have parameters related to pore-scale diffusion
length scales or reactive mineral surface area; such parameters may never be knownwith precision for realistic spatially heterogeneous soils. The impact of reaction parameter uncertainty upon predictions of multicomponent transport models has
been studied recently by Tebes-Stevens and Valocchi (1998). Bethke (1996) presents a general discussion of sources of uncertainty in geochemical modeling.
In spite of the reservations noted above, multicomponent transport models areuseful and valuable tools for analyzing complex environmental pollution problems
and developing practical management strategies. The proper role of reactive transport modeling was assessed by Steefel and Van Cappellen (1998) who stated,“While reactive transport modeling is a fundamentally quantitative approach, it canbe used to test the full range of implications of qualitative or conceptual models. It
can help guide field observations by identifying which parameters and processes
actually control the system behavior.” Elaborating upon the assessment by Steefel
and Van Cappellen (1998), we can identify several specific key ways that such models can be used:
Process Coupling and Interactions. While laboratory experimentation is thebest tool for conducting fundamental studies of geochemical and microbiologicalprocesses under controlled conditions, mathematical models are the best tools forinvestigating the impacts of multiple coupled biogeochemical reactions in the presence of complex flow fields and spatial heterogeneity. These models also enable
extrapolation to environmentally relevant temporal and spatial scales.
Interpretation of Laboratory and Field Data. Multicomponent transportmodels provide a useful framework for interpreting experimental results. The
model serves as a tool for understanding qualitative and quantitative trends and relationships present in the data. In many cases properly applied modeling will identify specific hypothesis-driven experiments that will lead to more effective quantitative understanding of underlying biogeochemical processes.
Sensitivity Analysis. One of the most powerful uses of multicomponent transport models is conducting sensitivity analyses. These analyses permit the systematic evaluation of the impact of model parameters (both reactive and hydrogeological), initial conditions, and boundary conditions on the model output. The resultsof a sensitivity analysis provide insight into the relative importance of individual
reactions within a complex biogeochemical system; thus, these results can help identify the most important parameters and processes and thereby help guide allocationof resources for laboratory and field investigations.
Integration and Synthesis. The transport model is a tool for integrating allof the knowledge obtained from simulation, sensitivity analyses, and laboratory andMISCIBLE SOLUTE TRANSPORT

1531

field experimentation. This integration will lead to more coherent and rigorous conceptual models for the underlying coupled transport and reactions processes, andto models that are validated through their agreement with this diverse body of knowledge. As multicomponent transport models become increasingly capable of describing a broad range of natural phenomena, they become increasingly powerfultools in the fields of environmental and soil science.

6.9.8 References
Abrams, R.H., K. Loague, and D.B. Kent. 1998. Development and testing of a compartmentalized reaction network model for redox zones in contaminated aquifers. Water Resour. Res.
34:1531–1541.
Allison, J.D., D.S. Brown, and K.J. Novo-Gradac. 1991. MINTEQA2/PRODEFA2. A geochemical assessment model for environmental systems: Version 3.0 user’s manual. EPA/600/3-91/021.
USEPA, Washington, DC.
Andre, C., M. Sardin, P. Vitorge, and M.-H. Fauré. 1998. Analysis of breakthrough curves of Np(V) in
clayey sand packed column in terms of mass transfer kinetics. J. Contam. Hydrol. 35:161–173.
Appelo, C.A.J., and D. Postma. 1993. Geochemistry, groundwater and pollution. A.A. Balkema, Rotterdam, the Netherlands.
Bahr, J.M. 1990. Kinetically influenced terms for solute transport affected by heterogeneous and homogeneous classical reactions. Water Resour. Res. 26:21–34.
Bethke, C.M. 1994. The Geochemist’s Workbench. Version 2.0, A Users Guide to Rxn, Act2, Tact, React,
and Gtplot. Hydrogeology Program, University of Illinois.
Bethke, C.M. 1996. Geochemical reaction modeling—Concepts and applications. Oxford Univ. Press,
New York, NY.
Brookins, D.G. 1984. Geochemical aspects of radioactive waste disposal. Springler-Verlag, New York,NY.
Brusseau, M.L., and P.S.C. Rao. 1989. The influence of sorbate-organic matter interactions on sorption
nonequilibrium. Chemosphere 18:1691–1706.
Bryant, S.L., R.S. Schechter, and L.W. Lake. 1986. Interactions of precipitation/dissolution waves andion exchange in flow through permeable media. AIChE J. 32:751–764.
Cameron, D.A., and A. Klute. 1977. Convective–dispersive solute transport with a combined equilibrium and kinetic adsorption model. Water Resour. Res. 13:183–188.
Cederberg, G.A., R.L. Street, and J.O. Leckie. 1985. A groundwater mass transport and equilibrium chemistry model for multicomponent systems. Water Resour. Res. 21:1095–1104.
Chilakapati, A. 1995. RAFT: A simulator for reactive flow and transport of groundwater contaminants.
PNL-10636. Pacific Northwest Lab., Ritchland, WA.
Chilakapati, A., T. Ginn, and J. Szecsody. 1998. An analysis of complex reaction networks in groundwater modeling. Water Resour. Res. 34:1767–1780.
Dane, J.H. 1977. Calculation of hydraulic conductivity decreases in the presence of mixed NaCl-CaCl2
solutions. Can. J. Soil Sci. 58:145–152.
Dane, J.H., and A. Klute. 1977. Salt effects on the hydraulic properties of a swelling soil. Soil Sci. Soc.
Am. J. 41:1043–1049.
de Marsily, G., P. Combes, and P. Goblet. 1992. Ground-water models cannot be validated—Comment.
Adv. Water Resour. 15:367–369.
Dutt, G.R., R.W. Terkeltoub, and R.S. Rauschkolb. 1972. Prediction of gypsum and leaching requirements for sodium-affected soils. Soil Sci. 114:93–103.
Dzombak, D.A., and F.M. Morel. 1990. Surface complexation modeling. John Wiley, New York, NY.
Engesgaard, P., and K.L. Kipp. 1992. A geochemical transport model for redox-controlled movementof mineral fronts in groundwater flow systems—A case of nitrate removal by oxidation of pyrite.
Water Resour. Res. 28:2829–2843.
Essaid, H.I., B.A. Bekins, E.M. Godsy, E. Warren, M.J. Baedecker, and I.M. Cozzarelli. 1995. Simulation of aerobic and anaerobic biodegradation processes at a crude oil spill site. Water Resour.
Res. 31:3309–3327.
Felmy, A.R. 1990. GMIN: A computerized chemical equilibrium model using a constrained minimizationof the Gibbs free energy. PNL-7281. Pacific Northwest Laboratory, Ritchland, WA.

1532

CHAPTER 6

Felmy, A.R., D.C. Girvin, and E.A. Jenne. 1984. MINTEQ A computer program for calculating aqueous geochemical equilibria. NTIS PB84-15718. EPA-600/3-84-032. Natl. Tech. Inform. Serv.,
Springfield, VA.
Felmy, A.R., and J.H. Weare. 1986. The prediction of borate mineral equilibria in natural waters: Application to Searles Lake, California. Geochim. Cosmochim. Acta 50:2771–2783.
Fryar, A.E., and F.W. Schwartz. 1994. Modeling the removal of metals from groundwater by a reactive
barrier: Experimental results. Water Resour. Res. 30:3455–3469.
Gerke, H.H., J.W. Molson, and E.O. Frind. 1998. Modelling the effect of chemical heterogeneity on acidification and solute leaching in overburden mine spoils. J. Hydrol. 209:166–185.
Giles, C.H., T.H. MacEwan, S.N. Nakhwa, and D. Smith. 1960. Studies in adsorption. XI. A study ofclassification of solution adsorption isotherms, and its use in diagnosis of adsorption mechanisms
and in measurement of specific surface areas of solids. J. Chem. Soc. Oct. p. 3973–3980.
Gureghian, A.B. 1981. A two-dimensional finite-element solution for the simultaneous transport of waterand multi solutes through a nonhomogeneous aquifer under transient saturated–unsaturated flow
conditions. Sci. Total Environ. 21:329–337.
Harmon, T.C., W.P. Ball, and P.V. Roberts. 1989. Nonequilibrium transport of organic compounds ingroundwater. p. 405–437. In B.L. Sawhney and K. Brown (ed.) Reactions and movement of organic chemicals in soils. Spec. Publ. 22. SSSA, Madison, WI.
Hassett, J.J., and W.L. Banwart. 1989. The sorption of nonpolar organics by soils and sediments. p.
243–270. In B.L. Sawhney and K. Brown (ed.) Reactions and movement of organic chemicals
in soils. Spec. Publ. 22. SSSA, Madison, WI.
Hunter, K.S., Y. Wang, P. Van Cappellen. 1998. Kinetic modeling of microbially-driven redox chemistry of subsurface environments: Coupling transport, microbial metabolism and geochemistry.
J. Hydrol. 209:53–80.
Huyakorn, P.S., J.B. Kool, and Y.S. Wu. 1991. VAM2D—Variably saturated analysis model in two dimensions. Version 5.2 with hysteresis and chained decay transport. Documentation and user’s
guide. NUREG/CR-5352, Rev. 1. U.S. Nuclear Regulatory Commission, Washington, DC.
Inskeep, W.P., and P.R. Bloom. 1985. An evaluation of rate equations for calcite precipitation kinetics
at pCO2 less than 0.01 atm and pH greater than 8. Geochim. Cosmochim. Acta 49:2165–2180.
Jennings, A. A. 1987. Critical chemical reaction rates for multicomponent groundwater contamination
models. Water Resour. Res. 23:1775–1784.
Jennings, A.A., D.J. Kirkner, and T.L Theis. 1982. Multicomponent equilibrium chemistry in groundwater quality models. Water Resour. Res. 18:1089–1096.
Jury, W.A., H. Frenkel, and L.H. Stolzy. 1978. Transient changes in the soil-water system from irrigation with saline water: I. Theory. Soil Sci. Soc. Am. J. 42:579–585.
Karickhoff, S.W., D.S. Brown, and T.A. Scott. 1979. Sorption of hydrophobic pollutants on natural sediments. Water Res. 13: 241–248.
Keating, E.H., and J.M. Bahr. 1998. Reactive transport modeling of redox geochemistry: Approaches
to chemical disequilibrium and reaction rate estimation at a site in northern Wisconsin. Water
Resour. Res. 43:3573–3584.
Kirkner, D.J., and H. Reeves. 1988. Multicomponent mass transport with homogeneous and heterogenous chemical reactions: Effect of the chemistry on the choice of numerical algorithm. 1. Theory. Water Resour. Res. 24:1719–1729.
Knapp, R.B. 1989. Spatial and temporal scales of local equilibrium in dynamic fluid-rock systems.
Geochim. Cosmochim. Acta. 53:1955–64.
Konikow, L.F., and J.D. Bredehoeft. 1992. Ground-water models cannot be validated. Adv. Water Resour. 15:75–83.
Lagerwerff, J.V., F.S. Nakayama, and M.H. Frere. 1969. Hydraulic conductivity related to porosity andswelling of soil. Soil Sci. Soc. Am. Proc. 33:3–11.
Lichtner, P.C. 1985. Continuum model for simultaneous chemical reactions and mass transport in hydrothermal systems. Geochim. Cosmochim. Acta 49:779–800.
Lichtner, P.C. 1993. Scaling properties of time-space kinetic mass transport equations and the local equilibrium limit. Am. J. Sci. 293:257–296.
Lichtner, P.C. 1996. Continuum formulation of multicomponent-multiphase reactive transport. p. 1–81.
In P.C. Lichtner et al. (ed.) Reactive transport in porous media. Reviews in Mineralogy. Vol. 34.
Mineralogical Society of America, Washington, DC.
Lichtner, P.C. 1998. Modeling reactive flow and transport in natural systems. p. 3–72. In L. Martini and
G. Ottonello (ed.) Proceedings of the Rome Seminar on Environmental Geochemistry. Dipartimento di Scienze dell Terra-Universita di Genova, Italy.
Lichtner, P.C. 2000. FLOWTRAN User’s manual. Los Alamos National Laboratory Document, NM.

MISCIBLE SOLUTE TRANSPORT

1533

Lichtner, P.C., and M.S. Seth. 1996. User’s manual for MULTIFLO: Part II MULTIFLOW 1.0 and GEM
1.0. Multicomponent-multiphase reactive transport model. CNWRA 96-010. Center for Nuclear
Waste Regulatory Analyses, San Antonio, TX.
Liu, C.W., and T.N. Narasimhan. 1989a. Redox-controlled multiple-species reactive chemical transport.
1. Model development. Water Resour. Res. 25:869–882.
Liu, C.W., and T.N. Narasimhan. 1989b. Redox-controlled multiple-species reactive chemical transport.
2. Verification and application. Water Resour. Res. 25:883–910.
Mangold, D.C., and C.-F. Tsang. 1991. A summary of subsurface hydrological and hydrochemical models. Rev. Geophys. 29:51–79.
Mattigod, S.V., and J.M. Zachara. 1996. Equilibrium modeling in soil chemistry. p. 1309–1358. Methods of soil analysis. SSSA Book Ser. 5. SSSA, Madison, WI.
McNeal, B.L. 1968. Prediction of the effect of mixed-salt solutions on soil hydraulic conductivity. Soil
Sci. Soc. Am. Proc. 32:190–193.
McNeal, B.L. 1974. Soil salts and their effects on water movement. p. 409–431. In J. van Schilfgaarde(ed.) Drainage for agriculture. Agron. Monogr. 17. ASA, Madison, WI.
Miller, C.W., and L.V. Benson. 1983. Simulation of solute transport in a chemically reactive heterogeneous system: Model development and application. Water Resour. Res. 19:381–391.
Miller, C.T., M.M. Poirier-McNeill, and A.S. Mayer. 1990. Dissolution of trapped nonaqueous phaseliquids: Mass transfer characteristics. Water Resour. Res. 26:2783–96.
Narasimhan, T.N., A.F. White, and T. Tokunaga. 1986. Groundwater contamination from an inactive
uranium mill tailings pile. 2. Application of a dynamic mixing model. Water Resour. Res.
22:1820–1834.
Oreskes, N., K. Shrader-Frechette, and K. Belitz. 1994. Verification, validation, and confirmation of numerical models in the earth sciences. Science 263:641–646.
Papelis, Ch., K.F. Hayes, and J.O. Leckie. 1988. HYDRAQL: A program for the computation of chemical equilibrium composition of aqueous batch systems including surface-complexation modeling of ion adsorption at the oxide/solution interface. Technical Rep. 306. Dept. of Civil Engineer., Stanford University, Stanford, CA.
Parkhurst, D.L., and C.A.J. Appelo. 1999. User’s guide to PHREEQC (Version 2)—A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. Water Resources Investigation Report 99-4259. U.S. Geol. Surv., Denver, CO.
Parkhurst, D.L., D.C. Thorstenson, and L.N. Plummer. 1980. PHREEQE—A computer program for geochemical calculations. U.S. Geol. Survey, Water Resources Investigations, Report 80-96.
Pignatello, J.J. 1989. Sorption dynamics of organics by soils and sediments. p. 243–270. In B.L. Sawhney and K. Brown (ed.) Reactions and movement of organic chemicals in soils. SSSA Spec. Publ.
22. SSSA and ASA, Madison, WI.
Pitzer, K.S. 1979. Activity coefficients in electrolyte solutions. p. 157–208. In R.M. Picowicz (ed.) Activity coefficients in electrolyte solutions. CRC Press, Boca Raton, FL.
Plummer, L.N., and E. Busenberg. 1982. The solubilities of calcite, aragonite and vaterite in CO2–H2O
solutions between 0 and 90°C, and an evaluation of the aqueous model for the system
CaCO2–CO2–H2O. Geochim. Cosmochim. Acta 46:1011–1040.
Plummer, L.N., T.M. Wigley, and D.L. Parkhurst. 1978. The kinetics of calcite dissolution in CO2 systems at 5° to 60°C and 0.0 to 1.0 atm CO2. Am. J. Sci. 278:179–216.
Powers, S.E., C.O. Loureiro, L.M. Abriola, and W.J. Weber. 1991. Theoretical study of the significanceof nonequilibrium dissolution on nonaqueous phase liquids in subsurface systems. Water Resour. Res. 27:463–477.
Reiniger, P., and G.H. Bolt. 1972. Theory of chromatography and its application to cation exchange insoils. Neth. J. Agric. Sci. 20:301–313.
Rhoades, J.D. 1997. Sustainability of irrigation: An overview of salinity problems and control strategies. p. 1–42. CWRA 1997 Annual Conf. Footprints of Humanity: Reflections on Fifty Yearsof Water Resource Developments” Lethbridge, Alberta, Canada. 3–6 June 1997.
Rittmann, B.E., and J.M. VanBriesen. 1996. Microbiological processes in reactive modeling. p. 311–334.
In P.C. Lichtner et al. (ed.) Reactive transport in porous media. Reviews in Mineralogy Vol. 34.
Mineralogical Society of America, Washington, DC.
Robbins, C.W., R.J. Wagenet, and J.J. Jurinak. 1980. A combined salt transport-chemical equilibrium
model for calcareous and gypsiferous soils. Soil Sci. Soc. Am. J. 44:1191–1194.
Rubin, J. 1983. Transport of reacting solutes in porous media: Relation between mathematical natureof problem formulation and chemical nature of reactions. Water Resour. Res. 19:1231–1252.
Rubin, J., and R.V. James. 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:1332–1356.

1534

CHAPTER 6

Russo, D. 1988. Numerical analysis of the nonsteady transport of interacting solutes through unsaturated soil: 1. Homogeneous systems. Water Resour. Res. 24:271–284.
Russo, D., and E. Bresler. 1977. Analysis of the saturated-unsaturated hydraulic conductivity in amixed Na–Ca soil system. Soil Sci. Soc. Am. J. 41:706–710.
Sardin, M., D. Schweich, F.J. Leij, and M.Th. van Genuchten. 1991. Modeling the nonequilibrium transport of linearly interacting solutes in porous media—A review. Water Resour. Res. 27:2287–2307.
Selim, H.M., J.M. Davidson, and R.S. Mansell. 1976. Evaluation of a two-site adsorption–desorption
model for describing solute transport in soils. p. 444–448. In Proc. Summer Computer Simulation Conference. American Institute of Chemical Engineers. Washington, DC.
Shainberg, I., and G.J. Levy. 1992. Physico-chemical effects of salts upon infiltration and water movement in soils. p. 37–93. In R.J. Wagenet et al. (ed.) Interacting processes in soil science. Lewis
Publishers, CRC Press, Boca Raton, FL.
Šimçnek, J., and D.L. Suarez. 1993a. The UNSATCHEM-2D code for simulating two-dimensional variably saturated water flow, heat transport, carbon dioxide transport, and solute transport with major
ion equilibrium and kinetic chemistry. Version 1.1. Research Report No. 128. U.S. Salinity Laboratory, USDA, ARS, Riverside, CA.
Šimçnek, J., and D.L. Suarez. 1993b. Modeling of carbon dioxide transport and production in soil: 1.
Model development. Water Resour. Res. 29:487–497.
Šimçnek, J., and D.L. Suarez. 1994. Two-dimensional transport model for variably saturated porous
media with major ion chemistry. Water Resour. Res. 30:1115–1133.
Šimçnek, J., and D.L. Suarez. 1997. Sodic soil reclamation using multicomponent transport modeling.
ASCE J. Irrig. Drain. Eng. 123:367–376.
Šimçnek, J., D.L. Suarez, and M. Šejna. 1996. The UNSATCHEM software package for simulating onedimensional variably saturated water flow, heat transport, carbon dioxide production and transport, and multicomponent solute transport with major ion equilibrium and kinetic chemistry. Version 2.0. Research Report No. 141. U.S. Salinity Laboratory, USDA, ARS, Riverside, CA.
Šimçnek, J., and M.Th. van Genuchten. 1994. The CHAIN_2D code for simulating two-dimensional
movement of water flow, heat, and multiple solutes in variably-saturated porous media. Version
1.1. Research Report No. 136. U.S. Salinity laboratory, USDA, ARS, Riverside, CA.
Sparks, D.L. 1989. Kinetics of soil chemical processes. Academic Press, New York, NY.
Sposito, G., and S.V. Mattigod. 1980. GEOCHEM: A computer program for the calculations of chemical equilibria in soil solutions and other natural systems. Dep. Soil Environ. Sci., Univ. California, Riverside, CA.
Sposito, G. 1981. The thermodynamics of soil solutions. Oxford Univ. Press, New York, NY.
Sposito. G. 1984. The surface chemistry of soils. Oxford Univ. Press. New York, NY.
Steefel, C.I., and K.T.B. MacQuarrie. 1996. Approaches to modeling of reactive transport in porous
media. p. 83–129. In P.C. Lichtner et al. (ed.) Reactive transport in porous media. Reviews in
Mineralogy Vol. 34. Mineralogical Society of America, Madison, WI.
Steefel, C.I., and P. Van Cappellen. 1998. Reactive transport modeling of natural systems. J. Hydrol.
209:1–7.
Steefel, C.I., and S.B. Yabusaki. 1996. OS3D/GIMRT, software for multicomponent-multidimensional
reactive transport. User Manual and Programmer’s Guide, PNL-11166. Pacific Northwest National Laboratory, Richland, WA.
Stokes, R.H. 1979. Thermodynamics of solutions. p. 1–28. In R.M. Pitkowicz (ed.) Activity coefficientsin electrolyte solutions. CRC Press, Boca Raton, FL.
Stumm, W., and J.J. Morgan. 1996. Aquatic chemistry: Chemical equilibria and rates in natural waters.
3rd ed. Wiley, New York, NY.
Suarez, D.L., and J. Šimçnek. 1993. Modeling of carbon dioxide transport and production in soil: 1.
Parameter selection, sensitivity analysis, and comparison of model predictions to field data. Water
Resour. Res. 29:499–513.
Tanji, K.K., L.D. Doneen, G.V. Ferry, and R.S. Ayers. 1972. Computer simulation analysis on reclamation of salt-affected soil in San Joaquin Valley, California. Soil Sci. Soc. Am. Proc. 36:127–133.
Tebes-Stevens, C., and A.J. Valocchi. 1998. Numerical solution techniques for reaction parameter sensitivity coefficients in multicomponent subsurface transport models. Hydraulic Engineering Series No. 59. UILI-ENG-98-2009. Department of Civil and Environmental Engineering, University
of Illinois at Urbana-Champaign.
Tebes-Stevens, C, A.J. Valocchi, J.M. VanBriesen, and B.E. Rittmann. 1998. Multicomponent transportwith coupled geochemical and microbiological reactions: model description and example simulations. J. Hydrol. 209:8–26.
Truesdell, A.H., and B.F. Jones. 1974. Wateq, a computer program for calculating chemical equilibria
of natural waters. J. Res. U.S. Geol. Surv. 2:233–248.

MISCIBLE SOLUTE TRANSPORT

1535

Valocchi, A.J. 1985. Validity of the local equilibrium assumption for modeling sorbing solute transport
through homogeneous soils. Water Resour. Res. 21:808–820.
Valocchi, A.J. 1988. Theoretical analysis of deviations from local equilibrium during sorbing solute transport through idealized stratified aquifers. J. Contam. Hydrol. 2:191–207.
Valocchi, A.J., and M. Malmstead. 1992. Accuracy of operator splitting for advection–dispersion-reaction problems. Water Resour. Res. 28:1471–1476.
Valocchi, A.J., R.L. Street, and P.V. Roberts. 1981. Transport of ion-exchanging solutes in groundwater: Chromatographic theory and field simulation. Water Resour. Res. 17:1517–1527.
VanBriesen, J. 1998. Modeling coupled biogeochemical processes in mixed waste systems. Ph.D. thesis. Civil Engineering. Northwestern University, Evanston, IL.
van der Zee, S.E.A.T.M., and G. Destouni. 1992. Transport of inorganic solutes in soil. p. 95–145. In
R.J. Wagenet et al. (ed.) Interacting processes in soil science. Lewis Publishers, Boca Raton, FL.
van der Zee, S.E.A.T.M., and W.H. Riemsdijk. 1994. Transport of reactive solutes in soils. p. 1–103. In
M.Y. Carapcioglu (ed.) Advances in porous media. Vol. 2. Elsevier, Amsterdam, the Netherlands.
van Genuchten, M.Th., and R.W. Cleary. 1979. Movement of solutes in soil: Computer simulation andlaboratory results. p. 349–386. In G.H. Bolt (ed.) Soil chemistry B. Physico-chemical models.
Elsevier, New York, NY.
van Genuchten, M.Th., and J. Šimçnek. 1996. Evaluation of pollutant transport in the unsaturated zone.
p. 139–172. In P.E. Rijtema and V. Eliáš (ed.) Regional approaches to water pollution in the environment. NATO ASI Series: 2. Environment. Kluwer, Dordrecht, The Netherlands.
van Genuchten, M.Th., and R.J. Wagenet. 1989. Two-site/two-region models for pesticide transport anddegradation: Theoretical development and analytical solutions. Soil Sci. Soc. Am. J.
53:1303–1310.
Viswanathan, H.S., B.A. Robinsons, A.J. Valocchi, and I.R. Triay. 1998. A reactive transport model of
neptunium migration from potential repository at Yucca Mountain. J. Hydrol. 209:251–280.
Wagenet, R.J., and J.L. Hutson. 1987. LEACHM: Leaching Estimation And CHemistry Model, A
process-based model of water and solute movement, transformations, plant uptake and chemical reactions in the unsaturated zone. Continuum 2. Water Resour. Inst., Cornell University, Ithaca,
NY.
Walsh, M.P., S. L. Bryant, R. S. Schechter, and L. W. Lake. 1984. Precipitation and dissolution of solids
attending flow through porous media. AIChE J. 30:317–328.
Walter, A.L., E.O. Frind, D.W. Blowes, C.J. Ptacek, and J.W. Molson. 1994. Modeling of multicomponent reactive transport in groundwater 1. Model development and evaluation. Water. Resour.
Res. 30:3137–3148.
Walter, A.L., E.O. Frind, D.W. Blowes, C.J. Ptacek, and J.W. Molson. 1994. Modeling of multicomponent reactive transport in groundwater 2. Metal mobility in aquifers impacted by acidic minetailings discharge. Water. Resour. Res. 30:3149–3158.
Weber, W.J., P.M. McGinley, and L.E. Katz. 1991. Sorption phenomena in subsurface systems: Concepts, models, and effects on contaminant fate and transport. Water Res. 25:429–528.
Westall, J.C. FITEQL. 1982. A computer program for determination of equilibrium constants from experimental data. Version 2.0, Rep. 82-02. Dept. Chem., Oregon State Univ., Corvallis, OR.
Westall, J.C., J.L. Zachary, and F.M.M. Morel. 1976. MINEQL: A computer program for the calculation of chemical equilibrium composition of aqueous systems. Tech. Note 18. Dep. of Civ. Eng.,
Mass. Inst. of Technol., Cambridge, MA.
White, N., and L.W. Zelazny. 1986. Charge properties in soil colloids. p. 39–81. In D. L. Sparks (ed.)
Soil physical chemistry. CRC Press, Boca Raton, FL.
Wolery, T.J. 1992. EQ3NR, a computer program for geochemical aqueous speciation-solubility calculations: Theoretical manual, user’s guide, and related documentation. Version 7. Publ. UCRLMA-110662, PT III. Lawrence Livermore National Laboratory, Livermore, CA.
Yaron, B., and G.W. Thomas. 1968. Soil hydraulic conductivity as affected by sodic water. Water Resour. Res. 4:545–552.
Yeh, G.T., and D.D. Huff. 1985. FEMA: A finite element model of material transport through aquifers.
Rep. ORNL-6063. Oak Ridge Natl. Lab., Oak Ridge, TN.
Yeh, G.T., K.M. Salvage, J.P. Gwo, J.M. Zachara, and J.E. Szecsody. 1998. HYDROBIOGEOCHEM:
A coupled model of hydrological transport and mixed biochemical kinetic/equilibrium reactions
in saturated-unsaturated media. Rep. ORNL/TM-13668. Oak Ridge Natl. Lab., Oak Ridge, TN.
Yeh, G.T., and V.S. Tripathi. 1989. A critical evaluation of recent developments in hydrogeochemicaltransport models of reactive multichemical components. Water Resour. Res. 25:93–108.
Yeh, G.T., and V.S. Tripathi. 1990. HYDROGEOCHEM: A coupled model of hydrological and geochemical equilibrium of multicomponent systems. Rep. ORNL-6371. Oak Ridge Natl. Lab., Oak
Ridge, TN.

1536

CHAPTER 6

Yeh, G.T., and V.S. Tripathi. 1991. A model for simulating transport of reactive multispecies components: Model development and demonstration. Water Resour. Res. 27:3075–3094.
Zysset, A., F. Stauffer, and T. Dracos. 1994a. Modeling of chemically reactive groundwater transport.
Water Resour. Res. 30:2217–2228.
Zysset, A., F. Stauffer, and T. Dracos. 1994b. Modeling of reactive ground water transport governed bybiodegradation. Water Resour. Res. 30:2423–2434.
Zyvoloski, G.A., B.A. Robinson, Z.V. Dash, and L.L. Trease. 1997. Summary of the models and methods for the FEHM application—A finite element heat- and mass-transfer code. Rep. LA-13307MS. Los Alamos National Laboratory, Los Alamos, NM.