Development of a geochemical model to predict leachate water quality associated with coal mining practices KN van Zweel 20706391 Dissertation submitted in fulfilment of the requirements for the degree Magister Scientiae in Environmental Sciences (specialising in Hydrology and Geohydrology) at the Potchefstroom Campus of the North-West University Supervisor: Dr SR Dennis May 2015 Page 2 of 77 Summary South Afric a mines coal to supply in the growing energy demands of the country. A majority of these mines are opencast result ing in back filled pits and above ground disposal facilities. Leachate emanating from these disposal sites are saline and in most cases highly acidic. Currently the standard testing procedure to quantify expected leachate qualities include Acid Base Accounting (ABA), Net -acid Generating test (NAG), static -and kinetic leaching. The aim of this study is to model standard humidity cell leach tests performed using the PHREEQC code. This model can then be scaled up to field conditions to model 1D reactive transport. It is commonly accepted that the rate of pyrite oxidation in backfil led pits and waste storage facilities is governed by the rate of oxygen ingress and that no pyrite oxidation take place in the saturated zone. This is not the case for humidity cells, as sufficient oxygen is available for reaction. Pyrite reactions rates i n humidity cells is expected to be governed by a combination of available reaction surface and ash layer resistance. This is modelled in PHREEQC (Parkhurst & Appelo, 2003 ) using the KINETIC block. Leachate composition is then modelled in the column by maki ng use of the TRANSPORT block. The experimental data is fitted by using the reactive surface and ash layer diffusion coefficient as a fitting parameter. PHREEQC does not have a gas transport module to model oxygen diffusion through the column. Due to this shortfall of PHREEQC, the influence of oxygen ingress in the system can not be directly modelled under kinetic conditions . Davis and Ritchie (1986 a) proposed an anlylitical solution in which the integrated sulphate production rate can be calculated for a waste heap dump. This rate takes into account the influence of oxygen ingress and the development of an ash layer resistance to the pyrite oxydation rate. This intergrated rate can then be defined in a RATES block in PHREEQC. The Aproximate Analytical So lution (AAS) model proposed by Davis and Ritchie (1986 a) is used to scale up the model used for the humidity leach cell experiment . It was found from the modelling results and comparison with PYROX that the model under predicts the integrated sulphate production rate in the initial stages of the reacting waste heap dump. It does however show results that are in close agreement Page 3 of 77 with the results obtained from PYROX in later stages of the lifespan of the waste heap dump. This highlight limitations to the AAS model’s applicability on geochemical problems. The model can only be applied to describe waste heap dumps where the particles at the top of the heap are fully oxidized. Keywords: PHREEQC, Humidity leach cell, Pyrite, Kinetic mod elling Page 4 of 77 List of Abbreviations HLC Humidity Leach Cell ABA Acid Base Accounting PHEEEQC A Computer Program for Speciation, Batch -Reaction, One - Dimensional Transport, and Inverse Geochemical Calculations NAG Net Acid Generation ARD Acid Rock Drainage CSTR Continuous Stirred Tank Reactor SK Surface Kinetic ASTM American Society for Testing and Materials PYROX Geochemical model developed by Wuderley and Blowes (1997) ABATE Acid -Base: Accounting, Techniques and Evaluation 1D One Dimensional 2D Two Dimensional 3D Three Dimensional XRD X -Ray Diffraction AAS Approximate Analytical Solution MATLAB Mathematical program developed by MathWorks ® SCM Shrinking Core Model Ash Layer Layer forming around the reactive core of a particle Page 5 of 77 Table of Contents Summary ............................................................................................................... 2 List of Abbreviations ............................................................................................ 4 Table of Contents .................................................................................................. 5 List of Figures ....................................................................................................... 7 List of Tables ......................................................................................................... 9 1. Introduction ..................................................................................................... 10 1.1 BACKGROUND .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2 PROBLEM STATEMENT .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3 AIMS AND OBJECTIVES .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2. Literature Survey ............................................................................................ 12 2.1 INTRODUCTION .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 PYRITE OXIDATION .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 REACTION KINETICS .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 Surface reaction kinetics .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.2 Shrinking Core Model .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2. 4 PREDICTIVE METHODS .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.4.1 Static test .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4.2 Kinetic tests .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5 MODE LLING OF ACID MINE D RAINAGE .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.6 MODELLING SPOILS HEA PS .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3. Development methodology and model code ............................................... 29 3.1 INTRODUCTION .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2 MODELLIN G HCT WITH THE SURFA CE KINETIC MODEL .... . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2.1 Developing a methodology to describe a CSTR in PHREEQC .... . . . . . . . . . . . . . . 33 3.3 MODEL SETUP IN PHREE QC .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3.1 Defining the reaction kinetics .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3.2 Equilibrium phases .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.3. 3 Calculating composition from mineralogy .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.3. 4 Defining the reaction surface .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Page 6 of 77 3.4 UP-SCALING TO FIELD CON DITIONS .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.4.1 O xygen transport .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4.2 The Approximate Analytical solution .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.4. 3 Calculating the diffusion coefficients .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.4. 4 Model setup .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 9 4. Model Results and Discussion ...................................................................... 51 4.1 INTRODUCTION .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 SURFACE KINETIC MODEL .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.3 UP-SCALING TO FIELD CON DITIONS .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3.1 O xygen transport .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5. Conclusions and Recommendations ............................................................ 58 6. References ...................................................................................................... 60 Appendix A: Mineralogy ..................................................................................... 66 Appendix B: Static Procedures ......................................................................... 69 Appendix C: Leach Test Data ............................................................................ 71 Appendix D: PHREEQC Code ............................................................................ 74 Page 7 of 77 List of Figures Figure 2.1: Schematically drawing of the three leach models discussed by Davis and Ritchie, (1986a) ........................................................................................... 20 Figure 2.2: Flow diagram of the ABATE strategy (Usher et al., 2002) ................. 23 Figure 2.3: Schematic representation of an oxidizing particle (Davis, 1983) ....... 26 Figure 3.1: Basic schematically drawing of a HCT after Usher et al. (2002) ....... 30 Figure 3.2: Basic schematically drawing of a HCT after Andre (2009) ................ 31 Figure 3.3: PHREEQC code that calculated the molar amount of dissolved oxygen.................................................................................................................. 33 Figure 3.4: Material balance for an element of volume of the reactor (Levenspiel, p84) ...................................................................................................................... 33 Figure 3.5: Fragment of a PHREEQC code used to model a CSTR by making use of the RATES block (Tiruta-Barna, 2008) ...................................................... 35 Figure 3.6: Conceptual model of a column modelled in PHREEQC by making use of the 1D transport capabilities of the software (Nicholson et al., 2003) 36 Figure 3.7: Conceptual diagram of a waste heap dump depicting the basic approach by Davis and Ritchie (1986a) ........................................................... 42 Figure 3.8: Fragment of the PHREEQC code used to calculate the oxidation rate of pyrite using the SCM kinetics....................................................................... 46 Figure 3.9: Fraction of PHREEQC code used to calculate the rate of sulphate production ........................................................................................................... 48 Figure 3.10: Conceptual illustration of AAS model setup ..................................... 50 Figure 4.1: Sulphate leach rate calculated by the SK model compared to HLC data ...................................................................................................................... 52 Figure 4.2: Calcium leach concentrations calculated by the SK model compared to HLC data ........................................................................................................ 53 Figure 4.3: Magnesium leach concentration calculated by the SK model compared to HLC data ...................................................................................... 54 Figure 4.4: Planar front movement through the heap as a function of time: AAS model compared to the PYROX model ............................................................ 55 Figure 4.5: Oxygen concentration profiles for various time intervals throughout the lifetime of the heap ...................................................................................... 56 Page 8 of 77 Figure 4.6: The integrated sulphate production rate calculated by the AAS model compared to the results obtained from PYROX ................................. 57 Page 9 of 77 List of Tables Table 2.1: List of possible acid-generating salts (Younger, 2000) ......................... 19 Table 3.1: Calculated molar amounts of HLC ........................................................ 39 Table 3.2: Input parameters for the AAS model and PYROX simulation ........... 49 Page 10 of 77 1. Introduction 1.1 Background Several commercial geochemical models are available that has similar capabilities to that of PHREEQC (Parkhurst & Appelo, 2003) , for example Geochemist’s Workbench (Bethke, 1998) and POLYMIN (Molson et al. , 2005). Geochemist Workbench has a friendly user interface and graphical output. The problem with commercial models are the associated cost which in certain cases restrict access to these software packages . It is shown in this dissertation that a reasonable model to describe HLC data can be developed using PHREEQC, which is an open source program available from the U.S. Geological Survey . It is further shown in literature that PHREEQC has the capabilities to model 1D reactive transport ( Appelo & Postma, 2005). The use of PHREEQC to model column experiments is by no means novel; it is a general approach to model batch and flow reactors ( Tiruta -barna, 2008 ). Currently the ABATE strategy (Usher et al. , 2002) or variations thereof is applied in South-Africa to develop a geochemical model to describe the impact of waste heap dumps or backfilling scenarios. Incorporating data obtained from these studies into a site specific geochemical numerical model is in most cases difficult. It is therefore important to model the data obtained from HCL tests and PHREEQC is used for this purpose. The following step is to up-scale the model results to represent field conditions to predict possible future seepage qualities that can be expected from these systems. It is shown from literature that the rate determining step governing the rate of pyrite oxidation can be attributed to diffusive oxygen transport into the heap (Davis, 198 3 ). Modelling this behavior require s solving complex differential equations. Davis (198 3 ) developed an analytical approximation to describe this behavior and can be implemented using PHREEQC. 1.2 Problem Statement The majority of South Africa n coal mines are opencast mines. The mining method results in back filled pits and above ground disposal facilities. Leachate emanating from these disposal sites is saline and could be highly acidic. Currently the standard testing procedure to quantify expected leachate qualities include Acid Base Accou nting, Ne t-acid test, static -and kinetic leaching as outlined in the ABATE stategy (Usher et al. , 2002). Page 11 of 77 It is required to incorporate data gathered from these tests in the development of a geochemical model that can be used to predict future water qualities expec ted from either waste heap dumps or back filled pits. The PHREEQC model will be applied to this problem. 1.3 Aims and Objectives The aim of this dissertation is to investigate and develop methodologies to model geochemical problems applicable to waste heap d umps or discard dumps using PHREEQC (Parkhurst & Appelo, 2003 ). The objectives of the dissertation are listed below: ? Research and develop a methodology to model humidity leach cell behaviour using the PHREEQC code; ? Research methodologies that can be used to describe the physicochemical property governing leaching behaviour of waste heap dumps or discard dumps; ? Validate this methodology by comparing results obtained from this model to PYROX (Wunderley & Blowes, 1997); The objectives listed above forms part of the strategy to develop a model in PHREEQC that utilis es site specific data to model long term leaching behaviour of waste heap dumps. Page 12 of 77 2. Literature Survey 2.1 Introduction Certain mining activities, including gold, ni ckel, copper and coal are known to promote the formation of acid mine drainage (AMD). AMD is commonly caused when sulfide-bearing materials are exposed to oxygen and water (Akcil & Koldas, 2006; Banwart & Malmström, 200 1). Mining activities allows for oxyg en to be introduced in deep geological environments containing minerals that are in a reduced state. The oxidation of iron disulfide (FeS 2), which is ubiquitous in most metal sulfides, is a complex series of reactions ( Appelo & Postma, 2005 ). Section 2.2 includes a simplified set of reactions that is generally accepted to describe the formation of AMD. 2.2 Pyrite Oxidation Pyrite oxidation is explained in this section by means of the following equations: Equation 2. 1 FeS2 + 7 2? O2 + H2O ? Fe 2+ + 2SO4 2? + 2H+ Equation 2. 2 Fe2+ + 1 4? O2 + H + ? Fe3+ + 1 2? H2O Equation 2. 3 ??3+ + 3?2? ? ??(??)3 + 3? + Equation 2. 4 ???2 + 14?? 3+ + 8?2? ? 15?? 2+ + 2??4 2? + 16?+ Page 13 of 77 Younger (2000) suggests that a general reaction equation for the oxidation of pyrite can be written as Equation 2.5. Equation 2. 5 ???2 + 15 4 ?2 + 7 2 ?2? ? ??(??)3 + 2??4 2? + 4?+ In the first reaction, iron disulfide is oxidized, causing the formation of acid and the release of sulfate in the water. The ferrous iron (Fe 2+ ) can further be oxidized to Ferric iron (Fe 3+ ). The ferric iron can then either be hydrolyz ed to form ferric hydroxide (Fe (OH) 3 ), also known as “yellow boy ” , or the ferric iron can act as a catalyst in oxidation of iron disulfide ( Skousen et al., 2002). As discussed later in the literature, the main contributor to acidity in post mining water is metal sulfides. There are several reactions contributing to the buffering of a low pH resulting from the oxidation of pyrite. The most obvious of which is the dissolution of calcite (CaCO 3 ) (Bain et al., 2001) and other carbonate rock forming minerals ( dolomite and aragonite) (Appelo & Postma, 2005). The presence of buffering minerals affects the mobility of dissolved metals, for example the precipitation of gypsum may occur due to the increase of Ca2+ at the calcite dissolution front, thereby providing an upper limit on dissolved Ca2+ and SO 4 2- concentrations (Bain et al. , 2001) It is well known from literature , (Appelo & Postma, 2005; Evangelou, 1995) , that systems exposed to the atmosphere and those isolated from the atmosphere behave differently with regard to the dissolution and precipitation of carbonate minerals. In closed systems the sum of dissolved carbonates species are considered constant and the distribution of the dissolved species is a function of pH. Equation 2. 6 shows the relation of Ca 2+ concentration and CO 2 pressure in the dissolution of CaCO 3 in pure water (Appelo & Postma, 2005). Equation 2. 6 ???2+ = ?10 ?6[???2]/4 3 Page 14 of 77 where m refer to the to the concentration of Ca2+ and P denotes the partial pressure of carbon dioxide. From previous research (Morin et al. , 198 7 ) on areas affected by AMD it is shown that with the lowering of pH , secondary carbonates and hydroxide minerals may dissolve in a consistent sequence (Bain et al. , 2001). With the initial flux of acid water , the first minerals to dissolve near the source are calcite, which will buffer t he water close to a pH of 6. The fo llowing mineral is secondary siderite, buffering the pH near 5. After the depletion of siderite, gibbsite will start to dissolve , that results in a pH close to 4. Lastly any ferrihydrite will go into solution at a pH in the vicinity of 3. Buffering effects of minerals in the system is a critical control in the mobility of metals in AMD discharge (Banwart & Malmström, 2001). Silicate weather ing also contributes to a buffering effect in post-mining water. This is a very slow process and the contribution will be a slow and gradual process (Appelo & Postma, 2005). Page 15 of 77 2.3 Reaction Kinetics 2.3.1 Surface reaction kinetics Pyrite can either dissolve and then be oxidized or can be directly oxidized by either ferric iron or oxygen. The oxidation of pyrite by oxygen is shown in Equation 2. 1. The rate is generally reported to be very slow (Appelo & Postma, 2005). Equation 2. 4 shows the oxidation of pyrite by ferric iron. Reactions 1 and 4 depicted in Equation 2. 1 and Equation 2. 4 respectively, progress simultaneously as long as there is a reaction converting Fe 2+ to Fe 3+ . Equation 2. 2 shows the conversion of Fe 2+ to Fe 3+ . Based on Equations 2. 1 and 2. 4 the rate law for the oxidation of pyrite can be written as follows (Evangelou, 1995) : Equation 2. 7 ? ?[???2] ?? = [(?1(?2) ?1 + ?2(?? 3+)?2](?) O2 and Fe 3+ refer to the partial pressure and concentration of O 2 and Fe 3+ respectively. S is the surface area. The general rate law for the change in solute concentration due to the dissolution or precipitation of a mineral can be written as : Equation 2. 8 ???????? = ? ?0 ? ( ? ?0 ) ? ?(?) where rmineral denotes the overall reaction rate with the units [mol/l/s]. The constant k is the specific rate constant with the units [mol/m 2/s]. A0 is the initial reaction surface area V is the reaction solution volume, m0 the initial moles of solid mineral and m denotes the current moles of solid mineral in the system. The constant n is an empirical value relating to the change in reactive surface area. The function g(C) accounts for the effect of the solution composition on the rate equation. Typical examples are the pH of the solution or the distance from equilibrium. The factor (m/m0) n accounts for the change in surface area where m0 is the initial moles of the solid and m the moles at a given time (Appelo & Postma, 2005). Page 16 of 77 Kinetics reported by Williams and Rimstidt (1994) and also referenced by Appe lo and Postma (2005) for the oxidation of pyrite by O 2 is documented as follows: Equation 2. 9 ???????1 = 10 ?8.19??2 0.5??+ ?0.11 w here rpyrite1 is the overall pyrite reaction rate [mol/m 2/s]. It can be seen that the rate has a square root dependency on the concentration of O 2. This results in a significant effect on the rate for low oxygen concentrations and it has a small effect on the rate for high oxygen concentrations. This effect can be explained by the pyrite surface becoming saturated with O 2 (Appelo & Postma, 2005 ). The rate of pyrite oxidation by ferric iron is described in Equation 2. 10 and the overall rate in Equation 2. 11. Equation 2. 10 ???????2 = 10 ?6.07???3+ 0.3 ???2+ ?0.32 Equation 2. 11 ???????3 = 10 ?8.58???3+ 0.3 ???2+ ?0.47??+ ?0.32 Fe 2+ and H+ compete with Fe 3+ for sorption on the cathodic site of the pyrite, as can be seen by the negative affect the Fe 2+ and H+ concentration has on the rate law. The rates will increase infinitely with the decrease of Fe 2+ and H+ concentrations (Appelo & Postma, 2005 ). To avoid this infinite increase in a numerical model, an inhibitor factor is introduced that approach 1 when approaching a limiting concentration. The equation for the inhibitor factor is shown below: Equation 2. 12 (1 + ???2+ ???? )? The exponent n is an empirical coefficient and ???? is a small limiting concentration. Using this inhibitor factor, Appelo and Postma (2005) propose the following rate equation for the oxidation of pyrite by oxygen and Fe 3+ : Page 17 of 77 Equation 2. 13 ???????2 = 6.3 × 10 ?4???3+ 0.92 (1 + ???2+/10 ?6)?0.43 In the absence of oxygen Equation 2.13 is expressed as: Equation 2. 14 ???????3 = 1.9 × 10 ?6???3+ 0.28 (1 + ???2+ 10?6 )?0.52??+ ?0.3 In many cases the equilibrium concept can give adequate explanation for what is observed in the system. This approach can however not account for the temporal nature of AMD. Appelo and Postma (2005) shows that the rate law f or the oxidation of ferrous iron can be written as follows: Equation 2. 15 ???????4 = ? ????2+ ?? = ? ? ???2+ ? [?? ?]2 ? ??2 The oxidation of pyrite at low pH is very slow, as can be inferred from t he rate law by looking at the quadratic dependence of the OH - concentration expressed in the rate law. For the dissolution of a mineral depicted in Equation 2. 16 , Collon et al., (2006 ) proposed a kinetic rate law for the mode lling of geochemical reactions. This reaction is shown in in Equation 2. 17 : Equation 2. 16 ??? ? ?? + + ?? In Equation 2. 16 , A and B denotes the reactant concentrations and n denotes the stoichiometric coefficient and the ion charge on the right hand side of the equation. Page 18 of 77 Equation 2. 17 ?(?) = ????? ? [?(???) ? (? +)?(???) In Equation 2. 17 , r(t) denotes the reaction rate, ? is a constant that is calculated in Equation 2. 19 , ???? ? represent the reaction constant for precipitation and K is the equilibrium constant. This rate law is based on an approach where two parameter s are used to account for experimental results observed in closed reactor tests. With this approach the reaction rate is the difference between the rate of dissolution and the rate of preci pitation, and can be written as: Equation 2. 18 ?(?) = ????? ? ???? T o account for the variance in accessibility of minerals A nB in natural rock and the induced variation of the apparent kinetic constants, K takes into account the relationship between the equilibrium and kinetic constant ( ? = ????? ? /???? ? ). The specific surface is described by ?. This approach is similar to that expressed in Equation 2. 8 (Appelo & Postma, 2005) and is presented in the equation below: Equation 2. 19 ? = ??[????(?)] 2 3? w here m denotes the mass of the mineral. Fg expresses the influence of the specific surface of a rock volume on the reaction rate and can be calculated by the followin g equation: Equation 2. 20 ?? = 0.532 6?? ???? Page 19 of 77 In the equation above MR represent the mass of rock per litre of water, GR represent the mass fraction of the mineral in the rock and ?? denotes the density of the rock. It is commonly observed that groundwater samples of completely flooded deep mines are much more polluted than when mining operations were still ongoing (Younger, 2000). The reason for this is the formation of ‘acid-generating salts’ according to Younger (2000). Pyrite oxidations start during mining operations in unsaturated areas in the mines. Due to the partially wetted conditions, Equation 2. 5 does not run to completion, causing the formation of intermediate solid phases. These intermediate products include the following minerals as shown in Table 2. 1. Table 2.1: List of possible acid-generating salts (Younger, 2000) Mineral Empirical formula melanterite FeSO 4 ?7H2O römerite Fe 2+ Fe 3+ (SO 4 )?14H2O coquimbite Fe 2(SO 4 ) 3 ?9H2O copiapite Fe 2+ Fe 3+ (SO 4 ) 6 (OH) 2?20H2O potassium jarosite KFe 3+ (OH) 6 (SO 4 )2 2.3.2 Shrinking Core Model A number of mathematical models to describe the kinetics of a heterogeneous reaction have been developed over the years. These models were developed for describing leaching behavior of copper heaps. Davis and Ritchie (1985 a) refer to three primary types that all describe in-situ leaching behavior by focusing on the leaching behavior of the individual particles in the heap (Figure 2. 1). All three th ese modelling approaches assume the particle is spherical and that the reaction front moves radially towards the center of the particle (Davis & Ritchie, 1985 a; Mayer et al. , 2002). Page 20 of 77 Figure 2.1: Schematically drawing of the three leach models discussed by Davis and Ritchie, (1986a) According to Davis and Ri tchie (198 6a ), the three models differ in the assumptions made about the reaction rate compared to the reactant transport rate of the reaction sites. The most well -known and widely used of these models is the Shrinking core model (SCM). This mathematical model describes the kinetics for a solid -gas and solid- liquid system (Levenspiel, 1999 ). The model finds it application in geochemistry by describing the oxidation of pyrite in heap leaching scenarios. For a heterogeneous reaction in which a fluid reacts with a solid the following reaction applies: Equation 2. 21 ?(?????) + ??(?????) ? ????? ???????? + ????? ??????? In the equation above the parameter b represent the stoichiometric coefficient of reactant B. The rate of reaction can be controlled by three mechanisms, th e first of which is diffusion through the gas film. This is not discussed further in this dissertation, as it is assumed that the reaction of the particles occur at the water -solid Page 21 of 77 interface (Mayer et al. , 2002). The second rate determining step is the diffusion of the reactant through the ash layer and can be mathematically described as follows: Equation 2. 22 ? dNA dt = rmineral = 4?r 2D2 dCA dr w here r is the radius of the unreacted core [m] and D2 is the diffusion coefficient of the reactant through the ash layer [m 2/s]. The ash layer or alteration rim refers to the reacted crust that forms around the unreacted core (Figure 2. 1). The third rate determining step is chemical rate controls, and can be mathematically expressed as follows: Equation 2. 23 ? 1 4??2 ??? ?? = ???????? = ?? 4??2 ??? ?? = ??????? w here k’’ is the first order rate constant for the surface reaction and CAg is the gas phase concentration of the reactant diffusing through the ash layer [moles/l]. Mayer et al. (2002) only distinguish between transport -controlled and surface-controlled reactions. Levenspiel ( 1999 ) shows that the abovementioned rate controlling steps can be combined to account for the influence of all three rate controlling steps: Equation 2. 24 ? ?? ?? = ???????? = ??? ??? (???)? ??2 + 1 ??? w here R is the radius of the initial unreacted particle and the density of the particle is denoted by ?B [kg/m 3 ]. Ardejani et al. (2008) used this approach to describe the kinetic oxidation of pyrite . Page 22 of 77 2.4 Predictive Methods In order to implement a successful environmental strategy, the source of pollution risk must be identified and there must be a good understanding of how this risk might evolve at the site in question (Berger et al. , 2000; Linklater et al. , 200 5). According to Nordstrom (2010), several points must be raised when looking at geochemical modelling. Firstly when referring to a model, it does not necessary imply that the model is a computer code or a mathematical approach. It can consist of a well constrained logical proposition with the goal to improve or refine a conceptual model of the problem. Usher et al. (2002) proposed the ABATE strategy as an all -inclusive approach of methodologies to develop a geochemical model as shown in Figure 2. 2. A geochemical model can provide invaluable insight with regards to u nderstanding environmental risk and can assist in the interpretation of field data (Linklater et al., 2005). Page 23 of 77 Figure 2.2: Flow diagram of the ABATE strategy (Usher et al., 2002) 2.4.1 Static test ABA is generally used as a first order classification procedure in acid rock drainage (ARD) prediction methodologies. ABA results provide information about the potential of a sample to generate ARD (Price, 1997); it is a valuable and widely used test that determines both the acid-generating and acid neutralizing potential of a sample (Skousen et al., 2002) . Net Acid Generating (NAG) tests are also performed to aid in a better interpretation of the ABA results (Miller et al. , 1997). The test procedure is not time consuming, is low cost and the analytical procedure is relatively simple (Usher et al. , 2002). Page 24 of 77 There is however limitations to this test procedure. The most important of which is the fact that ABA do es not address the transient behavior of ARD and this is also the case for all the static type tests. Usher et al. (2002) notes that the ABA test procedure assumes instant availability of reactants, simple reaction stoichiometry and account for size effects of solid reactants. 2.4.2 Kinetic tests In order to a ddress the shortfalls of static testing Usher et al. (2002) recommends a kinetic leach test methodology adapted from the D5744 -96 Standard Test Method (ASTM, 2001). The minimum suggested duration of this leach test is twenty weeks and is in some cases not enough time for sulphate production to stabilize. Humidity cell tests (HCT) supplement static tests in this regard, and can assist in obtaining reaction rates and an indication of the potential mineral leaching behavior of rock samples (Price, 1997). Curre ntly one of the biggest challenges faced by environmental practitioners is to integrate site specific data obtained from HCT ’s into predictive geochemical models (Sunkavalli, 2014). 2.5 Modelling of acid mine drainage Mathematical models have become an integral part of the development of a geochemical model and are the last step in the ABATE strategy. Modelling the evolution of water chemistry as it leaches through a waste heap dump can become very complex , as there are multiple processes influencing the change in leachate composition. Geochemical reactions that may occur include aqueous speciation, precipitation and dissolution of minerals, ion exchange and redox reaction (Wunderley et al., 1995; Linklater et al., 200 5) Both equilibrium and reaction path geochemical models can be used to describe the formation of AMD. The y can be linked with a combination of different transport models for CSTR, 1 D, 2D, 3D transport and dynamic or steady state conditions (Andre, 2009). It is seen from literatur e that a purely equilibrium based model is not sufficient to describe AMD (Appelo & Postma, 2005; Andre, 2009). This is due to the slow reaction kinetics of pyrite oxidation. Wunderley et al. (1995) states that when reactions between solid mineral s and reactants in the aqueous phase are rapid it is valid to adopt an equilibrium mass-action approach to the geochemical model. F or this assumption to be valid, the rate of reaction between the solid minerals and the aqueous phase must be infinitely faster than the governing flow phase (Brown et al., 1999). For a heap leach scenario this assumption is not valid, reactions are Page 25 of 77 described kinetically. In these systems there are several rate determining factors. The general approach to the modelling of AMD is assuming that oxygen diffusion is the rate determining step in the oxidation of pyrite. In the well -known work done by Davis and Ritchie (1986a) and later also by Molson et al. (2005) , the shrinking core model is used to describe the reaction kinetics of pyrite oxidation in a waste heap. In this approach it is assumed that the rate of oxygen diffusion of the reaction sites is the rate limiting step. This is also discussed as an important factor in the oxidation of pyrite soils by Appelo and Postma (2005). It is noted that oxygen transport is only important in the unsaturated zone. In the saturated zone , the approximate amount of oxygen available for pyrite oxidation is 0.33 mM O 2, suggesting that the unsaturated zone is the most reactive and the greatest source o f sulphate oxidation. Davis and Ritchie (198 5a) further assume that the rate of unreacted core shrinkage is much slower than the oxygen diffusion rate within the particle. It is further assumed that the reaction at the surface of the particle is instantaneous and that the ash layer diffusion is the rate determining step. In the model developed by Davis and Ritchie (1986a) it is assumed that diffusion is the only means of oxygen transport. Fick’s second law (Equation 2.25 ) with a consumption term is used to describe the diffusion process mathematically: Equation 2. 25 ??? ?? = ?1 ?2?? ??2 ? ?(?, ?) The model describes a two -stage diffusion process to the reaction sites. The first stage consists of ox ygen diffusion through the air filled pores. D1 is the diffusion coefficient for this stage. The consumption term q(z,t) represents the change in volume due to the oxidation of pyrite, and can be expressed as Equation 2.26 (Molson et al. , 2005). Equation 2. 26 ?(?, ?) = ?2 3(1 ? ?) ?2 ( ? ? ? ? ) [?2]? ? Page 26 of 77 w here H is Henry’s constant [kg/m 3 ] and [O2]a is the concentration of oxygen in the pore space. D2 represents the diffusion coefficient for the forming ash layer shielding the shrinking particle. The concentration of oxygen dissolving into the water film surrounding the particle is calculated using Henry’s law, where H represents Henry’s constant for a given temperature. The parameter ? represents the gas filled porosity of the material. Figure 2. 3 show a schematic drawing of the shrinking core mode l. Figure 2.3: Schematic representation of an oxidizing particle (Davis, 1983) It can be seen from Figure 2. 3 that a ash layer forms around the reacting particle core. Oxygen diffuses from surface of the particle through the ash layer to the unreacted core. From Equation 2.26 and Figure 2.3 it can be seen that the change in oxygen due to pyrite oxidation can be related to the change in the unreacted core radius (r): Equation 2. 27 ??? ?? = ?2(1 ? ?) ??? ? ?(? ? ?) [?2]? ? Page 27 of 77 w here ? is the ratio of oxygen to sulphur consumed on the basis of the reaction stoichiometry and rc represents the reaction rate. R represents the initial unreacted core radius and r is the transient unreacted core radius. The density ?S is calculated from the bulk density and the weight fraction of pyrite present. The majority of the most recent geochemical models used to model the oxidation of pyrite in wa ste heap dumps are based on the approach outlined above (Linklater et al. , 200 5; Gerke et al. , 2001; Molson et al. , 2008; Bain et al., 2001; Brown et al. , 1999) . Page 28 of 77 2.6 Modelling spoils heaps Substantial work has been done on the importance of chemical, physical and microbiological factors influencing the weathering rate of pyrite (Linklate r et al. , 200 5, Molson et al., 2008) . Leaching of metals contained in the constituents of mining ore or the leached product or spoils heaps and tailings dams rely on the oxidation of insoluble sul phide bearing minerals to soluble sulphates (Davis & Ritchie, 1986 a). Water transports the soluble oxidation products through the gangue w here further oxidation, dissolution and precipitation can take place (Brown et al., 1999). Thus the concentrat ion of chemical constituents in the effluent of a spoils heap is space and time dependent, which in turn depend on the pyrite oxidation rate, water infiltration rate and the chemical oxidation rates (Pantelis et al., 2001) . Several mathematical models describing the leaching process in spoils heaps, h ave been developed over the last couple of years (Davies & Ritchie, 1986 ; Molson et al., 2005; Linklater et al. , 200 5; Pantelis et al., 2001). In the work done by Davies and Ritchie (1986a) it is assumed that oxygen is transported into the heap by means of diffusion. In the work conducted by Lefebvre et al. , (2001) it is shown that the convective flow of oxygen caused by heat production of pyrite oxidation can also determine the rate of oxygen supply. Pantelis et al. (2001) modelled the spoils heap as a three-phase system consisting of gas and water phases flowing through a rigid solid porous phase. To simplify the mode i t was assumed that the sulphide containing minerals consist entirely of pyrite. There was also no compensat ion for the temperature dependence of the diffusion coefficients, viscosities and thermal conductivities. This was done under the assumption that physical conditions in the spoils heap reach pseudo steady state which can last in most cases for decades or longer. The pH of the spoils heap effluent is controlled by several factors. Strömber and Banwart (1999) identified two main controls, namely sulphide oxidation with calcite dissolution that can sustain a neutral pH and simultaneous oxidation of sulphide and weathering of primary sili cates. The latter is characteris ed by effluent water with a pH between 3 and 4 . Page 29 of 77 3. Development methodology and model code 3.1 Introduction In the previous chapter the approach to develop ing a geochemical model was discussed. It was noted that there are challenges to integrate site specific data obtained from HCT in to a geochemical model. I n this chapter the development of a model will be discussed to address this problem. The model is based on surface kinetics (Surface Kinetic Model). This model will be used to simulate a standard humidity cell leach test (ASTM, 2001). The next step is to scale the model up to field conditions by making use of the of the Approximate A nalytical Solution or AAS method as developed by Davis and Ritchie (1986a) . The model wi ll be simulated using PHREEQC (Parkhurst & Appello, 2003) , and compared to results obtained from PYROX (Wunderley & Blowes, 1997). The mathematical governing equations are discussed and the derivations shown which are used to describe the geochemical system. 3.2 Modelling HCT with the Surface Kinetic Model Several rates that describe the oxidation of pyrite and the factors that influence the rate of reaction were discussed in the previous chapter. Following the steps set out in the ABATE strategy; the next step , following data collection , is the development of a geochemical model. Si te specific data were obtained by means of a HCT. The HCT data were model led using PHREEQC code. PHREEQC has the ability to model dynamic processes such as reaction kinetics, 1D reac tive transport through the use of the RATES block and TRANSPORT block respectively. Simulating this process using PHRREEQC is complex and several approximations are required. T he basic physiochemical parameters governing the system will be discussed first to define the system. Air is introduced into the system by pumping ‘dry air’ into the cell; this is followed by air that has been passed through a humidifier. This usually consists of a three day period for each air cycle. The two air cycles are then followed b y leaching the cell with deionized water. Figure 3 . 1 show s a simplified schematic drawing of the system in question. Page 30 of 77 Figure 3.1: Basic schematically drawing of a HCT after Usher et al. (2002) T he objective of circulating dry and moist air through the system is to accelerate the weathering rate of the sample. The final leach step is then used to remo ve the weathered products from the system. In order to describe the behavior of a HCT by making use of geochemical modelling , several assumptions are required to sufficiently simplify the system. The following main assumptions that were applied are as follows : ? oxygen is assumed to be in excess; ? there are no stagnant zones in the leach column and all the particles are exposed to the deionized water for t he same period of time; and ? the residence time of the leachate in the cell is equivalent to one leach cycle. The assumptions listed above allow the HCT to be approximated as a continuous stirred tank reactor (CSTR). Equation 3.1 describes the residence time of the leachate with in the cell or reactor: Page 31 of 77 Equation 3 . 1 ? = ?????/(?? ? ??) The parameter ? represents residence time and Vtank the leachate volume of the system. Figure 3 . 2 depicts the CSTR approac h to modelling a HCT . The leachate movement through the particles comprising the sample is described as a well -mixed reactor. It is assumed that all the particles in the cell come into contact with the leachate in the form of a thin liquid film forming around the particles. The assumption that oxygen is available in excess is valid as air is pumped into the ce ll and is circulated through the sample material. It is further assumed that the liquid film surrounding the particle is thin, and is negligible compared to the area of the particle exposed to air , which leads to the assumption that there is no rate limiting effects in terms of oxygen transport to the reactive zone. Figure 3.2: Basic schematically drawing of a HCT after Andre (2009) The concentration of O 2 in solution can then be calculated using Henry’s law. Usually the Henry’s law constants are reported at a reference temperature of 25°C and 1 atm. The temperature can be evaluated through the use of Equation 3.2 and Equation 3.3 (Koretsky, 200 3 ): Page 32 of 77 Equation 3 . 2 ( ????? ?? ) = ?? ? ? ?? ? ? Where R is the gas constant, ?? ? are the partial molar enthalpy and ?? ? the pure species enthalpy. A similar equation describes the correction for pressure, but it is beyond the scope of this dissertation as the pressure will always be assumed to be at 1 atm. Assuming the mole fraction of oxygen in ambient conditions is 0.21, the partial pressure of O 2 is then calculated using the following equation: Equation 3 . 3 ??2 = (0.21)? T he concentration of oxygen in solution is calculated through the use of Henry’s coefficient for ox ygen in water ( ??2 = 44252.9 ??? @ 25? ??? 1 ???), and is expressed in the following equation: Equation 3 . 4 ??2 = ??2? ??2 = 0.21[???] 44253.9[???] = 4.75 × 10?6 T he concentration of oxygen [O 2] in molality is expressed as follows : Equation 3 . 5 [?2] = ??2 ? w here V is 1 liter of water and ??2 is the molar amount of O 2 dissolved. The concentration of dissolved oxygen can then be calculated through Equation 3.6. Equation 3 . 6 [?2] = ??2 ??2 ? = 2.63 × 10?4[?] Page 33 of 77 This can be simulated in PHREEQC by making use of the EQUILIBRIUM_PHASES block. The liquid is equilibrated with air at a defined temperature. Figure 3 . 3 shows an example of PHREEQC code that calculates the molar amount of oxygen that dissolves in water at 25°Celcius. The amount of oxygen is expressed as the log value of the partial pressure of oxygen, in this case 0.21 atm. Figure 3.3: PHREEQC code that calculated the molar amount of dissolved oxygen 3.2.1 Developing a methodology to describe a CSTR in PHREEQC Following the methodology described by Tiruta -Barna (2008), a mass balance is conduct over a finite element as shown in Figure 3.4. Figure 3.4: Material balance for an element of volume of the reactor (Levenspiel, p84) The mass balance depicted in Figure 3 . 4 is achieved through the following equation: Page 34 of 77 Equation 3 . 7 ( ??? ?? ) = [ ? ????????? (??? ? ? ??)] 1 + [?????]2 where c’ denotes the concentration of the reactant c, Vleachate refer to the volume of solute in the system. The term numbered with subscript 1 is the convection term and the term 2 is the reaction kinetics. The CSTR equation is modelled in PHREEQC by solving the differential equation in the RATES block. This is accomplished through the use of PHREEQC`s Runge -Kutta integrator. A fragment of the code is depicted in Figure 3 . 5 to solve the following differential equation (Tiruta -Barna, 2008): Equation 3 . 8 ?? , ?? = ? ????????? (??? ? ? ??) The analytical solution to Equation 3.8 is given as: Equation 3 . 9 ?? = ?0 ???? (? ?? ????????? ) where c’ is the current concentration of reactant in the reactor, Q is the volumetric inflow rate of liq uid into and out of the reactor and Vleachate is the volume of the reactor. The solution obtained from PHREEQC is compared to the analytical solution in Figure 3.5 . Page 35 of 77 Figure 3.5: Fragment of a PHREEQC code used to model a CSTR by making use of the RATES block (Tiruta-Barna, 2008) Now that a solution for the convection term of the CSTR exists , the reaction kinetic term of Equation 3 . 7 needs to be addressed. It is proposed to model the convection term using the TRANSPORT block in PHREEQC and model the reaction kinetics using the RATE block. The TRANSPORT block has the capability to call up the calculations of the RATE block to simulate the rate kinetics. The TRANSPORT block also has the capability to simulate multicomponent diffusion in an aqueous solution. All chemical processes in PHREEQC can be included in the transport s imulation (Parkhurst & Appelo, 2003 ). In the work conducted by Appelo et al. (1997) this methodology is used to model the evolution of leachate through a column containing pyrite and marine sediment. Nicholson et al. (2003) also used the 1D transport capab ilities to model waste rock piles containing uranium for closure purposes. In the abovementioned studies the column is modelled by dividing the column into a number of nodes. Figure 3 . 6 depicts a conceptual drawing of the modelling approach used in PHREEQC to model 1D reactive transport. Page 36 of 77 Figure 3.6: Conceptual model of a column modelled in PHREEQC by making use of the 1D transport capabilities of the software (Nicholson et al., 2003) Stagnant cells can be added to the ‘active cells’ to model stagnant zone or radial diffusion. This function ality can be used to implement a CSTR. In the TRANSPORT block stagnant cells can be defined that can be used to interact with ‘active’ transport cells using the MIX block t o model the Q value defined in E quation 3.7 . Reaction kinetics and solution composition can then be defined for in dividual ‘active cells’ or nodes. B oundary conditions can also be implemented (Parkhurst & Appelo, 2003 ). The first is the Dirichlet boundary condition where the concentration is constant and known: Equation 3 . 10 ?(????, ?) = ?0 The second boundary condition is the no flux boundary condition or Neumann boundary condition: Page 37 of 77 Equation 3 . 11 ??(???? , ?) ?? = 0 The third boundary condition is known as the Cauchy boundary condition or flux boundary condition: Equation 3 . 12 ?(????, ?) = ?0 + ?? ? ??(???? , ?) ?? where DL is the dispersion coefficient [ m2/s ] . The Dirichlet - and Cauchy boundary conditions are then used in an ‘active cell’ to mix with a stagnant cell where the flow direction defined in the TRANSPORT block is diff usion only. The flow rate into and out of the reactor (the ‘active cell’) can then be set by the shifts input. The shifts input define the rate at which the leachate moves through each cell . This methodology was adopted from a similar approach developed by Tiruta -Barna et al. (2006). 3.3 Model setup in PHREEQC 3.3.1 Defining the reaction kinetics Following the methodology proposed by Appelo and Postma (2005) the reaction rate for the dissolution/precipitation can be modeled by combining Equation 2. 8 and Equation 2. 9 resulting in the following equation : Equation 3 . 13 ? = ?0 ( ? ?0 ) ? (10?8.19??2 0.5??+ ?0.11)(1 ? ???/?) The ion activity product, IAP, has an influence on the rate of reaction and has an important influence on the rate as the solution approach equilibrium. For the dissolution of gypsum the same methodology was followed, rendering the following rate law for the dissolution of gypsum Page 38 of 77 Equation 3 . 14 ?[??] ?? = 3.1 ? 10?8?0(1 ? ???/?) The kinetic dissolution of dolomite can be described by the following rate equ ation (Appelo & Postma, 2005) : Equation 3 . 15 ? = ??? ?? = ??? ? ln ( ?? ???? ) This is a very elementary rate equation and several more complex rate equations were developed that describe the dissolution of dolomite under different condition (Brown et al. , 2003; Zhang et al. , 2007; Pokrovsky et al., 2009) . Partial pressure of carbon dioxide, temperature and concentration of magnesium can influence the dissolution rate of dolomite. Several kinetic rate laws are defined in the default PHREEQC database. These rates can be used in the simulations through the KINETICS block. 3.3.2 Equilibrium phases Several minerals dissolve or precipitate fast relative to the oxidation of pyrite or the dissolution of silicate minerals. The reactions of these minerals can then be described as equilibrium controlled as they happen instantaneous compared to the time scale of minerals like pyrite to react. Common equilibrium reactions found in waste heap s that is likely to occur, is the dissolution and precipitation of gypsum (Nicholson et al., 2003). In PHREEQC the EQUILIBRIUM_PHASES block is used to define these minerals. It is assumed that these reactions occur instantaneous an d is controlled by the solute composition and can be described by the law of mass action : Equation 3 . 16 ?? + ?? ? ?? + ?? Page 39 of 77 Equation 3 . 17 ? = [?]?[?]? [?]?[?]? Secondary minerals that may precipitate are also defined in the EQUILIBRIUM_PHASES block. 3.3.3 Calculating composition from mineralogy The composition of the HLC column is determined by making use of the mineralogy data obtained from the XRD analysis that can be found in Appendix A. A mass of 1000 grams was used in the HLC. The XRD data provides a mass percentage for the minerals present in the sample. From this the molar amounts of each mineral present can be calculated as shown in Table 3 . 1. Table 3.1: Calculated molar amounts of HLC Composition Mineral Amount (weight %) moles Anatase 5.45 0.682 Calcite 6 .22 0.621 Dolomite 15.17 0.823 Kaolinite 61.6 2.386 Pyrite 5.43 0.453 Quartz 6 .13 1.020 It is assumed that quartz will not react in the HLC as the reaction of this mineral is very slow relati ve to the duration of the test and therefore is not considered in the model. 3.3.4 Defining the reaction surface The reactive surface area of the minerals in the HLC is unknown and will be used as a fitting parameter. The reaction surface can be approximated using the following equation: Equation 3 . 18 ? = ( 6 ? ) ? Page 40 of 77 where d is the diameter of a spherical particle and ? is surface roughness factor. The latter is defined as the ratio of the true reactive surface area to the equivalent geometric surface area. 3.4 Up-Scaling to Field conditions The second problem that needs to be addressed in is the integrati on of data obtained from static and kinetic tests into a geochemical model that describe field conditions. For this purpose , methods for modelling physical and chemical processes described in literature must be implemented in PHREEQC. The discussion b elow describes the conceptualisation and development of theory that are widely used today in geochemical models to describe waste heap dumps ( Pantelis et al. , 200 1). Ritchie (1977) developed a mathematical model to better understand the rate limiting mechanisms of pyrite oxidation in waste heap dumps. The model also aimed to address the various parameters governing the active transport behavior in the heap. The model considere d the heap to be a homogeneous semi-infinite slab (Davies, 1983 ). It was assumed that oxygen transport to the reaction sites in the heap was the main rate determining mechanism. All reactants besides oxygen were assumed to be available, including bacteria catalyzing the oxidation reaction of pyrite. It was seen that the catalyzing effects of bacteria had little effects on the model as the rate of reaction was already assumed to be controlled by the availability oxygen in the heap. For this reason Davis and Ritchie (1986a) adopted the planar moving boundary approach to model the reactions in the waste heap dump . It was assumed that the diffusion of oxygen through pore space of the heap was the rate controlling step and therefore it implied that the oxygen ar riving at the reaction sites reacted instantly with the pyrite, resulting in a zero oxygen concentration at this point. The mass flux of oxygen through the reaction boundary in the model developed by Ritchie (1997) was assumed to be the velocity of the moving front, depending explicitly on the density of pyrite in waste dump, the ratio of oxygen consumed to pyrite reacted in the oxidation reaction and the diffusion coefficient of oxygen within the pores pace of the dump (Davies, 1983 ). Page 41 of 77 Davies (1983 ) extend ed this model by implementing a two stage model. The first stage is the diffusional transport of oxygen from the atmosphere, through the pore space, to the particles in the heap. The second stage describes the transport of oxygen to the reaction sites thro ugh the ash layer. The model proposed by Davies (1983 ) and later published in as a series of three very insightful articles (Davis & Ritchie, 1986a, b, c ) will be followed in an attempt to model waste heap dumps and the backfilled pits by incorporating data obtained from static and kinetic tests as outlined in the ABATE strategy described in C hapter 2. 3.4.1 Oxygen transport In Chapter 2 the mass balance in the pore space of a dump governed by diffusive transport was discussed . In order to use this equation, it must be assumed that all the particles are spherical and identical in size. The mass balance for oxygen in the spherical particle is given by the following equation (Molson et al., 2005): Equation 3 . 19 ?? ?? = ?2 3(1 ? ??) ?2 ( ? ? ? ? ) [?2]? ? where D2 is the effective diffusion coefficient, incorporating the diffusion properties of oxygen in water through the ash layer. The last term in the equation denotes the concentration of oxygen in pore space, where H represent Henry’s coefficient. This term calculates the molar concentration of oxygen in the liquid film surrounding the particle. R represents the particle ’s unreacted radius at time t = t 0. Page 42 of 77 3.4.2 The Approximate Analytical solution The conceptual model of a waste heap dump is shown in Figure 3 . 7 . Figure 3.7: Conceptual diagram of a waste heap dump depicting the basic approach by Davis and Ritchie (1986a) In order to model the waste heap dump as described by Davis and Ritchie (1986a) a pseudo-steady state approximation for the particles in the dump is assumed, rendering a set of two equations: Equation 3 . 20 ?1 ?? ?? = ?2? ??2 ? 3? ?? [1 ? ?] ??? 0 < ? < 1 Equation 3 . 21 ?? ?? = ? ?? ?(1 ? ?) ??? 0 < ? < 1 The symbol u represent the oxygen concentration, R the position of reaction front within the particle [m]. T he following boundary and initial conditions apply to the above equations: ?(0, ?) = 1 Page 43 of 77 ?? ?? = (1, ?) = 0 ?(?, 0) = 0 ?(?, 0) = 1 This set of equations needs to be solved numerically. Davies and Ritchie, (1986 a) proposed an analytical solution to this system. The following analytical equations form part of the solution to the above set of differential equations: Equation 3 . 22 ?? = 1 6? This equation is used to calculate the time it takes for the reaction front X(t) to reach the top of the waste heap dump, i.e. the time it takes for the particles on the surface of the heap to be fully oxidized as described by the shrinking core model. The constant k is calculated as follows: Equation 3 . 23 ? = ?1 3 The constant k1 is calculated by the following equation: Equation 3 . 24 ?1 = 3??2(1 ? ?)? 2 (?1?2) where the ? denotes a proportionality constant that include both the Henry’s law and gas law. The term D2 refers to the diffusion coefficient of oxygen through the ash layer, L denotes the total height of the heap. The term ? denotes the porosity of the Page 44 of 77 waste heap dump. D2 is the bulk diffusion coefficient for oxygen transport through the pore space of the heap and symbolises the particle size. The movement of the reaction front X(t) through the waste heap dump as a function o f time is calculated with the following equation: Equation 3 . 25 ?(?) = ?2? ? ?? ???? This now allows for the calculation of the oxygen profile as a function of time and can be expressed as follows: Equation 3 . 26 ?(?, ?) = 1 ? ? ?? [1 + ???] ??? 0 ? ? ? ?(?) w here ? = 6? and x denotes the location in the heap above the moving reaction front X(t). The oxygen profile below the reaction front X(t) is calculated by the following equation: Equation 3 . 27 ?(?, ?) = exp (???(? ? ?)) [1 + ???] ??? ?(?) ? ? ? 1 It can be shown that Equation 3 . 19 can be solved in PHREEQC by making use of the RATES block as discussed in the previous section. For illustration purposes a fragment of the code to model the kinetic oxidation of pyrite as de scribed by Mayer (1999) is depicted in Figure 3 . 8 . The following equation from the study by Mayer (1999) are modeled: Page 45 of 77 Equation 3 . 28 ?? ?? = ?103? ?2 3.5 ? (? ? ?)? [?2]?? where D2 is the diffusion coefficient of oxygen in water, r is the unreacted core radius, R is the particle radius (unreacted core plus ash layer) and S is the reactive surface and can be calculated as follows: Equation 3 . 29 ? = 3 ? 10?3 ?? ? w here ?? is the volumetric fraction of mineral at time t =t i defined as [ m 3 mineral m-3 bulk]. This parameter changes with time, compensating for the reduction in reaction surface as the oxidation of pyrite progresses . Figure 3 . 8 shows t he PHREEQC and MATLAB solutions to the equations above. The figure also shows a fragment of the PHREEQC code used to solve the equations. Page 46 of 77 Figure 3.8: Fragment of the PHREEQC code used to calculate the oxidation rate of pyrite using the SCM kinetics Page 47 of 77 The dissolved oxygen is defined usi ng the EQUILIBRIUM_PHASES block in PHREEQC due to the fact that PHREEQC do not have a gas transport module. For this reason it is proposed to make use of the integrated sulfide production rate proposed by Davis and Ritchie (1986a) : Equation 3 . 30 ?(?) = ????3 ? 1 ???[? ? (?(?)2/2] ? ?? 6?2(1 ? ?)2(1 + 2?) [(1 ? ?)2(1 + 2?)2 ? ?????2??(1 ? ?(?))] ?? 0 ?? where L is the total height of the waste heap, ?S is the ratio of mass of oxygen consumed per mas of sulfide generated in the pyrite oxidation reaction. The constant k3 is calculated as follows: Equation 3 . 31 ?3 = 3???? ?4 where ?s denotes the density of pyrite in the heap and ?4 is the characteristic time for transport. F or the calculation of this constant, the reader is referred to the work done by Davis and Ritchie (1986a) . For small particles it can be shown that the rate expressed in Equation 3.37 can be simplified as follows: Equation 3 . 32 ?(?) = ????? ?4 1 ?2? ? ?? This render s an expression that is only time dependent and can be solved with ease in PHREEQC, eliminating the problem of compensating for the oxygen transport when the shrinking core model is used. The RATES block in PHREEQC will automatically integrate the rate expression when the standard procedure is used to Page 48 of 77 input a rate expression in PHREEQC. For this reason the keyword total_time is used to recall the current time for each time step. This value then defines the time step t in Equation 3 . 32 , al lows the RATES block to simulate a rate expression that is only time dependent. Figure 3 . 9 shows the input to the rate block used to simulate the rate discussed above. Figure 3.9: Fraction of PHREEQC code used to calculate the rate of sulphate production 3.4.3 Calculating the diffusion coefficients The bulk diffusion coefficient D1, governing the transport of oxygen through the pore space of the heap is a function of void space, degree of saturation of the heap and temperature. The following equation can be used to calculate the bulk diffusion coefficient (Reardon & Moddle, 1985) : Equation 3 . 33 ?1 = 3.98 × 10 ?9 ? [ ?? ? 0.05 0.95 ] ? ?1.5 w here ?? is the air filled porosity and T denotes temperature in Kelvin. Page 49 of 77 3.4.4 Model setup A test scenario is set up to compare results obtained from the AAS mo del with results obtained from PYROX. A homogenous slap with a height of 1 meter is modelled. It is assumed that the top surface area is sufficiently wide that it can be assumed that oxygen is transported from the top surface b y means of diffusion as described in the previous chapters. It is assumed the particles are al l spherical and identical in size with a radius of 0.001515 meter. The remaining input parameters is listed in Table 3 . 2, where D 2 is calculated using Equation 3 . 33 . Table 3.2: Input parameters for the AAS model and PYROX simulation Parameter Unit Description L 1 m Depth of heap ? 0.38 porosity D1 3.31E -06 m2 /s Diffusion Coefficient of oxygen in pore space of dump D2 2.60E -09 m2 /s Diffusion Coefficient of oxygen in water a 0.001515 m radius of particle v 1.746 mass of oxygen used per mass of pyrite oxidized in the oxidation reaction rs 45 kg/m 3 Density of sulphur within the dump ? 0.03 Constant U0 0.265 kg/m 3 Oxygen concentration at the surface of the dump ? 22460 J/kg Heat produced from oxidation reaction per mass of sulphur oxidation (Davis, 1983 ) ?s 1.616 mass ratio, not mole ratio - error The model calculates the sulphate production on an area basis of 1 m2 and can be conceptualized as depicted in Figure 3 . 10. Page 50 of 77 Figure 3.10: Conceptual illustration of AAS model setup Page 51 of 77 4. Model Results and Discussion 4.1 Introduction In this chapter the results of the developed model will be presented and discussed. 4.2 Surface Kinetic model The data obtained from the HLC were modelled by making use of the SK model. In Figure 4 . 1 the release of sulphates were modelled. It can be seen from the data that there is initially a very high concentration of sulphate in the leachate, this decrease s exponentially ove r the course of the twenty week period the HCL test was conducted . The initial high value of sulphates observed can be attributed to the dissolution of gypsum and not so much to the oxidation of pyrite, for this reason gypsum was defined as a kinetic species and not as an equilibrium species. The amount of gypsum initially present was unknown, and was used as a calibration parameter to calculate the SO 4 2- release rate. The dissolution of gypsum is responsible for the release of SO 4 2- and Ca2+ . This complicates the calibration of the SO 4 2- leach rate because adjusting the amount of gypsum present also affects the Ca2+ release rate. Dolomite and calcite were the only two bicarbonate minerals present in the material. Both dolomite and calcite release Ca2+ when they dissolve. The dissolution rate of dolomite wa s calibrated using the Mg2+ , and in turn the calcite and gypsum release rates we re calibrated with the SO 4 2- and Ca2+ leach data. A reasonable model fit is obtained in the first 13 weeks, but over pr edicts the sulphate release thereafter. This can be attributed to the gypsum rate under predicting the speed of dissolution in the early stages of the leaching procedure as it is anticipated that the initial gypsum present in the system would have leached out in the first two to three weeks. Page 52 of 77 Figure 4.1: Sulphate leach rate calculated by the SK model compared to HLC data Figure 4 . 2 show the leach rat e of Ca2+ predicted by the surface kinetic model compared to HLC data. It can be seen that the model does not adequately predict the release of calcium in the first 13 weeks . The over prediction from the model in the first 8 weeks may be attributed to the contribution of gypsum. Figure 4 . 3 show the release rate of Mg2+ calculated by PHREEQ compared to HLC data. The figure show a similar trend to that of Figure 4 . 2 suggesting that the dolomite rate law may not be adequa te in describing the dissolution of dolomite. This may be contributed to the over prediction of calcium in the system. Page 53 of 77 Figure 4.2: Calcium leach concentrations calculated by the SK model compared to HLC data The almost sinusoidal oscillation of the leach data observed in all the parameters described above can be attributed to the leach test procedure with the cycling of dry and moist air and leaching with water only once a week. The PHREEQC model does not take this behaviour into account and thus explains the deviation of the model with the HLC data at weeks 6 -7 and weeks 12 -13. Page 54 of 77 Figure 4.3: Magnesium leach concentration calculated by the SK model compared to HLC data 4.3 Up-Scaling to field conditions 4.3.1 Oxygen transport Figure 4 . 4 below depicts the movement of the planar reaction from the top surface of the heap through the length of the heap as a function of time. The Approximate Analytical S olution is used to calculate the movement through the heap, the results from this model is compared to the results obtained from PYROX for the same conditions. Page 55 of 77 Figure 4.4: Planar front movement through the heap as a function of time: AAS model compared to the PYROX model Figure 4 . 5 shows the oxygen profiles as a function of time as calculated by Equation 3 . 25 and Equation 3 . 26 . This result is compared to results obtained from PYROX. Page 56 of 77 Figure 4.5: Oxygen concentration profiles for various time intervals throughout the lifetime of the heap Figure 4 . 6 shows the integrated sulphate production rate of sulphate compared to the integrated sulphate production rate calculated by PYROX for the same conditions. It can be seen from Figure 4 . 4 that in the initial stage of the simulation the AAS model over predicts the pore space oxygen concentration. This effect can also be observed in Figure 4 . 4 . The over prediction of oxygen is directly linked to the under prediction in sulphate production rate, as the pyrite is not oxidizing as fast in the initial stages of the AAS model compared to the PYROX model. The AAS model predicts the time for the heap to be fully oxidized to take around one and a half years, the same as PYROX. Page 57 of 77 Figure 4.6: The integrated sulphate production rate calculated by the AAS model compared to the results obtained from PYROX Page 58 of 77 5. Conclusions and Recommendations It was shown in this study that PHREEQC can be used to develop a reasonable model to describe the leach behaviour of a HLC. The results obtained from the model were compared with data obtained from a standard humidity leach cell as outlined by Usher et al. (2002). Reaction kinetics obtained from literature was used to describe the reactions anticipated in the HLC an d was found to describe the trends observed in the HCL data to a reasonable extent. Deviations of the model from observed data was attributed to the possible presence of gypsum in the system. This caused complications in describing the release of magnesium and calcium. The kinetic rate law used to describe the dissolution of dolomite was also simplified. This may have introducing a possible error into the model. It was seen from literature that diffusive oxygen transport into a waste heap dump is in most cases the rate determining step in the acid generation of a waste heap dump. Several mathematical models were identified that describe this behaviour. An analytical approach (AAS model) developed by Davis (198 3 ) was used to model the oxygen transport to th e reaction sites in the heap and the subsequent production of sulphates. The results obtained from this model were compared to the results obtained from PYROX. It was shown that the integrated sulphate production rate obtained from the AAS model can be mod elled in PHREEQC. This provides then a usable model that can be applied to model waste heap dump scenarios. It was found from the modelling results and comparison with PYROX that the model under predicts the integrated sulphate production rate in the initial stages of the reacting waste heap dump. It does however show results that are in close agreement with the results obtained from PYROX in later stages of the lifespan of the waste heap dump. This provides then limitations to the AAS model’s applicability on geochemical problems. The model can only be applied to describe waste heap dumps where the particles at the top of the heap are fully oxidized. Page 59 of 77 Several shortfalls were identified during the development of the model; the most prominent problems will b e listed below with recommendations to address these issues: ? Reaction rate laws found in literature may not adequately describe the kinetic dissolution of minerals. This may introduce errors in the model when the kinetic rate laws are calibrated to fit the data. It is therefore recommended additional XRD tests be conducted of the HLC sample on completion of the leach test in order to perform a mass balance to validate the calibrated reaction kinetic rate laws in terms of convers ion of the dissolution reaction. ? It is assumed in this work that a HLC can be approximated by a CSTR reactor model. This may be an over simplification as the CSTR approach assumes the system is well mixed and do not account for stagnant zones in the system. This may be rectified by ma king uses of a reactor in series approach to account for stagnant zones. Page 60 of 77 6. References AKCIL, A. & KOLDAS, S. 2006. Acid Mine Drainage (AMD): causes, treatment and case studies. Journal of Cleaner Production 14:1139 -1145 . ANDRE, B.J. 2009. Generation of acid mine drainage: reactive transport models incorporating geochemical and microbial kinetics, PhD thesis: University Of Colorado At Boulder . APPELO, C. A . , VERWEIJ, J.E. & SCHAFER, H. 1998. A hydrogeochemical transport model for an oxidation experiment with pyrite/calcite/exchangers/organic matter containing sand. Appl. Geochem., 13:257 –268 . APPELO, C.A.J. & POSTMA, D. 2005. Geochemistry, groundwater and pollution. Boca Raton: CRC Press. ARDEJANI, F. D. , SHOKRI, B.J. , MORADZADEH, A. , SOLEIMANI, E. & JAFARI, M.A. 2008. A combined mathematical geophysical model for prediction of pyrite oxidation and pollutant leaching associated with a coal washing waste dump. Int. J. Environ. Sci. Tech., 5 (4):517 -526 . AST M. 2001. Standard test method for accelerated weathering of solid materials using a modified Humidity Cell, D 5744 -96ASTM (2007) Standard test method for laboratory weathering of solid materials using a humidity cell, D5744 -07, pp. 1 –19. BAIN, J.G. , MAYER , K.U. , BLOWES, D.W. , FRIND, E.O. , MOLSON, J.W.H. , KAHNT, R. & JENK, U. 2001. Modelling the closure -related geochemical evolution of groundwater at a former uranium mine. Journal of Contaminant Hydrology 52,109 – 135 . BANWART, S.A. & MALMSTR?M, M.E. 2001. Hydrochemical modelling of preliminary assessment of mine water pollution. Journal of Geochemical Exploration 74, 73 -97 Page 61 of 77 BERGER, A.C. , BETHKE, C.M. & KRUMHANSL, J.L. 2000. A process model of natural attenuation in drainage from a historic mining district. A pplied Geochemistry 15:655 -666 . BETHKE C.M., 2008. Geochemical and Biogeochemical Reaction Modeling. Cambridge University Press, 547 pp. BROWN, P. , LUO, X.L. , MOONEY, J. & PANTELIS, G. 1999. The modelling of flow and chemical reactions in waste piles. In: 2nd Internat. Conf. CFD in the Minerals and Process Industries. CSIRO, Melbourne, Australia, December 6 –8, pp. 169 -174. BROWN, J.G. , PIERRE D. & GLYNN, P.D. 2003. Kinetic dissolution of carbonates and Mn oxides in acidic water: measurement of in situ field rates and reactive transport modelling. Applied Geochemistry 18:1225 –1239 . COLLON, P. , FABRIOL, R. & BUÈS, M. 2006. Modelling the evolution of water quality in abandoned mines of the Lorraine Iron Basin. Journal of Hydrology, vol. 328, pp 620 – 634 . DAVIS, G.B. & RITCHIE , A.I.M. 1986(a). A model of oxidation in pyritic mine waste: part 1: equations and approximate solution, Applied Mathematical. Modelling, vol. 10, pp. 3 14 –322. DAVIS, G.B. & RITCHIE , A.I.M. 1986(b). A model of oxidation in pyritic mine waste: part 2: comparison of numerical and approximate solution, Applied Mathematical. Modelling, vol. 10, pp. 314 –322. DAVIS, G.B. & RITCHIE, A.I.M. 1986(c). A model of oxidation in pyritic mine waste: part 3: import of particle size distribution, Applied Mathematical. Modelling, vol. 10, pp. 314 –322. DAVIS, G.B. 1983. Mathematical modelling of rate -limiting mechanisms of pyritic oxidation in overburden dumps, PhD thesis: The University of Wollongong . Page 62 of 77 EVANGELOU. V.P. 1995. Pyrite Oxidation and it’s Control: solution chemistry, surface chemistry, acid mine drainage (AMD), molecular oxidation mechanisms, microbial role, kinetics, control, ameliorates and limitations , microencapsulation. 1st ed. CRC Press. GERKE, H. H., MOLSON, J.W . & FRIND, E.O. 2001. Modelling the impact of physical and chemical heterogeneity on solute leaching in pyritic overburden mine spoils. Ecological Engineering 17: 91 –101. KORETSKY, M.D. 20 03. Engineering and Chemical Thermodynamics, USA: John Wiley & Sons . LEFEBVRE, R. , HOCLLEY, D. , SMOLENSKY, J. & GÉLINAS, P. 2001. Multiple transfer processes in waste rock piles producing acid mine drainage 1: Conceptual model and system characterization. Journal of Contaminant Hydrology 52: 137 -164. LEVENSPIEL, O. 1999. Chemical Reaction Engineering. 3d ed. New York: Wiley. LINKLATER, C.M. , SINCLAIR, D.J. & BROWN, P.L. 2005. Coupled chemistry and transport modelling of sulphidic waste rock dumps at t he Aitik mine site, Sweden. Applied Geochemistry 20:275 –293 . MAYER, U. K. 1999. A numerical model for multicomponent reactive transport in variably saturated porous media. Ph.D. Thesis, Dept. of Earth Sciences, University of Waterloo, Waterloo, Ontario, Canada. MAYER , U.K. , FRIND, E. O. & BLOWES, D.W. 2002. Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions. Water Resources Research, vol. 38 . MILLER, S., ROBERTSON, A. & DONAHUE, T. 1997. Advances in acid drainage prediction using the net acid generation (NAG) 4th international conference on acid rock drainage. Page 63 of 77 MOLSON, J.W. , FALA, O. , AUBERTIN, M. & BUSSIÈRE, B. 2005. Numerical simulations of pyrite oxidation and acid mine drainage in unsaturated waste rock piles. Journal of Contaminant Hydrology 78:343 -371 . MOLSON, J., FALA, O., AUBERTIN, M. & BUSSIÈRE, B. 2008. Geochemical transport modelling of acid mine drainage within heterogeneous waste rock piles. In Proceedings of the 61 st Canadian Geotechnical Conference and 9th Joint CSG/IAHCNC Groundwater Specialty Conference, Edmonton, Alta., 21 -24 September 2008. BiTech Publishers Ltd., Richmond, B.C. pp. 1586 -1593 . MORIN, K.A. , CHERRY, J.A. , DAV E, N.K. , LIM, T.P. & VIVYURKA, A.J. 1987. Migration of acidic groundwater seepage from uranium -tailings impoundments: 1. Field study and conceptual hydrogeochemical model. Journal of Contaminant Hydrology 2: 217 –303 . NICHOLSON, R.V. , RINKER, M.J. , ACOTT, G. & VENHUIS, M.A. 2003. Integration of field data and a geochemical transport model to assess mitigation strategies for an acid-generating mine rock pile at a uranium mine. Proceedings, Sudbury: Mining and the Environment. NORDSTROM, D.K. 2010. Advances in the Hydrogeochemistry and Microbiology of Acid Mine Waters. International Geology Review 42:6, 499 -515. PANTELIS, G. , RITCHIE, A.I.M. & STEPANYANTS, Y.A. 2001. A conceptual model for the description of oxidation and tra nsport processes in sulphidic waste rock dumps. Applied Mathematical Modelling 26: 751 -770 . PARKHURST, D.L. & APPELO, C.A.J. 2003. User’s guide to PHREEQC (version 2)— a computer program for speciation, batch -reaction, one -dimensional transport, and inverse geochemical calculations. Report 99 -4259. US Geological Survey Water-Resources Investigations. Page 64 of 77 POKROVSKY, O.S. , GOLUBEV, S.V. , SCHOTT, J. & CASTILLO, A. 2009. Calcite, dolomite and magnesite dissolution kinetics in aqueous solutions at acid to circumneutral pH, 25 to 150 °C and 1 to 55 atm pCO2: New constraints on CO2 sequestration in sedimentary basins. Chem. Geo., 265:20 –32 . PRICE, W.A. 1997. Guidelines and recommended methods for the prediction of metal leaching and acid rock drainage at minesi tes in British C olumbia (DRAFT). British Columbia Ministry of Employment and Investment, Energy and Minerals Division, Smithers: BC, pp.143. REARDON, E.J. & MODDLE, P.M., 1985. Gas diffusion coefficient measurements on uranium mill tailings: Implications to cover layer design. Uranium 2: 111 –131. SKOUSEN, J. , SIMMONS, J. , MCDONALD, L. M. & ZIEMKIEWICZ, P. 2002. Acid – Base Accounting to Predict Post -Mining Drainage Quality on Surface Mines. J. Environ. Qual. 31:2034 –2044. STRÖMBERG, B. & BANWART, S.1999. Weathering kinetics of waste rock from the Aitik copper mine, Sweden: scale dependent rate factors and pH controls in large column experiments. Journal of Contaminant Hydrology 39: 59 –89. SUNKAVALLI , S.P. 2014. Gap between Humidity Cell Testing Data and Geochemical Modeling of Acid Rock Drainage. ISSN: 2157 -7587 HYCR, an open access journal . TIRUTA -BARNA, L., RAKOT O ARISOA, Z. & MEHU, J., 2006. Assessment of the multi-scale leaching behaviour of compacted coal fly ash. J. Hazard. Mater. 137:1466 –147 . TI RUTA -BARNA, L. 2008. Using PHREEQC for modelling and simulation of dynamic leaching tests and scenarios. Journal of Hazardous Materials 157:525 –533 . Page 65 of 77 USHER, B.H., CRUYWAGEN, L -M., DE NECKER, E. & HODGSON, F.D.I. 2002. On - site and laboratory investigations of spoil in opencast collieries and the development of acid-base accounting procedures, Water Research Commission. WUNDERLEY, M.D. , BLOWES, D.W. , FRIND, E.O. , PTACEK. C.J. & AL, T.A. 1995. A Multicomponent reactive transport model incorporating kinetically controlled pyrite oxidation. Convergence on Mining and the Environment, Sudbury, Ontario. WUNDERLEY, M.D. & BLOWES, D.W. 1997. PYROX version 1.1 User Manual. Institute for Groundwater Research, Unvi. Of Waterloo, Waterloo, Ontario. YOUNGER , P. 2000. Predicting Temporal changes in total iron concentrations in groundwater flowing from abandoned deep mines: a fi rst approximation. Journal of Contaminant Hydrology 44 : 47 –69 . ZHANG, R. , HU, S. , ZHANG, X. & YU, W. 2007. Dissolution Kinetics of Dolomite in Water at Elevated Temperatures. Aquat Geochem 13:309 –338 . Page 66 of 77 Appendix A: Mineralogy INTRODUCTION The mineralogical composition of the rock samples was determined by means of X - ray Diffraction (XRD). The XRD was performed at Waterlab (PTY) Ltd, Pretoria. The total element analyses were performed by means of X -ray Fluorescence (XRF) at the Council for Geoscience, Pretoria. METHODOLOGY The material was prepared for XRD analysis using a backloading preparation method. It was analysed with a PANalytical Empyrean diffractometer with PIXcel detector and fixed slits with Fe filtered Co -K? radiation. The phases were identified using X’Pert Highscore plus software; The relative phase amounts (weight%) were estimated using the Rietveld method. Errors are on the 3 sigma level in the column to the right of the amount (in weight per cent). RESULTS The results of the XRD and XRF are listed below in Table A -1, Table A -2 and Table A -3 respectively. Table A-1: XRD results (weight % with error) Composition (%) [s] Sample 1 Mineral Amount Error (weight %) Anatase 5.45 0.45 Calcite 6.22 0.72 Dolomite 15.17 1.17 Kaolinite 61.6 1.35 Pyrite 5.43 0.6 Quartz 6.13 0.6 Page 67 of 77 Table A-2: XRF major elements (weight %) Trace Elements Trace Element Concentration Sample 1 U 2.36 V 54 W 1.39 Y 22.5 Yb 3.7 Zn 37 Zr 118 As <5.00 Ba 88 Bi <5.00 Br <1.00 Cd <5.00 Ce <5.00 Cl 327 Co 17.1 Cs <1.00 Cu 36.1 Ga 14.1 Ge <1.00 Hf 2.81 Hg <5.00 La 17.9 Lu <1.00 Mo 3.32 Nb 9 Nd 21.1 Ni 37 Pb <5.00 Rb 18.5 Sb <5.00 Sc 6.2 Se <1.00 Sm 3.15 Sn <5.00 Sr 156 Ta 1.23 Te 2.96 Th 11.1 Tl <1.00 Page 68 of 77 Table A-3: XRF trace element concentration (ppm) Trace Elements Trace Element Concentration Sample 1 U 2.36 V 54 W 1.39 Y 22.5 Yb 3.7 Zn 37 Zr 118 Page 69 of 77 Appendix B: Static Procedures INTRODUCTION Acid base accounting (ABA) is a valuable and widely used static test that determines both the acid-generating and acid neutralizing potential of a sample. A description of the different ABA components is given below: Acid potential (AP) is determined by the % S with a factor of 31.2 5. The units of AP is kg CaCO3/ t rock and indicates the theoretical amount of calcite that would be required to neutralize the acid produced; and The Neutralization Potential (NP) is determined by reacting the sample with a known excess of standardized hydrochloric or sulphuric acid. The paste is then back titrated with standardized sodium hydroxide. This is done to determine the unconsumed acid. NP is also expressed as kg CaCO3/t rock. Net Neutralization Potential (NNP) is determined by subtracting AP from NP. NNP is used to screen the ABA results in terms of their ARD potential. The criteria are as follows : If NNP = NP – AP < 0 The sample has the potential to generate acid If NN = NP – AP > 0 The sample has the potential to neutralise acid Research have shown that a range of values for NNP between -20 to 20 kg/t CO 3 can either become acidic or neutral, and is therefore uncertain (Usher et al., 20 02). Neutralising Potential Ratio (NPR) is the NP: AP ratio. Guidelines for screening criteria based on the ABA results are listed below in Table B -1. Table B-1: Screening methods using NPR (Price, 1997) Potential to generate Acid Rock Drainage NPR Comments Rock Type I: Likely <1:1 Likely to generate ARD Rock Type II: Possibly 1:1-2:1 Possible ARD generating Rock Type III: Low 2:1-4:1 Possibly AMD generating if NP is insufficiently reactive or is depleted at a faster rate than sulphides. Rock Type IV: None >4:1 No further AMD testing required unless materials are to be used as a source of alkalinity. Page 70 of 77 % S and NPR From T able B -1 and %S from the ABA (Leco) results a graph can be plotted to visually represent the ABA results. Samples with %S content of less than 0.3% is regarded as having insufficient oxidisable %S and will not sustain acid generation. Samples with a NRP value above 4 is considered to have enough neutralisation capacity to be classified to have a very low probability to produce acid. Samples with a NRP value of lower than 1 and a %S content that exceed 0.3% are regarded as highly probable of producing acid. NAG test terminology and screening methods The Net acid Generating (NAG) test consists of rapidly oxidizing sulphide minerals in a sample using hydrogen peroxide (H 2O 2). This value can serve as a rough guideline to assess the potential of a material to produce acid after a period of weathering. It should be noted that this represents a worst case scenario where all the sulphide minerals are reacted completely. This will not be the case under field conditions, as sulphide minerals are shielded in the rock matrix. Several other factor also control the rate of sulphide mineral oxidation, causing a slower release of sulphide oxidation products than what is observed in peroxide extraction. It is assumed with these tests that the elevation of SO 4 is a direct result of the sulphide mineral oxidation. The methodology is described by Usher et al. (2002). The potential of the sample to generate acidic drainage is calculated as the difference Table B-2: Screening criteria for NAG test Final pH in NAG Test Acid Generating Potential >5.5 Non-acid generating 3.5 -5.5 Low risk acid generating <3.5 High risk acid generating Page 71 of 77 Appendix C: Leach Test Data Page 72 of 77 Chemical Parameter (mg/l) WEEK 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 pH 7.35 7.6 9 7.7 4 8.24 8.11 8.4 9 8.51 8.31 8.04 8.07 8.06 7.8 7 7.8 8 8.04 8.05 8.00 8.00 8.01 7.9 6 8.02 Electrical Conductivity (mS/m) 97. 9 54. 3 42.0 30.6 31.0 33. 6 32.0 26.4 27.1 22.9 21.1 25.3 24.0 17.1 14.4 12.5 12.9 11.0 11.9 8.9 Total Dissolved Solids 669 365 291 210 214 229 218 182 186 159 146 172 166 120 100 88 92 74 80 62 Total Alkalinity as CaCO 3 120 140 148 120 108 136 136 108 104 100 108 136 138 83 65 60 60 53 55 46 Ammonia -N 0.3 <0.2 0.2 <0.2 <0.2 <0.2 <0.2 <0.2 0.2 <0.2 <0.2 <0.2 <0.2 < 0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 Nitrate-N <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 Chloride-Cl <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 <0.2 < 0.2 Sulphate-SO 4 443 143 81 41 45 66 63 31 45 30 24 33 33 11 8 7 8 6 8 <5 Fluoride -F <0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 <0.2 <0.2 <0.2 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 ICP -MS Scan (mg/l) WEEK 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Silver-Ag 0.001 0.000 0.000 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 Aluminium -Al 0.134 0.108 0.099 0.091 0.091 0.086 0.088 0.082 0.102 0.106 0.089 0.104 0.106 0.062 0.063 0.056 0.054 0.050 0.055 0.047 Arsenic -As 0.000 0.000 0.000 0.001 0.002 0.001 0.000 0.000 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 Boron -B 0.383 0.350 0.297 0.228 0.199 0.178 0.169 0.172 0.208 0.169 0.105 0.092 0.087 0.051 0.047 0.034 0.034 0.029 0.027 0.019 Barium -Ba 0.250 0.202 0.140 0.157 0.155 0.141 0.136 0.131 0.131 0.124 0.116 0.061 0.062 0.133 0.210 0.173 0.222 0.220 0.252 0.231 Beryllium -Be 0.000 0.001 0.000 0.000 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.000 0.000 0.000 0.001 0.000 0.001 0.001 0.000 Bismuth -Bi 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.002 0.003 0.003 0.000 0.002 0.000 0.000 0.001 0.003 Calcium-Ca 126.6 68.03 52.60 38.51 39. 60 40.30 37. 40 34. 80 51.04 50.56 37.59 59.99 61.12 23.71 22.37 17.13 17.38 15.89 17.28 13.14 Cadmium-Cd 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.000 0.001 0.001 0.000 Cobalt-Co 0.004 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.000 0.000 0.000 0.000 Chromium-Cr 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 Copper-Cu 0.003 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 Iron -Fe 0.076 0.021 0.021 0.010 0.061 0.044 0.036 0.029 0.012 0.009 0.006 0.005 0.024 0.010 0.002 0.002 0.004 0.003 0.007 0.003 Potassium-K 9.524 6.8 39 5.107 3.9 60 3.4 70 3.286 3.140 3.076 3.9 80 3.8 68 2.374 3.314 3.3 37 1.198 0.852 0.682 0.653 0.563 0.520 0.429 Lithium -Li 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.018 0.018 0.000 0.017 0.017 0.018 0.019 0.019 0.000 Magnesium-Mg 39. 85 20.71 15.97 10.53 10.24 9.8 60 8.3 40 8.070 14.67 13. 88 10.42 20.73 20.74 6.7 49 6.013 4.574 4.626 3.952 4.299 3.084 Manganese-Mn 0.255 0.104 0.077 0.061 0.066 0.071 0.063 0.059 0.082 0.081 0.058 0.050 0.049 0.023 0.020 0.012 0.010 0.008 0.008 0.005 Page 73 of 77 ICP -MS Scan (mg/l) WEEK 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Molybdenum-Mo 0.003 0.004 0.004 0.004 0.004 0.002 0.003 0.004 0.003 0.002 0.002 0.001 0.001 0.001 0.000 0.000 0.001 0.000 0.001 0.001 Sodium-Na 7.512 5.085 4.101 2.698 2.314 2.106 2.077 1.894 3.3 30 3.082 1.791 2.955 3.068 0.911 0.726 0.547 0.635 0.495 0.633 0.379 Nickel -Ni 0.046 0.071 0.027 0.013 0.015 0.011 0.009 0.012 0.025 0.037 0.025 0.034 0.036 0.006 0.007 0.004 0.003 0.002 0.001 0.001 Phosphorous-P 0.025 0.012 0.003 0.007 0.018 0.012 0.017 0.011 0.000 0.003 0.005 0.007 0.004 0.000 0.009 0.000 0.004 0.000 0.000 0.006 Lead -Pb 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.000 0.001 0.006 0.003 0.005 0.003 0.002 0.004 0.002 0.002 0.002 Sulphur-S 162.9 54. 37 28. 44 14. 49 16. 61 19.48 18. 82 14.60 53.12 51.65 31.29 40.11 40.14 5.761 6.310 4.404 4.6 61 3.9 79 4.7 33 2.952 Antimony -Sb 0.004 0.010 0.008 0.002 0.006 0.004 0.002 0.002 0.000 0.000 0.000 0.001 0.002 0.002 0.002 0.000 0.000 0.000 0.000 0.000 Selenium-Se 0.002 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.000 0.000 0.000 0.001 0.001 0.002 0.002 0.001 0.002 0.002 0.001 Silicon-Si 1.293 1.357 1.254 0.960 0.949 0.881 0.869 0.749 1.084 1.028 0.643 0.810 0.790 0.290 0.295 0.209 0.209 0.180 0.199 0.144 Tin -Sn 0.000 0.000 0.000 0.006 0.000 0.001 0.000 0.001 0.000 0.000 0.000 0.010 0.012 0.013 0.003 0.010 0.009 0.017 0.019 0.015 Stronsium-Sr 0.724 0.392 0.377 0.313 0.320 0.334 0.306 0.298 0.367 0.368 0.278 0.329 0.322 0.156 0.162 0.128 0.137 0.124 0.132 0.098 Titanium -Ti 0.032 0.012 0.002 0.004 0.019 0.011 0.077 0.010 0.017 0.028 0.016 0.004 0.000 0.000 0.007 0.015 0.013 0.009 0.005 0.002 Thallium -Tl 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.002 0.000 Vanadium -V 0.000 0.001 0.000 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 Wolfram-W 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.003 0.002 0.000 0.000 0.003 0.002 0.001 0.003 0.002 Zinc -Zn 0.064 0.054 0.037 0.022 0.054 0.048 0.042 0.040 0.039 0.019 0.088 0.033 0.033 0.017 0.024 0.013 0.012 0.011 0.010 0.008 Zirconium -Zr 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 Page 74 of 77 Appendix D: PHREEQC Code Page 75 of 77 CSTR function: TRANSPORT -cells 1; -boundary_conditions constant closed; -flow_direction diffusion_only -stagnant 1 # 1 stagnant layer, but more are possible, for modeling bad mixing in the tank -time_step 144 hour # each time_step, the MIX is performed -shifts 20 # number of time_steps, total time is 100 hours -punch_cells 3 # only graph the tank solution -punch_frequency 1 # sample every 5 hours Page 76 of 77 Rates: RATES Pyrite -start 1 A = 2* m0 # initial moles of the kinetic reactant ( initial surface per mole of pyrite (m2/mol pyrite)) 2 V = 0.75 10 if SI("Pyrite")> 0 then goto 100 20 fH = mol("H+") 30 fFe2 = (1 + tot("Fe(2)")/1e-6) # 40 if mol("O2") < 1e-6 then goto 80 50 rO2 = 10^-8.3 * mol("O2")^0.5 * fH^-0.11 90 rate = A/V * (m/m0)^ 5* rO2 * (1-SR("Pyrite")) # The shrinking reactive surface seems to control the reaction rate, as the rate is much faster when published values for n is used. A very large exponential factor for (m/m0) is requiered to fit the data 100 save rate * time -end # Gypsum # -start # 1 A = 4 *m0 # 2 V = 0.75 # 10 rate = (A/V)*(m/m0)^6* 5.8*10^-7.6*(1-SR("Gypsum")) # 20 save rate*time # -end Dolomite Page 77 of 77 -start 1 A = 20*m0 2 V = 0.75 10 rate = (10^-10 )* (1 - SR("Dolomite"))*((m/m0)^30)*((A/V)) 20 save rate * time -end Kinetics: KINETICS 3 Pyrite; -m0 0.045 # Gypsum; -m0 0.018 Dolomite;-m0 0.82