Published online 8 March 2006
Published in Vadose Zone J 5:430-444 (2006)
DOI: 10.2136/vzj2005.0039
© 2006 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Variably Saturated Reactive Transport of Arsenic in Heap-Leach Facilities
David L. Deckera,*,
Jirka 
im
nekb,
Scott W. Tylerc,
Charalambos Papelisd and
Mark J. Logsdone
a Desert Research Institute, 2215 Raggio Parkway, Reno, NV 89512 USA
b Dep. of Environmental Sciences, A135 Bourns Hall, Univ. of California, 900 University Ave., Riverside, CA 92521 USA
c Dep. of Natural Resources and Environmental Sciences and Dep. of Geological Sciences and Engineering, MS 175, Univ. of Nevada, Reno, NV 89557-0180 USA
d 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
e Geochimica Inc., 9045 Soquel Drive, Suite 2, Aptos, CA 95003 USA
* Corresponding author (dave.decker{at}dri.edu)
Received 14 March 2005.
 |
ABSTRACT
|
|---|
A 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 and western Utah, are characterized in part by the occurrence of As-bearing minerals. The large-scale mining methodology used on the Carlin Trend includes the use of large-scale heap-leach facilities in which 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 associated with other mining provinces around the world, controlling the flux of As 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.
|
INTRODUCTION
|
|---|
THE SUBSURFACE TRANSPORT of metals in variably saturated media is of interest because of the increasing demands on subsurface water supplies for drinking, industrial processes, and agriculture. Particular interest in arsenic, selenium, and chromium is driven by the biological toxicity of compounds formed with these elements. Arsenic toxicity has induced a reduction in the federal drinking water standard from 50 to 10 µg L–1 in the United States. The financial impacts of treating drinking water for As has led to renewed interest in understanding the processes controlling the fate and transport of As compounds in the subsurface.
Arsenic is prevalent in the known groundwater resources within the Great Basin of Utah, Nevada, and California. Much of this As is thought to be the result of dissolution of As-bearing, broadly distributed, primary rock and mineral formations (Welch et al., 2000). Arsenic bearing minerals are also present in the form of mineralized zones associated with the finely disseminated gold-bearing deposits of eastern and northeastern Nevada and western Utah. The co-occurrence of As with gold in these deposits is significant enough that As is a recommended gold-ore prospecting indicator. The resurgence in gold mining within the Great Basin since the discovery and initial development of the Carlin Trend has resulted in the most productive gold-producing region in North America (Driesner and Coyner, 2000). The large-scale mining techniques used on the Carlin Trend include the use of heap-leach mining technologies wherein the gold is extracted by leaching large volumes of rock that have 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 is leached with a sodium-cyanide solution. A byproduct of leaching is the dissolution of other metal and metalloid-bearing mineral phases. The dissolution of these minerals in turn results in significant aqueous concentrations of undesired constituents such as As compounds. Once the gold is extracted from the heap-leach facility, effective methods of either reducing the concentration or the flux of 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 in hydrologic regimes, necessitates the development of tools to assist in the design of effective heap-leach closure methodologies. While the Carlin Trend is not located near large metropolitan areas, it does encompass several large hydrographic 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 transport in heap-leach rock from mines located within the Carlin Trend. The physics of fluid transport in these highly heterogeneous, variably saturated environments were studied 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 the Carlin Trend in an effort to understand the mechanisms of As release and transport in a mine setting. As part of this work, empirically derived pH-dependent isotherm formulations were developed with the intention of incorporating pH-dependent sorption behavior into a variably saturated reactive flow and transport model 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.
|
MODELING REACTIVE ARSENIC OXYANION TRANSPORT
|
|---|
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 (imnek and Suarez, 1994; imnek 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 (imnek et al., 1998), the inclusion of carbonate geochemistry with the associated CO2 transport subroutines (imnek and Suarez, 1993), and the public-domain accessibility of the source code to enable the addition of As aqueous chemistry, sorption, and mineral dissolution processes.
Variably Saturated Flow
A successful variably saturated flow and reactive transport model must include a mathematical approach describing 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 (imnek et al., 1996, 1998). The governing equation for the isothermal Darcian flow of water in a variably saturated rigid porous medium using the assumptions that the air phase plays an insignificant role in the liquid flow process and that water flow due to thermal gradients can be neglected is given by Richards' equation as
 | [1] |
where 
w is the volumetric water content, h is the matric potential or pressure head, K(h) is the hydraulic conductivity 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
 | [2] |
The variably saturated hydraulic properties for soil are often described using van Genuchten's (1980) equations for the retention and unsaturated hydraulic conductivity functions as
 | [3] |
and
 | [4] |
where
 | [5] |
 | [6] |
 | [7] |
in which r and s 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 
and 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 to macropore 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 hydraulic conductivity function given above are reported in the literature, 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 specific problem of heap-leach flow modeling is given in Decker and Tyler (1999b).
Multicomponent Solute Transport
The partial differential equation describing one-dimensional advective–dispersive solute transport under transient water flow conditions in partially saturated porous media is taken as (imnek et al., 1996)
 | [8] |
where Ck is the total dissolved concentration of the aqueous species k, Sk is the total sorbed species concentration of the aqueous component k, Mk is the total solid phase concentration of the aqueous component k, 
is the bulk density of the medium, D is the dispersion coefficient, qw is the volumetric flux, and Nc is the number of primary aqueous species. The volumetric flux is calculated with Darcy's Law
 | [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, and solid components of the aqueous species k for all Nc species as
 | [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 as
 | [11] |
while the third type of boundary condition, Cauchy, describes the concentration flux as
 | [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 the soil profile.
|
GEOCHEMICAL REACTIONS IN UNSATCHEM
|
|---|
UNSATCHEM is a variably saturated reactive transport code that was designed to simulate major ion transport in agricultural soils (imnek et al., 1996). The code numerically solves the Richards equation for variably 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 Ca2+, Mg2+, Na+, K+, SO42–, Cl–, alkalinity, and CO2. The code accounts for temperature-dependent equilibrium chemical reactions between these components through complexation, cation exchange, and precipitation–dissolution. The mineral phases include calcite, dolomite, gypsum, nesquehonite, hydromagnesite, and sepiolite. The precipitation–dissolution reactions for calcite can be modeled using either equilibrium or kinetic expressions, while the dissolution of dolomite is always considered as a kinetic process. The model is capable of simulating saline water movement, and modified Debye-Hü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 the aqueous 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 of CO2 with a mass-balance relationship that includes diffusive and convective transport in both gas and fluid phases. A compilation of the aqueous complexes and mineral species included in UNSATCHEM is given in Table 1, while thermodynamic data, kinetic reaction constants, and other details are provided in imnek et al. (1996).
|
MODIFICATIONS TO UNSATCHEM
|
|---|
The UNSATCHEM code has been modified in this work to include the geochemistry associated with As in carbonate rock. The addition of As transport to the original 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 oxygen-dependent 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) [As3+] and As(V) [As5+] (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 mineral source of As is dependent on the specific mineralogical composition of each mine site. However, a common As-bearing mineral phase found in Carlin Trend ore bodies is As-bearing pyrite. UNSATCHEM has been modified to include three forms of this mineral including pyrite (FeS2), arsenopyrite (FeAsS), and arsenian pyrite (FeAsxS2–x), where x is the fraction of As. The dissolution of pyrite, or As-bearing pyrite, releases aqueous Fe. UNSATCHEM has been modified to include the precipitation reaction for iron hydroxide to remove the aqueous Fe mass from solution. In addition, the dissolution reactions of pyrite, As-bearing pyrite, and the precipitation reaction of iron hydroxide, 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 and pH 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 the modifications to UNSATCHEM that were undertaken to provide As transport capability are given in subsequent sections.
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 mass balance, two additional equations for the two As oxidation states considered are
 | [13] |
 | [14] |
where variables with the subscript T represent the total analytical concentrations in solution of the particular species, and where the brackets refer to molalities (mol kg–1). In addition to these mass-balance equations, the overall charge-balance equation, including As complexes, is expressed as
 | [15] |
Precipitation–Dissolution Reactions
Inorganic pyrite dissolution is described by a series of reactions that encompass a wide pH range. These reactions are commonly written as (Langmuir, 1997)
 | [16] |
 | [17] |
 | [18] |
The Carlin Trend heaps represent a somewhat unique environment to describe pyrite oxidation due to the presence of significant amounts of carbonate minerals relative to the low mass percentages of sulfide minerals. The equilibrium pH for the solution in equilibrium with both gypsum and calcite minerals varies between 7.8 and 8.3 (Krauskopf and Bird, 1995). Because of this high pH environment, the dissolution of pyrite by Fe3+ is minimized as a result of the rapid precipitation of Fe3+ to Fe(OH)3 by the following reaction:
 | [19] |
which, in the flow-through environment of the heap-leach facility, is assumed to be irreversible. Therefore, the reactions given in Eq. [17] and [18] are ignored, and the 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, the amount 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 a controlled laboratory environment, and is given by (Williamson and Rimstidt, 1994) as
 | [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 the Williamson and Rimstidt (1994) rate expression is incorporated into UNSATCHEM.
At high pH, Fe(OH)3 is relatively insoluble, and precipitates by Eq. [19]. An empirical rate law for this reaction is given by
 | [21] |
for which the rate of iron hydroxide precipitation is expressed in moles per second at 20°C (Langmuir, 1997).
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
 | [22] |
where x is the fraction of As present in pyrite. Note that setting x = 0 results in the pyrite oxidation reaction given in Eq. [16], and setting x = 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 invoked by 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 and sorbed phases is accomplished using the two pH-dependent isotherm expressions described in Decker et al. (2006). The first expression is a Sips-type, or Langmuir–Freundlich, isotherm that is pH dependent.
 | [23] |
This expression encompasses two analytically measured parameters, Ce and pH, and one calculated parameter, qe. There are three isotherm-fitting parameters from the Sips formulation, qm, ns, and as, and incorporating pH dependence adds two additional fitting parameters, 
and 
. The simulations included in this study simplified this model to a Langmuir form by assigning ns = 1. In addition, a Keren-type pH-dependent isotherm is also incorporated into UNSATCHEM
 | [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 was taken. An approximation of the redox potential can be made 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 controlling the oxidation state of As. In practice, this is accomplished by specifying the critical oxygen concentration that defines the threshold between oxidizing and reducing conditions as an input variable. This variable in turn is used to determine 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 to control As redox chemistry presumes that the geochemistry associated with Fe, O2, and As is coupled. In the presence of sulfide mineralogy, this assumption is probably reasonable since the coupled interaction of O2 and Fe drives the redox chemistry. However, in the absence of sulfide mineralogy, the use of the partial pressure for O2 may not be a satisfactory redox surrogate for As, and explicit As redox chemistry may need to be incorporated into UNSATCHEM to adequately represent this type of domain. The Carlin Trend ore deposits are associated with sulfide mineralogy, and therefore the assumption is made that the use of O2 is a reasonable surrogate representation of the complex redox chemistry for As.
The oxygen transport modeling approach is similar to that used for carbon dioxide transport (imnek and Suarez, 1993). The gas transport mass conservation equation is
 | [25] |
where Jda describes the O2 flux caused by diffusion in the gas phase, Jdw is the O2 flux caused by dispersion in the dissolved phase, Jca is the O2 flux caused by convection in the gas phase, and Jcw is the O2 flux caused by convection in the dissolved phase. The term CT is the total volumetric concentration of O2, and P is a source–sink term. These terms are defined as
 | [26] |
where Cw and Ca are the volumetric concentrations of O2 in the dissolved and gas phases, respectively, Da is the effective soil matrix diffusion coefficient of O2 in the 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 a 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 = 0). Another possibility is to assume that because of the much lower viscosity of air in comparison with water, a relatively 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 the gas 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 the total oxygen flux from convection in the air phase is 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 contribute to the overall oxygen mass budget as a result of the exothermic heat of reaction produced from the dissolution of sulfides (Ritchie, 1994). However, it is unlikely that convective O2 flux is a major feature of the Carlin Trend heap-leach facilities for three reasons. First, the sulfide 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 is highly likely that a temperature rise of <2°C could be expected. Second, the heaps are constructed in a way that limits the likelihood of vertically continuous high-conductivity 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 the heat generated during pyrite dissolution is included in this version of UNSATCHEM.
The term P in Eq. [25] represents a sink term for oxygen resulting from the consumption of oxygen by chemical reactions such as those given in Eq. [16], [19], and [22]. No sources of oxygen other than the atmosphere are included in this formulation. Henry's Law governs the distribution of gas between aqueous and gaseous phases
 | [27] |
where KO2 is the Henry's Law constant (1.26 x 10–3 mol bar–1 (Langmuir, 1997)), R is the universal gas constant (8.314 kg m2 s–2 K–1 mol–1), and T is the absolute temperature (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 of a water film. Therefore, Henry's Law is used to move gaseous 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 are calculated from the molecular diffusion coefficients for oxygen 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 as
 | [28] |
 | [29] |
where Dws = 2.10 x 10–9 m2 s–1 and Das = 2.01 x 10–5 m2 s–1 are the molecular diffusion coefficients for O2 in the dissolved and gas phases from Glinski and Stepniewski (1985), 
a and w are the tortuosity factors in both phases, and 
w is the dispersivity in the dissolved phase (imnek et al., 1996).
Numerical Solution Strategy
The numerical solution of Richards' equation is accomplished with the use of a mass-lumped, finite-element method that is reduced to a finite-difference scheme (imnek 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 the solute, heat, and gas transport equations are accomplished with the Galerkin finite element method. The flow and transport models are inclusive of a full suite of boundary conditions including specified head and concentration (Dirichlet), or specified flux (Neumann).
An iterative solution is used to solve the chemical system and is discussed in the UNSATCHEM manual (imnek et al., 1996). The iteration process must be accomplished at each node, for each time step, and between the solute transport and equilibrium components of the model. The solution process at each time step proceeds by iteratively solving the flow equation and solving the transport equations for temperature, CO2, and then multicomponent solute transport. The geochemical module is invoked at each time step, for each node, to redistribute mass between solid, aqueous, and sorbed phases.
The governing solute transport equation (Eq. [8]) includes three time-derivative components: total aqueous concentration, 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 in heap-leach systems are presented to illustrate the importance of incorporating reactive geochemistry into transport simulations that attempt to estimate As fluxes at long time-scales. Because the model examples deal with variably saturated media, gas-phase transport of oxygen and carbon dioxide is an important process controlling 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 is dependent on mineralogy, geochemistry, and the physics of flow in porous media. The examples are intended to demonstrate 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 representations of 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 to Carlin Trend ores. These examples include the following heap facility closure scenarios: no surface treatment, an engineered 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 is shown 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 are taken 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 with UNSATCHEM and HYDRUS-1D. The steady-state, upper flux boundary conditions are used to simplify the interpretation of the reactive transport results. The lower flux boundary condition is free drainage, or zero pressure 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 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 is deeper than 30 m.
The geochemical conditions within the heap and the covers 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 interstitial water of the heap as well as the geochemistry of the recharging 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 the experimental 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 to 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 for both gases. Thus, gas transport into the heap profile is limited to the top surface.
View this table:
[in this window]
[in a new window]
|
Table 7. pH-dependent isotherm surface model parameters. Units for the calculated isotherm surface plot are: Ce, µg L–1; qe, µg g–1; pH = –log[H+].
|
|
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 this broad 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 of long-term management of both the remaining mineral resource in the facility, and in ensuring effective environmental stewardship.
Mine operations are required to close heap-leach facilities such that a number of environmental regulations are met. This is typically accomplished through a regimen of rinsing leachate fluids from the heap with fresh 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 of such magnitude that designing an effective closure plan that includes understanding the impacts of engineering solutions to limit the precipitation flux into the top of the heap facility on the chemistry of the discharge water is key to minimizing costs and maximizing the effectiveness of the management plan. The three examples presented illustrate the long-term effects on As concentration and mass 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 yr–1 that is representative of a high-elevation 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 to represent a heap that is recovering from a period of gold production wherein the pH is maintained at very high levels (in excess of pH 10) to prevent the loss of cyanide (Bartlett, 1992) (Fig. 2e). This pH rapidly decreases to lower 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 the dissolution 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 rate that 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).
View larger version (46K):
[in this window]
[in a new window]
|
Fig. 2. Selected results for an uncovered heap-leach facility when the Sips-type isotherm for As sorption was used.
|
|
Arsenic is initially distributed as As(III) between the aqueous and sorbed phases in a proportion determined by 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 with the oxygen transport and arsenian pyrite oxidation fronts. In Fig. 2i and 2j, the presence of the oxidation front is graphically illustrated by the horizontal lines running from a nonzero As concentration value to the y axis. The As(III) results illustrate the position of the oxidation front at 25, 50, and 100 yr (Fig. 2i). At 50 yr, the oxidation front has reached a depth of 15 m. From the surface to 15 m the available As exists as As(V) as shown in Fig. 2j, while below 15 m reducing conditions exist and the available As exists as As(III) (Fig. 2i). The graphical representation of the oxidation front shown in Fig. 2i and 2j is expressed again as a continuous time series at the bottom of the heap profile (Fig. 2k). From the onset of the simulation until 115 yr, the discharging As from the heap profile is the reduced form, As(III). At 115 yr, the oxidation front reaches the bottom of the profile as a result of the complete oxidation of arsenian pyrite. As a result, the redox conditions at the bottom of the profile transition from reducing to oxidizing, and the discharging As is then present as As(V). The dramatic decrease in discharge As concentration following the arrival of the oxidation front is a direct result of the pH-dependent sorption characteristics of As(III) and As(V). The rock sorption capacity is much larger for As(V) than As(III), leading to a larger mass of As(V) bound to the rock mineral surfaces and a proportionally smaller aqueous concentration. This characteristic is apparent by the very 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 on oxidation 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 in a more dynamic aqueous concentration response both within the profile and at the lower boundary. This mobility characteristic has been observed for As(V) and As(III) in saturated environments (Welch and Stollenwerk, 2003).
Example 2: Heap Facility with Fluid Flux Cover
An approach to limit fluid flux into a heap-leach facility 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 or most of this fluid (Albright et al., 2004). The second example infers the performance of a surface augmentation or evapotranspiration cover system by limiting the fluid flux into the top of the simulated heap facility profile to 50 mm yr–1. The reduction in infiltrating flux by a factor of six is a reasonable performance for a non-RCRA 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 to a 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 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 iron hydroxide (Fig. 3d). The dissolution of calcite prevents the acidification of the profile during pyrite dissolution, maintaining a near-neutral pH profile (Fig. 3e).
View larger version (45K):
[in this window]
[in a new window]
|
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.
|
|
Arsenic is initially distributed as As(III) between the aqueous and sorbed phases in a proportion determined by the Keren-type isotherm (24) and for the parameter values 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 rate is the conversion from As(III) to As(V) at the bottom of 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 top of the profile is substantially reduced compared to the first example as a consequence of reducing the meteoric water flux by a factor of six (Fig. 3j).
Example 3: Heap Facility with Fluid and Gas Flux Cover
The third example is representative of a cover system that is designed to limit both fluid and gas flux into the heap facility through the top surface. While a typical carbonate-hosted heap facility has more than enough carbonate mineral mass to neutralize the smaller mass fraction 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 neutralize the acid generated from sulfide mineral oxidation. In either case, a cover system may be warranted to control the fluid and chemical flux at the facility discharge.
The incorporation of sophisticated earthen cover designs into heap facilities for the purposes of reducing the flux 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 system that reduces the gas-phase flux by maintaining a high water content in the cover, thereby reducing the air-phase content and producing a corresponding reduction in 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 this design 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-thick layer of clay loam incorporated into the top of the 30-m profile, while the fluid flux cover is inferred as in the previous two examples by applying a steady-state recharge flux of 50 mm yr–1. The hydraulic and geochemical properties of the heap facility rock are identical to the 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).
View larger version (44K):
[in this window]
[in a new window]
|
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.
|
|
The reduced gas flux cover significantly reduces the flux 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 a result, 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 such that As(III) speciation prevails for the duration of the simulation (Fig. 4k).
|
SUMMARY
|
|---|
The three examples are intended to provide a range of possible 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 the Sips-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 concentration in Fig. 5 at 50 and 115 yr is a manifestation of the oxidation 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 cause of 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 and water–rock geochemistry.
View larger version (25K):
[in this window]
[in a new window]
|
Fig. 5. Time series of aqueous As concentrations at the lower boundary discharge point for all examples and both isotherm types.
|
|
View larger version (17K):
[in this window]
[in a new window]
|
Fig. 6. Time series of cumulative aqueous As mass at the lower boundary discharge point for the three examples using the Sips-type isotherm.
|
|
The uncovered example, with a fluid flux rate six times higher than the other examples, produces the largest cumulative As mass at the simulated outfall located at the bottom of the modeled domain. As a result 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 the bottom 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 the second 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 fact that both examples are simulated with the same infiltrating fluid flux. Because the gas flux cover prevents the oxidation of As(III) to As(V), and because As(III) is substantially more mobile than As(V), the gas flux cover allows 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 concentration of 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 the incorporation of two pH-dependent isotherm formulations 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-specific secondary mineralogy.
Three examples were provided that are intended to demonstrate the code capabilities and the importance of including 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, then these 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 of constructing 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 of primary importance and can be done with a higher degree of 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 water balance of landfill final covers. J. Environ. Qual. 33:2317–2332.[Abstract/Free Full Text]
- Appelo, C.A.J., and D. Postma. 1993. Geochemistry, groundwater and pollution. A.A. Balkema, Brookfield, VT.
- Bartlett, R.W. 1992. Solution mining: Leaching and fluid recovery of materials. Gordon and Breach Science Publishers, Philadelphia.
- Decker, D.L., C. Papelis, S.W. Tyler, M. Logsdon, and J. imnek. 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 solute transport parameters for heap-leach recovery materials. J. Environ. Qual. 28:543–555.
TOP
ABSTRACT
INTRODUCTION
