Comput Geosci (2015) 19:445–478
DOI 10.1007/s10596-014-9443-x

ORIGINAL PAPER

Reactive transport codes for subsurface environmental
simulation
C. I. Steefel · C. A. J. Appelo · B. Arora · D. Jacques · T. Kalbacher · O. Kolditz ·V. Lagneau · P. C. Lichtner · K. U. Mayer · J. C. L. Meeussen · S. Molins ·D. Moulton · H. Shao · J. Šimůnek · N. Spycher · S. B. Yabusaki · G. T. Yeh

Received: 24 April 2014 / Accepted: 26 August 2014 / Published online: 26 September 2014
© Springer International Publishing Switzerland 2014

Abstract A general description of the mathematical andnumerical formulations used in modern numerical reactive
transport codes relevant for subsurface environmental simulations is presented. The formulations are followed by short
descriptions of commonly used and available subsurface
simulators that consider continuum representations of flow,transport, and reactions in porous media. These formulations are applicable to most of the subsurface environmental
benchmark problems included in this special issue. The listof codes described briefly here includes PHREEQC, HPx,
PHT3D, OpenGeoSys (OGS), HYTEC, ORCHESTRA,TOUGHREACT, eSTOMP, HYDROGEOCHEM, CrunchFlow, MIN3P, and PFLOTRAN. The descriptions include a

C. I. Steefel () · B. Arora · S. Molins · N. Spycher
Earth Sciences Division, Lawrence Berkeley National Laboratory,
1 Cyclotron Road, Berkeley, CA 94720, USAe-mail: CISteefel@lbl.gov
C. A. J. Appelo
Hydrochemical Consultant, Valeriusstraat 11, 1071 MB
Amsterdam, the Netherlands
e-mail: appt@hydrochemistry.eu
D. JacquesInstitute for Environment, Health, and Safety, Belgian
Nuclear Research Center, Mol, Belgium
T. Kalbacher · O. Kolditz · H. Shao
Department of Environmental Informatics, Helmholtz Centre
for Environmental Research, Leipzig, Germany
V. Lagneau
University, Centre de Géosciences, 35 rue Saint Honoré,
77305 Fontainebleau Cedex, France
P. C. Lichtner
Lichtner OFM Research, Santa Fe, NM, USA

high-level list of capabilities for each of the codes, alongwith a selective list of applications that highlight theircapabilities and historical development.
Keywords Reactive transport · Modeling · Environmental
simulation · Computer software · Code benchmark

1 Introduction
In this paper, we present a relatively general description of the mathematical and numerical formulations usedin modern numerical reactive transport codes based on

K. U. Mayer
Earth and Ocean Sciences, University of British Columbia,
Vancouver, Canada
J. C. L. Meeussen
WU Environmental Sciences, University of Wageningen,Wageningen, the Netherlands
D. Moulton
Mathematical Modeling and Analysis, Los Alamos National
Laboratory, Los Alamos, NM, USA
J. Šimůnek
Department of Environmental Sciences, UC Riverside,
Riverside, CA, USA
S. B. Yabusaki
Earth Systems Science Division, Pacific Northwest National
Laboratory, Richland, WA, USA
G. T. Yeh
National Central University, Jhongli, Taiwan

446

continuum representations relevant for subsurface environmental simulations. These formulations should be applicable to most of the benchmark problems included in thisspecial issue, although the diversity of formulations andoptions to be found in the simulators is beyond the scope ofa single manuscript. It is felt that the presentation of mathematical and numerical formulations is important, however,so as to avoid the use of either the simulators or the benchmarks as “black boxes”. Moreover, a presentation of theformulations gives readers and future users a clearer idea ofthe limits of the benchmark problems and the model comparisons shown. For example, one expects different results
for the benchmark problems using a pore scale rather than acontinuum representation, which is the basis for the benchmarks included in this special issue [1–3]. Similarly, the useof a different set of rate laws (e.g., an atomistic rather thana transition state theory formulation) will also most likely
produce different results. Thus, the assumption is that thebenchmarks presented in the manuscripts that follow areintended for continuum reactive transport simulators thatutilize in large part the mathematical and numerical formulations described below. Where important differences in
the formulations occur, these are detailed in the individual
manuscripts.
This is followed with brief descriptions of the codes usedin the benchmarking problems included in this special issue.
Again, these descriptions are not exhaustive—the interested
reader is referred to the publications cited for the individual
codes, or to user’s manuals where publication has not beenpossible. The focus of this special issue is on the presentation of benchmark problems and not the codes themselves.
However, we feel it is useful to present at least a brief summary of the capabilities and histories of the various codes
now in common use in the environmental simulation community. To first order, this should provide the reader with
an idea of which simulators they might want to explore
further.
A master set of tables of the subsurface environmental
simulation capabilities for each of the codes is also included.
Table 1 summarizes the flow and transport capabilities of
the codes, while Table 2 summarizes their geochemical
and microbial process capabilities. Table 3 provides a shortlist of the numerical features, including code parallelization, if implemented. The list of continuum-based codes
considered here includes PHREEQC, HPx, PHT3D, OpenGeoSys (OGS), HYTEC, ORCHESTRA, TOUGHREACT,eSTOMP, HYDROGEOCHEM, CrunchFlow, MIN3P, and
PFLOTRAN. These are many of the most commonly used,widely available reactive transport codes, although others
exist. Perhaps more importantly in this context, these are thecodes that were used in the benchmark problems included inthis special issue. An effort has been made to present a highlevel list of capabilities for each of the codes, along with a

Comput Geosci (2015) 19:445–478

summary of the historical development that highlights code
capabilities first and foremost.

2 Governing equations
2.1 Continuum approach
A rock or sediment or soil mass consisting of aggregates of mineral grains and pore spaces or voids isreferred to as a porous medium. An actual porous medium
is a highly heterogeneous structure containing physical discontinuities marked by the boundaries of porewalls which separate the solid framework from the void
space. Although it is possible, in principle at least, todescribe this system at the microscale of individual pores(e.g., [1–3]), such a description rapidly becomes difficult as the size of the system increases and many
pore volumes become involved. It is therefore necessary
to approximate the system by a more manageable one.
One quantitative description of the transport of fluids andtheir interaction with rocks is based on a mathematical
idealization of the real physical system referred to as a continuum. In this theory, the actual discrete physical system,
consisting of aggregates of mineral grains, interstitial porespaces, and fractures, is replaced by a system in which
physical variables describing the system vary continuously
in space. Allowance is made for the possibility of a discrete set of surfaces across which discontinuous changes inphysical properties may occur. In this fictitious representation of the real physical system, solids and fluids coexist
simultaneously at each point in space. This is the basis
for all of the benchmarks presented in this special issue.
Benchmarks for pore scale representations of flow andtransport have been presented recently in Oostrom et al. [4].
A representative selection of the most important governing equations for continuum subsurface environmental
reactive transport simulators are given below.
2.2 Reaction processes
A range of equilibrium and kinetically controlled biogeochemical reaction processes are typically consideredin modern reactive transport simulators, including (1)mineral dissolution and precipitation reactions, (2) homogeneous (aqueous phase) reactions, (3) multisite surface
complexation diffuse double layer or triple layer models,with or without an electrostatic correction for the energy
associated with charging of the mineral surface, (4) multisite ion exchange reactions, (5) gas–aqueous phaseexchange, and (6) microbially mediated homogeneous andheterogeneous
reactions.

Comput Geosci (2015) 19:445–478

447

2.2.1 Homogeneous phase reactions
Homogeneous or intra-aqueous phase reactions include bothclassical aqueous complexation reactions, often assumed
to be at equilibrium, and kinetically controlled aqueous
phase reactions (e.g., microbially mediated sulfate reduction, radioactive decay).
Equilibrium reactions If we assume that the various aqueous species are in chemical equilibrium, it is possible toreduce the number of independent concentrations, that is,
the number that actually need to be solved for. Mathematically, this means that in a system containing Ntot aqueous
species, the number of independent chemical components inthe system Nc is reduced from the total number of species by
the Nx linearly dependent chemical reactions between them
[5–12]. This leads to a natural partitioning of the system
into Nc primary or basis species, with concentrations designated here as Cj , and the Nx secondary species, referred toas Ci [8–10]. The equilibrium chemical reactions between
the primary and secondary species take the form
Ai 

Nc


(1)

νij Ajj =1

where Aj and Ai are the chemical formulas of the primaryand secondary species, respectively, and νij is the numberof moles of primary species j in one mole of secondary
species i. The equilibrium reactions provide an algebraic
link between the primary and secondary species via the lawof mass action for each reaction
Ci = Ki−1 γi−1

Ncγj C jνij(i = 1, ..., Nx )

(2)

j =1

where γj and γi are the activity coefficients for the primary and secondary species, respectively, and Ki is theequilibrium constant for the reaction, written here as theconsumption of one mole of the secondary species i. Here,the assumption is that the activities of the species are theproduct of the species concentrations and their activity
coefficients:
a j = γj C j(3)This partitioning between primary and secondary species isalso used as the basis for the definition of total concentrations, j , which are the actual conserved quantities in mostof the reactive transport simulators:
j = Cj +Nxνij Ci .

(4)

i=1

Kinetic reactions When the aqueous phase reactions are
kinetically controlled, then the reactions must be written interms of primary species alone. Treatment of redox disequilibrium is a good example. In this case, the reactions can be
written with a transition state theory (TST) type rate law asQs  nRs = k s 1 −
ai
(5)
Keq
where ks is the rate constant for the reaction, Qs is the ion
activity product defined by
Qs = ai−1Ncν

aj ij ,

(6)

j =1

where aj and ai are the activities of the primary and secondary species in the reaction, and the assumption is thatthe secondary species, i, has a stoichiometric coefficient of
1 on the left side of the reaction. In addition, the product ofvarious species (the number of which depends on the particular system) affecting the rate far from equilibrium (e.g.,H+ ) is included as the last term in Eq. 5, with n representing the exponential dependence on the catalytic or inhibitory
species. Alternatively, the reaction can be treated as completely irreversible, in which case there is no “affinity” orfree energy term (1 − Qs /Keq ):

ain
(7)
Rs = k s
This more general form of irreversible reactions is not limited by any thermodynamic equilibrium constraints and canbe used, for example, to model decay and ingrowth ofradioactive daughter products.
2.2.2 Gas reactions
Reactions with a discrete gas phase can be included witheither an equilibrium or kinetic treatment. In the equilibrium
case, mass action equations are similar in form to those forthe aqueous complexes, except insofar as the definitions ofconcentration and activity corrections change. The fugacityof a gas species can be written as
fl = Kl−1

Nc


ν
γj Cj lj ,

(8)

j =1

where fl is the fugacity of the lth gas. The fugacity of thegas is related to the partial pressure of the gas through thefugacity coefficient, ϕl :
Pl =

fl
.
φl

(9)

In some cases, the assumption of an ideal gas is valid, inwhich the concentration of the gas species is calculated fromthe ideal gas law:
nl
Pl
=
V
RT

(10)

448

Comput Geosci (2015) 19:445–478

where nl and Pl are the number of moles and partialpressures of the lth gas.
2.2.3 Surface complexation
The treatment of surface complexation models in modern reactive transport simulators is now standard, essential
in the view of many because it allows for the treatmentof sorption as influenced by variable chemical conditions[13, 14]. In this approach, the sorbing sediment surfaces
are considered to possess surface functional groups thatcan form complexes analogous to the formation of aqueous complexes in solution. These surface reactions include
proton exchange, cation binding, and anion binding via ligand exchange at surface hydroxyl sites (represented here
as XOH). For example, the sorption of a metal could be
represented as
X•OH + Mz+  X•OMz+-1 + H+ ,(11)

according to



zF ψ0
Keq = Kint exp
RT


(14)

where Kint is the intrinsic equilibrium constant that does notdepend on the surface charge.
Kinetic treatments of surface complexation are becoming
more common, particularly as it is observed that the equilibrium treatments cannot capture the slow release of metals
or contaminants noted in natural soils and aquifers [18].
2.2.4 Ion exchange
Ion exchange reactions can be described via a massaction expression with an associated equilibrium constant[19–22]. The exchange reaction can be written in genericform, assuming here a chloride solution, asvAClu (aq) + uBXν (s) ↔ uBClν (aq) + vAXu (s)(15)

where z+ is the valence of the metal. The mass action
equation for this reaction is given by

where X refers to the exchange site occupied by the cations
Au+ and Bv+ . The equilibrium constant or selectivity coefficient, Keq , for this reaction can be written as [18][X•OM z+ − 1 ][H+ ]Keq =.
[X•OH ][M z+ ]

Keq =

(12)

where the square brackets, [], refer to activities. In Eq. 12,
Keq is the equilibrium constant, which in the case of theelectrostatic models described below is corrected further for
the free energy of surface ionization [15], while the concentrations of the surface complexes are typically written interms of mole fractions or coverage fractions [14, 16].
It is possible to consider either an electrostatic modelin which Coulombic forces associated with charging of thesurface are used to correct the equilibrium constant, or a
non-electrostatic model in which surface charge (whether
present or not) is not used. Non-electrostatic surface complexation is usually implemented based on the generalized composite model where it is assumed that it can be
described in terms of “generic” surface functional groups[17]. In the case of the electrostatic surface complexation
model in which the electrical double layer (EDL) consistingof an electrical potential that decreases rapidly away fromthe charged mineral surface is accounted for, the free energy
associated with the electrical potential, Gcoul , is added tothe intrinsic free energy associated with chemical bonding,GintrGads = Gintr + Gcoul= Gintr − zF ψ0(13)

and ψ0 is the mean surface potential (V ). The equilibrium
constant thus can be adjusted for the electrostatic effects[BClν ]u [AXu ]ν
[AClu ]ν [BXν ]u(16)

where the brackets, [], refer to the thermodynamic activities.
Several activity conventions are commonly incorporated inmodern reactive transport simulators. One possibility is theGaines–Thomas activity convention, which assumes a reaction stoichiometry of the following form [21], written here
for the case of Na+ and Ca2+ exchange:(17)
Na+ + 1 2CaX2 ↔ NaX + 1 2Ca2+
where X refers to the exchange site. In the Gaines-Thomas
convention, each exchange site, X, has a charge of −1
and the activities of adsorbed species correspond to thecharge equivalent fractions. The Gapon activity convention is obtained by writing the reactions in every case witha single exchanger [21]. Alternatively, the Vanselow convention [19] describes the exchanger activity with mole
fractions.
2.2.5 Mineral dissolution and precipitation
Minerals may be treated as at equilibrium in a similar fashion to the aqueous complexes, although unless colloidal,
they are considered as immobile. Including the mineral
concentrations, Cm , in the total component concentrations
yields:j = Cj +

Nx

i=1

νij Ci +Nm
k=1

νkj Cm(18)

Comput Geosci (2015) 19:445–478

449

Alternatively, mineral dissolution and precipitation can betreated kinetically, which is often an acceptable practicesince mineral–water rates are generally much slower thanintra-aqueous reaction rates. Local equilibrium behavior
is possible by using sufficiently large rate constants, ofcourse—in all cases, the kinetic treatment of reactions is themore general [12]. The types of mineral reactions can bedivided into (1) TST type reactions (linear or nonlinear) in
which dissolution and precipitation are treated as reversibleat equilibrium such that a dependence on Gibbs energy or
saturation state is explicitly included, (2) irreversible reactions that include a dependence on Gibbs energy or saturation state, (3) irreversible reactions in which a Gibbs energy
or saturation state dependence is not included, although adependence on concentration can be, and (4) microbiallymediated reactions treated with a combination of thermodynamic and kinetic (Monod) terms.
Transition state theory rate laws A general (linear andnonlinear) form for the TST type rate laws for mineraldissolution and precipitation [23, 24] is given by
Rm = sgn [ ] Am km= sgn [ ] Am km





m
η G
−1
RT
η

m
Qm
−1
Keq


an
an

exp

(19)

  
where sgn [ ] = sgn log Qm Km gives the sign ofthe reaction term, negative for dissolution, positive for precipitation, Am is the reactive surface area of the reacting
km is the rate constantmineral (m2 m−3 porous medium),(mol m−2 s−1 ), and the term a n is product of all of the
far from equilibrium effects on the reaction (e.g., pH) that
may be catalytic or inhibitory (this is the product over all
species, with the assumption of a zero exponent if there is noeffect on the reaction rate). Qm and Keq are the ion activity
product and equilibrium constant, respectively. The parameter n can be interpreted as the reaction order, but it may bea non-integer value determined experimentally because of
the complicated reaction mechanisms in play. The last termon the right hand side is referred to as the “affinity term”
and gives the dependence of the rate on Gibbs energy or
saturation state, with η and m representing experimentallyor theoretically determined exponential dependences. For
example, η can be thought of as the inverse of the Temkincoefficient, which relates the stoichiometric coefficient forthe species in the rate-limiting reaction to its stoichiometriccoefficient in the mineral.
Dissolution-only or precipitation-only mineral rate laws In
Eq. 19, the assumption is that the reaction can be reversed
at equilibrium, with a continuous rate between dissolution
and precipitation. Another possibility is to use the same rateequation, but as applied to dissolution alone or precipitation
alone. This could be the approach in the case of a dissolution
rate law that is not continuous on either side of equilibrium(supersaturated or undersaturated) point, that is, the form of
the dissolution and precipitation rate laws are not the same.
Or it may be possible that there is no or an extremely slow
reverse reaction, as appears to be the case for the mineralquartz at 25 °C. In this case, the rates are given by



m

η G
n
Rm = sgn [ ] Am km
a
exp
−1 ;
RT
G < 0Rm = 0;G > 0.

(20)This kind of rate law or combination of rate laws canpresent problems for the reactive transport simulators, since
the derivative of the reaction rate that goes into the Jacobian (matrix of function derivatives) is not continuously
differentiable.
Irreversible mineral rate laws Alternatively, it is possible toconsider irreversible reactions (i.e., those with no affinity or
G dependence), with or without a dependence on species
concentration or activity

(21)
an ,
Rm = Am km
although the thermodynamic basis for such rate laws is lessclear, and implementation may result in reactions proceeding past their thermodynamic equilibrium point.
2.2.6 Microbially mediated mineral rate laws
Microorganisms play a critical role in many subsurface
environments, in particular, where bioremediation is considered for immobilization of redox-sensitive contaminantsby amendment with organic substrates. Microorganisms
obtain energy for growth and maintenance from redox reactions. Microbially mediated redox reactions involve twopathways, a catabolic pathway used for energy production
and an anabolic pathway used for cell synthesis, the sumof which is the overall stoichiometry of the microbiallymediated reaction [25]:
0 =

Nc


fs0 νjs k + fe0 νjek Ak

(22)

k

in which νjs k and νjek refer to the stoichiometric coefficientsfor the anabolic (cell synthesis) and catabolic (energy) pathways, respectively, and Ak refer to the species involved in
the reactions. The portion of electron equivalents used forcell synthesis (0 < fs0 < 1) is an input parameter thatdefines the growth yield of the reaction. The remaining portion of electron equivalents (fe0 = 1−fs0 ) defines its energy
yield.

450

Comput Geosci (2015) 19:445–478

The rate of microbially mediated reactions can be written
as a function of its maximum intrinsic rate (μmax , in unitsof mol per biomass per time), the biomass concentration(B), kinetic limitations (0 < F K < 1) and thermodynamic
limitations (0 < FT < 1)
Rm = μmax BFK FT(23)

Kinetic limitations are obtained using Monod terms forelectron donor (ED ) and acceptor (EA )



[ED ]
[EA ]
FK =
(24)
KED + [ED ]KEA + [EA ]Thermodynamic limitations affect the overall rate at low
substrate concentrations or for low-energy yielding reactions [26]

 

Gr + mGATP(25)
FT = max 0, 1 − expχRTwith mGATP being a minimum energy barrier, or energy
quantum, that organisms must overcome to catalyze thereaction.

only to isopycnic (“iso-density”) flow:
h = z +

Written in terms of the hydraulic conductivity and head,
Darcy’s Law is given by
q = −K∇h.

∂ φρf
∂t

= ∇ · ρf q + ρf Qa(26)

where ∇· is the divergence operator, q is the volumetricflux (or Darcy flux) of water [m3 H2 O m−2 medium s−1 ],φ is the porosity [m3 void m−3 medium], ρf is the fluid
density [kg m−3 H2 O], Qa is volumetric source term
[m3 H2 O m−3 medium s−1 ], and t is time [s]. This is a general formulation that is applicable to variable density flow
involving a single phase. Single phase flow in permeable
subsurface porous media is typically described with Darcy’s
law [27, 28]. The volumetric or Darcy flux is defined interms of the fluid pressure gradient:q = φv = −ksat∇P − ρf gez ,μ(27)

where v is the pore velocity [m s−1 ], z is the depth [m],∇P [kg m−1 s−2 , or Pa m−1 ] is the fluid pressure gradient, g is the acceleration due to gravity [m s−2 ], andksat [m2 ], ez is a unit vector in the vertical direction, and[m2 ] μ [kg s−1 m−1 ] are the permeability tensor for fully
saturated conditions and dynamic viscosity, respectively.
Alternatively, one can recast the equations in terms of thehydraulic head, h [m], although this is strictly applicable(30)

2.3.2 Variably saturated flow
Variably saturated flow is typically described either withthe full set of multiphase flow and conservation equations,or with the use of the Richards equation, which assumes apassive air phase.
Richards equation Variably saturated flow can be treatedusing Richards equation [29, 30] under the assumption of apassive air phase and neglecting hysteresis:
Sa Ss

The equation for fluid continuity in the case of single phaseflow is given by [27]:(28)

One can also define the hydraulic conductivity tensor, K[m s−1 ], asksat ρf g
K =
.
(29)
μ

2.3 Fluid flow
2.3.1 Single phase flow

P
.
ρf g

∂Sa
∂h
+ φ
= ∇ · [kra K∇h] + Qa ,
∂t
∂t

(31)

where Sa [m3 H2 O m−3 void] defines the saturation of
the aqueous phase, Ss [m−1 ] is the specific storage coefficient and kra [−] is the relative permeability. The equivalent
expression written in terms of pressure is∂ Sa φρf
kra ksat ∇P − ρf gez + ρf Qa .
= ∇· −ρf∂tμ(32)Multiphase flow The component conservation equations forthe most general case of multiphase flow can be written foreach component j in phase α:∂ φ ρ α S α Yj α
α

∂t


kra ksatρ α Yj α=∇· −(∇P − ραgez ) + ρα Qj
μ
α


(33)
where Yj α is the mass fraction of component j in phase α
and the variables are as defined above, except for the moregeneral case of phase α. Thus, a component like CO2 or
H2 O can be tracked between the various phases, gas, liquid,and solid.
Relative permeability and saturation formulations The relative permeability and aqueous phase saturation can beComput Geosci (2015) 19:445–478

451

calculated using the soil hydraulic functions based on van
Genuchten [31]:1 − SraSa = Sra +(34)(1 + αψan )m

2
1/m m
l
1
−
1
−
S
kra = SSea
(35)
eaa − Sra(36)
Sea =1 − Srawhere Sra [−] defines the residual saturation of the aqueous
phase, ψa is the pressure head [m] α[m−1 ], and n, m, andl are soil hydraulic function parameters, with m = 1 −1 . Sea is the effective saturation of the aqueous phase. Theneffective saturation of the aqueous phase can also be related
to the capillary pressure, Pc [Pa] [31]:Sea = 1 + (α |Pc |)n

−m

.

(37)

Other formulations for capillary pressure and aqueous phasesaturation are available, e.g., the Brooks–Corey curve [32]:|Pc | λ
(38)
Sea =
Pe
where Pe is the entry pressure [Pa] and λ is a fitting parameter. For the gas phase, the following formulation for relative
permeability has been suggested:krg = 1 − kra ,where Seg is defined in terms of the liquid saturation, Sa ,and the residual gas saturation, Srg , asSa.
1 − Srg(41)

2.4 Molecular diffusion
The simplest treatment of molecular diffusion is in termsof Fick’s first law, which states that the diffusive flux ofa species in solution is proportional to the concentration
gradient:
Jj = −Dj ∇Cj ,Jj = −uj Cj

Dj Cj ∂μj
∂μj
= −
∂x
RT ∂x

uj =

Dj
RT

(45)

.

To derive the Nernst–Planck equation, we can use a definition of the chemical potential that includes the explicit
contribution of the electrochemical potential:
μj = μ0j + RT ln aj + zj F ψ= μ0j + RT ln γj Cj + zj F ψ(46)where μ0j is the standard chemical potential, F is Faraday’s
constant, zj is the species charge, and ψ is the electricalpotential. By substituting Eq. 46 into Eq. 44, we obtain thefollowing:
Jj = −Dj ∇Cj −zi F
Di Ci ∇ψ − Dj Cj ∇ ln γj
RT

(47)

where T is the absolute temperature, and γj is the activity coefficient of the j th species. The first term on the righthand side of Eq. 47 is the Fickian diffusion term. The second
term on the right hand side is the gradient in electrostatic
potential (sometimes referred to as the electrochemical
migration term). The electrical potential may develop gradients due to a number of processes (e.g., the application of anelectrical current), but as captured in geochemical transport
codes at the present time, this effect typically arises due tothe diffusion of charged species at different rates (the diffusion potential). The last term in Eq. 47 is the contribution
from gradients in activity coefficients, which is presently
treated in only PHREEQC Nernst–Planck equation implemented in modern subsurface reactive transport simulators:(42)

where Dj is the diffusion coefficient, which is specific tothe chemical species considered as indicated by the subscript j . In the case of diffusion in porous media, it isnormally necessary to include a tortuosity correction as well(discussed further below). Integrating the diffusive fluxes
over a control volume leads to Fick’s second law [27]:∂Cj= −∇ · Jj = ∇ · Dj ∇Cj ,(43)∂twhich provides an expression for the change in concentration in terms of the divergence of the diffusive flux.

(44)

where uj is the species mobility [28] defined by thefollowing:(39)

but a formulation proposed by Luckner et al. [33] is often
preferred:1 3 

1 m 2m(40)
krg = 1 − Seg / 1 − Seg/

Seg =Fick’s law is a strictly phenomenological relationshipthat is more rigorously treated with the Nernst–Planck equation. The full version of the Nernst–Planck equation can bederived from considering the gradient in chemical potential,µj , as the driving force for diffusion:

Jj = −Dj ∇Cj −zi F
Di Ci ∇ψ.
RT(48)

2.5 Tortuosity
Tortuosity is defined as the square of the ratio of the path
length the solute would follow in water alone, L, relative tothe tortuous path length, it would follow in porous media,
Le [27]:
TL =

2
L
Le

(49)

452

Comput Geosci (2015) 19:445–478

With this formulation, the diffusion coefficient in porous
media is given by

2.7 Reactive transport equations

Di∗ = TL Di

Combining the transport and biogeochemical reaction network terms results in the reactive transport equation(50)

The diffusive flux, then, is given by
Jidiff = −φDi TL ∇Ci = −φDi∗ ∇Ci

(51)

Nr



∂ (φSL Ci )∗= ∇ · φSL Di ∇Ci − ∇ · (qCi ) −νir Rr∂tr=1where Di∗ is the diffusion coefficient for species i in porous
media.
Various approaches for calculating diffusion coefficientsin porous medium, Di∗ , are in use, with a formulation basedon Archie’s Law employing an exponential dependence on
the porosity, being perhaps the most common:

−

Nm


Ng

νim Rm −

m=1



(56)

νil Rll=1

where m is the cementation exponent. The pore aqueousand gas phase effective diffusion coefficients, Dea and Deg ,are also defined locally based on porosity and saturation
following the relationship given respectively by Millington[34]:

where the term on the left hand side is the accumulation
term, with the porosity and liquid saturation, SL , multipliedby the concentration yielding units of mole per cubic meter
medium per second. The reaction terms are partitioned
between aqueous phase reactions, Rr , mineral reactions,
Rm , and gas reactions, Rl , all considered as kinetically controlled. Where equilibrium aqueous complexation reactions
are present (the usual case), it is possible to define a total
concentration, as in Eq. 4, and write the reactive transport
equation as

Dea = φ 1/3 Sa / D0a73
Deg = φ 1/3 Sg / D0g

Nr



∂ (φSL i )= ∇ · φSL Di∗ ∇i − ∇ · (qi ) −νir Rr∂tDi∗ = φ m Di

(52)

73

(53)

r=1

where D0a and D0g are the free phase diffusion coefficientsin the aqueous and gas phases, respectively.

−

Nm


Ng

νim Rm −m=1

2.6 Mechanical dispersion
Mechanical dispersion is defined as the product of the fluid
velocity and dispersivity, α, with longitudinal and transverse
components
DL = αL Vi
DT = αT Vi(54)where Vi refers to the average velocity in the principal
direction of flow, and the subscripts L(longitudinal) andT (transverse) refer to the dispersion coefficient paralleland perpendicular to the principal direction of flow, respectively [27]. In tensorial index notation, the dispersion tensor
becomes
Dij = αT V δij + (αL − αT ) Vi Vj V
(55)
where V is the average velocity magnitude, δij is the Kronecker delta, and Vi and Vj are the velocity magnitudes
in the I and J coordinate directions, respectively [27].
Only a few of the reactive transport simulators have considered a full dispersion tensor, since the expressions becomecomplicated when flow is not parallel to the principal coordinate directions of the numerical grid. In addition, crossterms can only be considered in finite element formulations, which have the disadvantage of not being locally massconservative.

νil Rl(57)

l=1

The Fickian treatment of dispersion can be added to Eqs. 56
and 57 in a similar way as Fickian diffusion. The Fickian
treatment of diffusion in Eq. 57 can be replaced with the
Nernst–Planck equation.
Equations 56 and 57 can be expanded to consider transport in phases other that the liquid by adding the corresponding accumulation and transport terms [35, 36].
2.8 Reaction-induced porosity and permeability change
2.8.1 Porosity changes
Porosity changes in matrix and fractures are directly tied
to the volume changes as a result of mineral precipitationand dissolution. The molar volumes of minerals createdby hydrolysis reactions (i.e., anhydrous phases, such as
feldspars, reacting with aqueous fluids to form hydrous minerals such as zeolites or clays) are often larger than thoseof the primary reactant minerals; therefore, constant molar
dissolution–precipitation reactions may lead to porosity
reduction.
The porosity, φ, of the medium (fracture or matrix) can
be calculated fromφ = 1 −

Nm

m=1

φm

(58)

Comput Geosci (2015) 19:445–478

where Nm is the number of minerals, and φm is the volume
fraction of mineral m in the rock (m3 mineral m−3 medium).
As the volume fraction of each mineral changes due to mineral reactions, the porosity can be recalculated at each timestep.
2.8.2 Fracture permeability changes
Fracture permeability changes can be approximated using
the porosity change and an assumption of plane parallel
fractures of uniform aperture [12]. The modified permeability, k, is then given by 3
φk = k0(59)
φ0where k0 and φ0 are the initial permeability and porosity, respectively. This law yields zero permeability only
under the condition of zero fracture porosity. In most experimental and natural systems, permeability reductions tovalues near zero occur at porosities significantly greaterthan zero. This generally is the result of mineral precipitation in the narrower interconnecting fracture apertures. This can be described potentially with a higherorder dependence of the permeability on the porosity, orwith a percolation threshold model for the permeability[37].
2.8.3 Matrix permeability changes
Matrix permeability changes can be calculated fromchanges in porosity using ratios of permeabilities calculated from the Carman–Kozeny relation [27] and ignoringchanges in grain size, tortuosity, and specific surface area asfollows: (1 − φ0 )2 φ 3(60)
k = k0(1 − φ)2 φ0
The simple cubic law (59) and the Kozeny–Carman (60)porosity–permeability equations may not reflect the complex relationship of porosity and permeability in geologic
media that depends on an interplay of many factors, suchas pore size distribution, pore shapes, and connectivity [38].
Laboratory experiments have shown that modest reductions in porosity from mineral precipitation can cause large
reductions in permeability [39].
2.8.4 Effects of permeability and porosity changes
on capillary pressures
Permeability and porosity changes will likely result in modifications to the unsaturated flow properties of the rock.
Changes to unsaturated flow properties are approximated bymodification of the calculated capillary pressure (Pc ) using453

the Leverett scaling relation [40] to obtain a scaled Pc asfollows:
Pc = Pc

k0 φ
kφ0

(61)

2.9 Multicontinuum representations
Reliable predictions of flow and transport in fractured
porous media are important to address several issues,
including groundwater contamination, hydrothermal venting, degradation of concrete, carbon sequestration, and mining and mineralization activities. The existence of fracturesand heterogeneity in porous media, which are often hierarchical in nature, has led to the multicontinuum class ofmodels. Some of the commonly used continuum-scale models include (1) equivalent continuum model (ECM) [41];(2) dual permeability model (DPM), dual or multiporosity model [42, 43], multiple interacting continua approach(MINC) [44], and (3) discrete fracture and matrix model[45]. Among these three commonly used approaches, thedual-continuum approach has been extensively applied indifferent subsurface environments [46–48]. Broadly speaking, the dual-continuum approach considers two interacting regions, one associated with the less permeable intraaggregate pore region, or the soil matrix domain, and theother associated with the inter-aggregate pore region, or thefracture domain. In this approach, the REV is partitioned
into sub-volumes of each domain such that
VREV = Vf + Vm

(62)

where VREV refers to the volume of the porous medium,
Vf and Vm refer to the volume of the fracture and matrix
domains, and subscripts f and m refer to the fracture andmatrix domains, respectively. Thus, the porosity of eachdomain is given by
φf =

φf
φm
, φm =
VREV
VREV

(63)

As expected, these relations can be easily extended toinclude multiple interacting domains. The dual or multidomain conceptualizations can be distinguished by thedifferent formulations of the governing equations of flow in
the fracture domain and/or the different approaches throughwhich the exchange between the fracture and rock matrix(or multiple domains) is established. A review of the different multicontinuum approaches including the governing
equations and exchange functions are available elsewhere[49–52].

454

Comput Geosci (2015) 19:445–478

2.10 Thermal processes
Although thermal processes are not the focus of the benchmarks presented in this special issue, we present the governing equation used in most of the non-isothermal codes that
actually solve for the energy balance. One form of the governing equation for energy conservation in a porous medium
is given by



∂
φ
Sα ρα Uα + (1 − φ) ρr cr T∂tα+ ∇ · qα ρα Hα − κ∇T = Qe .

(64)

In this equation, T refers to temperature and the subscript α
represents a fluid phase with Darcy velocity α , density ρα ,saturation Sα , internal energy α , and enthalpy Hα . The coefficient κ denotes the thermal conductivity of the medium
and cr and ρr refer to the specific heat and density of theporous medium. The quantity Qe denotes a source/sink
term. Internal energy and enthalpy are related by the equatio

reaction equations. As described by Steefel and MacQuarrie [54], this could consist of 100 % kinetic reactions, inwhich case the number of conservation equations equals
the number of species in the system. Alternatively, whenequilibrium can be assumed for some of the aqueous or surface complexes, or for a mineral phase, the reaction rate
can be eliminated and replaced with a mass action equation (2). The algebraic equations can be included together
with the Nc conservation equations for evolution of thetotal primary or component species concentrations (4) as afunction of time (t), in which case, the number of equations (differential and algebraic) to be solved still equals
the number of species, complexes, and minerals in thesystem.
n+1

Ns
Nm


1
Cj +
νij Ci +νmj Cm
t
i=1

m=1

Ns


Nm



− Cj +

νij Ci +i=1

n= 0 j = 1, · · · Nc

νmj Cm

m=1

(67)
pα
.
Uα = Hα −ρα(65)

Thermal conductivity is often described by the phenomenological relationship given by Somerton et al. [53] κ = κdry + Sl κsat − κdry
(66)
where κsat and κdry are the fully saturated and dry thermal
conductivities and Sl is the liquid saturation.

log Cin + 1 =Nsνij log γj Cjn+1− log γin + 1

i=1− log Ki i = 1, · · · Ns
n+1
log Cm=Nmνmj log γjn + 1 Cjn + 1

m=1− log Km m = 1, · · · Nm
3 Numerical formulations
A complete review of the possible numerical formulations
is well beyond the scope of this paper. However, we present
here a brief summary of some of the more commonlyused approaches, which rely in most cases on Newton–
Raphson methods for the nonlinear reaction networks foundin modern environmental simulators. Although not considered here, similar approaches to coupling flow, solutetransport, and energy equations could be included as well.
Typically, these are handled via operator splitting, but global
implicit (one-step) approaches are possible.
3.1 Numerical approaches for biogeochemical reactions
The general Newton–Raphson approach is summarized in
Steefel and MacQuarrie [54]: a series of typically coupled conservation equations for the primary or component
species (4) are solved together with kinetic and equilibrium(68)

where n and n + 1 refer to the present and future timelevel, respectively, the γı ’s and γj ’s are the activity coefficients for the secondary and primary species, respectively,
Ns is the number of secondary species or complexes in equilibrium with the primary or component species, and Nm
is the number of minerals treated as equilibrium phases.
A similar approach can be taken to include equilibrium
surface complexes and gases. In addition, the system canbe further partitioned between species or mineral phases
linked to equilibrium or kinetically controlled reactions. As
shown here, the algebraic mass action equations can be linearized by treating the logarithms of the concentrations as
the primary unknowns.
An alternative method is referred to as the direct substitution approach (DSA), which can be used to eliminate
the algebraic equations corresponding to the mass action
equations. In this case, we obtain a nonlinear set of ordinary differential equations for the unknown concentrationsof the primary or component species, written here to include

Comput Geosci (2015) 19:445–478

455

explicitly a set of kinetic reactions involving both aqueous species (intra-aqueous kinetic reactions) and mineralphases:⎤n+1 ⎡
⎤n
⎡
Nc
Nc
Ns
Ns




1 ⎣
νij γi−1 Ki−1
(γj Cj )νij ⎦− ⎣Cj +
νij γi−1 Ki−1
(γj Cj )νij ⎦
Cj +
t
i=1

−

−

⎧
Ns
⎨
⎩

s=1

⎧
Nm
⎨
⎩

m=1

j =1

j =1

i=1

⎞η

⎛

m ⎫n+1
⎬

Nc




sgn log Qs Ks νis ks ⎝(γj Cj )νmj Ks−1 ⎠ −1
⎭
j =1

⎞η
⎛
m ⎫n+1
Nc
⎬




−1 ⎠
sgn log Qm Km νim Am km ⎝(γj Cj )νmj Km−1= 0 j = 1, · · · Nc .
⎭

(69)

j =1

Note that in this formulation, an implicit solved at the n + 1
time level is carried out for both equilibrium and kinetic
reactions.
Other approaches are possible in which kinetic and equilibrium reactions are decoupled. This approach allows forthe use of free energy minimization routines, for example[54–56].
3.2 Approaches for the reactive transport problem
Modern subsurface environmental reactive transport simulators typically use as many as three different approaches
to handle the reactive transport equations numerically. The
most common approach is operator splitting of the reaction and transport terms, usually in a sequential non-iterative
manner. In this approach, the transport of the conserved
quantities (typically the total concentrations) is solved,
and this is followed by a node by node reaction solve.
In the sequential iterative approach (SIA), this sequential solve of the transport and reaction is repeated untilconvergence is achieved. In the sequential non-iterative
approach (SNIA), however, only a single solve of thetransport and reaction terms is carried out. The otherapproach, referred to as the global implicit or one-stepmethod, reaction and transport are solved simultaneously[54, 57].
The chief advantage of the global implicit approach isthat, unlike the operator splitting method for the reactive
transport problem, one is not restricted to time steps lessthan the Courant condition. The Courant condition would
require that mass not be transported more than a single gridcell in any one time step. Even in the case of a Courant
number, Cr = 1 operator splitting errors have been reported[58].

3.2.1 Global implicit approach
In the one-step or global implicit approach, the transportand reaction terms are solved simultaneously. One approach
would be to solve directly for the total concentrations (thej ’s) in any one iteration and then follow this with a
distribution of species calculation to determine the concentrations of the individual species. Another method isbased on the DSA and consists of solving simultaneously
for the complexation (which are assumed to be at equilibrium), the heterogeneous reactions, and transport terms.
This means that the primary species (the Cj ’s) rather than
the total concentrations are the unknowns. Following thenotation of Lichtner [59], we can introduce the differential
operatorL(j ) =∂∗(φSL ) + ∇ • q − D ∇ j ,
∂t

(70)

and write the governing differential equations in terms ofthe total component concentrations as
aq

L(j ) = Rjmin + Rj (j = 1, ..., Nc ),(71)

where the reaction term Rjmin includes only kinetic reactions. Note that in this formulation, the diffusion coefficients are assumed to be the same for all of the aqueous
species (both primary and secondary) using a Fickian treatment of diffusion.
We proceed by using Eqs. 70 and 2 to obtain an
expression for the total soluble concentration in terms of the456Comput Geosci (2015) 19:445–478

primary species alone. Substituting this result into Eq. 71
gives

the nodal points. The Jacobian matrix in the case of onedimensional transport takes a block tridiagonal form⎡⎡
L ⎣Cj +Ns

νij γi−1 Ki−1

sgn [m ] νim Am km

⎝

Nc


⎞η

m

−1 ⎠
(γj Cj )νmj Km
−1

j =1

m=1

+

(γj Cj )νij ⎦

⎛

Nm


Ns


⎤

j =1

i=1

+

Nc


⎛
sgn [

s ] νis ks

⎝

Nc


⎞η
νmj

(γj Cj )Ks−1 ⎠

m

−1

=0

j =1

s=1

j = 1, · · · Nc ,(72)

where sgn [ m ] and sgn [ s ] give the sign of the reaction
term for the kinetic mineral and aqueous complex reactions,
respectively.
3.2.2 Nonlinear solver method for global implicit
Using a global implicit approach becomes more difficult in the case of a multicomponent, multispecies system, both because the coupling of species via reactions
enlarges the size of the coefficient matrix and because
it typically results in sets of nonlinear equations whichmust be solved. The size of the Jacobian matrix whichmust be constructed and solved becomes larger, since eachfunction will include contributions from the concentrations both at that nodal point itself and from neighboring
nodal points which are used in the finite difference discretization. For example, in the case of one-dimensionaltransport and Nc unknown concentrations at each nodal
point, the form of the Newton equations to be solved
is [12, 54]⎤
A1,1 A1,2 0 · · ·⎢ A2,1 A2,2 A2,3 · · ·
⎥
⎢
⎥
⎢
⎥
···
⎢
⎥
⎣
· · · AN−1,N−2 AN−1,N−1 AN−1,N ⎦···0
AN,N−1 AN,N
⎡
⎡
⎤
⎤
δC1
f1
⎢ δC2 ⎥
⎢ f2 ⎥
⎢
⎢
⎥
⎥
⎢
⎥
⎥
×⎢
(74)
⎢ ··· ⎥ = −⎢ ··· ⎥⎣ δCN−1 ⎦⎣ fN−1 ⎦δCN
fN
where the entries in the Jacobian matrix (the A) are submatrices of dimension Nc by Nc in the case where there are
Nc unknowns per nodal point. The entries δCi refer here
to the entire vector of unknown concentration corrections at
any particular nodal point, while the functions fi include theentire vector of equations for the unknown concentrations at
each nodal point. In practice, the logarithms of the primary
species concentrations are solved for rather than the concentrations themselves, which introduces a natural scaling of
the numerical approach that greatly improves convergenceproperties.
3.2.3 Operator splitting approach
Operator splitting of reaction and transport is an alternative
to the global implicit approach [12, 54]. In this approach,a single time step consists of a transport step followed bya reaction step using the transported concentrations. Physically, one can think of the method as equivalent to taking
part of the contents of one beaker (i.e., one node) andadding it to a second beaker downstream. After this physical process is completed, reactions occur which modify
the chemical composition of each of the beakers in thesequence. Mathematically, the method can be represented asa two step sequential process consisting of a transport steptransp

Ci
Nc

∂fi,j
k=1

∂Ck,j

+

Nc

k=1

δCk,j +

Nc

k=1

− Cin

t

∂fi,j
δCk,j + 1
∂Ck,j +1

∂fi,j
δCk,j − 1 = −fi,j∂Ck,j − 1

(i = 1, ..., Ntot ),(75)= Rin + 1 (i = 1, ..., Ntot )(76)= L(Ci )n ,followed by a reaction step
Cin + 1 − Citranspt(73)

where i refers to the component number, k is the unknown
component species number, and j , j + 1, and j − 1 are

where L is the spatial operator. The problem can also be
formulated using the total concentrations, thus reducing thenumber of equations which must be solved. In the casewhere the Courant number is 1, then the entire contents of
one grid cell are transferred to the next one downstream,
and the approach is equivalent to that of a mixing cell model

Comput Geosci (2015) 19:445–478

when explicit time weighting and a formulation without
numerical dispersion is used.

4 Code descriptions
The diversity of process models (Section 3) and numericalapproaches (Section 4) used in the simulation of reactive
transport in subsurface environments makes it nearly impossible for any code to include all potential capabilities or
to run on all computer architectures, from workstations to
distributed memory parallel systems. Further, code development has typically occurred incrementally over several
years, in some cases over a decade or more, driven by theneeds of research projects. The reactive transport codes used
for this benchmarking effort reflect the heterogeneity resulting from this process; even when re-engineered from theground up, the codes are the legacy of years of development
and research. In some instances, development of native geochemical capabilities from scratch has been circumvented
by coupling flow and transport codes, typically using an
operator splitting approach, to existing geochemical capabilities. In the code overview that follows, we start with
PHREEQC, a geochemical code that is used as geochemical engine in HPx, PHT3D, and OpenGeoSys. HYTEC uses
the same concept to bring in the geochemical capabilities
of CHESS. We follow with codes that have native geochemical modules. ORCHESTRA is distinctive in that allmodel definitions are separate from the calculation engine
or equation solver. TOUGHREACT and eSTOMP use theoperator splitting approach, while HYDROGEOCHEM,
CrunchFlow, MIN3P, and PFLOTRAN allow for implicit
coupling of transport and biogeochemistry.
4.1 PHREEQC
PHREEQC is a computer program developed by Parkhurst
and Appelo designed to perform a variety of aqueous
geochemical calculations [60, 61]. Originally a batchreaction model, the current capabilities of PHREEQC
include inverse modeling and one-dimensional transport
calculations with a wide range of physical and chemical
processes. PHREEQC can simulate a range of equilibrium and kinetic reactions. Several aqueous models are
implemented in the databases available for PHREEQC,
including ion-association models that use Debye–Hückel
formulations (phreeqc.dat, wateq4f.dat, llnl.dat, minteq.dat,
minteq.v4.dat, and iso.dat); the Pitzer specific-ion interaction model (pitzer.dat), and the specific-ion interaction
theory (SIT) model (sit.dat). All databases include temperature dependence of activity-coefficient constants and van’t
Hoff or analytical expressions for equilibrium constants
[60].

457

Heterogeneous equilibria that can be imposed include
surface complexation, ion exchange, fixed-composition
gases and minerals, solid solutions, and a multicomponent gas phase. Surface complexation can be calculated byusing the diffuse-double layer model or the CD-MUSICmodel. In addition, it is possible to calculate explicitly thecomposition of the diffuse layer that balances the charge
accumulated by surface complexation. Solid solutions may
be ideal-multicomponent or nonideal-binary solid solutions.
The pitzer.dat and phreeqc.dat databases have additionalfeatures. Density and specific conductance are calculated
as functions of the solution composition. Molar volumesof aqueous species are calculated as functions of pressure,temperature, and ionic strength, which allow calculationof the pressure dependence of equilibrium constants [62].
These databases use the Peng–Robinson formulation for
nonideal gas fugacities and can provide accurate simulations of mineral and gas solubilities up to 200 °C and1,000 bar for major elements [62]. It may be noted thatthe calculation of density-driven flow requires composition, temperature, and pressure-dependent solution densities
such as calculated by PHREEQC. The iso.dat database
implements an approach to isotopes where each isotope
is considered to be a separate thermodynamic component [63]. The database allows simulation of kinetic andequilibrium fractionation of isotopes of carbon, hydrogen,nitrogen, oxygen, and sulfur and can be extended to otherisotopes.
Homogeneous or heterogeneous kinetic reactions can be
defined with arbitrary rate expressions that are functions ofsolution composition and availability of solid reactants byusing Basic programs that are evaluated with an embedded
Basic interpreter. The general kinetic reaction capability of
PHREEQC can be used for any biological reactions. It is
possible to formulate any rate expression, including Monodor kinetics of any order, inhibition factors, rates that dependon free energy, and rate variation as a function of available electron acceptors. It is also possible to consider theeffects of a consortium of bacteria that have specific substrate requirements and individual rates of consumption,
growth, and decay.
Solute-transport calculations may consist of advection
and dispersion or pure diffusion. Advection is simulatedwith the maximal time step allowed by the Courant criterion, while dispersion is simulated with sub-time steps whenrequired by the Von Neumann criterion. Reactions are calculated after each advective and dispersive step. The method
can accurately simulate advective–dispersive systems andcan propagate a square wave with no numerical dispersion
for unretarded solutes. In diffusive calculations, a single
diffusion coefficient can be used for all aqueous species,
or individual diffusion coefficients can be used for eachaqueous species subject to the constraint of no net transport

Yes
Yes

No
Diagonal
Yes
Yes
No
Yes

Yes
Diagonal
No
No
Yes
Yes

No
No
No

No
No
No

Yes
Yes

Yes
Yes

Yes
No

No
Diagonal
No
No
No
No

Yes
Yes

No
Yes
No

Yes
No

No
Diagonal
Yes
Yes
Yes
Yes

Yes
Yes

No
Yes
Yes

Yes
Yes

1,2,3D

Yes
Diagonal
Yes
Yes
Yes
Yes

Yes
Yes

No
No
No

Yes
Yes

1D

No
No
Yes
Yes
Yes
Yes

Yes
Yes

Yes
Yes
Yes

Yes
Yes

1,2,3DNo
Diagonal
Yes
Yes
No
No

Yes
Yes

Yes
Yes
Yes

Yes
Yes

1,2,3D

No
Full
Yes
Yes
Yes
No

Yes
Yes

Yes
Yes
Yes

Yes
Yes

1,2,3D

Yes
Diagonal
Yes
Yes
No
Yesa

Yes
Yes

No
Yes
No

Yes
No

1,2,3D

Yes
Diagonal
Yes
Yes
No
Yesb

Yes
Yes

No
Yes
Yes

Yes
Yes

1,2,3D

ORCHESTRA TOUGHREACT eSTOMP HYDROGEOCHEM CrunchFlow MIN3P

b Not available in MIN3P-THCm, but dual porosity model is available in MIN3P-Dual [203]No
Diagonal
Yes
Yes
Yes
Yes

Yes
Yes

Yes
Yes
No

Yes
Yes

1,2,3D

OpenGeoSys HYTEC

a Not available in CrunchFlow, but available in CrunchFlowEDL

Transport
Advection
Molecular diffusion
Electrochemical
migration
Dispersion tensor
Gas phase advection
Gas phase diffusion
Colloids
Multiple continua

Flow
Saturated flow
Richards equation
Multiphase–
multicomponent flow
Variable density flow
Non-isothermal flow

1,2,3D

1D

Dimensions

1,2,3D

PHREEQC HP1/HPx PHT3D

Capabilities/features

Table 1 A comparison of the key flow and transport features of reactive transport modeling codes

No
Diagonal
Yes
Yes
Yes
Yes

Yes
Yes

Yes
No
Yes

Yes
Yes

1,2,3DPFLOTRAN

458
Comput Geosci (2015) 19:445–478

Yes

Yes

DDL,
Non-edl
Yes
Yes

Yes

Yes

Yes

Yes

Yes

Yes
Yes

Yes

Yes

DDL,
Non-edl
Yes
Yes

Yes

Yes

Yes

Yes

Yes

Yes
Yes

Yes
Yes
Yes

Yes

Yes

Yes
Yes
Yes

Yes

Yes

Yes
Yes
Yes

Yes
Yes

Yes

Yes

Yes

Yes

Yes

DDL,
Non-edl
Yes
Yes

Yes

No

Yes

Yes

PHREEQC HP1/HPx PHT3D

Yes
Yes
Yes

Yes
Yes

No

Yes

Yes

Yes

Yes

DDL,
Non-edl
Yes
Yes

Yes
Yes
Yes

Yes
Yes

No

Yes

No

No

Yes

DDL,
CC
Yes
Yes

Yes

No

Yesa
Yes

Yes

Yes

Yes

Yes

Yes
Yes
Yes

Yes
Yes

Yes

Yes

Yes

Yes

Yes

Yes
Yes

DDL

Yes

No

Yes

No

Yes
Yes
Yes

Yes
Yes

Yes

Yes

Yes

No

Yes

DDL,
CC
Yes
Yes

Yes

No

No

Yes

Yes
Yes
Yes

Yes
Yes

No

No

Yes

No

Yes

DDL,
Non-edl
Yes
Yes

Yes

Yes

Yes

Yes

Yes
Yes
Yes

Yes
Yes

Yes

Yes

Yes

No

Yes

DDL,
Non-edl
Yes
Yes

Yes

No

Yes

Yes

Yes
Yes
Yes

Yes
Yes

Yes

Yes

Yes

Yes

Yes

DDL,
Non-edl
Yes
Yes

Yes

No

Yes

Yes

Yes

No

No

Yes

PFLOTRAN

Yes
Yes
Yes

Yes
No
Yes

Yes
Yes

No

Yesc
Yes
Yes

No

No

No

Yes

Yes
Yes

No

No

No

Yes

Yes
Yes

Non-edl Non-edl

Yes

Yes

Yes

Yes

OpenGeoSys HYTEC ORCHESTRA TOUGHREACT eSTOMP HYDROGEOCHEM CrunchFlow MIN3P

c Not available in MIN3P-THCm, but available in a customized version of MIN3P [227]b DDL denotes Diffuse Double Layer, non-edl denotes Non-Electrostatic model, and CC denotes Constant Capacitancea Not available in OGS, but OGS has built-in interfaces to different geochemical and biogeochemical solvers (GEM, PHREEQC, ChemApp, BRNS) that can implement Pitzer activity model

Microbial
Monod kinetics
Thermodynamic
Biomass growthGeochemistry
Extended DebyeHückel
Davies activity
model
Pitzer activity
model
Non-isothermalgeochemistrySurfacecomplexationb
Ion exchange
Aqueous–gas
exchange
Kinetic mineral
precipitation–
dissolution
Mineralnucleation
Mineralsolid-solutions
Equilibrium
isotopefractionation
Kinetic isotope
fractionation
Aqueous kinetics
Radioactive
decay chains

Capabilities/
features

Table 2 A comparison of the key geochemical and microbial features of reactive transport modeling codes

Comput Geosci (2015) 19:445–478
459

No
Yes
Yes
Yes

Computational
Parallelization
Inverse estimation
Graphical user interface
Open source
NofWith PESTg
Yes
No
Yes
With PESTg
Yes
Yes

FVM, MMC FEM
TW
BE

Yes
No
Yes (FCT)c

No

Yes

Yes
Yes

Yes

Yes
No
No
No

FVM
TW

Yes
No
No

No

Yes

No
No

Yes

Partial
Yes
No
No

FVM
BE/CN

Yes
No
No

No

Yes

Yes
C-E

Yes

Yes
No
No
No

FVM
BE

Yes
No
Yes

No

Yes

No
No

Yes

Yes
Yes
Yes

Yes

Yes

No
C-E

Yes

No
No
No
Yes

No
With PESTg
No
No

FEM, MMC FVM
TW,Mid
BE

Yes
No
Yes

Yes

Yes

Yes
EGM

Yes

Yes
With PESTg
No
No

FVM
BE

No
Yes
Yes

No

No

No
Nob

Yes

Yes
With PESTg
No
Yes

FVM
BE

Yes
Yes
Yes

No

Yes

Yes
Yes

Yes

PFLOTRAN

and RK denotes Runge-Kutta
f An MPI-based parallel version of PHT3D has been developed and is in the testing phase [228]g PEST denotes Parameter Estimation Software Toolkit [229]e BE denotes Backwards Euler, TW denotes time weighted difference including forward difference, backward difference, and Crank-Nicolson (CN) Central Difference, Mid denotes Mid difference,d MC denotes mixing cell, FEM denotes Finite Element Method, FVM denotes Finite Volume Method and MMC denotes Modified Method of Characteristic or backward particle tracking

c FCT refers to linearized algebraic flux corrected transport

b Simplified 1-D approach following Neuzil [214]

Yes
No
Yes
No

MC
FE/ RK

Yes
No
Yes

Yes

Yes

No
No

Yes

OpenGeoSys HYTEC ORCHESTRA TOUGHREACT eSTOMP HYDROCrunchFlow MIN3P
GEOCHEMa C-E denotes compaction-erosion, and EGM denotes Equilibrium Geo-Mechanics ModelNo
With PESTg
Yes
No

FEM
TW

MC
TW

Yes
No
Yes

No

No

Yes
No
No

Yes

Yes

Yes
No
Yes

No
No

No
No

Numerical Scheme
Operator splitting
Global implicit
High Peclet number
transport
Spatial discretizationd
Time discretizatione

No

PHT3D

Yes

PHREEQC HP1

Coupling
Reaction-inducedNo
porosity–permeability
Coupled heat of reaction No
Coupled deformation/
No
compactiona
Reaction consumptionYesof phase (H2 O, CO2 )Electrical double layer Yes
transport

Capabilities/features

Table 3 A comparison of the key numerical capabilities of reactive transport modeling codes

460
Comput Geosci (2015) 19:445–478

Comput Geosci (2015) 19:445–478

of charge, essentially solving the Nernst–Planck equation
given in Eq. 48. The porosity of a cell can be updated in
multicomponent diffusion simulations to account for dissolution or precipitation of minerals. Colloids, which can reactas surface-complexation units, can be transported, coagulated, and resuspended depending on system properties andsolute composition. Diffusion of interlayer (exchangeable)
cations and of solutes in the diffuse layer of surfaces can be
simulated.
Nearly any set of geochemical parameters and physicalproperties can be output to tab-delimited files or graphedwith a built-in plotting capability that can produce multiple
real-time plots simultaneously. Input for PHREEQC can be
defined and run with a specialized module for Notepad++,or with one of two graphical user interfaces (PhreeqcI andPhreeqc for Windows; the latter uses PHREEQC version
2.18).
The wide range of geochemical capabilities offered by
PHREEQC, its open-source code, as well as continued support and development by the U.S. Geological Survey havemade it attractive for many researchers to couple it as a geochemical module to existing flow and transport codes thatexceed PHREEQC’s native transport capabilities. Specifically for this purpose, PHREEQC has been encapsulated
in the module IPhreeqc [64] to allow all PHREEQC functions to be called from C, C++, Fortran, or from scripting
languages (IPhreeqcCOM), such as Excel macros, Python,or Matlab. Coupling is performed with a sequential noniterative approach (SNIA), i.e., the same method as usedin PHREEQC to couple transport and reactions. Some ofthese couplings are used to simulate unsaturated zone processes such as HPx [65, 66], saturated zone environments
such as PHAST [64, 67] and PHT3D [68], and many morethat are beyond the scope of this overview. At the time ofthis writing (April, 2014), PHREEQC is in its version 3
release. When the coupling to PHREEQC is based on theearlier release (version 2), it is noted as PHREEQC-2 in thedescriptions below. All PHREEQC-related programs canbe found at wwwbrr.cr.usgs.gov/projects/GWC coupled orhydrochemistry.eu/downl.html.
4.2 HPx (HP1, HP2, HP3)
HPx, an acronym for HYDRUS [69] and PHREEQC [60],is a versatile coupled reactive transport code for variably
saturated flow conditions. Originally released as a onedimensional flow and transport model about a decade ago[70, 71], it is now also available for multidimensional flow
and transport problems [65, 66]. Additionally, HPx codes
are fully integrated in the well-designed graphical userinterfaces of the HYDRUS codes [51]. Combining the capabilities of the individual HYDRUS and PHREEQC codes,
HPx can be used to simulate many different processes in461

variably saturated porous media: (1) transient water flow, (2)transport of multiple aqueous and gaseous components, (3)heat transport, and (4) a broad range of low-temperature biogeochemical reactions, including interactions of the liquid
phase with minerals, gases, exchangers, and sorption surfaces, described using thermodynamic equilibrium, kinetic,or mixed equilibrium–kinetic formulations.
Variably saturated water flow is described using the
Richards equation, with the possibility of considering either
uniform or dual-porosity (i.e., with water exchange between
mobile and immobile zones) water flow models [69]. As
HPx is specifically developed for soil systems, it includes(1) several single- and dual-porosity models describing
soil hydraulic properties, with an option to consider hysteresis; (2) specific boundary conditions (BCs), such asatmospheric BCs with time-variable precipitation, evaporation and transpiration rates, free drainage, seepage face,and drain BCs; (3) an option to describe macro and small
scale variability using different material types and linear
scaling factors, respectively; and (iv) a compensated root
water and solute uptake model [72], linked to a staticor dynamic root distribution. Heat transport is describedusing an advection–diffusion–dispersion type equation forvariably saturated conditions. Also, solute transport is simulated using the advection–diffusion–dispersion equation,with options to consider uniform, mobile/immobile water,and dual-porosity transport models [69]. As in the caseof water flow, some features are included specifically forsoil systems, such as passive and active uptake (or exudation) of solutes by plant roots. Diffusion in the gas phaseis accounted for by explicitly taking into account variable gas contents. For two- and three-dimensional systems,
anisotropy in the hydraulic conductivity and dispersivity canbe defined.
PHREEQC-2 is currently used as the geochemical solver
within HPx, preserving most functionalities of PHREEQC2 (described in Section 4.1). Within the HPx framework,
the BASIC scripting tool can also be used for a number ofadvanced options. One can define passive and active root
solute uptake, where uptake parameters can depend on solution composition. The BASIC scripting tool can also be used
to define the links between geochemical state variables andphysical transport parameters. Relations between geochemical state variables and physical parameters, such as porosity,
saturated hydraulic conductivity, dispersivity, tortuosity of
the aqueous and gaseous phase (in calculations of effective
diffusion coefficients), heat capacity, heat conductivity, andheat dispersivity, are not predefined but are all adjustable byusers. This means, for example, that porosity–permeability–
diffusion updates are not restricted to the Kozeny–Carmen
approach or Archie’s law, but can consider more complex
relations accounting for the multiscale nature of some typesof porous media [73].

462

Coupling between the transport and geochemical solvers(HYDRUS and PHREEQC-2, respectively) is by a sequential non-iterative approach. Operator-splitting errors canbe controlled using an appropriate spatial and temporal
discretization [71]. The flow and transport equations are
solved numerically using Galerkin-type linear finite element schemes. The fully implicit temporal scheme is used
to solve the flow equation. The water content term is evaluated using the mass conservative method proposed by
Celia et al. [74]. A Crank–Nicholson or an implicit temporal scheme is used to solve the transport equations. The
available spatial weighting schemes include, in addition to
the Galerkin scheme, an upstream weighting scheme and ascheme with artificial dispersion. Additional measures aretaken to improve solution efficiency for transient problems,
including automatic time step adjustment and adherence topreset ranges of the Courant and Peclet numbers. Nevertheless, numerical schemes in HPx are applicable only for
relatively low Peclet numbers, which seems appropriate formost flow and transport problems in soil systems. Timesteps are adapted to ensure convergence and stability whensolving the transport and flow equations.
HPx codes have been used to model a range of different processes. The versatility of HP1 was for example
demonstrated for long-term transient flow and transport ofmajor cations and heavy metals in a soil profile involving cation exchange processes [75], long-term U-migration
in a soil profile following inorganic P-fertilization [76],inverse optimization of cation exchange parameters duringtransient flow [77], fate and transport of mercury in soils[78–80], chemical degradation of cement-based materials
and alterations of transport parameters using a microstructural approach [73, 81], and carbon fluxes in soil systems[82].
4.3 PHT3D
PHT3D is a three-dimensional multicomponent reactive transport code specifically designed for simulating
flow and reactive transport processes in saturated porous
media. PHT3D couples the multispecies transport simulator MT3DMS [83] with the geochemical reaction model
PHREEQC-2 [60] via sequential operator-splitting. The
coupling to MT3DMS provides several options of solution techniques for transport: (1) method of characteristics (MOC), (2) modified method of characteristics(MMOC), (3) hybrid method of characteristics (HMOC),(4) upstream finite difference method, and (5) thirdorder TVD scheme [84]. The solution technique should
be selected by the user based on the nature of thesystem being modeled (such as presence of sharp concentration gradients, redox reactions) and the selecteddiscretization.

Comput Geosci (2015) 19:445–478

Through the flexible, generic nature of PHREEQC-2,
PHT3D can handle a broad range of equilibrium andkinetic reactive processes, including aqueous complexation,mineral precipitation/dissolution, multisite cation exchange,and surface complexation reactions. The PHT3D code also
has the capability to simulate kinetically reacting immobile
components. Typical reactive processes such as microbialactivity, NAPL dissolution, and isotopic fractionation canalso be simulated. The structure of PHT3D code is suchthat these and other processes can be included througha PHREEQC-2 database files without modification of the
PHT3D source code.
PHT3D was developed in the 1990s and has been
applied to diverse groundwater problems. In the original
version, transport was calculated during a time step, an input
file was written for PHREEQC for geochemical calculations, and these steps were repeated for subsequent time
steps until finished [85]. A second version was released
in 2010 that tightened this reaction-transport coupling in
the code to reduce computation time. In this version,
PHREEQC is used for storing the chemical data of themodel, including the chemical activities and the composition of surface complexes from the previous time step.
A review of the developments in PHT3D version 2 andits potential applications were summarized by Appelo andRolle [86].
Previous modeling studies using PHT3D targeted natural and enhanced remediation of organic (e.g., BTEX,
CHCs) and inorganic (e.g., metals) contaminants and thegeochemical changes occurring during their degradation.
Applications of PHT3D also included modeling the fate
and transport of pesticides, reactive transport of arsenic anduranium, trace metal remediation by in situ reactive zones,
temperature-dependent pyrite oxidation, and isotopic fractionation during natural attenuation of organic pollutants[84]. Recently, PHT3D has been used to model interactionsof seasonal aquifer thermal energy storage and groundwater contamination plumes [87] and radially symmetric flow
and solute transport around point sources and sinks [88].
Recent enhancements of the code include implementation ofelectrokinetic transport to deliver a remediation compound
within heterogeneous and low-permeability sediments forin situ chemical oxidation [89]. Further information regarding installation, user documentation, and benchmark test
problems can be found at http://www.pht3d.org/.
PHT3D has been integrated into the PMWIN preor postprocessing package so as to enhance its graphical user interface (GUI) and visualization capabilities. This PMWIN coupling provides users the flexibility to select from existing reaction packages/databasefiles or develop and integrate newly developed packages [90]. Another GUI available to PHT3D users is
Visual MODFLOW.

Comput Geosci (2015) 19:445–478

4.4 OpenGeoSys
OpenGeoSys (OGS) is a scientific open-source
initiative for numerical simulation of thermo–hydro–
mechanical/chemical processes in porous and fractured
media. The basic concept is to provide a flexible
numerical framework for solving coupled multifieldproblems. The OGS code is targeting primarily applications in environmental geoscience, e.g., in the fields
of contaminant hydrology, water resources management, waste deposits, or geothermal systems, but it
has also been applied to new topics in energy storage[91].
The software is written with an object-oriented finite
element method concept, including a broad spectrum ofinterfaces for pre- and postprocessing [92]. The OGS developer environment provides an open platform that includes
professional software-engineering tools such as platformindependent compiling and automated result testing tools.
To ensure the quality of the code and facilitate communications among different developers, more than 150 benchmarks are recently validated before any changes in thesource code are released. In addition, the OGS community has been active in several international benchmarking
initiatives, including e.g. DEVOVALEX (with applications
mainly in waste repositories, Rutqvist [93], CO2BENCH
and Sim-SEQ (in CO2 sequestration, [94, 95]). The bookcontaining over 100 wellverified benchmarks has been published recently [96].
In OGS, heterogeneous or porous-fractured media canbe handled by dual continua or discrete approaches,i.e., by coupling elements of different dimensions. OGS
has a built-in random-walk particle tracking method
for Euler-Lagrange simulations [97] and can simulate
incompressible and compressible flow [98], densitydependent flow [99], unsaturated flow [100], two-phaseflow [101, 102], and overland flow [103]. All flow
model components are available for non-isothermal conditions and for deformable porous media. Equations ofstate (EoS) for several fluids (gases, supercritical fluids) are available to module temperature-driven dynamics[104]. Multicomponent mass transport and (equilibrium andnonequilibrium) heat transport [105] modules are part of
OGS.
The OGS Tool-Box [106] promotes functionality extension through code coupling interfaces for example coupling external geochemical and biogeochemical simulation
tools such as PHREEQC/IPhreeqc, GEMS, and BRNS,
in addition to the underlying internal process couplings.
This flexibility has made it possible for OGS to simulate
geochemical processes such as ion exchange, mineral dissolution/precipitation, and equilibrium reactions in partly saturated bentonites (OGS#PHREEQC-2 [107]), formation of463

solid solutions and radionuclide retardation (OGS#GEMS[108]), and thermochemistry for technical energy storage[91], as well as microbiological reactions such as
Michaelis–Menten degradation kinetics to estimating in situ
degradation rates [109] and organic carbon degradation via
multiple microbial degradation pathways, dispersive mixing controlled bio-reactive transport with different reaction
kinetics (OGS#BRNS [110]).
The OpenGeoSys community has been in development
since the mid-1980s and the numerical simulator has been
continuously re-engineered through Fortran, C, and C++
implementations and in the meantime the code has been
optimized for massively parallel machines, e.g., JUQUEEN.
The OGS community also developed continuous workflowsthat include various interfaces for pre- and postprocessing.
Additionally, visual data integration became an important
tool for establishing and validating data driven models[111] since three-dimensional visualization facilities, suchas the UFZ-VISLAB (Leipzig, Germany), can be used forcomprehensive systems analysis, computational steering,and demonstration purposes for stakeholder and authorities. Detailed information regarding code, benchmarks, andcommunity is available through: opengeosys.org
4.5 HYTEC
HYTEC is a reactive transport code that integrates a wide
variety of features and options that have evolved, after
more than a decade of development, to a widely usedand versatile simulation tool [112]. Solution capabilities
for biogeochemistry are provided by the code CHESS(http://chess.ensmp.fr). The model accounts for many commonly encountered processes, including interface reactions(surface complexation with electrostatic correction andcation exchange), precipitation and dissolution of solid
phases (minerals and colloids), organic complexation, andredox and microbial reactions. All reactions can be modeled using a full equilibrium, a full kinetic, or a mixed
equilibrium–kinetic approach. All reactions can be modeled using a full equilibrium, a full kinetic, or a mixed
equilibrium–kinetic approach. Thermodynamic data aretaken from international databases, e.g., EQ 3/6 [113].
The hydrodynamic module of HYTEC is adapted forhydrodynamic conditions commonly encountered in the laboratory or in the field. The code allows for dual porosity,
saturated, unsaturated or two-phase flow, variable boundary
conditions, sinks, and sources [114]. HYTEC searches foran accurate solution to the multicomponent transport problem using an iterative, sequential, so-called strong coupling
scheme. Strong coupling permits variable hydrodynamic
parameters as a function of the local chemistry [115].
For example, the porosity of a porous medium decreases
after massive precipitation of newly formed mineral phases,

464

which modifies the water flow paths and transport parameters, e.g., diffusion coefficients. HYTEC solves thisinterdependency accurately, which makes the tool particularly useful for, e.g., cement alteration at long timescales(e.g., waste storage and performance assessment).
HYTEC has been used at various scales from micrometer to kilometer and within various geosciences fields:
i.e., cement degradation, cement/clay interaction, polluted
soils, waste dumps, radioactive waste disposal (performanceassessment, near- and far-field processes [116, 117]), geological storage of acid gases [118, 119], and copper anduranium in situ recovery processes. Other applications concern the evolution and degradation of (geo)materials suchas ashes, concrete, and cements; the latter often being simulated by a typical CEM-I cement, but more sophisticatedmodels for cements can be used including sorption on
primary or secondary calcium silicate hydrate phases, carbonation, and sulfate attack on the material [120, 121]. The
strong coupling approach as outlined above makes HYTEC
particularly useful for the modeling of long-term leachingof solidified wastes [122–124].
Efforts to develop, test, and validate the HYTEC model
largely exceed the scope of a single laboratory. The Reactive Transport Consortium (Pôle Géochimie-Transport,
http://pgt.geosciences.ensmp.fr) provides a national
research framework with the objective of creating a longterm framework for the development of reactive transportmodels, reference studies, and new application domains.
Already operational for several years, the collaborative
efforts within the PGT have made it possible to makeconsiderable progress in the domain of reactive transportmodeling. Recent improvements include new solid solution, gas, and isotope modules. HYTEC is now able todeal with complex gas mixtures and to model accurately
their interactions with water–salt–rock systems, and thisability has been tested against laboratory batch experiments [125]. The new isotope module also allows for the
management of complex radioactive decay chains andchemical isotopic fractionation within the core of the code,
allowing strong coupling between chemical and transportmodules.
4.6 ORCHESTRAORCHESTRA is a platform developed by Meeussen [126]for modeling equilibrium chemistry, with the option ofincluding kinetics and/or transport processes. ORCHESTRA does not enforce a specific model, but is rather ageneric toolbox that allows users to create their own models
using an object representation, which is stored in an object
database and is fully user-definable.
In ORCHESTRA, reactive transport processes are implemented by a mixing-cell concept. This implies that transport

Comput Geosci (2015) 19:445–478

systems can be composed of (well-mixed) cells and connections between these cells. The cells contain the informationon the local physical and chemical composition, while theconnections between the cells contain the mass transport
equations (diffusion, convection, etc.). Both the connections
between the cells, but also the literal equations involvedin mass transport, are defined in text input files. Thus,
ORCHESTRA provides the flexibility to the users to adapt
the model structure and equations.
Effectively, the mixing-cell concept results in a finite difference scheme. In ORCHESTRA, it is possible to choosea sequential iteration (SIA) approach (predictor corrector
method) or a non-iterative approach (SNIA) solved with afourth order Runge–Kutta method. With SNIA, the kinetic
reactions are solved using the same time step as the transport
processes.
ORCHESTRA is coded in Java, which results in an executable code that runs on different operating systems (e.g.,
Windows, Linux, Unix, Solaris, OSX). Java also contains
standard language constructs for parallel programming,
which facilitates parallelization and makes efficient use
of multiprocessor hardware. In contrast with the IPhreeqc
method for parallelization that uses domain decomposition with nodes/cells predefined for each processor, the
ORCHESTRA approach uses a single node pool, fromwhich nodes are taken by each available processor. The
advantage of the IPhreeqc method is that it is more suitablefor distribution work over different physical machines. The
advantage of the ORCHESTRA method is better load balancing for systems with strongly varying calculation times
per node.
In ORCHESTRA, it is possible to obtain real-timegraphical output for any chemical/physical parameter (e.g., concentration or pH profiles) during acalculation.
Historically, ORCHESTRA was first developed at the
Macaulay Institute in Aberdeen, UK. Earliest applicationsof the code targeted chromium remediation and releasefrom chromite ore-processing residue [127, 128]. Modeldevelopment and verification then shifted to Wageningen
University and Research Center, the Netherlands. Here,
ORCHESTRA was used to model interactions between
humic and fulvic acids, ions, and mineral surfaces [129–131]. These studies led to the implementation of CDMUSIC, NICA-Donnan, Donnan ion exchange, and Ligand
and Charge Distribution (LCD) model in ORCHESTRA for
multi-surface adsorption modeling of cations and anions
on organic matter, clay, and Fe oxides and the determination/interpretation of model parameters. ORCHESTRA
has also been coupled to BIOCHEM for evaluating chemical speciation and eco-toxicological impacts of heavy metal
uptake by biota [132]. ORCHESTRA has capabilities tomodel both saturated and unsaturated zones and has beenComput Geosci (2015) 19:445–478

extensively applied to model data from leaching tests, landfill leachate, pilot studies as well as groundwater data. For
example, it has been used to implement a simple rhizospheremodel by simulating diffusion and solute uptake towards asingle cylindrical root [133]. Since 2010, code development
has continued at the Nuclear Research and Consultancy
Group. Recent applications of the model include couplingof ORCHESTRA with LeachXS database to improve understanding and prediction of the long-term performance of
cementitious barriers used in nuclear applications. This
thermo–hydro–mechanical (THM) coupling has also beenused to assess long-term release behavior of sulfidic mining wastes [134], carbonation-induced corrosion in highlevel waste tanks [135], and uncertainty quantification ofthermodynamic parameters [136]. More information oncode development and example applications is available at
http://www.macaulay.ac.uk/ORCHESTRA/.
4.7 TOUGHREACTTOUGHREACT [137–139] was developed by introducing reactive geochemistry into the framework of TOUGH2
V2 [140, 141]. TOUGH2 is a general purpose multiphaseflow simulator that relies on a suite of equation of state(EOS) modules to compute fluid phase partitioning, andflow/transport of various components (e.g., water, CO2 ,
salt, air, tracers, radionuclides) in liquid, gas, and nonaqueous phases. Fluid flow by advection is computed using themultiphase extension of Darcy’s Law, and diffusive masstransport is considered in all phases. The flow of heat byconduction and convection is also computed, taking into
account temperature effects from phase changes. The evolution of each phase, including appearance and disappearance,
and the partitioning of components between the different
phases, are computed dynamically for each model grid
block assuming local thermodynamic equilibrium for given
conditions of temperature, pressure, and composition. The
mass and energy balance equations are solved implicitly by
Newton–Raphson iterations. An unstructured finite-volume(integrated finite-difference) numerical approach is implemented for the flow/transport computations. This allows
modeling of porous and/or fractured systems and other typesof systems with zones of different component mobility,
through the use of multiple interacting continua (MINC)[44, 142].
By building on TOUGH2, TOUGHREACT simulates
non-isothermal, multicomponent reactive transport of aqueous and gaseous components in variably saturated media.
Reactive transport is solved by an operator-splitting
approach that can be either iterative or non-iterative. Reactive processes considered include aqueous and surface complexation, ion exchange, mineral precipitation/dissolution,and gas exsolution/dissolution. Reactions among aqueous

465

species and between aqueous species and minerals proceed under equilibrium and/or kinetic constraints. Other
reactions operate under local thermodynamic equilibrium.
The precipitation and dissolution of minerals is optionallycoupled to porosity, permeability, and capillary pressureusing various correlations that feed back into the multiphaseflow computations. Aqueous and gaseous components aretransported via advection and/or diffusion, following the
calculation of multiphase fluxes. The transport step is followed with a fully implicit solve of the combined kinetic
and equilibrium geochemistry.
TOUGHREACT has been applied to simulate a wide
range of subsurface hydrological and biogeochemical
environments. Applications have included, for example,
geothermal systems [143–145], nuclear waste repositories[146–149], geologic carbon sequestration [150–154], andenvironmental remediation [155–157]. The development of
TOUGHREACT has been ongoing since the late 1990s. Version V2 [137] was used in the benchmarks described in thisspecial issue. More recent developments of TOUGHREACT include the coupling of geomechanical effects [158]and hybrid MPI-OpenMP parallelization of the reactive
transport routines.
4.8 eSTOMPeSTOMP is the parallel processing version of the Subsurface Transport Over Multiple Phases (STOMP) simulator[159]. The eSTOMP simulator can address variably saturated subsurface flow, transport, and reactions in heterogeneous porous media on the most powerful computers. The
scalability of the formulation has been demonstrated on up
to 131,000 processor cores for large complex biogeochemical reaction networks in highly resolved physically andgeochemically heterogeneous sediments driven by dynamic
variably saturated flow conditions.
Transient or steady-state flow can be confined/unconfined and/or variably saturated. Spatiallyvariable material properties include anisotropic permeability/hydraulic conductivity (which can besaturation-dependent); diffusive porosity, longitudinal andtransverse dispersivity, saturation function parameters(van Genuchten, Brooks-Corey, tabular data; Mualem,
Burdine, Haverkamp), tortuosity, specific storativity, andcompressibility.
Transport algorithms address advection, anisotropic dispersion, and interphase partitioning. High Peclet number
total variation diminishing (TVD) transport schemes minimize numerical dispersion.
Arbitrarily complex reaction and rate law forms canbe addressed in the biogeochemical reaction solver,
as well as standard equilibrium and kinetic reactions
for aqueous complexation, surface complexation, ion

466

exchange, precipitation–dissolution reactions, and aqueous–gas exchange. Microbially mediated reactions include
species-specific biomass production and decay, Monod-type
rate laws with biomass dependencies, and thermodynamic
controls. The modular reaction design currently provides
the user with a range of solver capabilities including Pitzer,
Debye–Hückel, and Davies activity models through tworeaction solvers, ECKEChem [160] and BIOGEOCHEM[161].
A unique feature of the eSTOMP simulation capability is the option to couple highly detailed, mechanistic,
genome-scale metabolic in silico models into the reaction module. In the constraint-based flux balance approach,hundreds of intracellular and environmental exchange reactions fundamentally account for the cell-scale physiology
controlling cell maintenance, biosynthesis, and respiration.
Under a maximization principle (e.g., growth), influxesof nutrients, electron donors, and acceptors drive a system of reactions with associated fluxes that can be determined using linear programming solvers. In this formulation, the biologically mediated reaction stoichiometry
and rates are a function of the metabolic status of theorganism, providing a more accurate depiction of theseprocesses.
The eSTOMP simulator solves the partial-differential
equations that describe the conservation of mass or energy
quantities by employing integrated-volume finite-difference
discretization to the physical domain and backward
Euler discretization to the time domain. An operator splitting solution method is used for flow, transport, andreactions, each of which employs implicit time-stepping
schemes. The resulting equations are nonlinear coupled
algebraic equations, which are solved using Newton–
Raphson iteration. Flow and solute transport are solved
using parallel iterative conjugate gradient techniques as
implemented in the Portable Extensible Toolkit for Scientific Computation (PETSc) libraries [162]. Differential algebraic equation solver techniques are used forthe standard biogeochemical reaction network and the
IBM ILOG CPLEX optimization solver package is used
for the constraint-based system of metabolic pathway reactions. Input is directed through semi-formatted text files
and output is available through a variety of user-directed
formats. The simulator recognizes a number of boundary condition types and allows their specification bothinternally and externally to the computational domain. Gridding options include structured Cartesian and boundary
fitted.
The memory and computational requirements for simulating coupled variably saturated and multicomponent
biogeochemical reactive transport in highly resolved threedimensional porous media typically exceeds the capabilitiesof standard desktop computing. Accordingly, eSTOMP was

Comput Geosci (2015) 19:445–478

designed for execution on the full range of parallel processing computer architectures. eSTOMP uses one-sided
communication and a global shared memory programming
paradigm from the Global Array Toolkit (GA) library [163]for scalability, performance, and extensibility on massively
parallel processing computers. The approach is compatible
with the more commonly used Message Passing Interface(MPI) that is used by the PETSc global implicit solver in
eSTOMP [162]. The principal benefit to this approach isthat simulations with highly resolved process and property
detail can be performed on a massively parallel computerin hours instead of weeks on a desktop computer. These
efficiencies enable comprehensive exploration of alternative
conceptual process models, hypothesis testing, and uncertainty quantification that represent a much larger expenditure of computational effort than the final production
runs.
The theoretical and numerical approaches applied in thesimulator have been documented in published theory guides[159, 160]. The simulator has undergone a rigorous testingprocedure (http://stomp.pnnl.gov/QA/stomp-w-r.stm) andis maintained under version control procedures that meet
NQA-1-2000 software requirements, as well as the requirements specified under DOE Order 414.1C for Safety
Software.
4.9 HYDROGEOCHEMHYDROGEOCHEM has evolved from a reactive chemical transport model developed by Yeh and co-workers[164] into a generic model that couples thermal (T),hydrology (H), geomechanics (M), and chemical (C)(THMC) processes (HC) [165], (THC) [166, 167]; (THMC)[168, 169]). It includes many numerical options toenhance its robustness and efficiency. The latest versionof HYDROGEOCHEM involves four fully coupled modules: fluid flow, thermal transport, reactive biogeochemicaltransport, and geomechanic displacement and deformation[168, 169].
The pressure-based multiphase flow equations of arbitrary L fluid phases are transformed into a set of a totalpressure equation and (L-1) equations governing the evolutions of (L-1) accumulated degree of saturations. This setof equations are explicitly coupled via the storage coefficients and diffusion coefficients of all primary variables inevery equation [151], which make the system much less stiff
for numerical approximations. The transport equations governing the evolution of M species are transformed into a
set of M equations governing the transport of M reaction
extents. This effectively decouples the transport of components and kinetic variables from those for the rates of
equilibrium reactions. As a result, it facilitates the iterative approach between a much-reduced number of partial

Comput Geosci (2015) 19:445–478

differential equations of component and kinetic-variabletransport and nonlinear algebraic equations that implicitly define the rates of equilibrium reactions [168, 170,171]. In the geomechanics module, the media are treatedas consisting of nonlinear viscous-elastic materials. Thisconsideration enables the model to deal with large deformation and prevents a complete recovery after the stresses
are released, which is more realistic than an elastic material
treatment.
Similar types of governing equations to those describedin Section 3 were implemented in HYDROGEOCHEM up
to version 2. However, a different approach has been takenin the later versions. This alleviates the necessity of writing
reactions in the canonical form for many types of reactions.
Starting with a reaction network of N reactions involving M
species, each reaction can involve any number of species,
one can write M species transport equations based on theprinciple of chemical kinetics [172]. Performing Gauss–
Jordan reduction of the reaction network [173–175], one
obtains a system of M transport equations for M reaction extents, each of which is a linear combination of the
M species concentrations [170, 171]. A reaction extent iscalled a component when there is no reaction term appearing in its transport equation. A reaction extent is called akinetic-variable when there is one and only one independent kinetic reaction appearing in its transport equation.
Similarly a reaction extent is called an equilibrium variable
when there is one and only one independent equilibrium
reaction appearing in its transport equation. This approacheffectively decouples the computation of all N reaction
rates from the simulation of a much reduced set of transport equations of components and kinetic variables as statedin the previous paragraph. It also allows the flexibility of
programming any type of reaction rates into the code.
HYDROGEOCHEM can include a wide range ofbiogeochemical reactions such as aqueous complexation, adsorption/desorption, ion-exchange, precipitation/
dissolution, redox, acid-base reactions, and microbialmediated reactions. Multi-adsorbing sites and/or multi-ion
exchange sites options are available for sorption reactions.
Twenty-nine types of kinetic rate equations and four typesof equilibrium rate equations have been implemented in
the latest version of HYDROGEOCHEM [176, 177]. Other
types of users’ defined kinetic reactions and equilibrium
reactions (e.g., the hydration reaction which evolves thesite rather than conserves it) can be accommodated. Fourschemes are available for coupling hydrologic transportand biogeochemical reactions: (1) fully implicit iteration approach [178], (2) operator-splitting approach, (3)predictor–corrector approach, and (4) mixed predictor–
corrector and operator-splitting approach [169]. To enable
the model to deal with complex media geometry, unstructured finite element meshes are used. For two-dimensional

467

problems, both triangular elements and quadrilateral elements or the mixture of two types of elements can be
employed to approximate complicated global geometry as
well as irregular media interfaces. For three-dimensional
problems, tetrahedral, hexahedral, triangular prism elements, or the mixture of these three types of elements can beselected. To make the model robust and efficient, Galerkinfinite element is the choice for numerical approximation of
the governing equations in multiphase phase flow module.
For scalar transport modules of thermal and reactive transport, the Galerkin finite element, Petrov–Galerkin finite
element, and/or Lagrangian–Eulerian finite element methods can be used to approximate the governing equations.
For the geomechanic module, the linear successive approximation algorithm [179] is employed to approximate thegoverning equation of displacement.
HYDROGEOCHEM has been applied to many problems for research investigations as well as for practicalsimulations of radioactive waste disposal and redox zone
migration. For example, it was employed to assess the adequacy of Kd-based models [180], used to investigate theperformance of waste trench designs [181], and applied to
remediation of contaminated lands [182]. It was applied
to a sandy site to assess the effect of redox zones on
the migration of metals [183], adapted to investigate themulti-segment reactions in contaminated aquifers [184],
and employed to assess the redox geochemistry at shallow
groundwaters [185]. The design capacity of fully coupled
THMC processes is demonstrated in Yeh et al. [170].
4.10 CrunchFlow
CrunchFlow is a software package for simulating reactive
transport developed by Steefel and co-workers and applied
since 1988 to a variety of problems in the earth and environmental sciences. CrunchFlow is based on a finite volume
discretization of the governing coupled partial differential
equations that link flow, solute transport, and multicomponent equilibrium and kinetic reactions in porous and/or
fluid media. Two approaches for coupling reactions andtransport are available at runtime, (1) a global implicit
approach that solves transport and reactions simultaneously(up to two dimensional), and (2) a time or operator splitting approach based on SNIA (up to three dimensional).
The global implicit approach is based on a backwards Euler
time discretization, with a global solution of the coupled
reactive transport equations using Newton’s method. The
global implicit approach in CrunchFlow is implemented
with the direct substitution approach (DSA) (72). An individual Newton iteration is solved using iterative conjugate
gradient techniques as implemented in the PETSc libraries[162]. The global implicit approach in CrunchFlow can be
run at arbitrarily large Courant numbers since there is no468

associated operator-splitting error, although other forms oferror may still occur. Currently, the global implicit approach
uses upstream weighting to simulate advective transport.
The operator splitting approach, in contrast, requires time
steps that are everywhere less than or equal to the Courant(for advection) and von Neumann (for diffusion) criteria.
Advective transport is solved with a third-order accurate
TVD approach that minimizes numerical dispersion while
maintaining monotonicity. The use of the FORTRAN 90
language allows for dynamic allocation of memory forarrays, making it possible to choose the reaction-transport
coupling scheme at runtime. Fully saturated flow is coupled
sequentially to the reactive transport calculations—partiallysaturated flow and transport (including a gas phase) canbe included, but currently requires calculation of the flow
and liquid saturation field with separate software. Feedbacks on transport (flow and diffusion) rates as a result ofreaction-induced porosity and permeability changes can beconsidered as well.
Diffusion can be modeled with either Fick’s law or with
the Nernst–Planck equation (currently without accounting
for gradients in activity coefficients), making it possible toinclude differing diffusion coefficients for charged species
while maintaining electroneutrality (48). CrunchFlow canalso simulate accumulation and transport of ions within adiscrete electrical double layer (EDL), with dynamic balancing of the surface charge on the mineral and in the
Stern layer calculated with a surface complexation model.
Dynamic treatment of the EDL thickness as a function ofionic strength can be used to simulate the changes in bulk
versus EDL porosity. The treatment is based on a Donnan equilibrium or mean electrostatic approach to the EDL[186]. In the absence of charged mineral surfaces, thisoption becomes a more conventional dual porosity treatmentof transport. Gas phase diffusion is handled by explicitly accounting for the composition of the gas phase andassociated gradients.
The CrunchFlow software makes use of the canonicalapproach to reaction networks in which a set of equilibrium secondary or non-component aqueous complexes aredefined in terms of primary or basis species that constitutethe actual unknowns in the system, thus reducing substantially the number of governing equations and unknowns
that need to be solved for once the secondary species are
eliminated algebraically. Aqueous complexes can also betreated as kinetically controlled, in which case they areretained as primary species, while mineral reactions aretreated in every case as kinetically controlled, although
local equilibrium behavior can be simulated by choosing suitably large reaction rate constants. Multisite surface
complexation reactions, with or without a formal treatment of electrostatic corrections, multisite ion exchange

Comput Geosci (2015) 19:445–478

reactions (including Gaines–Thomas, Vanselow, and Gapon
conventions), and gas–aqueous phase reactions are currently
treated as equilibrium processes. Using an automatic readof a thermodynamic and kinetic database, the code canbe used for reactive transport problems of arbitrary complexity and size (i.e., there is no a priori restriction on
the number of species or reactions considered). Microbiallymediated processes using classical Monod and inhibition
terms, with or without an explicit thermodynamic term [26],can be simulated with a treatment of dynamic microbial
populations.
The first use of predecessors of CrunchFlow was insimulations that took into account nucleation, precursors,
Ostwald ripening, and crystal growth [187]. Steefel and
Lasaga [188] made use of an early version of CrunchFlow
to investigate non-isothermal, multidimensional numerical
reactive transport model for hydrothermal systems, withkinetic treatment of mineral–water reactions. In a seriesof papers, Steefel and Lichtner [189–191] investigateddiffusion-controlled reactive transport bordering active fractures, with consideration of porosity and reactivity feedbacks as a result of geochemical reactions. Making useof new capabilities for modeling the Nernst–Planck equation (electrochemical migration), Giambalvo et al. [192]
investigated elemental fluxes from an off-axis hydrothermal system. Multisite, multicomponent simulation of ion
exchange processes affecting the transport of 137 Cs below
leaking tanks at the Hanford contaminated site was presented [22]. CrunchFlow was used as the basis of a diffusion and reaction + sediment burial model to explore
uranium isotopic exchange and alpha recoil effects andincluded a rigorous treatment of steady-state burial andcompaction [193]. In 2008, CrunchFlow was applied toreactive transport within a single pore, with a discussion ofscale effects on geochemical reaction rates [1]. In Maher et
al. [194], CrunchFlow was also used in a reactive transport
analysis of weathering rates in the Santa Cruz chronosequence where it was shown that laboratory and field rates
could be reconciled when nonlinear reaction kinetics were
included. Navarre-Sitchler et al. [195] used CrunchFlow
to simulate weathering rind formation based on the diffusion threshold model published earlier [196], which madeuse of CrunchFlow for three-dimensional diffusion simulations of the microtomography resolved pore network.
A series of studies making use of the CrunchFlow software focused on microbially mediated biogeochemical
reactions (Fe(III), U(VI), and sulfate reduction) associated with biostimulation at the Old Rifle site [197–201].
More recently, the capability for modeling of kinetic andequilibrium isotopic equilibration was added to the code.
More information on the code can be found at http://www.
csteefel.com/CrunchFlowIntroduction.html.

Comput Geosci (2015) 19:445–478

4.11 MIN3P
MIN3P is a multicomponent reactive transport code,
specifically designed for simulating flow and reactive
transport processes in variably saturated media. Unsaturated/saturated subsurface flows are described by Richards
equation and the code uses the global implicit approach
implemented using the direct substitution approach (DSA)
for solution of the multicomponent advection–dispersion
equation and biogeochemical reactions (72). Spatial discretization is performed based on the finite volume method,
facilitating simulations in one, two, and three spatial
dimensions. Advective transport terms can be describedby upstream weighting, centered spatial weighting, orusing a flux limiter technique to minimize numericaldispersion. Implicit time weighting is employed, whichallows to use large time steps for problems that arestrongly affected by water–rock interaction, without lossof accuracy. The highly nonlinear and coupled reactive transport equations are linearized using Newton’s
method.
The original MIN3P code developed by Mayer andco-workers [202] includes a generalized framework for
kinetically controlled reactions, which can be specified
through a database along with equilibrium processes.
The general kinetic formulation includes intra-aqueousand dissolution–precipitation reactions in addition to geochemical equilibrium expressions for hydrolysis, aqueous
complexation, oxidation–reduction, ion exchange, surface
complexation, and gas dissolution–exsolution reactions.
The generalized approach allows consideration of fractional order terms with respect to any dissolved species
and mineral phases, which facilitates the incorporationof a variety of experimentally derived rate expressions.
Monod and inhibition terms can be used to describe microbially mediated reactions. Dissolution–precipitation reactions can be described as surface-controlled or transportcontrolled. The formulation also facilitates the consideration of any number of parallel reaction pathways,
and reactions can be treated as irreversible or reversible
processes.
Several code developments have been performed over the
past decade to further improve the capabilities of MIN3P.
Specifically, MIN3P has been enhanced to include a dual
porosity formulation (MIN3P-DUAL [203]), allowing forreactive transport simulations involving mobile and immobile regions. This code was used in the assessment of thefate and transport of MTBE in a Chalk aquifer [204]. To
describe gas entrapment and release due to water table fluctuations, as well as formation and migration of a discontinuous gas phase, the MIN3P-BUBBLE model was developed[205, 206]. Simulations conducted with this version of the469

code were related to problems involving in situ gas production or consumption at hydrocarbon spill sites under naturalattenuation conditions. Investigations focused on oxygenand carbon transfer across the capillary fringe and the effectof these processes on contaminant degradation and groundwater plume migration. In addition, MIN3P was extended
by implementing the dusty gas model (DGM), which
describes multicomponent gas transport by advection anddiffusion in the vadose zone (MIN3P-DUSTY, [36]). The
code was used to simulate gas attenuation in landfill coversoils [207], methane production and oxidation in aquifers
contaminated by organic compounds (e.g. an oil spill site),and sulfide oxidation in mine tailings [36, 208]. MIN3P
was also enhanced by implementing the density coupling between flow and reactive transport (MIN3P-DENS,[209]) and was used to simulate permanganate-based
remediation of chlorinated solvents under free convection
conditions.
Recently, MIN3P-DENS was designated as the basis
for all future MIN3P developments. A thermo–hydro–
chemical (THC) formulation that extends the applicabilityof the code to reactive transport problems involving brines
was implemented. The revised code (MIN3P-THCm) also
accounts for the hydro-mechanical effects caused by surface
loading [210]. Specifically, code enhancements included
the implementation of the Harvie–Moller–Weare (HMW)model [211] based on Pitzer’s equations [212]; fluid density
calculations based on Pitzer’s equations [213]; a simplified hydro-mechanical coupling formulation to account for
the vertical compaction of the porous medium as presented
by Neuzil [214]; and the energy balance in the porous
medium. The MIN3P-THCm code was further enhanced byincluding water vapor transport and atmospheric boundary
conditions [215]. The resulting code was used to illustrate the long-term stability of a sedimentary basin during
a glaciation/deglaciation cycle [210] and to simulate carbon
sequestration in ultramafic mine tailings [215]. MIN3PTHCm is the code used to contribute to the solution of thebenchmarks presented in this special issue, unless otherwise
noted.
In addition to the capabilities outlined above,
MIN3PTHCm includes an option to account for permeability and tortuosity changes as a function of porosity
modifications due to mineral dissolution–precipitation
reactions. Current work involves the implementation ofmulticomponent diffusion and electrochemical migration
based on the Nernst–Planck equation, and the parallelization of the code for both shared and distributed memory
platforms.
Many applications of the MIN3P code family have utilized the capabilities to effectively simulate gas transport in
the unsaturated zone, coupled with reactive solute transport.

470

For example, the code has been used by a number of
researchers to simulate the generation and attenuation ofacid rock drainage in mine waste (e.g., [202]), sourcezone natural attenuation of hydrocarbon spills in variably
saturated media (e.g., [208, 216]), land application of foodprocessing waste water [217], and biogeochemical cyclingof silicon in forest soils [218].
The MIN3P-THCm code is maintained under fullversion control and is verified after each revision for a large
suite of verification examples. More information on thecode can be found at http://www.eos.ubc.ca/research/hydro/
research/min3p/reactive tran web.htm.
4.12 PFLOTRANPFLOTRAN is an open-source, massively parallel multiscale and multiphysics code for subsurface flow andtransport applications by Lichtner and co-workers [219]. A
system of nonlinear partial differential equations describing
non-isothermal multiphase flow, biogeochemical transportand geomechanics in porous media are solved. The code
is designed to run on machines ranging from laptops andworkstations to petascale leadership class computing architectures. Parallelization is achieved through domain decomposition using the Portable Extensible Toolkit for Scientific
Computation (PETSc) libraries for its parallel framework[162]. PFLOTRAN has been developed from the ground up
for parallel scalability and has been run on up to 218 processor cores with problem sizes of over 2 billion degrees of
freedom.
The code is written in object-oriented Fortran 2003.
The subsurface flow processes available in PFLOTRAN
include single-phase variably saturated flow (Richards
equation), non-isothermal two-phase water-supercritical
CO2 and air-water systems, thermal–hydrologic coupledheat and mass conservation, and three-phase ice–liquid–
vapor (water) flow. The flow modes can be coupled
to multicomponent reactive transport equations using
either a fully implicit Newto–Raphson algorithm or operator
splitting. The multicomponent reactive transport processes
include aqueous complexation, sorption (isotherm-based,
ion exchange, surface complexation), mineral precipitation/dissolution, and microbially mediated biodegradation.
In addition, the mean age of a tracer and water canbe calculated using Good’s formulation [220]. All theprocesses can be solved on either structured or unstructured grids. To model non-isothermal flow and transport
in fractured porous media, either a discrete fracture
network approach can be used, or a multiple interacting
continua formulation. A computationally efficient algorithm is used to rigorously decouple primary and secondary
continuum equations for the multiple interacting
continua method. Additional processes that can be

Comput Geosci (2015) 19:445–478

solved include surface flow (under development)
and geomechanics for elastic deformations. The surface flow mode can be coupled to the subsurface
Richards mode, while the geomechanics mode can besequentially coupled to any of the subsurface flowmodes.
Subsurface and surface flow equations are discretizedusing the finite volume method while the geomechanics
mode uses the finite element method. A novel feature of thecode is its ability to run multiple input files and multiple
realizations with heterogeneous permeability and porosity fields simultaneously on one or more processor cores
per run. This feature is especially powerful when conducting sensitivity analyses and quantifying model uncertainties
enabling computation of hundreds of realizations in theamount of time needed for a single realization. Since I/O isthe bottleneck for large simulations especially with multiple
realizations, PFLOTRAN uses the HDF5 API to read andwrite data in parallel. Obtaining better performance using
parallel I/O is currently an active area of investigation.
PFLOTRAN is a multiple DOE laboratory effort withcore developers from LANL, SNL, LBNL, and ORNL.
PFLOTRAN is currently being used for several subsurface
applications related to energy, climate, and national security. Some of these include nuclear waste disposal, CO2
sequestration, enhanced geothermal systems, groundwatercontamination, multiscale processes, and Arctic hydrology[221–225]. PFLOTRAN is released under GNU Lesser
General Public License. PFLOTRAN has users from allover the world including research laboratories, industries,
and universities. A regression suite of test problems is runnightly and with each new push of changes to the bitbucket repository for QA purposes. Information regarding
installation and user documentation can be found at http://
www.pflotran.org/, and can be freely downloaded fromhttps://bitbucket.org/pflotran/pflotran-dev.

5 Summary and path forward
Hopefully it is clear from the foregoing discussion that
subsurface reactive transport simulators have evolved dramatically in the last 10–15 years. It is now possible tosimulate a wide range of complex, coupled subsurface processes with a number of different simulators, and withcomparable results as demonstrated by the benchmark problems described in this volume. A range of coupling options
have emerged as different groups have sought modularity
of processes and applications in some cases, tightly coupledsolution procedures in others. In essentially every case, there
are tradeoffs between the various approaches and the choiceof one or the other should likely depend on the specificsubsurface environmental problem that needs to be solved.

Comput Geosci (2015) 19:445–478

We would argue, therefore, that the progress to date has
been impressive, to the point where reactive transport modeling might even be recognized in the future as a distinct discipline within the earth and environmental sciences, rather
than a mongrel activity wedged between hydrology and geochemistry. Nonetheless, there is much that remains to bedone in this field. Solid solutions need to be incorporated
into reactive transport models in a more extensive manner than achieved in current implementations. Coupling ofhydrology and geochemistry with mechanics, for example,
is still in the initial stages of development. More sophisticated (and inevitably complex) treatments of microbialcommunity function that are coupled to environmental conditions are also possible in the future. Treatments of multiscale phenomena are only now beginning to be explored,
in most cases through multiple continuum approaches. This
topic is likely to be a primary focus in the future [222,226]. Eventually, one expects to have subsurface environmental simulators that couple to surface water and takemore fully into account the interactions with plants and theatmosphere. In this regard, links to landforms are likely to
become important.

References
1. Li, L., Steefel, C.I., Yang, L.: Scale dependence of mineraldissolution rates within single pores and fractures. Geochim.
Cosmochim. Acta. 72, 360–377 (2008)
2. Molins, S., Trebotich, D., Steefel, C.I., Shen, C.: An investigationof the effect of pore scale flow on average geochemical reaction
rates using direct numerical simulation. Water Resour. Res. 48,
W03527 (2012)
3. Steefel, C.I., Molins, S., Trebotich, D.: Pore scale processes associated with subsurface CO2 injection and sequestration. Rev.
Mineral. Geochem. 77, 259–303 (2013)
4. Oostrom, M., Mehmani, Y., Romero-Gomez, P., Tang, Y., Liu,
H., Yoon, H., Zhang, C.: Pore-scale and continuum simulationsof solute transport micromodel benchmark experiments. Comput.
Geosci., 1–23 (2014)
5. Aris, R.: Prolegomena to the rational analysis of systems ofchemical reactions. Arch. Ration. Mech. Anal. 19, 81–99 (1965)
6. Bowen, R.M.: On the stoichiometry of chemically reacting materials. Arch. Ration. Mech. Anal. 29, 114–124(1968)
7. Hooyman, G.J.: On thermodynamic coupling of chemical reactions. Proc. Natl. Acad. Sci. 47, 1169–1173 (1961)
8. Kirkner, D.J., Reeves, H.: Multicomponent mass transport withhomogeneous and heterogeneous chemical reactions: Effect ofchemistry on the choice of numerical algorithm I. Theory. Water
Resour. Res. 24, 1719–1729 (1988)
9. Lichtner, P.C.: Continuum model for simultaneous chemical
reactions and mass transport in hydrothermal systems. Geochim.
Cosmochim. Acta 49, 779–800 (1985)
10. Reed, M.H.: Calculation of multicomponent chemical equilibria and reaction processes in systems involving minerals, gases,
and an aqueous phase. Geochim. Cosmochim. Acta 46, 513–528(1982)

471
11. Van Zeggeren, F., Storey, S.H.: The computation of chemical equilibria, p. 176. Cambridge University Press, Cambridge(1970)
12. Steefel, C.I., Lasaga, A.C.: A coupled model for transport ofmultiple chemical species and kinetic precipitation/dissolution
reactions with application to reactive flow in single phasehydrothermal systems. Am. J. Sci. 294, 529–592 (1994)
13. Davis, J.A., Meece, D.E., Kohler, M., Curtis, G.P.: Approaches
to surface complexation modeling of uranium (VI) adsorption on
aquifer sediments. Geochim. Cosmochim. Acta 68, 3621–3641(2004)
14. Kulik, D.A.: Thermodynamic concepts in modeling sorption at
the mineral-water interface. Rev. Mineral. Geochem. 70, 125–
180 (2009)15. Dzombak, D.A., Morel, F.M.M.: Surface complexation modeling: hydrous ferric oxide, p. 393. Wiley, New York (1990)
16. Wang, Z., Giammar, D.E.: Mass action expressions for bidentate adsorption in surface complexation modeling: Theory andpractice. Environ. Sci. Technol. 47, 3982–3996 (2013)
17. Davis, J.A., Coston, J.A., Kent, D.B., Fuller, C.C.: Application of the surface complexation concept to complex mineral
assemblages. Environ. Sci. Technol. 32(19), 2820–2828 (1998)
18. Liu, C., Zachara, J.M., Qafoku, N.P., Wang, Z.: Scale-dependent
desorption of uranium from contaminated subsurface sediments.
Water Resour. Res. 8, 44 (2008)
19. Vanselow, A.P.: Equilibria of the base-exchange reactions of bentonites, permutites, soil colloids, and zeolites. Soil Sci. 33, 95–
114 (1932)
20. Sposito, G.: The thermodynamics of soil solutions, p. 223.
Oxford University Press, Oxford (1981)
21. Appelo, C.A.J., Postma, D.: Geochemistry, groundwater, andpollution. A.A. Balkema, Rotterdam, p 649 (1993)
22. Steefel, C.I., Carroll, S., Zhao, P., Roberts, S.: Cesium migrationin Hanford sediment: A multi-site cation exchange model basedon laboratory transport experiments. J. Contam. Hydrol. 67, 219–
246 (2003)
23. Aagaard, P., Helgeson, H.C.: Thermodynamic and kinetic constraints on reaction rates among minerals and aqueous solutions; I, Theoretical considerations. Am. J. Sci. 282(3), 237–
285 (1982)
24. Lasaga, A.C.: Chemical kinetics of water-rock interactions. J.
Geophys. Res. Solid Earth (1978–2012) 89(B6), 4009–4025(1984)
25. Rittmann, B.E., McCarty, P.L.: Environmental biotechnology.
McGraw Hill, New York (2001)
26. Jin, Q., Bethke, C.M.: Predicting the rate of microbial respirationin geochemical environments. Geochim. Cosmochim. Acta 69,
1133–1143 (2005)
27. Bear, J.: Dynamics of fluids in porous media, p. 764. Dover
Publications (1972)
28. Steefel, C.I., Maher, K.: Fluid-rock interaction: A reactive transport approach. Rev. Mineral. Geochem. 70, 485–532 (2009)
29. Neuman, S.P.: Saturated-unsaturated seepage by finite elements.
J. Hydrol. Div. Am. Soc. Civ. Eng. 99(HY12), 2233–2250 (1973)
30. Panday, S., Huyakorn, P.S., Therrien, R., Nichols, R.L.:
Improved three-dimensional finite-element techniques for field
simulation of variably-saturated flow and transport. J. Contam.
Hydrol. 12, 3–33 (1993)
31. Van Genuchten, M.T.: A closed form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.
44, 892–898 (1980)
32. Brooks, R.H., Corey, A.T.: Hydraulic properties of porous media.
Colorado State Univ., Hydrology Paper No. 3, p. 27 (1964)
33. Luckner, L., Van Genuchten, M.T., Nielsen, D.R.: A consistent
set of parametric models for the two-phase flow of immiscible

472
fluids in the subsurface. Water Resour. Res. 25(10), 2187–2193(1989)
34. Millington, R.J.: Gas diffusion in porous media. Science
130(3367), 100–102 (1959)
35. Mayer, K.U., Blowes, D.W., Frind, E.O.: Reactive transportmodeling for the treatment of an in situ reactive barrier forthe treatment of hexavalent chromium and trichloroethylene ingroundwater. Water Resour. Res. 37, 3091–3103 (2001)
36. Molins, S., Mayer, K.U.: Coupling between geochemical reactions and multicomponent gas diffusion and advection – a reactive transport modeling study. Water Resour. Res. 43, W05435(2007)
37. Cheadle, M.J., Elliott, M.T., McKenzie, D.: Percolation thresholdand permeability of crystallizing igneous rocks: The importanceof textural equilibrium. Geology 32, 757–760 (2004)
38. Verma, A., Pruess, K.: Thermohydrological conditions and silica redistribution near high-level nuclear wastes emplaced insaturated geological formations. J. Geophys. Res. Solid Earth
93(B2), 1159–1173 (1988)
39. Vaughan, P.J.: Analysis of permeability reduction during
flow of heated aqueous fluid through Westerly Granite.
In: Tsang, C.-F. (ed.) Coupled processes associated with
nuclear waste repositories. Academic Press, New York(1989)
40. Slider, H.C.: Practical petroleum reservoir engineering methods, An Energy Conservation Science. Petroleum Publishing
Company, Tulsa (1976)
41. Wu, Y.S.: On the effective continuum method for modelingmultiphase flow, multicomponent transport, and heat transfer infractured rock. In: Faybishenk, B., Witherspoon, P.A., Benson,
S.M. (eds.) Dynamics of Fluids Fractured Rock, pp. 299–312(2000)
42. Barenblatt, G.I., Zheltov, I.P., Kochina, I.N.: Basic concepts inthe theory of seepage of homogeneous liquids in fissured rocks[strata]. J. Appl. Math. Mech. 24, 1286–1303 (1960)
43. Warren, J.E., Root, P.J.: The behavior of naturally fractured
reservoirs. SPE J. 3, 245–255 (1963)
44. Pruess, K., Narasimhan, T.N.: A practical method for modeling
fluid and heat flow in fractured porous media. Soc. Pet. Eng. J.
25, 14–26 (1985)
45. Snow, D.T.: A parallel plate model of fractured permeable media,
Ph.D. dissertation, p. 331. University of California, Berkeley(1965)
46. Arora, B., Mohanty, B.P., McGuire, J.T.: Inverse estimationof parameters for multidomain flow models in soil columns
with different macropore densities. Water Resour. Res. 47, 1–17
(2011)
47. Wu, Y.S., Di, Y., Kang, Z., Fakcharoenphol, P.: A multiplecontinuum model for simulating single-phase and multiphaseflow in naturally fractured vuggy reservoirs. J. Pet. Sci. Eng. 78,
13–22 (2011)
48. Aradóttir, E.S.P., Sigfússon, B., Sonnenthal, E.L., Björnsson, G.,Jónsson, H.: Dynamics of basaltic glass dissolution–capturing
microscopic effects in continuum scale models. Geochim. Cosmochim. Acta 121, 311–327 (2013)
49. Lichtner, P.C.: Critique of dual continuum formulations of multicomponent reactive transport in fractured porous media. In: Faybishenk, B., Witherspoon, P.A., Benson, S.M. (eds.) Dynamics
of Fluids Fractured Rock, pp. 281–298 (2000)
50. Berkowitz, B.: Characterizing flow and transport in fractured
geological media: A review. Adv. Water Resour. 25, 861–884(2002)
51. Šimu̇nek, J., Jarvis, N.J., Van Genuchten, M.T., Gärdenäs,A.: Review and comparison of models for describing nonequilibrium and preferential flow and transport in the vadoseComput Geosci (2015) 19:445–478
zone. J. Hydrol. 272, 14–35 (2003)
52. MacQuarrie, K.T., Mayer, K.U.: Reactive transport modeling in
fractured rock: A state-of-the-science review. Earth Sci. Rev. 72,
189–227 (2005)
53. Somerton, W.H., El-Shaarani, A.H., Mobarak, S.M.: High temperature behavior of rocks associated with geothermal type reservoirs. In: SPE California Regional Meeting. Society of Petroleum
Engineers (1974)
54. Steefel, C.I., MacQuarrie, K.T.B.: Approaches to modeling reactive transport in porous media. In: Lichtner, P.C., Steefel, C.I.,
Oelkers, E.H. (eds.) Reactive Transport in Porous Media, vol. 34,
pp. 83–125. Reviews in Mineralogy (1996)
55. Kulik, D.A., Wagner, T., Dmytrieva, S.V., Kosakowski,
G., Hingerl, F.F., Chudnenko, K.V., Berner, U.R.: GEMSelektor geochemical modeling package: Revised algorithm and
GEMS3K numerical kernel for coupled simulation codes. Comput. Geosci. 17(1), 1–24 (2013)
56. Leal, A.M.M., Blunt, M.J., LaForce, T.C.: Efficient chemical
equilibrium calculations for geochemical speciation and reactive transport modelling. Geochim. Cosmochim. 131, 301–322(2014)
57. Saaltink, M.W., Carrera, J., Ayora, C.: On the behavior of
approaches to simulate reactive transport. J. Contam. Hydrol. 48,
213–235 (2001)
58. Calderhead, A., Mayer, K.U.: Comparison of the suitabilityof the global implicit method and the sequential non-iterative
approach for multicomponent reactive transport modelling. In:
Proceedings of 5th Joint IAH-CNC/CGS Conference, Québec
City, Québec, Canada, pp. 24–28 (2004)
59. Lichtner, P.C.: Time-space continuum description of fluid/rockinteraction in permeable media. Water Resour. Res. 28, 3135–
3155 (1992)
60. Parkhurst, D.L., Appelo, C.A.J.: User’s guide to PHREEQ C(Version 2) – a computer program for speciation, batch-reaction,
one-dimensional transport and inverse geochemical calculations,
Water-Resources Investigations, Report 99–4259, Denver, Co,
USA, p. 312 (1999)
61. Parkhurst, D.L., Appelo, C.A.J.: Description of input and examples for PHREEQC version 3—a computer program for speciation, batch-reaction, one-dimensional transport, and inverse
geochemical calculations, U.S. Geological Survey Techniquesand Methods, book 6, chap. A43, 497 p., available only at http://
pubs.usgs.gov/tm/06/a43 (2013)
62. Appelo, C.A.J., Parkhurst, D.L., Post, V.E.A.: Equations forcalculating hydrogeochemical reactions of minerals and gases
such as CO2 at high pressures and temperatures. Geochim.
Cosmochim. Acta 125, 49–67 (2014)
63. Thorstenson, D.C., Parkhurst, D.L.: Calculation of individual isotope equilibrium constants for geochemical reactions. Geochim. Cosmochim. Acta
11(68), 2449–2465(2004)
64. Charlton, S.R., Parkhurst, D.L.: Phast4Windows: A 3D graphical user interface for the reactive-transport simulator PHAST.
Groundwater 4(51), 623–628 (2013)
65. Šimu̇nek, J., Jacques, D., Šejna, M., van Genuchten, M.T.:
The HP2 Program for HYDRUS (2D/3D): A coupled code forsimulating two-dimensional variably-saturated water flow, heattransport, and biogeochemistry in porous media, Version 1.0, PC
Progress, Prague, Czech Republic, p. 76 (2012)
66. Šimu̇nek, J., Jacques, D., Langergraber, G., Bradford, S.A.,
Šejna, M., van Genuchten, M.T.: Numerical modeling of contaminant transport using HYDRUS and its specialized modules,
vol. 93, pp. 265–284 (2013). ISSN: 0970-4140 Coden-JIISAD
67. Parkhurst, D.L., Kipp, K.L., Charlton, S.R.: PHAST Version 2—
a program for simulating groundwater flow, solute transport, andComput Geosci (2015) 19:445–478
multicomponent geochemical reactions: U.S. Geological Survey
Techniques and Methods 6–A35, p. 235 (2010)
68. Prommer, H., Post, V.E.A.: PHT3D, A Reactive Multicomponent
Transport Model for Saturated Porous Media. User’s Manual
v2.10 (2010). http://www.pht3d.org
69. Šiu̇nek, J., van Genuchten, M.T.: Modeling non-equilibrium flow
and transport processes using HYDRUS. Vadose Zone J. 7(2),782–797 (2008)
70. Jacques, D., Šimu̇nek, J.: User Manual of the Multicomponent Variably-Saturated Flow and Transport Model HP1,
Description, Verification and Examples, Version 1.0, SCK•CENBLG-998, Waste and Disposal, SCK•CEN, Mol, Belgium,
p. 79 (2005)
71. Jacques, D., Šimu̇nek, J., Mallants, D., van Genuchten, M.T.:
Operator-splitting errors in coupled reactive transport codes
for transient variably saturated flow and contaminant transport in layered soil profiles. J. Contam. Hydrol. 88, 197–218(2006)
72. Šimu̇nek, J., Hopmans, J.W.: Modeling compensated root
water and nutrient uptake. Ecological Modeling 220, 505–
521 (2009)73. Jacques, D., Šimu̇nek, J., Mallants, D., van Genuchten, M.T., Yu,
L.: A coupled reactive transport model for contaminant leaching from cementitious waste matrices accounting for solid phasealterations. In: Proceedings Sardinia 2011, Thirteenth International Waste Management and Landfill Symposium (2011)
74. Celia, M.A., Bououtas, E.T., Zarba, R.L.: A general massconservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26, 1483–1496 (1990)
75. Jacques, D., Šimu̇nek, J., Mallants, D., van Genuchten, M.T.:
Modelling coupled water flow, solute transport and geochemical reactions affection heavy metal migration in a Podzol soil.
Geoderma 145, 449–461 (2008)
76. Jacques, D., Šimu̇nek, J., Mallants, D., van Genuchten, M.T.:
Modeling coupled hydrogeologic and chemical processes: Longterm uranium transport following phosphorus fertilization.
Vadose Zone J. 7(2), 698–711 (2008)
77. Jacques, D., Smith, C., Šimu̇nek, J., Smiles, D.: Inverse optimization of hydraulic, solute transport, and cation exchange
parameters using HP1 and UCODE to simulate cation exchange.
J. Contam. Hydrol. 142–143, 109–125 (2012)78. Bessinger, B.A., Marks, C.D.: Treatment of mercurycontaminated soils with activated carbon: A laboratory, field,
and modeling study. Remediation (The journal of Environmental Cleanup Costs, Technologies, & Techniques), 21(1), 115–
135 (2010). doi:10.1002/rem.20275
79. Leterme, B., Jacques, D.: Modelling of mercury fate and transport in soil systems. In: Proceedings of the 12th International
UFZ-Deltares Conference on Groundwater-Soil-Systems andWater Resource Management - AQUACONSOIL 2013 (2013).
www.aquaconsoil.org/AquaConSoil2013/Proceedings.html
80. Leterme, B., Blanc, P., Jacques, D.: A reactive transport model
for mercury fate in soil – application to different anthropogenic pollution sources. Environmental Science and Pollution
Research, (under review) (2014)
81. Jacques, D., Maes, N., Perko, J., Seetharam, S.C., Phung, Q.T.,
Patel, R., Soto, A., Liu, S., Wang, L., De Schutter, G., Ye, G.,
van Breugel, K.: Concrete in engineered barriers for radioactive
waste disposal facilities - phenomenological study and assessment of long term performance. In: Proceedings of the ASME
2013 15th International Conference on Environmental Remediation and Radioactive Waste Management, ICEM2013-96282(2013)

473
82. Thaysen, E.M., Jacques, D., Jessen, S., Andersen, C.E., Laloy,
E., Ambus, P., Postma, D., Jakobsen, I.: Controls on carbon dioxide fluxes across the unsaturated zone of cropped and unplanted
soil mesocosms. Biogeosciences 11, 4251–4299 (2014)
83. Zheng, C., Wang, P.P.: MT3DMS: A modular three-dimensional
multispecies transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems: documentation and user’s guide. Contract Report
SERDP-99–1, U.S. Army Engineer Research and Development
Center, Vicksburg, Mississippi (available at http://hydro.geo.ua.
edu/mt3d) (1999)
84. Prommer, H., Post, V.E.A.: PHT3D, A Reactive multicomponent
transport model for saturated porous media. User’s Manual v2.10(2010). http://www.pht3d.org
85. Prommer, H., Barry, D.A., Zheng, C.: MODFLOW/MT3DMSbased reactive multicomponent transport modeling. Ground
Water 41, 247–257 (2003)
86. Appelo, C.A.J., Rolle, M.: PHT3D: A reactive multicomponent
transport model for saturated porous media. Ground Water 5(48),627–632 (2010)
87. Zuurbier, K.G., Hartog, N., Valstar, J., Post, V.E.A., van
Breukelen, B.M.: The impact of low-temperature seasonalaquifer thermal energy storage (SATES) systems on chlorinatedsolvent contaminated groundwater: modeling of spreading anddegradation. J. Contam. Hydrol. 147, 1–13 (2013)
88. Wallis, I., Prommer, H., Post, V.E.A., Vandenbohede, A.,
Simmons, C.T.: Simulating MODFLOW-based reactive transport
under radially symmetric flow conditions. Groundwater 51, 398–
413 (2013)
89. Wu, M.Z., Reynolds, D.A., Prommer, H., Fourie, A., Thomas,
D.G.: Numerical evaluation of voltage gradient constraints onelectrokinetic injection of amendments. Adv. Water Resour. 38,
60–69 (2012)
90. Chiang, W.H. 3D Groundwater modeling with PMWIN, 2nd ed.,p. 397. Springer-Verlag, The Netherlands (2005)
91. Nagel, T., Shao, H., Singh, A.K., Watanabe, N., Roßkopf, C.,
Linder, M., Wörner, A., Kolditz, O.: Non-equilibrium thermochemical heat storage in porous media: Part 1 – conceptualmodel. Energy 60, 254–270 (2013)
92. Kolditz, O., Bauer, S., Bilke, L., Böttcher, N., Delfs, J.O.,
Fischer, T., Görke, U.J., Kalbacher, T., Kosakowski, G.,
McDermott, C.I., Park, C.H., Radu, F., Rink, K., Shao, H., Shao,
H.B., Sun, F., Sun, Y.Y., Singh, A.K., Taron, J., Walther, M.,
Wang, W., Watanabe, N., Wu, N., Xie, M., Xu, W., Zehner, B.:
OpenGeoSys: An open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes
in porous media. Environ. Earth Sci. 67, 589–599 (2012)
93. Rutqvist, J., Barr, D., Birkholzer, J.T., Chijimatsu, M., Kolditz,
O., Liu, Q., Oda, Y., Wang, W., Zhang, C.: Results from an international simulation study on coupled thermal, hydrological, andmechanical processes near geological nuclear waste repositories.
J. Nucl. Technol. 163, 101–109 (2008)
94. Kolditz, O., Bauer, S., Beyer, C., Böttcher, N., Dietrich, P.,Görke, U.-J., Kalbacher, T., Park, C.-H., Sauer, U., Schütze, C.,
Shao, H.B., Singh, A.K., Taron, J., Wang, W., Watanabe, N.: A
systematic benchmarking approach for geologic CO2 injection
and storage. Environ. Earth Sci. 67, 613–632 (2012)
95. Mukhopadhyay, S., Birkholzer, J.T., Nicot, J.P., Hosseini, S.A.:
A single site multi-model comparative study for CO2 injection
field test: An introduction to Sim-SEQ. Environ. Earth Sci. 67,
601–611 (2012)
96. Kolditz, O., Shao, H., Görke, U.-J., Wang, W. (eds.): Thermohydro-mechanical-chemical processes in fractured porous media.
Lecture Notes in Computational Science and Engineering,
vol. 86. Springer, Heidelberg (2012)

474
97. Park, C.-H., Beyer, C., Bauer, S., Kolditz, O.: A study of preferential flow in heterogeneous media using random walk particle
tracking. Geosci. J. 12, 285–297 (2008)
98. Gräbe, A., Rink, K., Fischer, T., Sun, F., Wang, W., Rödiger, T.,
Siebert, C., Kolditz, O.: Numerical analysis of the groundwater
regime in the Western Dead Sea Escarpment. Environ. Earth Sci.
69, 571–586 (2013)
99. Singh, A.K., Baumann, G., Henninges, J., Görke, U.-J., Kolditz,
O.: Numerical analysis of thermal effects during carbon dioxide injection with enhanced gas recovery: A theoretical casestudy for the Altmark gas field. Environ. Earth Sci. 67, 497–
509 (2012)100. Wang, W., Rutqvist, J., Gorke, U.-J., Birkholzer, J.T., Kolditz,
O.: Non-isothermal flow in low permeable porous media: A comparison of unsaturated and two-phase flow approaches. Environ.
Earth Sci. 62, 1197–1207 (2011)
101. Kolditz, O., de Jonge, J.: Non-isothermal two-phase flow inporous media. Comput. Mech. 33, 345–364 (2004)
102. Park, C.-H., Boettcher, N., Wang, W., Kolditz, O.: Are upwind
techniques in multi-phase flow models necessary. J Comp. Phys.
230, 8304–8312 (2011)
103. Delfs, J.-O., Park, C.-H., Kolditz, O.: A sensitivity analysis of
Hortonian flow. Adv. Water Resour. 32, 1386–1395 (2009)
104. Böttcher, N., Kolditz, O., Liedl, R.: Evaluation of equations ofstate for CO2 in numerical simulations. Environ. Earth Sci. 67,
481–495 (2012)
105. Watanabe, N., McDermott, C., Wang, W., Taniguchi, T., Kolditz,
O.: Uncertainty analysis of thermo-hydro-mechanical processes
in heterogeneous porous media. Comput. Mech. 45, 263–
280 (2010)
106. Kalbacher, T., Delfs, J.O., Shao, H., Wang, W., Walther, M.,
Samaniego, L., Schneider, C., Musolff, A., Centler, F., Sun, F.,
Hildebrandt, A., Liedl, R., Borchardt, D., Krebs, P., Kolditz,
O.: The IWAS-ToolBox: Software coupling for an integrated
water resources management. Environ. Earth Sci. 65, 1367–
1380 (2012)107. Xie, M., Bauer, S., Kolditz, O., Nowak, T., Shao, H.: Numerical simulation of reactive processes in an experiment withpartially saturated bentonite. J Contam. Hydrol. 83, 122–
147 (2006)
108. Shao, H., Kolditz, O., Kulik, D.A., Pfingsten, W., Kosakowski,
G.: Reactive transport of multiple non-ideal solid solutions. Appl.
Geochem. 24, 1287–1300 (2009)109. Beyer, C., Bauer, S., Kolditz, O.: Uncertainty assessment of contaminant plume length estimates in heterogeneous aquifers. J.
Contam. Hydrol. 87, 73–95 (2006)
110. Centler, F., Shao, H., Park, C.-H., de Biase, C., Kolditz, O.,
Thullner, M.: GeoSysBRNS – a flexible multi-dimensional reactive transport model for simulating biogeochemical subsurface
processes. Comput. Geosci. 36, 397–405 (2010)111. Rink, K., Kalbacher, T., Kolditz, O.: Visual data management forhydrological analysis. Environ. Earth Sci. 65, 1395–1403 (2012)112. van der Lee, J., De Windt, L., Lagneau, V., Goblet, P.: Moduleoriented modeling of reactive transport with HYTEC. Comput.
Geosci. 29, 265–275 (2003)
113. Wolery, T.: EQ3/6. A software package for geochemical modelling of aqueous systems: Package overview and installation guide (version 7.0). Technical Report UCRL-MA-110662
PT I ed. Lawrence Livermore National Laboratory, USA(1992)
114. van der Lee, J., Lagneau, V.: Rigorous methods for reactive
transport in unsaturated porous medium coupled with chemistry and variable porosity. In: Miller, C.T., Farthing, M.W.,
Gray, W.G., Pinder, G.F. (eds.) Computational Methods in
Water Resources (CMWR XV), vol. 48, pp. 861–868. Elsevier(2004)

Comput Geosci (2015) 19:445–478
115. Lagneau, V., van der Lee, J.: Operator-splitting-based reactive
transport models in strong feedback of porosity change: The
contribution of analytical solutions for accuracy validation andestimator improvement. J. Contam. Hydrol. 112, 118–129 (2010)116. Debure, M., De Windt, L., Frugier, P., Gin, S.: HLW glass dissolution in the presence of magnesium carbonate: Diffusion cell
experiment and coupled modeling of diffusion and geochemicalinteractions. J. Nucl. Mater. 443, 507–521 (2013)
117. De Windt, L., Marsal, F., Corvisier, J., Pellegrini, D.: Modelingof oxygen gas diffusion and consumption during the oxic transient in a disposal cell of radioactive waste. Appl. Geochem. 41,
115–127 (2014)
118. Lagneau, V., Pipart, A., Catalette, H.: Reactive transport modelling and long term behaviour of CO2 sequestration in saline
aquifers. Oil Gas Sci. Technol. 60, 231–247 (2005)
119. Jacquemet, N., Pironon, J., Lagneau, V., Saint-Marc, J.: Armouring of well cement in H2 S-CO2 saturated brine by calcite coating– experiments and numerical modeling. Appl. Geochem. 27,
782–795 (2012)120. De Windt, L., Devillers, P.: Modeling the degradation of Portland
cement pastes by biogenic organic acids. Cem. Concr. Res. 40,
1165–1174 (2010)
121. De Windt, L., Deneele, D., Maubec, N.: Kinetics oflime/bentonite pozzolanic reactions at 20 and 50 °C: Batch tests
and modeling. Cem. Concr. Res. 59, 34–42 (2014)122. Dabo, D., Badreddine, R., De Windt, L., Drouadaine, I.: Tenyear chemical evolution of leachate and municipal solid waste
incineration bottom ash used in a test road site. J. Hazard. Mater.
172, 904–913 (2007)
123. De Windt, L., Badreddine, R., Lagneau, V.: Long-term reactive
transport modelling of stabilized/solidified waste: From dynamic
leaching tests to disposal scenarios. J. Hazard. Mater. 139, 529–
536 (2007)
124. De Windt, L., Badreddine, R.: Modelling of long-term dynamic
leaching tests applied to solidified/stabilised waste. Waste
Manag. 27, 1638–1647 (2007)125. Corvisier, J., Bonvalot, A.F., Lagneau, V., Chiquet, P., Renard, S.,
Sterpenich, J., Pironon, J.: Impact of co-injected gases on CO2
storage sites: Geochemical modeling of experimental results.
Proceedings of the International Conference on Greenhouse Gas
Technology 11, vol. 37, pp. 3699–3710. Energy Procedia, Kyoto(2013)
126. Meeussen, J.C.: ORCHESTRA: An object-oriented framework
for implementing chemical equilibrium models. Environ. Sci.
Technol. 37, 1175–1182 (2003)
127. Farmer, J.G., Graham, M.C., Thomas, R.P., Licona-Manzur,C., Paterson, E., Campbell, C.D., Geelhoed, J.S., Lumsdon,
D.G., Meeussen, J.C.L., Roe, M.J., Conner, A., Fallick, A.E.,
Bewley, R.J.F.: Assessment and modelling of the environmental
chemistry and potential for remediative treatment of chromiumcontaminated land. Environ. Geochem. Health 21, 331–33(1999)
128. Geelhoed, J.S., Meeussen, J.C.L., Hillier, S., Lumsdon, D.G.,
Thomas, R.P., Farmer, J.G., Paterson, E.: Identification andgeochemical modeling of processes controlling leaching of
Cr(VI) and other major elements from chromite ore processing residue. Geochim. Cosmochim. Acta 66, 3927–3942(2002)
129. Filius, J.D., Lumsdon, D.G., Meeussen, J.C.L., Hiemstra, T.,
Van Riemsdijk, W.H.: Adsorption of fulvic acid on goethite.
Geochim. Cosmochim. Acta 64, 51–60 (2000)
130. Weng, L., Hiemstra, T., Meeussen, J.C.L., Koopal, L., Van
Riemsdijk, W.H.: Interactions of calcium and fulvic acid at thegoethite-water interface. Geochim. Cosmochim. Acta 69, 325–
339 (2005)

Comput Geosci (2015) 19:445–478
131. Van Riemsdijk, W.H., Koopal, L.K., Kinniburgh, D.G.,
Benedetti, M.F., Weng, L.: Modeling the interactions between
humics, ions, and mineral surfaces. Environ. Sci. Technol. 40,
7473–7480 (2006)
132. Vink, J.P.M., Meeussen, J.C.L.: BIOCHEM–ORCHESTRA: A
tool for evaluating chemical speciation and ecotoxicologicalimpacts of heavy metals on river flood plain systems. Environ.
Pollut. 148, 833–841 (2007)
133. Nowack, B., Mayer, K.U., Oswald, S.E., Van Beinum, W.,
Appelo, C.A.J., Jacques, D., Seuntjens, P.: Verification and intercomparison of reactive transport codes to describe root-uptake.
Plant Soil 285, 305–321 (2006)
134. Van der Sloot, H.A., van Zomeren, A.: Characterisation leaching
tests and associated geochemical speciation modelling to assesslong term release behaviour from extractive wastes. Mine Water
Environ. 31, 92–103 (2012)135. Brown, K.G., Arnold, J., Sarkar, S., Flach, G., van der Sloot, H.,
Meeussen, J.C.L., Kosson, D.S.: Modeling carbonation of highlevel waste tank integrity and closure. In EPJ Web of Conferences
56, 05003, EDP Sciences (2013)136. Sarkar, S., Kosson, D.S., Mahadevan, S., Meeussen, J.C.L., vander Sloot, H., Arnold, J.R., Brown, K.G.: Bayesian calibration ofthermodynamic parameters for geochemical speciation modelingof cementitious materials. Cem. Concr. Res. 42, 889–902 (2012)137. Xu, T., Pruess, K.: Modeling multiphase non-isothermal fluid
flow and reactive geochemical transport in variably saturated
fractured rocks: 1, Methodology. Am. J. Sci. 301, 16–33 (2001)
138. Xu, T., Sonnenthal, E., Spycher, N., Pruess, K.: TOUGHREACT - a simulation program for non-isothermal multiphasereactive geochemical transport in variably saturated geologic
Media: Applications for geothermal injectivity and CO2 geologic
sequestration. Comput. Geosci. 32, 145–165 (2006)
139. Xu, T., Spycher, N., Sonnenthal, E., Zhang, G., Zheng, L.,
Pruess, K.: TOUGHREACT Version 2.0: A simulator for subsurface reactive transport under non-isothermal multiphase flow
conditions. Comput. Geosci. 37, 763–774 (2011)
140. Pruess, K., Oldenburg, C., Moridis, G.: TOUGH2 user’s guide,
version 2.0. Lawrence Berkeley Laboratory Report LBL-43134,
Berkeley, CA, p. 192 (1999)
141. Finsterle, S., Sonnenthal, E.L., Spycher, N.: Advances in subsurface modeling: The TOUGH suite of simulators. Comput.
Geosci. 65, 2–12 (2014)142. Pruess, K., Narasimhan, T.N.: On fluid reserves and the production of superheated steam from fractured, vapor-dominatedgeothermal reservoirs. J. Geophys. Res. 87, 9329–9339 (1982)
143. Dobson, P.F., Salah, S., Spycher, N., Sonnenthal, E.L.: Simulations of water rock interactions in the Yellowstone geothermal
system using TOUGHREACT. Geothermics 33, 493–502 (2004)
144. Xu, T., Ontoy, Y., Molling, P., Spycher, N., Parini, M., Pruess,
K.: Reactive transport modeling of injection well scaling andacidizing at the Tiwi field, Philippines. Geothermics 33, 477–491(2004)
145. Wanner, C., Peiffer, L., Sonnenthal, E.L., Spycher, N., Iovenitti,J., Kennedy, B.M.: Reactive transport modeling of the Dixie
Valley geothermal area: Insights on flow and geothermometry.
Geothermics 51, 130–141 (2014)146. Spycher, N., Sonnenthal, E.L., Apps, J.: Fluid flow and reactive
transport around potential nuclear waste emplacement tunnels at
Yucca Mountain, Nevada. J. Contam. Hydrol. 62–63, 653–674(2003)
147. Sonnenthal, E.L., Ito, A., Spycher, N., Yui, M., Apps, J., Sugita,
Y., Conrad, M., Kawakami, S.: Approaches to modeling coupled
thermal, hydrological, and chemical processes In the Drift Scale
Heater Test at Yucca Mountain. Int. J. Rock. Mech. Min. Sci. 42,
698–719 (2005)

475
148. Xu, T., Senger, R., Finsterle, S.: Corrosion-induced gas generation in a nuclear waste repository: Reactive geochemistry andmultiphase flow effects. Appl. Geochem. 23, 3423–3433 (2008)
149. Marty, N.C.M., Tournassat, C., Burnol, A., Giffaut, E., Gaucher,
E.C.: Influence of reaction kinetics and mesh refinement on
the numerical modelling of concrete/clay interactions. J. Hydrol.
364, 58–72 (2009)
150. Xu, T., Apps, J.A., Pruess, K.: Mineral sequestration of carbon
dioxide in a sandstone–shale system. Chem. Geol. 217, 295–318(2005)
151. Audigane, P., Gaus, I., Czernichowki-Lauriol, I., Pruess, K., Xu,
T.: Two-dimensional reactive transport modeling of CO2 injection in a saline aquifer at the Sleipner site. Am. J. Sci. 307, 974–
1008 (2007)
152. Aradóttir, E.S.P., Sonnenthal, E.L., Björnsson, G., Jónsson, H.:
Multidimensional reactive transport modeling of CO2 mineralsequestration in basalts at the Hellisheidi geothermal field, Iceland. Int. J. Greenh. Gas. Con. 9, 24–40 (2012)153. Zheng, L., Spycher, N., Birkholzer, J., Xu, T., Apps, J., Kharaka,
Y.: On modeling the potential impacts of CO2 sequestration on
shallow groundwater: Transport of organics and co-injected H2S
by supercritical CO2 to shallow aquifers. Int. J. Greenh. Gas.
Con. 14, 113–127 (2013)154. Dalkhaa, C., Shevalier, M., Nightingale, M.: Mayer, B.: 2-D reactive transport modeling of the fate of CO2 injected into a saline
aquifer in the Wabamun Lake Area, Alberta, Canada. Appl.
Geochem. 38, 10–23 (2013)155. Wu, Y., Ajo-Franklin, J.B., Spycher, N., Hubbard, S.S., Zhang,G., Williams, K.H., Taylor, J., Fujita, Y., Smith, R.: Geophysicalmonitoring and reactive transport modeling of ureolyticallydriven calcium carbonate precipitation. Geochem. Trans. 12(2011)
156. Meima, J.A., Graupner, T., Rammlmair, D.: Modeling the effectof stratification on cemented layer formation in sulfide-bearing
mine tailings. Appl. Geochem. 27, 124–137 (2012)157. Bea, S.A., Wainwright, H., Spycher, N., Faybishenko, B.,
Hubbard, S.S., Denham, M.E.: Identifying key controls on the
behavior of an acidic-U (VI) plume in the Savannah River Siteusing reactive transport modeling. J. Contam. Hydrol. 151, 34–
54 (2013)
158. Kim, J., Sonnenthal, E., Rutqvist, J.A.: A sequential implicit
algorithm of chemo-thermo-poro-mechanics for fractured
geothermal reservoirs. Computers and Geosciences (2014)159. White, M.D., Oostrom, M.: STOMP subsurface transport overmultiple phases version 4.0 user’s guide. Pacific Northwest
National Laboratory, Washington (2006)
160. White, M.D., McGrail, B.P.: STOMP subsurface transport
over multiple phases version 1.0, Addendum: ECKEChemequilibrium-conservation-kinetic equation chemistry and reactive transport. PNNL-15482. Pacific Northwest National Laboratory, Washington (2005)
161. Fang, Y.L., Yabusaki, S.B., Yeh, G.T.: A general simulator forreaction-based biogeochemical processes. Comput. Geosci. 32,
64–72 (2006)
162. Balay, S., Brown, J., Buschelman, K., Eijkhout, V., Gropp, W.D.,
Kaushik, D., Knepley, M.G., Curfman McInnes, L., Smith, B.F.,
Zhang, H.: PETSc Users Manual, ANL–95/11 – revision 3.4.
Argonne National Laboratory (2013)
163. Nieplocha, J., Palmer, B., Tipparaju, V., Krishnan, M., Trease,
H., Apra, E.: Advances, applications and performance of the
Global Arrays Shared Memory Programming Toolkit. Int. J.
High Perform. Comput. Appl. 20, 203–231 (2006)
164. Yeh, G.T., Tripathi, V.S.: HYDROGEOCHEM: A coupled modelof HYDROlogical transport and GEOCHEMical equilibrium ofmulti component systems, ORNL 6371, Oak Ridge National

476
Laboratory, p. 37831. Oak Ridge National Laboratory, Oak
Ridge (1990)
165. Yeh, G.T., Siegel, M.D., Li, M.H.: Numerical modeling of
coupled fluid Fflows and reactive transport including fast andslow chemical reactions. J. Contam. Hydrol. 379–390(2001), 47(2001)
166. Yeh, G.T., Li, Y., Jardine, P.M., Burgos, W.D., Fang, Y.L.,
Li, M.H., Siegel, M.D.: HYDROGEOCHEM 4.0: A coupledmodel of fluid flow, thermal transport, and HYDROGEOCHEMical transport through saturated unsaturated media Version 4.0.
ORNL/TM-2004/103, p. 37831. Oak Ridge National Laboratory,
Oak Ridge (2004)
167. Yeh, G.T., Sun, J.T., Jardine, P.M., Burgos, W.D., Fang, Y.L., Li,
M.H., Siegel, M.D.: HYDROGEOCHEM 5.0: A three dimensional model of coupled fluid flow, thermal transport, and
HYDROGEOCHEMical transport through variably saturated
conditions version 5.0. ORNL/TM-2004/107, p. 37831. Oak
Ridge National Laboratory, Oak Ridge (2004)
168. Yeh, G.T., Tsai, C.H.: HYDROGEOCHEM 6.0: A twodimensional model of coupled fluid flow, thermal transport,
geomechanics, and HYDROGEOCHEMical transport throughmultiple phase systems version 6.0 (FACTM2D: A Model for
multi-phase flow analysis and reactive chemical transport, thermal transport, and mechanics simulation, 2-dimensional version)- theoretical basis and numerical approximation. Graduate Institute of Applied Geology, National Central University, Jhongli(2013)
169. Yeh, G.T., Tsai, C.H., Ni, C.F.: HYDROGEOCHEM
6.0: A model to couple thermal-hydrology-mechanicschemical (THMC) processes user guide. Graduate Instituteof Applied Geology, National Central University, Jhongli(2013)
170. Yeh, G.T., Fang, Y.L., Zhang, F., Sun, J.T., Li, Y., Li, M.H.,
Siegel, M.D.: Numerical modeling of coupled fluid flow andthermal and reactive biogeochemical transport in porous andfractured media. Comput. Geosci. 14, 149–170 (2010)
171. Yeh, G.T., Tripathi, V.J., Gwo, J.P., Cheng, H.P., Cheng, R.J.,
Salvage, K.M., Li, M.H., Fang, Y.L., Li, Y., Sun, J.T., Zhang, F.,
Siegel, M.D.: Chapter 1: HYDROGEOCGEM: A coupled modelof variably saturated flow, thermal transport, and reactive biogeochemical transport. Groundwater reactive transport models.
In: Zhang, F., Yeh, G.T., Parker, J.C. (eds.) Bentham e-Books.
Bentham Science Publishers (2012). http://www.bentham.org
172. Smith, J.M.: Chemical engineering kinetics, p. 676. R. R.
Donnelley & Sons Company (1981)
173. Chilakapati, A., Ginn, T., Szecsody, J.: An analysis of complex
reaction networks in groundwater modeling. Water Resour. Res.
34, 1767–1780 (1998)
174. Fang, Y., Yeh, G.T., Burgos, W.D.: A New Paradigm to model
reaction-based biogeochemical processes. Water Resour. Res. 39,
1083–1108 (2003)
175. Kräutle, S., Knabner, P.: A reduction scheme for coupled multicomponent transport-reaction problems in porous
media: Generalization to problems with heterogeneous equilibrium reactions. Water Resour. Res. 43, W03429 (2007).
doi:10.1029/2005WR004465
176. Yeh, G.T., Tsai, C.H., Fang, Y., Yabusaki, S., Li, M.H.:
BIOGEOCHEM 1.5: A numerical model to simulate BIOGEOCHEMical reactions under multiple phase system. Graduate Institute of Applied Geology, National Central University,
Jhongli (2014)
177. Yeh, G.T., Tsai, C.H.: User’s manual for BIOGEOCHEM 1.5.
Graduate Institute of Applied Geology, National Central University, Jhongli (2014)

Comput Geosci (2015) 19:445–478
178. Yeh, G.T., Tripathi, V.S.: A critical evaluation of recent developments of hydrogeochemical transport models of reactive
multi-chemical components. Water Resour. Res. 25, 93–108(1989)
179. Liu, I.-S., Cipolatti, R.A., Rincon, M.A.: Successive linear
approximation for finite elasticity. Comput. Appl. Math. 29, 465–
478 (2010)180. Yeh, G.T., Tripathi, V.S.: A model for simulating transport of reactive multispecies components: Model development and demonstration. Water Resour. Res. 27, 3075–3094(1991)
181. Gwo, J.P., D’Azevedo, E.F., Frenzel, H., Mayes, M., Yeh, G.T.,
Jardine, P.M., Salvage, K.M., Hoffman, F.M.: HBGC123D: A
high performance computer model of coupled hydrogeologicaland bigeochemical processes. Comput. Geosci. 27, 1231–1242(2001)
182. Yeh, G.T., Gwo, J.P., Siegel, M.D., Li, M.H., Fang, Y.L., Zhang,F., Luo, W.S., Yabusaki, S.B.: Innovative mathematical modeling in environmental remediation. J. Environ. Radioact. (2011).
doi:10.1016/j.jenvrad.2011.06.010
183. Kent, D.B., Davis, J.A., Anderson, L.C.D., Rea, B.A., Waite,
T.D.: Transport of chromium and selenium in the suboxic zoneof a shallow aquifer: Influence of redox and adsorption reactions.
Water Resour. Res. 30, 1099–1114 (1994)
184. Abrams, R.H.: A compartmentalized approach to simulating
redox zones in contaminated aquifers. PhD Dissertation, Department of Geological and Environmental Sciences, Stanford University, Palo Alto, California (1999)
185. Bahr, J.M.: Keating, E.H.:. Redox geochemistry of shallow
groundwater discharging to Allequash Creek in northern Wisconsin. WRC GRR 95-08. Water Resources Center, Universityof Wisconsin-Madison, p. 41 (1995)
186. Tournassat, C., Appelo, C.A.J.: Modelling approaches for anionexclusion in compacted Na-bentonite. Geochim. Cosmochim.
Acta 75, 3698–3710 (2011)
187. Steefel, C.I., Van Cappellen, P.: A new kinetic approach to modeling water-rock interaction: The role of nucleation, precursors,
and Ostwald ripening. Geochim. Cosmochim. Acta 54, 2657–
2677 (1990)
188. Steefel, C.I., Lasaga, A.C.: A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow insingle phase hydrothermal systems. Am. J. Sci. 294, 529–592(1994)
189. Steefel, C.I., Lichtner, P.C.: Diffusion and reaction in rockmatrix bordering a hyperalkaline fluid-filled fracture. Geochim.
Cosmochim. Acta 58, 3592–3612 (1994)
190. Steefel, C.I., Lichtner, P.C.: Multicomponent reactive transport
in discrete fractures: I. Controls on reaction front geometry. J.
Hydrol. 209, 186–199 (1998)
191. Steefel, C.I., Lichtner, P.C.: Multicomponent reactive transport indiscrete fractures: II. Infiltration of hyperalkaline groundwater at
Maqarin, Jordan, a natural analogue site. J. Hydrol. 209, 200–224(1998)
192. Giambalvo, E.R., Steefel, C.I., Fisher, A.T., Rosenberg, N.D.,
Wheat, C.G.: Effect of fluid-sediment reaction on hydrothermal fluxes of major elements, eastern flank of the Juande Fuca Ridge. Geochim. Cosmochim. Acta 66, 1739–
1757 (2002)
193. Maher, K., Steefel, C.I., DePaolo, D., Viani, B.: The mineral dissolution rate conundrum: Insights from reactive transport modeling of U isotopes and pore fluid chemistry inmarine sediments. Geochim. Cosmochim. Acta 70, 337–
363 (2006)

Comput Geosci (2015) 19:445–478
194. Maher, K., Steefel, C.I., White, A.F., Stonestrom, D.A.: The
role of reaction affinity and secondary minerals in regulatingchemical weathering rates at the Santa Cruz Soil Chronosequence, California. Geochim. Cosmochim. Acta 73, 2804–2831(2009)
195. Navarre-Sitchler, A., Steefel, C.I., Sak, P.B., Brantley, S.L.: A
reactive transport model for weathering rind formation on basalt.
Geochim. Cosmochim. Acta 75, 7644–7667 (2011)196. Navarre-Sitchler, A., Steefel, C.I., Yang, L., Tomutsa, L.,
Brantley, S.L.: Evolution of porosity and diffusivity associatedwith chemical weathering of a basalt clast. J. Geophys. Res., 114(2009). doi:10.1029/2008JF001060
197. Li, L., Gawande, N., Kowalsky, M.B., Steefel, C.I., Hubbard,
S.S.: Physicochemical heterogeneity controls on uranium bioreduction rates at the field scale. Environ. Sci. Technol. 45, 9959–
9966 (2011)
198. Li, L., Steefel, C.I., Kowalsky, M.B., Englert, A., Hubbard, S.S.:
Effects of physical and geochemical heterogeneities on mineral transformation and biomass accumulation during uranium
bioremediation at Rifle, Colorado. J. Contam. Hydrol. 11, 45–63(2010)
199. Li, L., Steefel, C.I., Williams, K.H., Wilkins, M.J., Hubbard,
S.S.: Mineral transformation and biomass accumulation during
uranium bioremediation, Rifle, Colorado. Environ. Sci. Technol.
43, 5429–5435 (2009)
200. Druhan, J.L., Steefel, C.I., Molins, S., Williams, K.H., Conrad,
M.E., DePaolo, D.J.: Timing the onset of sulfate reduction overmultiple subsurface acetate amendments by measurement andmodeling of sulfur isotope fraction. Environ. Sci. Technol. 46,
8895–8902 (2012)201. Druhan, J.L., Steefel, C.I., Conrad, M.E., DePaolo, D.J.: A large
column analog experiment of stable isotope variations during
reactive transport: I. A comprehensive model of sulfur cyclingand δ 34S fractionation. Geochim. Cosmochim. Acta 124, 366–
393 (2014)202. Mayer, K.U., Frind, E.O., Blowes, D.W.: Multicomponent
reactive transport modeling in variably saturated porous
media using a generalized formulation for kinetically controlled reactions. Water Resour. Res. 38, 1174 (2002).
doi:10:1029/2001WR000862
203. Cheng, L.: Dual porosity reactive transport modeling, PhD thesis.
University of Sheffield, UK (2006)
204. Cheng, L., Lerner, D., Thornton, S., Mayer, K.U.: Managing
MTBE attenuation in a dual porosity chalk aquifer – field observations and modelling results. Proceedings of ModelCARE 2009
Groundwater and the Environment, p. 341. IAHS Publishing,
Wuhan (2009)
205. Amos, R.T., Mayer, K.U.: Investigating the role of gas bubble
formation and entrapment in contaminated aquifers: Reactive
transport modeling. J. Contam. Hydrol. 87, 123–154 (2006)
206. Amos, R.T., Mayer, K.U.: Investigating ebullition in a sandcolumn using dissolved gas analysis and reactive transport modelling, Environ. Sci. Technol. 40, 5361–5367 (2006)
207. Molins, S., Mayer, K.U., Scheutz, C., Kjeldsen, P.: Role of transport mechanisms in the attenuation of landfill gas in cover soils:
A multicomponent modelling study. J. Environ. Qual. 37, 459–
468 (2008)
208. Molins, S., Mayer, K.U., Amos, R.T., Bekins, B.A.: Vadose zone
attenuation of organic compounds at a crude oil spill site - Interactions between biogeochemical reactions and multicomponent
gas transport. J. Contam. Hydrol. 112, 15–29 (2010)
209. Henderson, T., Mayer, K.U., Parker, B., Al, T.: Threedimensional density-dependent flow and multicomponent reactive transport modeling of chlorinated solvent oxidation by potassium permanganate. J. Contam. Hydrol. 106, 195–211 (2009)

477
210. Bea, S.A., Mayer, K.U., MacQuarrie, K.T.B.: Modelling reactive transport in sedimentary rock environments - Phase II
MIN3P-THCm code enhancements and illustrative simulations
for a glaciation scenario. Technical report: NWMO TR-2011-13(2011)
211. Harvie, C.E., Moller, N., Weare, J.H.: The prediction of mineralsolubilities in natural waters: The Na-K-Mg-Ca-H-Cl-SO4 -OHHCO3 -CO3 -CO2 -H2 O system to high ionic strengths at 25 o C.
Geochim. Cosmochim. Acta 48, 723–751 (1984)
212. Pitzer, K.S.: Thermodynamics of electrolytes. I. Theoreticalbasis and general equations. J. Phys. Chem. 77, 268–277 (1973)
213. Monnin, C.: Density calculation and concentration scale conversions for natural waters. Comput. Geosci. 20, 1435–1445(1994)
214. Neuzil, C.: Hydromechanical coupling in geologic processes.
Hydrogeol. J. 11, 41–83 (2003)
215. Bea, S., Wilson, S., Mayer, K.U., Dipple, G., Power, I., Gamazo,P.: Reactive transport modeling of natural carbon sequestration
in ultramafic mine tailings. Vadose Zone J. 11, 1–17 (2012)
216. Sihota, N.J., Mayer, K.U.: Characterizing vadose zone hydrocarbon biodegradation using CO2 -effluxes, isotopes, andreactive transport modeling. Vadose Zone J., 11 (2012).
doi:10.2136/vzj2011.0204
217. Miller, G.R., Rubin, Y., Mayer, K.U., Benito, P.H.: Modeling vadose zone processes during land application of foodprocessing waste water in California’s Central Valley. J. Environ.
Qual. 37, 43–57 (2008)
218. Gérard, F., Tinsley, M., Mayer, K.U.: Preferential flow revealed
by hydrologic modeling based on predicted hydraulic properties
and intensive water content monitoring. Soil Sci. Soc. Am. J. 68,
1526–1538 (2004)
219. Lichtner, P.C., Hammond, G.E., Lu, C., Karra, S., Bisht, G.,
Andre, B., Mills, R.T., Kumar, J.: PFLOTRAN User manual:
A Massively Parallel Reactive Flow and Transport Model for
Describing Surface and Subsurface Processes (2013)
220. Goode, D.J.: Direct simulation of groundwater age. Water
Resour. Res. 32, 289–296 (1996)
221. Hammond, G.E., Lichtner, P.C.: Field-scale model for thenatural attenuation of uranium at the Hanford 300 area
using high performance computing. Wat. Res. Res 46, 1–
31 (2010)
222. Lichtner, P.C., Karra, S.: Modeling multiscale-multiphasemulticomponent reactive flows in porous media: Application
to CO2 sequestration and enhanced geothermal energy using
PFLOTRAN. In: Al-Khoury, R., Bundschuh, J. (eds.) Computational Models for CO2 Geo-sequestration & Compressed Air
Energy Storage, pp. 81–136. CRC Press (2014)
223. Karra, S., Painter, S.L., Lichtner, P.C.: Three-phase numericalmodel for subsurface hydrology in permafrost-affected regions.
Cryosphere Discuss. 8, 149–185 (2014)
224. de Vries, L.M., Molinero, J., Ebrahimi, H., Svensson, U.,
Lichtner, P.: High performance reactive transport simulation of
hyperalkaline plume migration in fractured rocks. Mineral. Mag.
77, 982 (2013)
225. Navarre-Sitchler, A.K., Maxwell, R.M., Siirila, E.R.,
Hammond, G.E., Lichtner, P.C.: Elucidating geochemical
response of shallow heterogeneous aquifers to CO2 leakage
using high-performance computing: Implications for monitoring of CO2 sequestration. Adv. Water Resour. 53, 45–
55 (2013)
226. Steefel, C.I., DePaolo, D., Lichtner, P.C.: Reactive transportmodeling: An essential tool and a new research approach for the
Earth sciences. Earth Planet. Sci. Lett. 240, 539–558 (2005)

478
227. Gibson, B.D., Amos, R.T., Blowes, D.W.: S-34/S-32 fractionation during sulfate reduction in groundwater treatment systems: Reactive transport modeling. Env. Sci. Technol. 45, 2863–
2870 (2011)
228. Greskowiak, J., Prommer, J.H., Liu, C., Post, V.E.A., Ma,
R., Zheng, C., Zheng, C., Zachara, J.M.: Comparison ofparameter sensitivities between a laboratory and field-scaleComput Geosci (2015) 19:445–478
model of uranium transport in a dual domain, distributed
rate reactive system. Water Resour. Res. 46, 9: W09509(2010)
229. Doherty, J.: PEST: A unique computer program for modelindependent parameter optimisation. Water Down Under
94: Groundwater/Surface Hydrology Common Interest Papers;
Preprints of Papers, p. 551 (1994)