Published online March 8, 2006

Variably Saturated Reactive Transport of Arsenic in Heap-Leach Facilities
David L. Decker,* Jirka Šimůnek, Scott W. Tyler, Charalambos Papelis, and Mark J. Logsdon
eralized zones associated with the finely disseminatedgold-bearing deposits of eastern and northeastern Nevada and western Utah. The co-occurrence of As withgold in these deposits is significant enough that As is arecommended gold-ore prospecting indicator. The resurgence in gold mining within the Great Basin since thediscovery and initial development of the Carlin Trend has
resulted in the most productive gold-producing region in
North America (Driesner and Coyner, 2000). The largescale mining techniques used on the Carlin Trend include
the use of heap-leach mining technologies wherein thegold is extracted by leaching large volumes of rock thathave been mined using open-pit methods and stacked
on an impervious plastic liner. The low-grade gold ore,coupled with the large scale of many modern gold mines,
produces heap-leach facilities that can commonly contain 100 million tonnes of rock. The gold within the ore isleached with a sodium-cyanide solution. A byproduct ofleaching is the dissolution of other metal and metalloidbearing mineral phases. The dissolution of these minerals
in turn results in significant aqueous concentrations ofundesired constituents such as As compounds. Once thegold is extracted from the heap-leach facility, effective
methods of either reducing the concentration or the fluxof dissolved metals from the heap-leach facility must be
designed and implemented to satisfy regulatory criteria.
The geological and mineralogical variability between
individual ore deposits, coupled with the large range inhydrologic regimes, necessitates the development of tools
to assist in the design of effective heap-leach closure
methodologies. While the Carlin Trend is not located nearlarge metropolitan areas, it does encompass several largehydrographic basins, which together represent a significant future groundwater and surface water resource. The
arid climate, coupled with unprecedented western population growth, adds to the importance of the water resources in these areas.
In an effort to understand the potential for As impacts
on subsurface water supplies from mining-related activities, we examined the processes controlling As transportin heap-leach rock from mines located within the Carlin
Trend. The physics of fluid transport in these highly
heterogeneous, variably saturated environments werestudied with a series of large-scale column experiments
using conservative tracers and time domain reflectometry to characterize advective–dispersive transport behavior (Decker and Tyler, 1999a). In a follow-up study,
the redox and pH-dependent sorption behavior of As
was examined for two distinct mine locations within theCarlin Trend in an effort to understand the mechanismsof As release and transport in a mine setting. As part ofthis work, empirically derived pH-dependent isotherm
formulations were developed with the intention ofincorporating pH-dependent sorption behavior into avariably saturated reactive flow and transport model

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACTA variably saturated reactive transport software package,
UNSATCHEM, has been modified for the purposes of simulating
As oxyanion transport in gold mine heap-leach facilities. Gold-ore
deposits in the Carlin Trend, located in northeastern Nevada andwestern Utah, are characterized in part by the occurrence of Asbearing minerals. The large-scale mining methodology used on theCarlin Trend includes the use of large-scale heap-leach facilities inwhich an extractive solution is used to dissolve the gold from the ore.
This solution also is capable of dissolving other minerals, including
those that contain As. While these facilities are generally environmentally benign compared with sulfide waste rock dumps associatedwith other mining provinces around the world, controlling the flux ofAs from the heap-leach facility after mining has ceased is of interest to
the mining and regulatory communities. We present three examples
that use the As transport capable version of UNSATCHEM to simulate the effect of three commonly applied heap-leach facility closure
methodologies. The results demonstrate the importance of including
geochemistry in transport simulations and illustrate the highly nonlinear behavior that is produced from a coupled flow and reactive
geochemistry numerical solver.

T

HE SUBSURFACE TRANSPORT of metals in variably satu-

rated media is of interest because of the increasing
demands on subsurface water supplies for drinking, industrial processes, and agriculture. Particular interest inarsenic, selenium, and chromium is driven by the biological toxicity of compounds formed with these elements. Arsenic toxicity has induced a reduction in thefederal drinking water standard from 50 to 10 mg L21 in
the United States. The financial impacts of treating
drinking water for As has led to renewed interest inunderstanding the processes controlling the fate andtransport of As compounds in the subsurface.
Arsenic is prevalent in the known groundwater
resources within the Great Basin of Utah, Nevada, andCalifornia. Much of this As is thought to be the result of
dissolution of As-bearing, broadly distributed, primaryrock and mineral formations (Welch et al., 2000). Arsenic
bearing minerals are also present in the form of minD.L. Decker, Desert Research Institute, 2215 Raggio Parkway, Reno,
NV 89512 USA; J. Šimůnek, Dep. of Environmental Sciences, A135
Bourns Hall, Univ. of California, 900 University Ave., Riverside, CA
92521 USA; S.W. Tyler, Dep. of Natural Resources and Environmental Sciences and Dep. of Geological Sciences and Engineering,MS 175, Univ. of Nevada, Reno, NV 89557-0180 USA; C. Papelis,
Desert Research Institute, 755 E. Flamingo Road, Las Vegas, NV
89119 and Water Resources Management Graduate Program, Univ.
of Nevada, 4505 Maryland Pkwy., Las Vegas, NV 89154-4029 USA;
M.J. Logsdon, Geochimica Inc., 9045 Soquel Drive, Suite 2, Aptos,
CA 95003 USA. Received 14 Mar. 2005. *Corresponding author(dave.decker@dri.edu).
Published in Vadose Zone Journal 5:430–444 (2006).
Original Researchdoi:10.2136/vzj2005.0039ª Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

430

431

www.vadosezonejournal.org

for As (Decker et al., 2006). This paper continues
these efforts with the inclusion of As geochemistry in
UNSATCHEM, a public domain, variably saturated, reactive transport solver.

and

5 KK K (h) hh $, 00

[4]

m 2
Kr 5 Se/2 [1 2 (1 2 S1/m
e ) ]

[5]

K(h) 5srs

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

where

MODELING REACTIVE ARSENICOXYANION TRANSPORT

1

There are several reactive transport codes available
that address a wide range of flow, transport, and chemical
modeling capabilities. The codes that are generally applicable to the reactive flow problem of heap-leach facilities are: TOUGH-REACT (Pruess, 1991), MULTIFLO(Lichtner and Seth, 1996), UNSATCHEM (Šimůnek andSuarez, 1994; Šimůnek et al., 1996), FEHM (Viswanathan
et al., 1998), MIN3P (Mayer et al., 2002), and FLOTRAN(Lichtner et al., 2004). The model chosen for modification,
UNSATCHEM, was selected because of the robust flow
and transport solver HYDRUS-1D (Šimůnek et al., 1998),
the inclusion of carbonate geochemistry with the associated CO2 transport subroutines (Šimůnek and Suarez,
1993), and the public-domain accessibility of the sourcecode to enable the addition of As aqueous chemistry,
sorption, and mineral dissolution processes.

Variably Saturated FlowSe 5

uw 2 ur
us 2 ur[6]

m 5 1 2 1/n, n . 1[7]

in which ur and us denote the residual and saturated
water contents, respectively, Se is the effective saturation, and Kr and Ks are the relative and saturated hydraulic conductivities, respectively. The parameters aand n are empirical fitting parameters (van Genuchten,1980). Richards’ equation and van Genuchten’s soil hydraulic functions are applicable to porous media
structures where Darcy’s Law is valid, that is, when a
low Reynold’s number characterizes fluid flow. However, the fluid flow regime in mine heap-leach facilities
may not necessarily exhibit Darcian behavior due tomacropore flow through the heterogeneous rock (Diodato and Parizek, 1994; Mallants et al., 1997; Nichol,
2002; Nichol et al., 2003). While efforts to develop
empirical modifications to the unsaturated hydraulicconductivity function given above are reported in theliterature, this analysis assumes Darcy’s Law to be applicable to the heap-leach facility flow regime. An overview of the applicability of this approach to the specificproblem of heap-leach flow modeling is given in Decker
and Tyler (1999b).

A successful variably saturated flow and reactive
transport model must include a mathematical approachdescribing fluid flow in variably saturated porous media. While the mathematics of variably saturated flow
are given in rigorous detail elsewhere (e.g., Hillel, 1998),
the flow equations are provided here as written into
UNSATCHEM (Šimůnek et al., 1996, 1998). The governing equation for the isothermal Darcian flow of waterin a variably saturated rigid porous medium using theassumptions that the air phase plays an insignificant rolein the liquid flow process and that water flow due tothermal gradients can be neglected is given by Richards’
equation as

The partial differential equation describing one-dimensional advective–dispersive solute transport under transient water flow conditions in partially saturated porous
media is taken as (Šimůnek et al., 1996)]uw]K
5 =  [K(h)=h] 1
2S
]t
]z

]uw Ck
]Sk
]Mk
]
]Ck
1r
1r
5
2 qw C k
uw D
]t
]t
]t
]z
]z

[1]

where uw is the volumetric water content, h is the matricpotential or pressure head, K(h) is the hydraulicconductivity function dependent on matric potential, z
is the spatial coordinate, t is time, and S is the source–
sink term. The present modeling effort is limited to one
dimension; thus, Richards’ equation can be reduced to


]uw
]
]h
[2]
5
K
11 2S
]t
]z
]z

ð

Þ

The variably saturated hydraulic properties for soil
are often described using van Genuchten’s (1980)
equations for the retention and unsaturated hydraulicconductivity functions as
uw (h) 5

5

ur 1
us

us 2 urn(1 1 |ah| )m

h,0
h$0

[3]

Multicomponent Solute Transportð

k 5 1, 2, …, NcÞ[8]where Ck is the total dissolved concentration of theaqueous species k, Sk is the total sorbed species
concentration of the aqueous component k, Mk is thetotal solid phase concentration of the aqueous component k, r is the bulk density of the medium, D is thedispersion coefficient, qw is the volumetric flux, and Nc is
the number of primary aqueous species. The volumetricflux is calculated with Darcy’s Law
qw 5 2K

ð ]h]z 1 1Þ

[9]

The solution of Eq. [8] requires the specification of an
initial condition and appropriate boundary conditions.
An appropriate initial condition for a heap-leach facility
simulation is a specification of the initial aqueous, sorbed,

432

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

and solid components of the aqueous species k for all Nc
species as

Component typeC(z, 0) 5 C(z)S(z, 0) 5 S(zÞReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Aqueous ionsAqueous complexest50

M(z, 0) 5 M(z)[10]Note that to simplify the notation, the index k in Eq.
[10] and below was dropped. Two types of boundary
conditions are applicable to heap-leach facility transport
problems. The first type of boundary condition, Dirichlet,
prescribes the concentration at boundaries asC(z, t) ¼ Co (z, t) at z 5 0 or z 5 L[11]while the third type of boundary condition, Cauchy, describes the concentration flux as
2uw D

Table 1. Carbonate geochemistry included in original UNSATCHEM
software (Šimůnek et al., 1996).

]c
1 qw C 5 qw Co at z 5 0 or z 5 L
]z

[12]

where Co is the concentration of the incoming fluid and L
is the depth of the soil profile. A zero concentration
gradient may also be prescribed at the bottom of thesoil profile.

GEOCHEMICAL REACTIONSIN UNSATCHEM
UNSATCHEM is a variably saturated reactive transport code that was designed to simulate major ion
transport in agricultural soils (Šimůnek et al., 1996). The
code numerically solves the Richards equation forvariably saturated flow and convection–dispersion-type
equations for heat, gas, and solute transport. The major
chemical variables considered in the original version of
the code are Ca21, Mg21, Na1, K1, SO422, Cl2, alkalinity,
and CO2. The code accounts for temperature-dependentequilibrium chemical reactions between these components through complexation, cation exchange, andprecipitation–dissolution. The mineral phases include
calcite, dolomite, gypsum, nesquehonite, hydromagnesite, and sepiolite. The precipitation–dissolution reactions for calcite can be modeled using either equilibriumor kinetic expressions, while the dissolution of dolomite
is always considered as a kinetic process. The modelis capable of simulating saline water movement, andmodified Debye-Hückel and Pitzer expressions are
utilized to calculate single-ion activity coefficients.
UNSATCHEM represents the geochemistry of the carbonate system including the precipitation–dissolution
reactions of the mineral phases described above, and theaqueous speciation of carbonate, bicarbonate, and carbonic acid. The reaction sequences in the carbonate
geochemical system are dependent on CO2 concentration. Therefore, UNSATCHEM models the transport ofCO2 with a mass-balance relationship that includes diffusive and convective transport in both gas and fluid
phases. A compilation of the aqueous complexes andmineral species included in UNSATCHEM is given in
Table 1, while thermodynamic data, kinetic reaction
constants, and other details are provided in Šimůnek
et al. (1996).

Mineral species
Sorbed species
CO2–H2O species

Components
Ca21, Mg21, Na1, K1, SO422, Cl2
1
1
CaCO°,
3 CaHCO3 , CaSO°,
4 MgCO°,
3 MgHCO3 ,2
2
2
MgSO°,
4 NaCO3 , NaHCO°,
3 NaSO4 , KSO4 ,2
22
H4SiO°,
4 H3SiO4 , H2SiO4
CaCO3 , CaSO4  2H2O, MgCO3  3H2O,
Mg5(CO3)4(OH)2  4H2O, Mg2Si3O7.5(OH)3 
H2O, CaMg(CO3)2
Ca, Mg, Na, K
PCO2, H2CO3, CO322, H1, OH2, H2O

MODIFICATIONS TO UNSATCHEM
The UNSATCHEM code has been modified in thiswork to include the geochemistry associated with As in
carbonate rock. The addition of As transport to theoriginal geochemistry included in UNSATCHEM required the inclusion of several new geochemical reaction
pathways. These included As aqueous complexation
reactions for two As oxidation states, pH- and oxygendependent mineral dissolution and precipitation reactions, and pH- and redox-dependent As sorption(Decker et al., 2006).
The two oxidation states for As that are stable in soils
and the subsurface in general are As(III) [As31] andAs(V) [As51] (Welch and Stollenwerk, 2003). These cations form aqueous oxyanion complexes with a net negative charge. The sorption behavior of these oxyanion
complexes to heap-leach facility rock is mineral, redox,
and pH dependent to such a degree that simple sorption
isotherm approaches are insufficient to adequately describe this behavior (Decker et al., 2006). The mineralsource of As is dependent on the specific mineralogicalcomposition of each mine site. However, a common Asbearing mineral phase found in Carlin Trend ore bodies is
As-bearing pyrite. UNSATCHEM has been modified toinclude three forms of this mineral including pyrite (FeS2),
arsenopyrite (FeAsS), and arsenian pyrite (FeAsxS22x),
where x is the fraction of As. The dissolution of pyrite, orAs-bearing pyrite, releases aqueous Fe. UNSATCHEM
has been modified to include the precipitation reaction foriron hydroxide to remove the aqueous Fe mass fromsolution. In addition, the dissolution reactions of pyrite,
As-bearing pyrite, and the precipitation reaction of ironhydroxide, require O2, and UNSATCHEM has been modified to include the transport of O2, using the same approach taken for CO2 transport in the original version of
the code.
Aqueous As oxyanion complexes exhibit significant
sorption behavior that is dependent on both redox andpH conditions. The redox and pH-dependent sorption
behavior was determined for two samples of Carlin
Trend heap facility rock. Two pH-dependent sorption
isotherm approaches were used to empirically describe
the sorption behavior as described in detail in Decker
et al. (2006), and both of these approaches have been
added to UNSATCHEM. A detailed description of themodifications to UNSATCHEM that were undertaken
to provide As transport capability are given in subsequent sections.

433

www.vadosezonejournal.org

Table 3. Arsenic aqueous complex thermodynamic values(Appelo and Postma, 1993; Welch and Stollenwerk, 2003).

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Aqueous Complexation Reactions
The set of As aqueous complexes that has been added
to the geochemical database built into UNSATCHEM
is given in Table 2; the values of the thermodynamic
parameters associated with these reactions are given in
Table 3. Since these complexes affect the overall massbalance, two additional equations for the two As oxidation states considered are
32
As(III)T 5 [H2 AsO23 ] 1 [HAsO22
3 ] 1 [AsO3 ]1 [H3 AsOo3 ]

[13]

32
As(V)T 5 [H2 AsO24 ] 1 [HAsO22
4 ] 1 [AsO4 ]1 [H3 AsOo4 ][14]

where variables with the subscript T represent the total
analytical concentrations in solution of the particularspecies, and where the brackets refer to molalities (molkg21). In addition to these mass-balance equations, theoverall charge-balance equation, including As complexes, is expressed as
2[Ca21 ] 1 2[Mg21 ] 1 [Na1 ] 1 [K1 ] 1 [CaHCO13 ]2
22
1 [MgHCO13 ] 1 [H1 ] 2 2[CO223 ] 2 [HCO3 ] 2 2[SO4 ]2 [Cl2 ] 2 [OH2 ] 2 [NaCO23 ] 2 [NaSO24 ] 2 [KSO24 ]32
2
2 [H2 AsO23 ] 2 2[HAsO22
3 ] 2 3[AsO3 ] 2 [H2 AsO4 ]32
2 2[HAsO22
4 ] 2 3[AsO4 ] 5 0[15]Precipitation–Dissolution Reactions
Inorganic pyrite dissolution is described by a series ofreactions that encompass a wide pH range. These reactions are commonly written as (Langmuir, 1997)
7
2

1
FeS2 1 O2 1 H2 O W Fe21 1 2SO22
4 1 2HFeS2 1 14Fe

31

1 8H2 O W 15Fe

21

1 2SO

22
4

[16]
1

1 16H

[17]
21

Fe

1
1
1 O2 1 H1 W Fe31 1 H2 O
4
2

[18]

The Carlin Trend heaps represent a somewhat uniqueenvironment to describe pyrite oxidation due to thepresence of significant amounts of carbonate minerals
relative to the low mass percentages of sulfide minerals.
The equilibrium pH for the solution in equilibrium withboth gypsum and calcite minerals varies between 7.8 andTable 2. Components added to UNSATCHEM to simulate As
geochemistry.
Component type

Components

Aqueous ionsAqueouscomplexes
Mineral species

As31, As51, Fe21
22
32
2
H2AsO23 , H3AsO°,
3 HAsO3 , AsO3 , H2AsO4 ,22
32
HAsO4 , AsO4 , H3AsO°4
Fe(OH)3
Pyrite minerals:
Pyrite–FeS2
Arsenopyrite–FeAsS
Arsenian Pyrite–FeAsxS2-x (0 , X , 1)
As(III), As(V)PO2

Sorbed species
O2 species

Complexation Reaction

Log K

Dhr (kcal mol21)

1
H3As(III)O3 W H2As(III)O2
3 1 H
H3As(III)O3 W HAs(III)O322 1 2H1
H3As(III)O3 W As(III)O332 1 3H1
1
H3As(V)O4 W H2As(V)O2
4 1 H
H3As(V)O4 W HAs(V)O422 1 2H1
H3As(V)O4 W As(V)O432 1 3H129.228221.33234.744
22.24329.1
220.597

6.56
14.199
20.25
21.69
20.92
3.43

8.3 (Krauskopf and Bird, 1995). Because of this highpH environment, the dissolution of pyrite by Fe31 is
minimized as a result of the rapid precipitation of Fe31
to Fe(OH)3 by the following reaction:
1
4

5
2

Fe31 1 O2 1 H2 O!Fe(OH)3 1 2H1[19]

which, in the flow-through environment of the heapleach facility, is assumed to be irreversible. Therefore,
the reactions given in Eq. [17] and [18] are ignored, andthe primary dissolution reaction for pyrite in Carlin
Trend type deposits, as given in Eq. [16], is controlled by
the availability of O2.
The acidity generated from the oxidation of pyrite(Eq. [16]) and the hydrolysis of ferric iron (Eq. [18]) is
buffered by the dissolution of calcite and the precipitation of gypsum. In Carlin Trend heap-leach facilities, theamount of carbonate mass far exceeds the total mass of
sulfides in the heap, and a high pH environment should
prevail at all times.
Simulating the oxidation of pyrite and the precipitation of ferric hydroxide requires the inclusion of kinetic
versions of these reactions. An experimental rate law
for Eq. [16] was determined in a stirred reactor in acontrolled laboratory environment, and is given by(Williamson and Rimstidt, 1994) asd(FeS2 )(mO2 )0:5
5 1028:10
dt(mH1 )0:11[20]

where the rate of pyrite dissolution is expressed in moles
per square meter per second. This rate law is particularly
useful to the present modeling effort, as it depends only
on pO2 and pH. An alternative rate law for pyrite
oxidation in moist air proposed by Jerz and Rimstidt(2004) is of interest because it is conceptually applicable
to a variably saturated environment. The Jerz and
Rimstidt (2004) experiment was performed under stagnant flow conditions and with initially unoxidized pyrite. These experimental conditions are not analogous
to a field-scale heap-leach facility, and therefore theWilliamson and Rimstidt (1994) rate expression is incorporated into UNSATCHEM.
At high pH, Fe(OH)3 is relatively insoluble, andprecipitates by Eq. [19]. An empirical rate law for thisreaction is given by
d[Fe(OH)3 ](Fe21 )
PO
5 21025:98
dt(H1 )2

2[21]

for which the rate of iron hydroxide precipitation is expressed in moles per second at 208C (Langmuir, 1997).

434

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

The occurrence of As as a minor constituent in pyrite
is a possibility and has been documented in precious
metals heaps around the world (Jambor, 1994). The oxidation of arsenian pyrite can be written as

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

FeAsx S22x 1

ð 72 22xÞO 1 H O!Fe 1 (22x)SO
12

2

22
4

2

1 2H1 1xAsaq[22]

where x is the fraction of As present in pyrite. Note thatsetting x 5 0 results in the pyrite oxidation reaction
given in Eq. [16], and setting x 5 1 results in a representation for the oxidation of the mineral arsenopyrite,
FeAsS. The generalized reaction in Eq. [22] is included
in UNSATCHEM as a master equation for pyrite, arsenopyrite, and arsenian pyrite oxidation and is invokedby specifying the As fraction as an input parameter. The
dissolution rate for this reaction is assumed to proceed
under the same rate law as pure-phase pyrite as given in
Eq. [20].

Arsenic Sorption
The partitioning of As mass between the aqueous andsorbed phases is accomplished using the two pHdependent isotherm expressions described in Decker
et al. (2006). The first expression is a Sips-type, orLangmuir–Freundlich, isotherm that is pH dependent.
1(aS Ce )n
1 1 de2epH 1 1 (aS Ce )nS

qe 5 qm

S[23]This expression encompasses two analytically measured parameters, Ce and pH, and one calculated parameter, qe. There are three isotherm-fitting parametersfrom the Sips formulation, qm, ns, and as, and incorporating pH dependence adds two additional fitting parameters, d and e. The simulations included in this study
simplified this model to a Langmuir form by assigningns 5 1. In addition, a Keren-type pH-dependent isothermis also incorporated into UNSATCHEM
KL Ce
1 1 KL Ce
A 1 B 3 10(pH214)
KL 5
1 1 C 3 10(pH214) 1 D 3 [10(pH214) ]2
qe 5 qm[24]

where A, B, C, and D are empirically determined
constants.

Oxygen Transport
A critical component to modeling As sorption is
including redox chemistry to enable speciating As in
both the aqueous and sorbed phases. The experimental
sorption results given in Decker et al. (2006) demonstrate that the sorption behavior of As(III) and As(V)exhibits significant pH dependence. The speciation of
As(III) and As(V) is controlled by redox chemistry.
Therefore, a reactive transport model for As should be
inclusive of redox processes. The inclusion of redox
chemistry is particularly important to a realistic inclusion of both As(III) and As(V) sorption behavior.

Water chemistry data for the two facilities were neither
available nor obtainable due to the anonymous nature of
the donations of rock material from the mining industry.
Therefore, a surrogate redox chemistry approach wastaken. An approximation of the redox potential can bemade through dissolved oxygen concentration in that
reducing conditions are characterized by low oxygen concentrations, and oxidizing conditions are characterized
by relatively higher oxygen concentrations. This generalized redox chemistry approach is incorporated into
UNSATCHEM for the singular purpose of controllingthe oxidation state of As. In practice, this is accomplishedby specifying the critical oxygen concentration that defines the threshold between oxidizing and reducing conditions as an input variable. This variable in turn is used todetermine the speciation of As on a node-by-node basis
at each time step. Arsenic speciation is allowed to proceed through this oxygen surrogate in a reversible manner. The use of oxygen and iron as a redox surrogate tocontrol As redox chemistry presumes that the geochemistry associated with Fe, O2, and As is coupled. In thepresence of sulfide mineralogy, this assumption isprobably reasonable since the coupled interaction of O2
and Fe drives the redox chemistry. However, in theabsence of sulfide mineralogy, the use of the partialpressure for O2 may not be a satisfactory redox surrogate
for As, and explicit As redox chemistry may need to beincorporated into UNSATCHEM to adequately represent this type of domain. The Carlin Trend ore deposits
are associated with sulfide mineralogy, and therefore theassumption is made that the use of O2 is a reasonablesurrogate representation of the complex redox chemistry
for As.
The oxygen transport modeling approach is similar tothat used for carbon dioxide transport (Šimůnek andSuarez, 1993). The gas transport mass conservation
equation is]CT]
5 2 (Jda 1 Jdw 1 Jca 1 Jcw ) 1 P
]t
]z

[25]

where Jda describes the O2 flux caused by diffusion in thegas phase, Jdw is the O2 flux caused by dispersion in the
dissolved phase, Jca is the O2 flux caused by convectionin the gas phase, and Jcw is the O2 flux caused byconvection in the dissolved phase. The term CT is thetotal volumetric concentration of O2, and P is a source–
sink term. These terms are defined asJda 5 2ua Da]Ca]z

Jdw 5 2uw Dw ]Cw]z
Jca 5 qa Ca
Jcw 5 qw Cw
CT 5 Ca ua 1 Cw ua[26]

where Cw and Ca are the volumetric concentrations of
O2 in the dissolved and gas phases, respectively, Da isthe effective soil matrix diffusion coefficient of O2 in the435

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

gas phase, Dw is the effective soil matrix dispersion
coefficient of O2 in the dissolved phase, qa is the soil air
flux, qw is the soil water flux, and ua is the volumetric air
content. Since UNSATCHEM does not consider coupled water and air movement, the flux of air, qa, is
unknown, and thus must be approximated somehow
using additional assumptions. One possibility is to assume that the advective transport of O2 in response to
the total pressure gradient is negligible compared with
O2 diffusion, and therefore to assume a stagnant gas
phase in which only diffusive transport occurs (qa 5 0).
Another possibility is to assume that because of themuch lower viscosity of air in comparison with water, arelatively small pressure gradient will lead to significant
gas flow. Changes in the total volume of water in the soil
profile caused by water flow are thus assumed to be
immediately matched by corresponding changes in thegas volume, leading to air fluxes. This simple model thus
takes into account “rainfall effects” on soil aeration,
which was estimated to contribute about 7 to 9% of normal aeration (Jury and Horton, 2004). Jury and Horton(2004) also stated that other factors, such as wind, barometric pressure changes, and temperature effects, do
not contribute to more than about 1% of normal aeration, and thus may be neglected in the majority of applications. Although the modified UNSATCHEM model
can consider simplified gas convection, the examples
given in this paper are run under steady-state, constant
fluid flux conditions. Therefore, the contribution to thetotal oxygen flux from convection in the air phaseis neglected.
Sulfide waste rock dumps in other mining districts
around the world can contain sulfide mass fractions of 2
to 30% (Ritchie, 1994). In some of these cases, convective oxygen flux is known to substantially contributeto the overall oxygen mass budget as a result of theexothermic heat of reaction produced from the dissolution of sulfides (Ritchie, 1994). However, it is unlikelythat convective O2 flux is a major feature of the Carlin
Trend heap-leach facilities for three reasons. First, thesulfide concentration of the ores, at 2.5% by weight, is so
low that significant heat generation is not likely to occur.
In fact, following the analysis given in Ritchie (1994), it ishighly likely that a temperature rise of ,28C could be
expected. Second, the heaps are constructed in a waythat limits the likelihood of vertically continuous highconductivity features such as commonly are found in
end-dumped waste-rock piles that provide preferential
flowpaths for gas transport from the toe of a facility to
the upper reaches of the dump. Third, the continuous
flux of downward-directed fluid applied at ambient
temperature at the surface would moderate the thermal
regime in such heap-leach facilities, thus damping the
upward convective flow of air and heat. Therefore,
neither this thermally driven convective flow nor theheat generated during pyrite dissolution is included inthis version of UNSATCHEM.
The term P in Eq. [25] represents a sink term foroxygen resulting from the consumption of oxygen bychemical reactions such as those given in Eq. [16], [19],and [22]. No sources of oxygen other than the at-mosphere are included in this formulation. Henry’s Law
governs the distribution of gas between aqueous andgaseous phases
Cw 5 KO RTCa

[27]

2

23

where KO is the Henry’s Law constant (1.26 3 10 mol
bar21 (Langmuir, 1997)), R is the universal gas constant(8.314 kg m2 s22 K21 mol21), and T is the absolutetemperature (K). The dissolution–precipitation reactions given in Eq. [16], [19], and [22] include water, and it
is assumed that these reactions occur in the presence ofa water film. Therefore, Henry’s Law is used to movegaseous oxygen to the aqueous phase such that the reactions in Eq. [16], [19], and [22] can occur.
The effective dispersion coefficients for O2 in soils arecalculated from the molecular diffusion coefficients foroxygen in free air and the dispersion coefficient in water.
UNSATCHEM calculates an effective diffusion coefficient for gas transport by multiplying the molecular
diffusion coefficient by a tortuosity factor. The resulting
effective dispersion coefficients, Dw, and diffusion coefficients, Da, are then written as2

Dw 5 Dws Hw 1 lw

7

qw
u /3
qw
5 Dws w2 1 lw
uw
uw
us

| |

| |

[28]

7

Da 5 Das Ha 5 Das

ua/3
u2s

[29]

where Dws 5 2.10 3 1029 m2 s21 and Das 5 2.01 3 1025
m2 s21 are the molecular diffusion coefficients for O2 inthe dissolved and gas phases from Glinski and Stepniewski (1985), Ha and Hw are the tortuosity factors inboth phases, and lw is the dispersivity in the dissolved
phase (Šimůnek et al., 1996).

Numerical Solution Strategy
The numerical solution of Richards’ equation isaccomplished with the use of a mass-lumped, finiteelement method that is reduced to a finite-difference
scheme (Šimůnek et al., 1996). This finite-difference
numerical approach is combined with a Picard iterative
solution technique to solve the time-derivative in
Richards’ equation. The numerical solutions to thesolute, heat, and gas transport equations are accomplished with the Galerkin finite element method. The
flow and transport models are inclusive of a full suiteof boundary conditions including specified head andconcentration (Dirichlet), or specified flux (Neumann).
An iterative solution is used to solve the chemicalsystem and is discussed in the UNSATCHEM manual(Šimůnek et al., 1996). The iteration process must be
accomplished at each node, for each time step, and between the solute transport and equilibrium componentsof the model. The solution process at each time step
proceeds by iteratively solving the flow equation andsolving the transport equations for temperature, CO2,
and then multicomponent solute transport. The geochemical module is invoked at each time step, for eachnode, to redistribute mass between solid, aqueous, andsorbed phases.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

436

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

The governing solute transport equation (Eq. [8]) includes three time-derivative components: total aqueousconcentration, sorbed concentration, and mineral-phase
concentration. The solute transport equation is highly
nonlinear because of these terms and must be solved
iteratively. Iteration continues until the differences between the calculated aqueous concentrations for the new
iteration and the previous iteration are less than a prescribed concentration tolerance.

MODEL DEMONSTRATION
Three examples of variably saturated transport inheap-leach systems are presented to illustrate theimportance of incorporating reactive geochemistry intotransport simulations that attempt to estimate As fluxes
at long time-scales. Because the model examples deal
with variably saturated media, gas-phase transport ofoxygen and carbon dioxide is an important processcontrolling redox chemistry, pH, and mineral dissolution.
Mineral dissolution and precipitation reactions control
redox chemistry, which in turn directly control As speciation. Since As sorption behavior is species and pH
dependent, the overall reactive transport of As in porous
media is a highly complex interactive process that isdependent on mineralogy, geochemistry, and the physics
of flow in porous media. The examples are intended todemonstrate that a more representative model of As
transport is obtained through the integration of flow
physics, geochemistry, and reactive transport than is possible with simple conservative or linear isotherm approaches. While these examples are not representationsof a specific heap-leach facility, they provide a range of
behavior that might be expected for three possible
closure scenarios that are believed to be applicable toCarlin Trend ores. These examples include the following
heap facility closure scenarios: no surface treatment, anengineered surface treatment to reduce the flux of infiltrating water, and a more sophisticated cover that reduces infiltrating water and gas flux.

Problem Definition
The problem definition for these three examples isshown in Fig. 1. The 30-m-tall, one-dimensional profiles
are the same in all three cases, with the exception of a
57-cm-thick clay-loam cover that is built into the 30-m
profile in Example 3 (Fig. 1c). The hydraulic properties
for the heap material are reported in Table 4, and aretaken from a study conducted on heap facility material
from the Carlin Trend (Decker and Tyler, 1999a). The
hydraulic properties for the clay-loam cover in the third
example are taken from the soils database included withUNSATCHEM and HYDRUS-1D. The steady-state,
upper flux boundary conditions are used to simplify theinterpretation of the reactive transport results. The
lower flux boundary condition is free drainage, or zeropressure head gradient. An alternative lower boundary
condition that can be used in UNSATCHEM is a seepage face, which allows for the development of a saturated zone. Incorporating this boundary condition into

Fig. 1. Schematics of the three model examples.

a heap facility model is appropriate when the effects of
the heap liner are of interest. For the examples shown
here, it is assumed that the location of the heap liner isdeeper than 30 m.
The geochemical conditions within the heap and thecovers are identical for all three examples and are summarized in Tables 5, 6, and 7. Table 5 includes the initial
aqueous geochemical conditions within the interstitialwater of the heap as well as the geochemistry of therecharging water from precipitation. The initial total As
mass (Table 6) is identical for all examples and both
isotherm formulation types. However, the initial aqueous As concentrations (Table 5) are different for the
Sips- and Keren-type isotherms as a result of the difference in fit between each of the isotherm types and theexperimental sorption data. The Sips- and Keren-type
isotherm parameters used in the examples are given in
Table 7 and are derived from As sorption experiments
conducted on heap facility rock (Decker et al., 2006).
The remaining minerals listed in Table 6 are intended toTable 4. Soil hydraulic parameters for the heap material and gasflux cover.
Parameter
ur
us
a
n
Ks
rb
DL

Unit

cm21
cm d21
g cm23cm

Heap material,
all examples

Gas-flux cover, Example 3

0.07
0.4
0.05
3.2
114.57
1.5
100

0.095
0.41
0.019
1.31
6.24
1.5
100

437

www.vadosezonejournal.org

Table 5. Initial and boundary liquid-phase concentrations andother geochemistry parameters.
Units

Initial heap
interstitial water

Recharge watermmol L21
mmol L21
mmol L21
mmol L21
mmol L21
mmol L21
mmol L21
mmol L21 (mg L21)

20
0.225
5
0.05
37.5
3
2
4.731 (354)

0.5
0.025
0.05
0.05
0.05
0.025
0.05
0

mmol L21 (mg L21)

3.330 (250)0mmol L21 (mg L21)

0 (0)
8.76
0
0

0 (0)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Parameter
Ca12
Mg12
Na11
K11
Alkalinity
SO22
4
Cl21
As(III) Sips
Isotherm
As(III) Keren
Isotherm
As(V)
pH
pO2
Tracer

atm

0.2
1.0

be representative of an average, carbonate-hosted,sulfide-bearing ore. In particular, the 20% by weight
carbonate and 2.5% by weight pyrite are characteristics
of a carbonate-hosted ore body. The inclusion of As in
the pyrite mass is meant to represent the mineral origin
of As that has been reported in several heap facilities.
The gas transport boundary conditions for all three
examples are modeled as constant partial pressures in
the atmosphere, with the top of the heap profile
represented with values of 0.2 atm for O2 and 10-3.5 atm
for CO2 and the bottom of the profile set to 0 atm forboth gases. Thus, gas transport into the heap profile islimited to the top surface.
The geographical region encompassing the Carlin
Trend is characterized by an arid climate with low
annual precipitation with correspondingly high evapotranspiration rates that can be several times larger than
the average annual precipitation. However, the large
range in heap facility site elevation complicates thisbroad categorization, with the possibility of significant
snowmelt recharge occurring at high-elevation sites in
the winter months (Kampf et al., 2002). Controlling recharge flux into heap facilities is a key component oflong-term management of both the remaining mineral
resource in the facility, and in ensuring effective environmental stewardship.
Mine operations are required to close heap-leachfacilities such that a number of environmental regulations are met. This is typically accomplished through aregimen of rinsing leachate fluids from the heap withfresh water and establishing a long-term fluids management plan to control water discharge for the facility. The

costs associated with the fluids management plan are ofsuch magnitude that designing an effective closure planthat includes understanding the impacts of engineeringsolutions to limit the precipitation flux into the top of theheap facility on the chemistry of the discharge water is
key to minimizing costs and maximizing the effectivenessof the management plan. The three examples presented
illustrate the long-term effects on As concentration andmass flux from precipitation and gas fluxes at the heap
facility surface.

Example 1: Uncovered Heap-Leach Facility
The first example represents an uncovered heap facility with a relatively large steady-state recharge flux
rate of 300 mm yr21 that is representative of a highelevation environment with significant snowmelt recharge. The calculated steady-state water content is
0.097, which corresponds to a pore volume of water of
291 cm3 and a time to infiltrate a single pore volume of
9.7 yr (Fig. 2h). The initial pH of the profile is set at 8.7 torepresent a heap that is recovering from a period of gold
production wherein the pH is maintained at very highlevels (in excess of pH 10) to prevent the loss of cyanide(Bartlett, 1992) (Fig. 2e). This pH rapidly decreases tolower pH values due to the influence of carbonate mineral dissolution, sulfide mineral dissolution, and gypsum
precipitation reactions. As arsenian pyrite is dissolved(Fig. 2b), the acid generated is buffered through thedissolution of calcite (Fig. 2a). The dissolution of arsenian pyrite is limited by the rate of oxygen transport
through the relationship given in Eq. [20]. Oxygen transport (Fig. 2g) in turn is limited by the effective dispersion
coefficient for oxygen. The linear profiles shown in
Fig. 2g are indicative of an arsenian pyrite oxidation ratethat is in excess of the oxygen transport rate. Sulfate is
produced through the dissolution of arsenian pyrite(Fig. 2f) and is precipitated as gypsum (Fig. 2c). The
aqueous iron generated from arsenian pyrite dissolution
is precipitated as iron hydroxide through the kinetic
relationship in Eq. [21] (Fig. 2d).
Arsenic is initially distributed as As(III) between the
aqueous and sorbed phases in a proportion determinedby the Sips-type isotherm, (Eq. [23]), and for the parameter values given in Table 7 for an initial pH of 8.76
(Fig. 2i). The speciation of As from As(III) to As(V)occurs along a moving redox front that is coincident withthe oxygen transport and arsenian pyrite oxidation
fronts. In Fig. 2i and 2j, the presence of the oxidation

Table 6. Initial conditions for mineral- and sorbed-phase concentrations.
Parameter
CalciteGypsumPyriteFe(OH)3
Pyrite surface area
Arsenian pyrite fraction
Sorbed As(III)
Sorbed As(V)
Total As(III) mass(aqueous and sorbed)

Unit
21

21

mol kg (g kg )mol kg21 (g kg21)
mol kg21 (g kg21)
mol kg21 (g kg21)
m2 g21mmol kg21 (mg kg21)
mmol kg21 (mg kg21)
mmol kg21 (mg kg21)

Keren isotherm

Sips isotherm

2 (200)0.005 (0.68)
0.208 (25)0
7.5 3 1028
0.00075
1666.663 (124.8691)
0 (0)1666.667 (124.8693)

2 (200)0.005 (0.68)
0.208 (25)0
7.5 3 1028
0.00075
1666.662 (124.8689)
0 (0)1666.667 (124.8693)

438

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Table 7. pH-dependent isotherm surface model parameters. Units
for the calculated isotherm surface plot are: Ce, mg L21; qe, mg
g21; pH 5 2log[H1].
ModelParameter

As(V)

As(III)

Sips typede
21
qm (mg kgrock )as (L mg21)
ns
A (L mg21)
B (L mg21)
CD
21
qm (mg kgrock )0.000004891
21.404960 000
0.001151
1.0
0.0007962
46.99
548 600
0.12501 200 000

4.925
0.4211
337 000
0.002011
1.0
0.0002831
6297
3 526 000
1.003
377 000

Keren type

front is graphically illustrated by the horizontal lines
running from a nonzero As concentration value to the yaxis. The As(III) results illustrate the position of theoxidation front at 25, 50, and 100 yr (Fig. 2i). At 50 yr, theoxidation front has reached a depth of 15 m. From thesurface to 15 m the available As exists as As(V) as shownin Fig. 2j, while below 15 m reducing conditions exist andthe available As exists as As(III) (Fig. 2i). The graphicalrepresentation of the oxidation front shown in Fig. 2i and2j is expressed again as a continuous time series at thebottom of the heap profile (Fig. 2k). From the onset ofthe simulation until 115 yr, the discharging As from theheap profile is the reduced form, As(III). At 115 yr, theoxidation front reaches the bottom of the profile as aresult of the complete oxidation of arsenian pyrite. As aresult, the redox conditions at the bottom of the profiletransition from reducing to oxidizing, and the discharging As is then present as As(V). The dramatic decrease indischarge As concentration following the arrival of theoxidation front is a direct result of the pH-dependent
sorption characteristics of As(III) and As(V). The rocksorption capacity is much larger for As(V) than As(III),
leading to a larger mass of As(V) bound to the rockmineral surfaces and a proportionally smaller aqueousconcentration. This characteristic is apparent by the very

Fig. 2. Selected results for an uncovered heap-leach facility when the Sips-type isotherm for As sorption was used.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

slow decline in aqueous As(V) concentrations throughout the heap profile with time (Fig. 2j).
The shapes of the As(III) and As(V) concentration
profiles are determined by the dissolution of arsenian
pyrite and the geochemistry that controls pH. The dissolution of arsenian pyrite releases As, which is subjected to the aqueous complexation reactions for As(III)or As(V) and is then distributed between the sorbed
and aqueous phases. This distribution is dependent onoxidation state and pH. The large As(V) sorption capacity of the rock results in an apparent relatively stable concentration profile following the early time redox
front propagation and related pH behavior. The substantially smaller sorption capacity for As(III) results ina more dynamic aqueous concentration response bothwithin the profile and at the lower boundary. This
mobility characteristic has been observed for As(V)and As(III) in saturated environments (Welch andStollenwerk, 2003).

439

Example 2: Heap Facility with Fluid Flux CoverAn approach to limit fluid flux into a heap-leachfacility is to augment the surface with low conductivity
materials or install a cover system that matches fluid flux
with a plant community that is capable of respiring all ormost of this fluid (Albright et al., 2004). The second
example infers the performance of a surface augmentation or evapotranspiration cover system by limiting thefluid flux into the top of the simulated heap facility profile to 50 mm yr21. The reduction in infiltrating flux bya factor of six is a reasonable performance for a nonRCRA cover using on-site borrowed materials. The results of this simulation are shown in Fig. 3. The calculated
steady-state water content is 0.076, which corresponds toa pore volume of 228 cm3 and a time to infiltrate a single
pore volume of 45.6 yr (Fig. 3h). The larger volumetric
air content associated with this reduced water content
results in larger air-phase diffusion-driven transport

Fig. 3. Selected results for a heap-leach facility with a surface treatment to limit fluid flux into the top of the heap when the Keren-type isotherm for
As sorption was used.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

440

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

rates of oxygen (Fig. 3g). As a consequence, the oxygen
limited arsenian pyrite oxidation front also moves faster(Fig. 3b). The corresponding mineral phases also exhibit
substantially faster dissolution or precipitation front velocities for calcite (Fig. 3a), gypsum (Fig. 3c), and ironhydroxide (Fig. 3d). The dissolution of calcite prevents
the acidification of the profile during pyrite dissolution,
maintaining a near-neutral pH profile (Fig. 3e).
Arsenic is initially distributed as As(III) between the
aqueous and sorbed phases in a proportion determinedby the Keren-type isotherm (24) and for the parametervalues given in Table 7 for an initial pH of 8.76 (Fig. 2i).
The As redox front velocity is identical to the pyrite
oxidation and oxygen front velocities (Fig. 3i and 3j). A
consequence of the higher arsenian pyrite oxidation rateis the conversion from As(III) to As(V) at the bottomof the profile after 50 yr (Fig. 3k), 65 yr earlier than in
the previous, uncovered example (Fig. 2k). The rate of

decline in As(V) aqueous concentrations near the topof the profile is substantially reduced compared to thefirst example as a consequence of reducing the meteoric
water flux by a factor of six (Fig. 3j).

Example 3: Heap Facility with Fluid andGas Flux Cover
The third example is representative of a cover systemthat is designed to limit both fluid and gas flux into theheap facility through the top surface. While a typicalcarbonate-hosted heap facility has more than enoughcarbonate mineral mass to neutralize the smaller massfraction of sulfide minerals, reducing the oxidation rate
could slow the release of As from As-bearing mineral phases. In addition, the development of large-scale
preferential flowpaths in a heap structure could reduce
the effectiveness of the carbonate mineral mass to

Fig. 4. Selected results for a heap-leach facility with a cover to limit fluid and gas fluxes into the top of the heap. Arsenic sorption is described using
the Sips-type isotherm.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

neutralize the acid generated from sulfide mineraloxidation. In either case, a cover system may be warranted to control the fluid and chemical flux at thefacility discharge.
The incorporation of sophisticated earthen cover designs into heap facilities for the purposes of reducing theflux of oxidizing gases into a sulfide heap or waste-rock
dump represents an interesting alternative to controlling
sulfide mineral oxidation rates (Kim and Benson, 2004).
This approach incorporates an earthen cover systemthat reduces the gas-phase flux by maintaining a highwater content in the cover, thereby reducing the airphase content and producing a corresponding reductionin the diffusion coefficient for oxygen as given in Eq.
[29]. While the concept of reducing the rate of sulfide
mineral oxidation to reduce the dissolution rate of metals bearing minerals is well founded, the impact of thisdesign on the geochemistry of a sulfide bearing heap
facility is not well understood. The coupled relationships
between oxygen concentration, sulfide mineral oxidation rate, redox- and pH-dependent As aqueous complexation, and water–rock sorption reactions produce
highly nonlinear geochemical behavior. To demonstrate
the efficacy of a gas-flux cover and the nonlinear geochemical response within a heap facility to the installation of a gas-flux cover, this final example is presented.
The gas-flux cover is simulated with an 80-cm-thicklayer of clay loam incorporated into the top of the 30-m

441

profile, while the fluid flux cover is inferred as in theprevious two examples by applying a steady-state recharge flux of 50 mm yr21. The hydraulic and geochemical properties of the heap facility rock are identical tothe previous examples (Tables 5 and 6), and the Sips-type
isotherm is used to describe As sorption (Table 7). The
results are shown in Fig. 4. The oxygen flux rate is
dramatically reduced as compared with the previous
examples (Fig. 4g). The partial pressure of oxygen only
approaches atmospheric values in the cover, while within the profile, the oxidation of arsenian pyrite proceeds
much faster than the oxygen flux rate, maintaining low
partial pressures of oxygen within the profile. The reduced oxygen flux rate slows the arsenian pyrite oxidation rate such that the oxidation front has only reached
20 m after 200 yr (Fig. 4b). The remaining mineralogy
that is dependent on arsenian pyrite oxidation chemistry
is similarly affected (Fig. 4a, 4c, and 4d).
The reduced gas flux cover significantly reduces theflux of oxygen through the top of the simulated heap
profile such that reducing conditions are maintained
within the profile for the entire 200-yr simulation. As aresult, the initial reducing conditions and As(III) speciation are maintained throughout the profile (Fig. 4i),
with the exception of the cover, wherein conditions are
such that As occurs as As(V) (Fig. 4j). The cover prevents the complete oxidation of the sulfide mineral mass
within the heap profile, and maintains redox conditions

Fig. 5. Time series of aqueous As concentrations at the lower boundary discharge point for all examples and both isotherm types.

442

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

such that As(III) speciation prevails for the duration of
the simulation (Fig. 4k).

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

SUMMARY
The three examples are intended to provide a range ofpossible behaviors that could be expected for As oxyanion transport for the most likely heap facility closure
options for a carbonate-rich, low-sulfide ore such as
those typical of the Carlin Trend operations. The As
breakthrough curves at the bottom of the simulation
profile for the three examples are presented in Fig. 5. The
cumulative As mass released from the bottom of the profile is plotted in Fig. 6 for all three examples using theSips-type isotherm. Each of the three examples was run
with both the Sips-type and Keren-type isotherm formulations (Fig. 5). The differences between the isotherm
formulations are a result of the differences in the mathematics of each and the overall fit to the experimental data
(Decker et al., 2006). The sharp drop in As concentrationin Fig. 5 at 50 and 115 yr is a manifestation of theoxidation front arriving at the bottom of the simulated
heap profile, and As transitioning from As(III) to As(V).
The arrival of the oxidation front is similarly the causeof the break in slope at the same times in Fig. 6. The
differences in release characteristics for As are a result of
the coupled effects of variably saturated transport andwater–rock geochemistry.
The uncovered example, with a fluid flux rate six
times higher than the other examples, produces thelargest cumulative As mass at the simulated outfall
located at the bottom of the modeled domain. As a re-sult of the high fluid flux rate, aqueous As concentration
drops rapidly as the sorbed As is rinsed from rock sorption sites. The second example, with a simple fluid flux
cover, produces a proportionally smaller cumulative As
mass, and an overall higher As concentration at thebottom of the simulated heap profile. The reduced flux
rate results in a low volumetric water content and a
proportionally larger air-phase volumetric content. This
produces a larger effective diffusion coefficient for oxygen, and a correspondingly faster sulfide mineral oxidation rate. The arrival of the oxidation front in thesecond example 65 yr earlier than the uncovered simulation is a direct result of modifying the rate of fluid flux
into the simulated profile.
The third example, with a fluid and gas flux cover,produces a larger cumulative As mass at the outfall than
the simple, fluid-control cover example, despite the factthat both examples are simulated with the same infiltrating fluid flux. Because the gas flux cover prevents theoxidation of As(III) to As(V), and because As(III) is
substantially more mobile than As(V), the gas flux coverallows the faster release of As for a given fluid flux rate.
A consequence of maintaining reducing conditions within
the simulated profile is the persistent high concentrationof As at the outfall (Fig. 5).

CONCLUSIONS
An existing variably saturated reactive transport code
was modified so that As transport in heap-leach facilities
could be simulated. The code modifications included theincorporation of two pH-dependent isotherm formula-Fig. 6. Time series of cumulative aqueous As mass at the lower boundary discharge point for the three examples using the Sips-type isotherm.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

tions for As, arsenian pyrite oxidation and related rock–water chemistry, and oxygen transport. The modified
code is capable of simulating the primary processes
controlling redox- and pH-dependent As transport. One
possible application of the code is to estimate the performance of heap-leach facility closure options that involve modifying the hydrologic or geochemical regime
to control As release.
The As transport capable version of UNSATCHEM
is useful as an initial estimator for As release behavior
with a modest cost in obtaining whole-rock As sorption
behavior in the laboratory. The rapid estimation of
As release behavior with a numerical solver that incorporates the primary physical and chemical processes
controlling transport is expected to be useful to the industrial, engineering, and regulatory communities. Additional improvements to UNSATCHEM could include
the incorporation of the Davis–Ritchie shrinking-core
pyrite oxidation model, and the inclusion of site-specificsecondary mineralogy.
Three examples were provided that are intended todemonstrate the code capabilities and the importance ofincluding the geochemistry that controls oxyanion behavior in transport simulations. The results of the simulations
provide a basis for a first round of engineering cost–benefit analysis using the solute (As in this case) cumulative mass and concentration behavior. If maintaining
low As concentration at the heap outfall is the primary
metric for heap facility management performance, thenthese simulations show that the heap should be oxidized
as quickly as possible to convert the mass of As to As(V)within the heap structure. Because the sorption capacity
for As(V) is higher than that for As(III), the release of As
will occur over a longer period of time, but at much lower
concentrations in an oxidized heap, as suggested by the
results from the second example. If the long-term costs ofconstructing and maintaining a discharge fluid-treatment
system are of primary concern, then optimizing the As
release behavior to achieve a particular As release function that is coupled to the treatment plant design is ofprimary importance and can be done with a higher degreeof realism by including the coupled physical and geochemical processes that control As release behavior.
REFERENCES
Albright, W.H., C.H. Benson, G.W. Gee, A.C. Roesler, T. Abichou, P.
Apiwantragoon, B.F. Lyles, and S.A. Rock. 2004. Field waterbalance of landfill final covers. J. Environ. Qual. 33:2317–2332.
Appelo, C.A.J., and D. Postma. 1993. Geochemistry, groundwater andpollution. A.A. Balkema, Brookfield, VT.
Bartlett, R.W. 1992. Solution mining: Leaching and fluid recovery ofmaterials. Gordon and Breach Science Publishers, Philadelphia.
Decker, D.L., C. Papelis, S.W. Tyler, M. Logsdon, and J. Šimůnek.
2006. Arsenate and arsenite sorption on carbonate-hosted precious
metals ore. Available at www.vadosezonejournal.org. Vadose Zone
J. 5:xx–xx (this issue).
Decker, D.L., and S.W. Tyler. 1999a. Evaluation of flow and solutetransport parameters for heap-leach recovery materials. J. Environ.
Qual. 28:543–555.
Decker, D.L., and S.W. Tyler. 1999b. Hydrodynamics and solutetransport in heap-leach mining. p. 1–13. In G.C. Miller and D.Y.
Kosich (ed.) Closure, remediation and management of precious
metals heap-leach facilities. Univ. of Nevada, Reno.

443

Diodato, D.M., and R.R. Parizek. 1994. Unsaturated hydrogeologic properties of reclaimed coal strip mines. Ground Water 32:
108–118.
Driesner, D., and A. Coyner. 2000. Major mines in Nevada 1999.
Nevada Bureau of Mines and Geology, Reno.
Glinski, J., and W. Stepniewski. 1985. Soil aeration and its role forplants. CRC Press, Boca Raton, FL.
Hillel, D. 1998. Environmental soil physics. Academic Press, New York.
Jambor, J.L. 1994. Mineralogy of sulfide-rich tailings and their oxidation products. p. 59–102. In J.L. Jambor and D.W. Blowes (ed.) The
environmental geochemistry of sulfide mine-wastes. Vol. 22. Mineralogical Association of Canada, Ottawa, Canada.
Jerz, J.K., and J.D. Rimstidt. 2004. Pyrite oxidation in moist air.
Geochim. Cosmochim. Acta 68:701–714.
Jury, W.A., and R. Horton. 2004. Soil physics. 6th ed. John Wiley andSons, New York.
Kampf, S., M. Salazar, and S. Tyler. 2002. Preliminary investigations ofeffluent drainage from mining heap-leach facilities. Available atwww.vadosezonejournal.org. Vadose Zone J. 1:186–196.
Kim, H., and C.H. Benson. 2004. Contributions of advective anddiffusive oxygen transport through multilayer composite caps overmine waste. J. Contaminant Hydrol. 71:193–218.
Krauskopf, K.B., and D.K. Bird. 1995. Introduction to Geochemistry.
3rd ed. McGraw-Hill, New York.
Langmuir, D. 1997. Aqueous environmental geochemistry. Prentice
Hall, Upper Saddle River, NJ.
Lichtner, P.C., and M.S. Seth. 1996. User’s manual for MULTIFLO:Part II MULTIFLO 1.0 and GEM 1.0, Multicomponent-multiphasereactive transport model. CNWRA 96–010.
Lichtner, P.C., S. Yabusaki, K. Pruess, and C.I. Steefel. 2004. Role ofcompetitive cation exchange on chromatographic displacementof cesium in the vadose zone beneath the Hanford S/SX Tank
Farm. Available at www.vadosezonejournal.org. Vadose Zone J. 3:
203–219.
Mallants, D., B. Mohanty, A. Vervoort, and J. Feyen. 1997. Spatial
analysis of saturated hydraulic conductivity in a soil with macropores. Soil Technol. 10:115–131.
Mayer, U., E.O. Frind, and D.W. Blowes. 2002. Multicomponent
reactive transport modeling in variably saturated porous mediausing a generalized formulation for kinetically controlled reactions.
Water Resour. Res. 38(9):1174. doi:10.1029/2001WR000862.
Nichol, C.F. 2002. Transient flow and transport in unsaturated
heterogeneous media; field experiments in mine waste rock.
Ph.D. diss. University of British Columbia, Vancouver, Canada.
Nichol, C., L. Smith, and R. Beckie. 2003. Long-term measurement of
matric suction using thermal conductivity sensors. Can. Geotech. J.
40:587–597.
Pruess, K. 1991. TOUGH2—A general purpose numerical simulator
for multiphase fluid and heat flow. LBL-29400. Lawrence Berkeley
Natl. Lab., Berkeley, CA.
Ritchie, A.I.M. 1994. Sulfide oxidation mechanisms: Controls andrates of oxygen transport. p. 201–246. In D.W. Blowes and J.L.
Jambor (ed.) The environmental geochemistry of sulfide minewastes. Vol. 22. Mineralogical Association of Canada, Ottawa.
Šimůnek, J., and D.L. Suarez. 1993. Modeling of carbon dioxidetransport and production in soil: 1. Model development. Water
Resour. Res. 29:487–497.
Šimůnek, J., and D.L. Suarez. 1994. Two-dimensional transport model
for variably saturated porous media with major ion chemistry.
Water Resour. Res. 30:1115–1133.
Šimůnek, J., D.L. Suarez, and M. Šejna. 1996. The UNSATCHEM
software package for simulating the one-dimensional variably
saturated water flow, heat transport, carbon dioxide production
and transport, and multicomponent solute transport with major ion
equilibrium and kinetic chemistry. Res. Rep. 141. U.S. Salinity
Laboratory, Riverside, CA.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1998. The
HYDRUS-1D software package for simulating the one-dimensional
movement of water, heat, and multiple solutes in variably-saturated
media. Version 2.0. IGWMC-TPS-70. International Ground Water
Modeling Center, Colorado School of Mines, Golden, CO.
van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:
892–898.

444

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Viswanathan, H.S., B.A. Robinson, A.J. Valocchi, and I.R. Triay. 1998.
A reactive transport model of neptunium migration from thepotential repository at Yucca Mountain. J. Hydrol. (Amsterdam)
209:251–280.
Welch, A.H., and K.G. Stollenwerk (ed.). 2003. Arsenic in groundwater. Kluwer Academic Publishers, Dordrecht, The Netherlands.

Welch, A.H., D.B. Westjohn, D.R. Helsel, and R.B. Wanty. 2000.
Arsenic in ground water of the United States: Occurrence andgeochemistry. Ground Water 38:589–604.
Williamson, M.A., and J.D. Rimstidt. 1994. The kinetics andelectrochemical rate-determining step of aqueous pyrite oxidation.
Geochim. Cosmochim. Acta 58:5443–5454.