UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Development and testing of a 2D axisymmetric water flow and solute transport model for heap leaching Afewu, Kodjo Isaac 2009

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2009_fall_afewu_kodjo.pdf [ 9.93MB ]
Metadata
JSON: 24-1.0067786.json
JSON-LD: 24-1.0067786-ld.json
RDF/XML (Pretty): 24-1.0067786-rdf.xml
RDF/JSON: 24-1.0067786-rdf.json
Turtle: 24-1.0067786-turtle.txt
N-Triples: 24-1.0067786-rdf-ntriples.txt
Original Record: 24-1.0067786-source.json
Full Text
24-1.0067786-fulltext.txt
Citation
24-1.0067786.ris

Full Text

DEVELOPMENT AND TESTING OF A 2D AXISYMMETRIC WATER FLOW AND SOLUTE TRANSPORT MODEL FOR HEAP LEACHING  by KODJO ISAAC AFEWU B.Sc. (Engineering), University of Science & Technology, Ghana, 1985 M.Sc. (Metallurgy), University of Zambia, 1990 M.Sc. (Engineering), University of the Witwatersrand, 1993  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY  in  THE FACULTY OF GRADUATE STUDIES (Materials Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  October 2009 © Kodjo Isaac Afewu, 2009  Abstract Heap leaching is an industrially relevant process for extracting metal values from marginal ores and/or pretreating ores for subsequent extraction of valuable metals.  Its advantages include low capital and operating costs, rapid  construction times, and relatively low environmental impact. The movement of water and solutes through a heap of ore has serious implications for efficiency of heap leach operations. Solutes (whether introduced into or generated within a heap) move with the water and tend to mix with the bulk solution by the mechanisms of advection and dispersion. Though it is well known that ore-water characteristics influence water flow and solute transport properties during heap leaching, these properties have not been adequately incorporated into the design and operation strategies of heap leaching. In this work, a mathematical model and the results of special laboratory experiments are used to study the reaction and/or transport of solutes during heap leaching. The model is 2D axisymmetric, comprising transport of solutes by advection and dispersion in a bed of broken ore under a point source of solution representing a single drip emitter.  The ore hydraulic parameters characterize the water  retention and permeability properties, while the transport and chemical parameters affect the distribution of water and solutes, and the rates and extents of reactions of solutes in the heap. The close agreement between the observed data and the predicted results lends credence to the modeling approach and its solution methodology.  It was  observed that, within the limited water content range relevant to safe and successful operation of a heap leach facility, dispersivity increased with increasing irrigation flux and hence water content.  The model captured the  trends in water and lixiviant distribution, and solute transport in heaps, clearly demonstrating the importance of incorporating ore hydraulic parameters in modeling heap leaching. It was also found that transport of solutes is slower than leaching and hence the former is the rate limiting step in the process. The  ii  model can be used to simulate various process enhancement strategies for efficient operation and optimal design, as well as economic evaluations and potential environmental impact assessment of heap leaching operations.  iii  Table of Contents Abstract................................................................................................................. ii Table of Contents................................................................................................. iv List of Tables ...................................................................................................... vii List of Figures .....................................................................................................viii Acknowledgements............................................................................................. xii Chapter 1 – Introduction .......................................................................................1 1.1 Heap Leaching – Overview.....................................................................1 1.2 Problem Definition ..................................................................................4 1.2.1 Modelling Water Flow and Solute Transport in Heap Leaching.......9 1.3 Objectives of this Investigation .............................................................12 1.4 Thesis Outline.......................................................................................14 Chapter 2 – Flow and Transport in Unsaturated Beds ........................................15 2.1 Heap Leaching: Current Developments and Applications.....................15 2.1.1 Heap Leaching of Gold Ores .........................................................16 2.1.2 Heap Leaching of Copper Ores.....................................................17 2.1.3 Heap Leaching of Uranium Ores ...................................................19 2.2 Flow and Transport in Unsaturated Porous Media ...............................20 2.2.1 Water Flow in Unsaturated Media .................................................21 2.2.2 Solute Transport ............................................................................31 2.3 Review of Current Methods of Modeling Heap Hydrology ....................40 2.3.1 Solution Flow in Heaps..................................................................41 2.3.2 Solute Transport in Heaps.............................................................42 2.4 Conclusions and Focus of this Study....................................................43 Chapter 3 – Model Development ........................................................................45 3.1 The Water Flow Model..........................................................................45 3.1.1 Constitutive Relations....................................................................47 3.1.2 Discretization of the Water Transport Model .................................51 3.1.3 Solution of the Water Transport Model ..........................................54 3.1.4 Testing the Water Transport Model ...............................................55 3.2 The Solute Transport Model .................................................................56 3.2.1 Formulation of the Solute Transport Model ...................................57 3.2.2 Discretization of the Solute Transport Model.................................59 3.2.3 Solution of the Solute Transport Model .........................................62 3.2.4 Testing the Solute Transport Model...............................................66 Chapter 4 – Experimental Methodologies ...........................................................68  iv  4.1 Introduction...........................................................................................68 4.2 Materials, Equipment and Procedures..................................................69 4.2.1 Ore ................................................................................................69 4.2.2 Reactor Geometries ......................................................................71 4.2.3 Tracer ............................................................................................71 4.2.4 Measuring Devices and Data Acquisition Systems .......................72 4.3 Test Procedures ...................................................................................76 4.3.1 Tracer Tests ..................................................................................77 4.3.2 Small Column Apparatus...............................................................77 4.3.3 The Large Column.........................................................................81 4.4 Performance Tests ...............................................................................83 4.4.1 Mini-Columns ................................................................................84 4.4.2 Acid Leaching in the Large Column...............................................84 Chapter 5 – Results and Discussion – Tracer Tests ...........................................86 5.1 The Small Column ................................................................................86 5.1.1 Determination of Hydraulic Parameters.........................................86 5.1.2 Estimation of Transport Parameters ..............................................96 5.2 Estimation of Parameters for the 3D Axisymmetric Transport Model .105 Chapter 6 – Results and Discussion – Performance Tests...............................117 6.1 Introduction.........................................................................................117 6.1.1 The Concept of Leaching Large Particles....................................117 6.2 Reactive Solute Transport Model .......................................................119 6.2.1 Dissolution and Transport of Copper ...........................................119 6.2.2 Kinetics of Leaching ....................................................................121 6.2.3 Numerical Solution ......................................................................123 6.2.4 Calibration of the Model...............................................................123 6.2.5 Estimation of Hydraulic and Dispersive Parameters....................124 6.2.6 Estimation of Chemical Parameters ............................................124 6.3 Results and Discussion of the Reactive Solute Transport Model .......131 Chapter 7 – Sensitivity of Copper Extraction to Heap and Solute Parameters .141 7.1 Introduction.........................................................................................141 7.2 Hydraulic Properties ...........................................................................144 7.2.1 Emitter Spacing ...........................................................................144 7.2.2 Heap Height ................................................................................147 7.2.3 Irrigation Flux...............................................................................149 7.2.4 Hydraulic Conductivity .................................................................152 7.2.5 Air Entry Capillary Head ..............................................................155 7.3 Transport Properties...........................................................................157 7.3.1 Dispersion Coefficients................................................................157 7.3.2 Diffusion Coefficient.....................................................................162 7.4 Chemical Reaction Parameters ..........................................................163 7.4.1 Particle Reaction Rate Constants................................................163 7.4.2 Acid Concentration ......................................................................168 v  7.5  Practical Applications of Sensitivity Analysis for Heap Leaching ........171  Chapter 8 – Conclusions and Recommendations.............................................173 8.1 8.2  Summary and Conclusions.................................................................173 Recommendations..............................................................................176  References .......................................................................................................178 Appendix A – Deriving the Equations of Continuity...........................................188 A.1 A.2  Deriving the Equation of Continuity for Water Flow ............................189 Deriving the Equation of Continuity for Solute Transport....................190  Appendix B – Time Domain Reflectometry .......................................................192 B.1 Overview.............................................................................................192 B.2 A Brief History of the Evolution of the TDR Technique .......................193 B.3 TDR Theory ........................................................................................194 B.3.1 Volumetric Water Content ...........................................................194 B.3.2 Electrical Conductivity .................................................................196 B.4 Orientation of Probe Installation .........................................................200 B.4.1 Vertical Installation ......................................................................201 B.4.2 Horizontal Installation ..................................................................201 B.5 Limitations of TDR ..............................................................................202  vi  List of Tables Table 4.1  Chemical composition of the ore (major metallic elements only)...70  Table 5.1  Moisture retention with irrigation flux (drying cycle) .......................87  Table 5.2  Moisture retention with irrigation flux (wetting cycle) .....................87  Table 5.3  Hydraulic conductivity parameters derived from the Small Column ......................................................................................................94  Table 5.4  Capillary pressure parameters derived from the Small Column ....94  Table 5.5  Variation of dispersivity with applied flux in the Small Column ......99  Table 5.6  Grid points queried for comparison of model and TDR results ....106  Table 5.7  Summary of concentration profiles of TDR probes......................109  Table 5.8  Optimal αT/αL ratios for various solution fluxes of the tracer data 109  Table 6.1  Kinetic copper oxide leaching parameters from the ore by bottle roll-tests.......................................................................................125  Table 6.2  Stoichiometric sulphuric acid required to leach copper oxides ....130  Table 6.3  Optimal values of kinetic parameters for copper and gangue leaching from the10-cm mini-columns and the Large Column.....132  Table 7.1  Parameter values used as the baseline case for simulations in the Large Column..............................................................................143  Table 7.2  Ratio of quantity of copper extracted from heaps of varying emitter spacing with time .........................................................................143  Table B.1  Dielectric constants of soil materials (20–25ºC) [Or et al., 2003].195  vii  List of Figures Figure 1.1  Schematic flow diagram of a typical copper heap leach and solvent extraction circuit ..............................................................................3  Figure 2.1  Typical relationship between permeability and degree of saturation for finely textured and coarsely textured materials [after Lal and Shukla, 2004]. ...............................................................................26  Figure 2.2  Preferential flow induced by irrigation rate and particle size [after O’Kane et al., 1999].......................................................................28  Figure 3.1  A schematic of a heap under irrigation from drip emitters .............46  Figure 3.2  Schematic diagram of a discretised 2D axisymmetric grid ............51  Figure 3.3  Predicted column wetting patterns at extreme values of capillary head in a 4 × 1 m column: a) hc,0 = 1.0 m; b) hc,0 = 0.001 m (all other parameters as determined in Chapter 5). .............................56  Figure 3.4  Predicted relative solute concentrations for a first-order chemical reaction (rate constant k = 1/day) at extreme values of dispersion parameters in a 4 × 1 m column: a) steady-state saturation profile used in the simulations (all hydrological parameters as determined in Chapter 5); b) Baseline dispersion parameters (as determined in Chapter 5); c) Dispersion parameters set to zero; d) Dispersion parameters set to 100 times the baseline......................................67  Figure 4.1  Particle size distribution curve of the as-received ore....................70  Figure 4.2  The four-electrode sensor circuit ...................................................75  Figure 4.3  A photograph of a mini-tensiometer with transducer .....................76  Figure 4.4  A photograph of the Small Column................................................78  Figure 4.5  Schematic diagram of the Small Column.......................................79  Figure 4.6  A photograph of the Large Column, with TDR probe layout ..........82  Figure 4.7  Schematic diagram of the Large Column and instrumentation, with TDR probe layout ..........................................................................83  Figure 5.1  Final drain-down curve for the Small Column................................88  Figure 5.2  Moisture retention curves with best power-law fits (drying) ...........89  Figure 5.3  Moisture retention curves with best power-law fits (wetting)..........90  viii  Figure 5.4  Capillary head vs effective saturation with best power-law fits ......91  Figure 5.5  Capillary pressure power-law fits based on average hydraulic parameters from the Small Column: capillary pressure .................95  Figure 5.6  Capillary pressure power-law fits based on average hydraulic parameters from the Small Column: derivative of capillary pressure ......................................................................................................96  Figure 5.7  Small Column normalized effluent tracer data and model fits........98  Figure 5.8  Linear correlation of volumetric water content with normalizationcorrected volume flux vw ..............................................................100  Figure 5.9  Linear correlation of longitudinal dispersivity with normalizationcorrected volume flux vw ..............................................................101  Figure 5.10 Diffusivity function vs volumetric water content from several studies [adapted from Lim et al., 1998]....................................................103 Figure 5.11 Small Column tracer curves at various depths and their corresponding model fits (at hc,0 = 0.05 m, vw = 16.2 L/m2/h).......104 Figure 5.12 Simulated (at 16 radial nodes) and observed tracer profiles in the Large Column at an irrigation flux of 6 L/m2/h..............................106 Figure 5.13 Se profile in the Large Column ....................................................110 Figure 5.14 Model and experimental effluent tracer responses from the Large Column at fluxes of 6, 9 and 12 L/m2/h assuming hc,0 = 0.082 m 112 Figure 5.15 Model and experimental effluent tracer responses from the Large Column at fluxes of 6, 9 and 12 L/m2/h assuming h0 = 0.05 m ....113 Figure 5.16 Axial velocity profiles for the single point irrigation in the Large Column tests predicted using the parameters determined from the Small Column tests......................................................................114 Figure 6.1  Schematic of leaching copper from a round particle of copper ore ....................................................................................................118  Figure 6.2  Copper extraction from bottle-roll test at 5 g/L initial acid ............126  Figure 6.3  Copper extraction from bottle-roll test at 10 g/L initial acid ..........126  Figure 6.4  Copper extraction from bottle-roll test at 20 g/L initial acid ..........127  Figure 6.5  Calculated rates of acid consumption during bottle roll leach tests ....................................................................................................128 ix  Figure 6.6  Rate of acid consumption by copper oxides ................................129  Figure 6.7  Effluent acid concentrations from the 10-cm mini-columns .........134  Figure 6.8  Effluent copper concentrations from the 10-cm mini-columns .....135  Figure 6.9  Effluent copper and acid concentrations from the Large Column 135  Figure 6.10 Predicted water content profile in the Large Column ...................136 Figure 6.11 Predicted aqueous Cu concentration profile in the Large Column ....................................................................................................137 Figure 6.12 Predicted acid concentration profile in the Large Column ...........139 Figure 7.1  Simulation of varying emitter spacing on copper extraction.........145  Figure 7.2  Acid profiles for copper extraction with varying emitter spacing ..145  Figure 7.3  Simulation of varying heap height on copper extraction ..............147  Figure 7.4  Acid profiles for copper extraction with varying heap height........148  Figure 7.5  Simulation of varying irrigation flux on copper extraction.............150  Figure 7.6  Acid profiles for copper extraction with varying irrigation flux ......150  Figure 7.7  Simulation of varying heap permeability on copper extraction.....153  Figure 7.8  Acid profiles for copper extraction with varying hydraulic conductivity..................................................................................153  Figure 7.9  Simulation of varying air entry head on copper extraction ...........155  Figure 7.10 Acid profiles for copper extraction with varying air-entry head ....156 Figure 7.11 Simulation of varying the slope of αL vs vw on copper extraction .158 Figure 7.12 Simulation of varying the transverse-to-longitudinal dispersivity ratio αL/αT on copper extraction ...........................................................158 Figure 7.13 Acid profiles for copper extraction with varying longitudinal dispersivity...................................................................................159 Figure 7.14 Acid profiles for copper extraction with varying transverse-tolongitudinal dispersivity ratio........................................................160 Figure 7.15 Simulation of varying molecular diffusivity coefficient on copper extraction.....................................................................................162  x  Figure 7.16 Simulation of varying gangue dissolution rate constant on copper extraction.....................................................................................163 Figure 7.17 Acid profiles for copper extraction with varying gangue acid consumption rate constant ..........................................................164 Figure 7.18 Simulation of varying the copper leaching rate constant on copper extraction.....................................................................................166 Figure 7.19 Acid profiles for copper extraction with varying copper leaching rate constant.......................................................................................167 Figure 7.20 Simulation of effects of varying acid concentration on copper extraction.....................................................................................169 Figure 7.21 Acid profiles for copper extraction with varying acid concentration ....................................................................................................170 Figure A.1  Schematic of a cylindrical differential element.............................188  Figure B.1  A schematic diagram of TDR wave showing points of measurement for water content and load resistance [Ward et al., 1994]............197  xi  Acknowledgements I wish to express my heartfelt gratitude to the following people who have played distinguished roles during the execution of this research: Professor David Dixon for his guidance, financial support and supervision throughout the course of this research. His knowledge and insight of heap leaching, as well as his keen and tireless devotion to improve and optimize heap leaching technology have been a valuable source of inspiration for this work. Members of my supervising team, namely Professors Ward Wilson and Leslie Smith. Amira P768 for providing the research funds. University of British Columbia for providing me with Phd Program Scholarship. Mantoverde Copper Mine for providing the copper oxide samples. Professor Michael Novak for lending me the TDR Cable Tester Tektronix 1502B and Dr Craig Nichol for helping me with the TDR setup. Professor David Dreisinger, Dr Winfred Assibey-Bonsu and Dr Stephen Lam for their encouragement. My colleagues, Berny Rivera, Darren Mayne, Gonzalo Viramonte, Jozsef Miskolcczi, Louis Kabwe for the help with the experimental setup. Workshop and store staff, namely, Ross McLeod, David Torok, Carl Ng, Serge Milaire and Glenn Smith for the equipment fabrication and ordering. Mineral Processing Laboratory staff, Sally Finora, Pius Lo and Keith Jellis, for providing space for the experimental setup.  xii  My brother Kofi and sister Sarah for their encouragement. My mother Cesia Afewu, who first taught me to count numbers. She could not live to see the completion of this work; may she rest in peace. My dear wife Gladys and daughters Kekeli, Makafui and Elorm, for their encouragement, support and patience. This work is dedicated to them.  xiii  Chapter 1 INTRODUCTION 1.1  Heap Leaching – Overview  Worldwide, heap leaching has become an acceptable method of extracting minerals from marginal ores and/or pretreating ores for subsequent extraction of valuable minerals, owing to its low capital and operating costs, short construction times and hence early return on capital, and environmental friendliness [Bartlett, 1997].  Many challenges confront the mineral extraction industry, including:  declining ore grades, increasing ore complexity, rising power and labour costs associated with conventional milling and slurry agitation or smelting and refining, ever-increasing environmental restrictions, and remaining competitive in shrinking markets due to materials substitution. These have prompted research and development of alternative mineral extraction processes that yield high extractions in a reasonable time frame at low capital and operating costs [Harrington et al., 1993; Montero et al., 1994; Bartlett, 1997]. This investigation is on heap hydrology, dealing with flow of solution and transport of solutes during heap leaching. As an introduction to the study, an overview of the heap leach process will be given. This will set the background required to explain the problem, the objectives and the relevance of this investigation. Heap leaching is probably one of the oldest methods of extracting metals from ores. It was practised for extraction of copper from cupriferous pyrite as early as 1752 at Rio Tinto in Spain [Gupta and Mukherjee, 1990; Bartlett, 1997], though the actual scientific inquiry into the manifold and intricate processes involved in heap leaching only started around the middle of the last century [Bartlett, 1997]. A modern heap is constructed on engineered and impervious ground as the base. The base is built with a slight slope [Gupta and Mukherjee, 1990;  Introduction  1  www.straits.com.au/technology] towards a series of solution drains. The surface of the ground is covered with an impervious pad (made from geosynthetic materials such as High Density Polyethylene (HDPE)) to further prevent solution leakage, which constitutes both a loss to production and a potential source of ground water pollution. In some cases, the ore is sent directly from the mine (run-of-mine ore) with little or no preparation. However, in many heap leach operations, the ore is crushed and agglomerated, prior to loading on the pad, to enhance heap permeability. Commercial heaps are huge reactors, extending tens or hundreds of meters in width and breadth and typically 6 to 15 meters in height [Dixon and Petersen, 2003]. The solution-piping network generally comprises main lines that run at ground level and header pipes along the heap surface, which distribute the solution over the surface of the heap through sprinklers or drip emitters [Orr and Vesselinov, 2002].  The use of drip emitters reduces agglomerate degradation and  evaporation losses of solution due to solution spraying. Emitters are therefore generally  preferred  to  sprinklers  www.straits.com.au/technology].  [Orr  and  Vesselinov,  2002,  However, since drip emitters were originally  designed to localize flow rather than to spread it uniformly, a dense grid is required to maximize leach uniformity and efficiency. They are also prone to clogging by solids and, in the case of copper leaching, micro-organisms which feed on kerosene residues in the leach solutions (from recycling Solvent Extraction (SX) raffinate for leaching). Adding a bio-filtration system to the leach solution circuit may minimize the effects of clogging [Orr and Vesselinov, 2002]. These additional modifications, however, increase the capital and operating costs of drip emitters. The drip emitters are commonly on a grid of 50 to 100 cm spacing [Petersen and Dixon, 2003]. The solution dissolves the target metal or mineral as it percolates through the pores and interstices of the ore particles in the heap onto the impervious pad and then drains by virtue of the sloping surface into the collection basins. The solution in the drains may be recycled (to the surface of the heap) or sent for further processing for metal extraction.  Introduction  2  Commercial operations of heap bioleaching, which is a variant of heap leaching for sulfide ores, began in 1980 [Brierley, 1980]. In heap bioleaching, certain microbes are used to catalyze oxidation reactions in the heap. The flow diagram of a typical copper heap leach with solvent extraction is given in Figure 1.1. Drip Emitters  Heap  Sloping Impervious surface Pregnant Solution  Raffinate Pregnant Solution Pond  Raffinate/Barren Pond Acid/(Bugs)  Loading  Loaded organic  Solvent Extraction  Stripped organic  Stripping Spent Electrolyte from Electrowinning  Figure 1.1  Advance Electrolyte to Electrowinning  Schematic flow diagram of a typical copper heap leach and solvent extraction circuit  Once the product metal concentration in the heap effluent has fallen below an economic cutoff point, production stops and the heap is decommissioned. Prior to heap decommissioning, the heap is rinsed to remove any excess reagents (such as cyanide or acid) and residual dissolved products (such as gold or base metals) to abate the potential for pollution of the environment, especially the groundwater. In the case of heap bioleaching of refractory gold ores, the spent heap is the valuable portion, since the liberated gold remains within the heap. In this case, the ore may be neutralized in place by lime for further heap leaching  Introduction  3  with cyanide or thiosulfate solutions, or reclaimed for grinding and feeding to a conventional Carbon-In-Leach (CIL) circuit [Bhakta and Arthur, 2002].  1.2  Problem Definition  Today, commercial heaps exist for the extraction of gold, copper and uranium from oxide and sulfide ores with varying degrees of modification and success [Gupta, 1990; Bartlett, 1997; Dixon, 2003; Petersen and Dixon, 2003; Marchbank et al., 2003]. With increasing use of this technology, it has become clear that a thorough insight into the effects of the driving forces, including chemical, physical, and microbial factors and their interactions, is the key to the successful implementation and improvement of the technology. To improve the economics of the process, concerted efforts are being put toward the development of this understanding by industry and academia. For example, the study of the thermodynamics, kinetics and mechanisms of the reactions in heap leaching (with and without microbial mediation) [Braun et al., 1974; Madsen et al., 1975; Madsen and Wadsworth, 1981] has provided much insight into the effects of limiting reactions and interpretation of heap leach reaction data. These have led to various strategies for optimizing the reactions. However, heap leaching technology still experiences some hydrological problems, resulting in relatively low rates and extents of extraction and necessitating long retention times to achieve a high degree of mineral extraction. For example, many copper heaps take more than a year to achieve 80% copper extraction [Petersen and Dixon, 2003]. The main hydrological problem in heap leaching is the inadequate and non-uniform distribution of lixiviant (water and solutes) to targeted zones in the heap under drip emitters [Schlitt and Nicolai, 1987; Petersen and Dixon, 2003; Dixon and Petersen, 2003].  Though the  importance of the effects of heap hydrology (which comprises flow and storage of solution, and transport of solutes and oxygen) on heap leach performance has been identified and documented since 1968 [Howard, 1968] and in later works [Roman et al., 1974, Woolery et al., 1978; Schlitt and Nicolai, 1987; Roman and  Introduction  4  Bhappu, 1993; Bartlett, 1997; Sanchez-Chacon and Lapidus, 1997; Orr and Vesselinov, 2002; Orr, 2002], little has been done to shed light on the underlying hydrological factors that contribute to heap leach inefficiency. Sanchez-Chacon and Lapidus [1997] undertook simulation of gold heap leaching and found that it only required about six days to dissolve the gold, but, due to the slow diffusion of the soluble gold complex, it required about seven weeks for the soluble gold to decrease to acceptable levels in the pores.  Thus, internal mass transport  seemed to be the most important resistance in the process [Sanchez-Chacon and Lapidus, 1997].  Compared to agitated systems, low rates of mineral  extraction from heap leaching are expected, because the particles of a heap are much coarser than those involved in agitated systems. This is exacerbated by the fact that a heap is an unconsolidated fixed bed that is being irrigated from point sources to maintain unsaturated conditions. The hydrological problems encountered during heap leaching are caused by three main factors: poor heap stacking, inappropriate irrigation schedule or method, and precipitation reactions within the heap. The effects of any one of these factors or a combination of them lead to one of the most important problems confronting heap operators: preferential flow. The details of the three factors and their contributions to the hydrology problem will be discussed. In the design of a heap, permeability to solution and gas is an important factor. It is common knowledge that fine materials generally leach faster than coarse materials. However, in heaps, fine materials tend to pack tightly, resulting in low heap permeability to the flow of solution and gas. The effective permeability will decrease more rapidly for uniform materials than for well-graded materials [Hillel, 1998; O’Kane et al., 1999]. During the stacking of a heap, segregation of coarse and fine materials, inclusion of a high proportion of clayey materials, and compaction due to the weight of loading trucks may result in low heap permeability, and in particular, areas of significantly different local permeability. Wide variations of permeability in the heap enhance preferential flow of solution, as the solution tends to flow through zones of higher permeability.  Introduction  5  Agglomeration of the crushed ore prior to stacking in the heap is the main method of achieving uniform heap permeability and rates of reactions when fine materials are in the ore. The primary purposes of irrigating the heap are to transport reacting solutes for leaching to the various parts of the heap, to transport the products of reaction out of the heap, and to maintain heap temperatures within the limits of activity of microbes in the case of heap bioleaching. Solution dripped onto the surface gets distributed throughout the depth and breadth of the heap by the combined effects of gravity and capillary action [Schlitt and Nicolai, 1987; Orr, 2002; Dixon, 2003]. However, due to random packing of the ore and/or uneven distribution of the solution on the surface from widely spaced drip emitters (typically 50 to 100 cm apart), the degree of saturation of the wetted zone tends to vary laterally. Solution infiltration occurs in the region along the axis of the emitter, which is small compared to the targeted volume of a particular emitter. The crosssectional area of the targeted volume under an emitter in a network with emitter spacing of 50 to 100 cm is roughly 0.25 to 1 m2. The flow from the emitter must spread over the entire cross-sectional area, but drips on only a very small area, resulting in a relatively high local flux, at least at the top of the heap. Laterally, the only mechanism for solution ingress into the zones between emitter axes is by capillary action [Schlitt and Nicolai, 1987; Orr, 2002; Dixon, 2003]. However, as governed by the laws of solution flow in porous media (discussed in Section 2.3.1), capillary forces tend to decrease dramatically as the relative wetting-phase permeability increases with increasing saturation. Furthermore, capillary forces are generally low in the coarse ores found in heaps [Orr, 2002; Dixon, 2003]. Hence, very little lateral flow of solution occurs away from the channels following the initial wetting under the emitters.  Dixon [2003] has  confirmed this by solving the Richards equation at steady state in axisymmetric coordinates, with solution being emitted from a point source at the origin of an isotropic porous medium.  He found that, although the entire bed attained a  significant degree of saturation, the actual downward flow of solution was  Introduction  6  restricted to within just a few centimeters of the vertical axis. The high local irrigation flux and little lateral spreading favour the formation and maintenance of vertical preferential flow channels along the axis of the emitter.  When this  phenomenon is particularly severe, it is referred to as ‘channelling’ [Orr and Vesselinov, 2002; Dixon, 2003]. Due to the non-uniform distribution of saturation and flow, aqueous solutes must diffuse over large distances between channels of flowing solution. The fact that diffusion is a slow process of transport means that the delivery of reagents to reaction sites, and the movement of dissolved species from reaction sites to advection channels and out of the heap, would also be slow. Rates and extents of metal extraction from heap leaching are, as a result, slow and unfavorable for the economics of the process. Of the many parameters that impact heap leaching processes, solution application rate is the one that can be mostly effectively manipulated by the operators. Other evidence of the heap hydrology problem is the absence of a scale-up procedure of laboratory data for full-scale heap design. Due to the usually small diameters of laboratory columns (typically 50 to 200 mm) used in heap leach amenability tests, the effect of inadequate and non-uniform distribution of the lixiviant is typically not detected. This will mislead full-scale heap design and analysis, if only column data is used for full-scale design. It has been realized that the rates and extents of leaching obtained in full-scale heaps nearly always fall below those obtained from column tests in the laboratory [Petersen and Dixon, 2003]. This is due to the different lengths of the diffusion distances in confined beds as found in laboratory columns and in unconfined beds as in commercial heaps. In a laboratory column, depending on the size of the column and the mode of irrigation, the length of these pathways may range from the radius of an ore particle or agglomerate to the radius of the column, times a suitable factor for tortuosity [Dixon and Petersen, 2003]. In a full scale heap under drip emitters, this distance may be as large as half the separation distance between drip points, again times a suitable factor for tortuosity [Dixon and Petersen, 2003].  Introduction  7  Lastly, solid reaction products can also cause hydrological problems in heap leaching operations. Heap leaching is a complex unit operation, comprising a series of interacting physico-chemical processes, including transport of solutes, solution and air through a porous unconsolidated pile of ore with a sequence of highly interdependent gas-liquid-solid reactions. What is more, the physical, chemical and biological conditions within a heap are continuously changing with time and space, resulting in a varying pore structure of the bed [Batarseh and Stiller, 1994].  For example, pore opening may result from dissolution and  product removal, while pore plugging may result from poor heap construction techniques and haulage methods, areas sealed off by migration of fines during leaching by the influence of high solution flows, precipitation of iron compounds such as jarosite (KFe3(SO4)2(OH)6) and scorodite (FeAsO4.2H2O), formation of elemental sulfur (in heap bioleaching), formation of lime scales during alkaline cyanide leaching of gold and silver, and decrepitation of the ore [Schlitt and Nicolai, 1987]. The evolution of hydrological problems results in considerable heterogeneity and, sometimes, wide variations of hydraulic properties within the heap. Typically, the porosity of the ore bed tends to vary, and so does the effective diffusivity, which eventually leads to variations in diffusion rates of solutes and gases, and to the formation of preferential solution flow channels through the heap [Eriksson and Destouni, 1997; Orr, 2002].  Preferential flow refers to the situation where  extreme local anisotropy results in channels of very high permeability which divert the leach solution, thus causing it to bypass large portions of the heap. Preferential flow leads to unusually long extraction times, less than expected extractions of the metal or mineral of interest and/or development of hot spots, which can lead to destruction of microbes in heap bioleaching. These problems adversely affect performance, predictability and profitability of the process. From the foregoing discussion of the problems of heap hydrology, it is clear that due to the enormous size and the long lifespan (which can be months to years) of heaps, evaluation of the effects of process design and changes on the  Introduction  8  performance of the heap can be a challenge. Mathematical modelling of the process incorporating all the significant parameters is a major tool that can be used for optimal design and various process enhancement strategies for efficient operation, as well as economic evaluations and potential environmental impact assessment of heap leaching.  1.2.1  Modelling Water Flow and Solute Transport in Heap Leaching  Heap leaching depends on a number of phenomena such as the liquid dynamics, mass transport of the lixiviant and product, mineral leaching, and auxiliary reactions (physical, chemical and biological). Solution flow and solute transport in a porous medium are dependent on the hydraulic properties of the medium including the degree of saturation and/or capillary pressure [Orr, 2002; Orr and Vesselinov, 2002; Dixon, 2003]. Thus, the concepts of unsaturated zone hydrology will generally provide the framework for understanding and optimizing the hydrological factors relevant to heap leaching [O’Kane et al., 1999]. Modelling of water flow and solute transport in unsaturated porous media is generally based on the Richards equation for soil water dynamics, and on the advection-diffusion equation (ADE) and/or the Dual Domain Model (DDM) for solute transport [Decker and Tyler, 1999; Buyuktas and Wallender, 2002].  In unsaturated porous media, the unsaturated hydraulic  conductivity and the pressure head are functions of the water content or effective saturation. Thus, for solving the models, the hydraulic functions of the ore are required. Unfortunately, at the field scale these hydraulic functions vary widely in space and time because of inherent spatial variability of the ore and heterogeneities introduced during ore handling and all the time-dependent physico-chemical factors mentioned above. A number of models have been developed over the years to simulate heap leach practice. All of these modeling efforts have much in common, accounting for air, water, heat and solute transport and consumption of acid or cyanide and oxygen. However, the equations and constitutive relations incorporated in these different  Introduction  9  models vary, with the variations reflecting emphasis on details of processes comprising the overall transport process. Some of them are briefly discussed below. Casas et al. [1993, 1998] developed a 3D model for copper bioheap leaching operations incorporating biological effects. Their model involves 3D transport of gas and heat, but ignores the effects of solution flow and solute transport. Montero et al. [1994] developed a 1D solute transport model, accounting for advection, dispersion and consumption of sulphuric acid, and for advection, dispersion and dissolution, adsorption and desorption of copper from saturated heap leach tests in columns of copper tailings. Muñoz et al. [1997] proposed a 3D solute transport model for leaching copper from tailings material under saturated conditions. Decker and Tyler [1996] analyzed flow and solute transport during cyanide heap leach decommissioning using 1D configuration of the advection-dispersion and dual porosity models. Thus they considered only axial dispersion of the solute. Pantelis et al. [2002] combined the Richards equation with the formalism of van Genuchten and Mualem [van Genuchten, 1981; Mualem, 1976] (as discussed in Section 3.2.1 as a common technique for modelling flow through unsaturated soils) to predict saturation as a function of depth and temperature. Leahy et al. [2005a,b; 2006] developed a three-phase CFD model for bioheap leaching chalcocite, accounting for air flow in 1D and 3D configurations in the heap. They assumed a homogeneous porous medium with 1D constant and uniform liquid saturation, thus accounting for only the axial dispersion, in the direction of flow. A group at the University of Swansea in Wales [Bennett et al., 2003a,b; McBride et al., 2005; Cross et al., 2006a,b] discussed a 3D comprehensive heap leach model for chalcocite. They excluded solute transport by dispersion (longitudinal and transverse) as being significant only for fully saturated flow conditions. They  Introduction  10  also failed to consider irrigation from point sources as in current industrial practice, but rather assumed uniform irrigation over a wide area. The Biohydrometallurgy Group at UBC (under the leadership of Dr David Dixon) has developed a heap leach model, called HeapSim [Dixon 2003; Dixon and Petersen, 2003], which incorporates such parameters as solution flow and solute transport, microbial growth and oxidation kinetics, mineral oxidation kinetics, oxygen absorption, and heat and gas transport [Dixon 2003; Dixon and Petersen, 2003].  Most parts of HeapSim have been developed from experimental and  theoretical work. The hydrology model within HeapSim assumes plug-flow of solution through the heap in discrete channels separated by stagnant regions accessible only by diffusion, as in the ‘Turner structure’ concept [Bouffard and Dixon, 2001; Dixon 2003; Dixon and Petersen, 2003]. The ‘Turner structure’ concept assumes vertical flow channels separated by stagnant diffusion pathways to describe advection and diffusion within heaps [Dixon and Petersen, 2003]. The models developed by Bouffard [2003] are essentially variants of the Turner structure-based models. Though good model fits using these models have been obtained from laboratory column tests [Bouffard and Dixon, 2001; Dixon 2003; Dixon and Petersen, 2003], these models are ‘non-parameterizable’ in the sense that the parameters involved have no meaningful physical analogues, and hence are only amenable to empirical treatment. The models that have 1D hydrology assume that the fluid is completely mixed in the lateral direction, and hence neglect variations in saturation with time and space.  However, 1D models generally cannot account for the effects of  geometry and the dynamics of water flow and solute transport during irrigation with drip emitters. In practice, radial mixing occurs and this should be accounted for, especially when using large diameter column test data and field situations, as required in this investigation. This is particularly important when using laboratory data for full-scale designs. In such situations, representing flow and transport in 2D, or even 3D, is necessary. For heap irrigation under a drip emitter grid, a 3D  Introduction  11  axisymmetric model will be adequate. It consists of a central axial flow, with axial and radial dispersion components. From the foregoing discussions, a comprehensive hydrology model of solution flow and solute transport in heap leaching operations, accounting for non-uniform degrees of saturation in the targeted volume of irrigation from a point source, is non-existent. The need, therefore, remains to model the transport processes in heaps based on the physics of water flow and solute transport through an unsaturated porous medium.  This will be accomplished by developing an  axisymmetric model for solution flow and solute transport in the heap under a point source of irrigation. Heap operators need such a tool to simulate multiple scenarios at high speed and low cost on computers, to aid heap design and operation decisions, and to diagnose problems in, and improve the performance of, existing operations. This is the focus of this investigation.  1.3  Objectives of this Investigation  The main objective of this work is to study and quantify the physics of solution flow and solute transport during heap leach operations under point sources of irrigation onto the surface of the heap. This main objective will be subdivided into the following: 1.  To identify the physical processes which occur during heap leaching operations under irrigation from a point source.  2.  To identify the relevant ore properties and variables that control water flow and solute transport during heap leaching.  3.  To investigate the scale and flow rate dependency of solute transport parameters and the occurrence of preferential flow in heaps.  4.  To formulate a 3D axisymmetric model for the evaluation of water flow and solute transport during heap leaching.  Introduction  12  5.  To undertake various laboratory tests to verify the model.  6.  To obtain heap leach performance data from tests at different scales for validation of the model.  7.  To demonstrate (using simulations) the application of the model to fullscale heap leach performance, focusing on the effects of emitter spacing, solution flow schedules and rates, and hydraulic properties of the ore.  The scope of this investigation will encompass the development of a hydrological model for the movement of water and solutes during heap leaching operations under isothermal and adiabatic conditions. This approach is not intended to imply that anisothermal or non-adiabatic conditions are not important.  The  anisothermal case is very important especially in heap bioleaching of sulphide ores. Inclusion of such conditions may identify other important phenomena such as the decrease in oxygen partial pressure with humidification and microbial cell growth inhibition at high temperatures [Bouffard and Dixon, 2003]. However, its inclusion at this stage would only increase the complexity of the problem. This work is intended as a first step for developing a hydrological model for heap leaching operations in general. Also the effects of gas flow, water generation and/or consumption by reactions within the heap, water evaporation from the surface of the heap and rainfall onto the heap will be considered negligible. The inclusion of the effects of temperature gradients and other conditions as mentioned above may form part of a subsequent study. Since the parameters estimated from the model will impact the performance of the heap leach process differently, a sensitivity analysis will be undertaken to identify parameters (such as solution rates, drip emitter grid spacing, heap height and ore hydraulic properties) that have dominant effects on heap performance for optimization of the process. Once the hydrological model is in place, operators will be able to use it to simulate multiple scenarios at high speed and low cost on a computer to aid design new heap leach operations and to improve the  Introduction  13  performance of existing heaps.  This will minimize the use of lengthy and  expensive laboratory and pilot-plant experimentation.  1.4  Thesis Outline  This thesis is divided into eight chapters and organized in the following manner: The introduction, motivation, objectives and scope of this thesis are described in Chapter 1. A literature review, comprising discussions of related research and existing knowledge, and applicable theory related to water flow and solute transport in porous media, is given in Chapter 2. The development of an axisymmetric heap hydrology model and techniques of its solution are given in Chapter 3. Chapter 4 describes the criteria for the design and procedures of the experimental work. The results and discussion of the experimental findings of the inert tracer tests and reactive solute transport tests are presented in Chapters 5 and 6, respectively. The results and discussion of simulations and sensitivity analysis of extraction of copper by heap leaching from the ore are presented in Chapter 7. Overall conclusions of the research and recommendations for future work are given in Chapter 8.  Introduction  14  Chapter 2 FLOW AND TRANSPORT IN UNSATURATED BEDS The objective of this chapter is to review heap leaching practice and the interdependent processes involved in solution flow and solute transport during such operations. While the focus of this study is on modeling solution flow and solute transport in heaps, it is not easy to understand the concepts without prior insight into the fundamental principles of flow and transport in unsaturated porous media.  Thus, in this chapter, the current applications of heap leach  technology and the principles of flow and transport in unsaturated porous media, and methods of evaluating them, will be reviewed.  2.1  Heap Leaching: Current Developments and Applications  The most important applications of heap leaching of oxidized ores include copper and gold, using sulphuric acid and sodium cyanide, respectively, as the leaching solutions. Heap bioleaching of refractory gold ores was also extensively piloted as a pre-treatment, required to liberate the gold occluded in sulphide mineral grains prior to cyanide or thiosulphate leaching.  Rapid expansion of heap  leaching operations worldwide coupled with improved pregnant solution recovery methods such as solvent extraction and electrowinning of copper and carbon adsorption of gold cyanide complex, has significantly decreased metal production costs, even in the face of declining ore grades. In the US, from 1995 to 1999, about one-third of gold production came from heap leach operations [Amey, 1999]. Worldwide, about 1.2 million tonnes of copper is produced annually by heap leaching and from mine waste dumps. Heap and dump leach account for about 30% of total new copper production in the US [Bartlett, 1997]. A summary of commercial applications of heap leaching for the extraction of specific minerals or metals is given below.  Flow and Transport in Unsaturated Beds  15  2.1.1  Heap Leaching of Gold Ores  Gold heap leach operations generally treat oxidized ores containing little or no sulphide mineralization. These ores are readily leached with cyanide solution. The dissolution of the gold grains is usually rapid due to their small size and the negligible consumption of cyanide from reacting with the gold. Thus the initial injection of cyanide into the heap is usually adequate to quickly leach all accessible gold [Bartlett, 1999]. As a result some researchers [Bartlett, 1999, Sanchez-Chacon and Lapidus, 1997] believe that the gold extraction rate during heap leaching from these ores can usually be quantitatively described by considering only the diffusion of dissolved gold out of the rock matrix, through solution-filled rock pores. This implies that cyanide concentration (especially in the absence of cyanicides) would be fairly uniform from the top to the bottom of the heap. Mount Poseidon mine in Australia successfully operated a Cu-Au heap leach facility. The ore was crushed and heap leached with water for a few weeks. The leached Cu was recovered from the solution by cementation, following which the leach solution was neutralized, clarified and reapplied to the heap. After rinsing, the heap was reclaimed, agglomerated with lime, and placed adjacent to its former position. The second heap was then leached by conventional methods for gold recovery [Marchbank et al., 2003]. For refractory gold ores, heap bio-oxidation is used as a pre-treatment step to liberate the gold occluded by sulphide minerals prior to conventional leaching with cyanide. Newmont Mining Corporation commissioned a 3.8 Mt/annum heap bioleach facility at its Carlin Operation in Nevada in December 1999 [Bhakta and Arthur, 2002]. The average gold grade was 2.5 g/t with only 23% of the gold being cyanide-soluble prior to bio-oxidation. The run-of-mine (ROM) ore was crushed to 80% passing 12.7 mm and heap leached with microbes operating between 32 and 50°C. After the heap bio-oxidation cycle, the ore was reclaimed, neutralized, ground (to 65% passing 75 µm), and leached with cyanide in a  Flow and Transport in Unsaturated Beds  16  conventional CIL circuit [Bhakta and Arthur, 2002]. The average gold recovery for the first three months of operation was between 55 and 60%, a little short of the target of 65% gold recovery. Newmont Mining is also conducting testwork on extracting Au from low-grade carbonaceous refractory ore by heap bio-oxidation prior to heap leaching using ammonium thiosulphate [Bhakta, 2003]. Gold dissolves during heap leaching in the same way as it does in aerated cyanide solution, according to Elsner’s equation [Dai et al, 2003]: 4 Au + 8 CN– + O2 + 2 H2O → 4 Au(CN)2– + 4 OH–  (R2.1)  Where gold is leached by copper (II) catalyzed ammonium thiosulphate solution, the reaction may be expressed as follows [Muir and Aylmore, 2004]: Au + 5 S2O32– + Cu(NH3)42+ → Au(S2O3)23– + Cu(S2O3)35– + 4 NH3  2.1.2  (R2.2)  Heap Leaching of Copper Ores  The copper industry has experienced more heap and dump leaching practice than any other. Approximately 25% of world copper production is estimated to come from that source [Marchbank et al., 2003]. Some operations are simple oxide heap leach operations, for example Manteverde Copper Mines in Chile [Zarate and Kelly, 1995]; others are dumps, and recently sulphide heaps are being operated successfully [Songrong et al., 2002, Marchbank et al., 2003]. Common copper oxide minerals include malachite [CuCO3.Cu(OH)2], azurite [2Cu(CO3).Cu(OH)2], cuprite [Cu2O], tenorite [CuO], chryosocolla [CuSiO3.2H2O] and brochantite [CuSO4.3Cu(OH)2]. The main leachant is sulphuric acid and to a lesser extent ferric iron as shown in the following reactions [Gupta and Mukherjee, 1990, www.straits.com.au/technology]. Cu2CO3(OH)2 + 2 H2SO4 → 2 CuSO4(aq) + 3 H2O + CO2  (R2.3)  Cu3(OH)2(CO3)2 + 3 H2SO4 → 3 CuSO4(aq) + 4 H2O + 2 CO2  (R2.4)  Flow and Transport in Unsaturated Beds  17  Cu2O + H2SO4 → Cu + CuSO4(aq) + H2O  (R2.5)  Cu2O + H2SO4 + Fe2(SO4)3 → 2 CuSO4(aq) + H2O + 2 FeSO4  (R2.6)  CuO + H2SO4 → CuSO4(aq) + H2O  (R2.7)  CuSiO3·2H2O + H2SO4 → CuSO4(aq) + SiO2 + 3 H2O  (R2.8)  Cu4(OH)6SO4 + 3 H2SO4 → 4 CuSO4(aq) + 6 H2O  (R2.9)  Sulphide heap leaching evolved from oxide heap leach operations.  Minera  Pudahuel developed the so-called “thin layer” leaching method for the Lo Aguirre property in Chile. This was followed by similar operations at Cerro Colorado and Quebrada Blanca [Marchbank et al., 2003]. In Australia, Girilambone and Nifty have operated extensive heap leach operations, pioneering the use of low-pressure air injection in the base of the heaps [Marchbank et al., 2003]. The first commercial copper heap bioleach facility in China was started in 1997 at the Dixing Copper Mine.  In 2001, China’s second heap bioleach plant was  started at the Zijin Mine, treating digenite and covellite for the extraction of copper [Songrong et al., 2002]. Oxidation of base metal sulphide minerals occurs in acidic ferric sulphate solutions, generating ferrous ions and base metal ions. For the case of a simple metal sulphide MeS, the reaction may be given as follows [Dixon, 2003]: MeS + Fe2(SO4)3 → MeSO4(aq) + 2 FeSO4 + S  (R2.10)  where Me is a divalent metal cation. Secondary copper sulphide minerals, including chalcocite (Cu2S), covellite (CuS) and bornite (Cu5FeS4), are the most amenable to heap bioleaching. Chalcocite is leached using microbially-mediated oxidation, using ferric sulphate and sulphuric acid, in a two-stage process as follows [Petersen and Dixon, 2003]:  Flow and Transport in Unsaturated Beds  18  Cu2S + 0.8 Fe2(SO4)3 → 0.8 CuSO4(aq) + Cu1.2S + 1.6 FeSO4  (R2.11)  Cu1.2S + 1.2 Fe2(SO4)3 → 1.2 CuSO4(aq) + 2.4 FeSO4 + S  (R2.12)  Natural covellite is leached as follows: CuS + Fe2(SO4)3 → CuSO4(aq) + 2 FeSO4 + S  (R2.13)  The ferric ions are re-generated by microbially-assisted oxidation of ferrous ions: 2 FeSO4 + 0.5 O2(a) + H2SO4 → Fe2(SO4)3 + H2O  (R2.14)  The aqueous oxygen is obtained from dissolution of oxygen gas [Dixon and Petersen, 2003]: O2(g) → O2(a)  (R2.15)  The elemental sulphur generated in reactions R2.10, R2.12 and R2.13 is further oxidized to sulphuric acid (in a parallel reaction to R2.14) by the catalytic action of microbes [Petersen and Dixon, 2003, Dixon and Petersen, 2003]: S + 1.5 O2(a) + H2O → H2SO4  (R2.16)  When pyrite is present in the ore, it may be bio-oxidized at high solution redox potentials as follows [Petersen and Dixon, 2003]: 2 FeS2 + 7.5 O2 + H2O → Fe2(SO4)3 + H2SO4  (R2.17)  Thus, during heap leaching of sulphide ores, the acid requirements can be met in three ways: as supplied in the feed solution, by the bio-oxidation of generated elemental sulphur and by the oxidation of pyrite.  2.1.3  Heap Leaching of Uranium Ores  Western Nuclear Inc and Union Carbide operated uranium heap leach operations around Gas Hills, Wyoming in the United States in the late 1960’s and early 1970’s [Gupta and Mukherjee, 1990].  Flow and Transport in Unsaturated Beds  An example of current application of  19  uranium heap leaching is at Denison Mines in Elliot Lake, Ontario, where underground heap bioleaching is successfully undertaken using acidified ferric sulphate solution produced from the bacterial oxidation of pyrite contained in the ore. Heap leaching contributes about 25 percent of overall production at this mine. The heap bioleaching cycle is about 18 months to obtain 70% uranium extraction from the ore [Marchbank et al., 2003].  2.2  Flow and Transport in Unsaturated Porous Media  A heap leach facility is essentially a packed-bed reactor under unsaturated flow. In hydrogeology, the term ‘flow’ refers to the movement of liquid and gaseous phases while ‘transport’ refers to the migration of solutes (for example acid, cyanide and metal ions) [Freeze and Cherry, 1979]. Solution flow through porous media has been studied in other scientific and engineering disciplines [Roman and Bhappu, 1993]. These include hydrology, chemical engineering, civil/soil engineering, agricultural engineering and petroleum engineering. The interests of hydrologists have been in ground water flow and ground water contamination by leakage from surface ponds. Chemical engineers have studied trickle- and fixed-bed reactors, which are similar in operation to leaching in columns. Soil/Civil engineers deal with flow of water in aquifers, the movement of moisture through and under engineered structures, and the associated stresses induced by moisture flow under the foundations of such structures. Agricultural engineers deal with movement of water and solutes in the root zone of soils. Petroleum engineers study the use of fluids to displace oil from oil reservoirs.  All these have contributed immensely to the  understanding of the flow of solution and the transport of solutes during heap leaching.  However, due to the wide particle size distribution of the ore, the  solution irrigation schedules in heap leaching and the scales of heap leach operations, the wealth of knowledge from these sources is still of limited use. The relevant concepts will therefore have to be adapted to experimental evidence.  Flow and Transport in Unsaturated Beds  20  Unsaturated flow is a special type of multiphase flow through porous media, with two phases, namely, air and water, in the pore channels of the solid particles in the porous medium [Freeze and Cherry, 1979; Fetter, 1993; Looney and Falta, 2000]. As water or solution infiltrates the unsaturated ore bed, the incoming water displaces the air present in the pores of the particles, increasing the moisture content of the ore bed. This process is governed by the equation of continuity, which states that the difference between the rates of inflow and outflow of water is equal to the change of water storage in the ore bed. Variations in the relative proportions of air and water in the pores cause variations in the hydraulic properties of the porous medium. These variations as well as the forces acting on both the air-water and the solid-water interfaces affect the flow regime in the medium. In terms of transport of solutes, the cross sectional area available for flow and solute movement under partially saturated conditions is limited to the part of the pore space that actually contains water. Soluble constituents are transferred from the solid phase to the liquid phase; precipitating and adsorbing constituents are transferred from solution onto the solid phase, while water vapour is carried by the gas phase [Zheng and Bennett, 2002].  2.2.1  Water Flow in Unsaturated Media  Water in an unsaturated soil is under a negative (sub-atmospheric) pressure, which results from the affinity of the water for soil/ore particle surfaces and capillary pores [Fetter, 1993; Hillel, 1998; Lal and Shukla, 2004]. This pressure is referred to in soil physics literature as capillary potential or matric potential. It is a function of the moisture content of the porous medium.  The lower the  moisture contents of the medium the higher the magnitude of the matric potential. The total medium moisture potential is the sum of the matric potential and gravitational potential, osmotic potential and electrochemical potential (Fetter, 1993; Lal and Shukla, 2004). The osmotic and electrochemical potentials are usually negligibly small and do not contribute significantly to the gradient of the  Flow and Transport in Unsaturated Beds  21  total potential [Fetter, 1993]. Thus the total moisture potential is reduced to the sum of the matric and gravitational potentials:  Φ = φm (θ ) + φz where  (2.1)  Φ = total moisture potential (Pa) φm = matric potential (Pa) θ = volumetric water content (m3 water/m3) φz = gravitational potential (Pa)  This equation can also be stated in terms of head: h=  where  Φ = hc (θ ) + z ρw g  (2.2)  h = total pressure head (m) ρw = density of water (~1000 kg/m3) g = gravitational acceleration (9.81 m/s2) hc = capillary head (m) z = height above a certain datum (m)  The gradient of the total moisture potential constitutes the driving force for flow of water in the unsaturated porous medium. In unsaturated flow, the volume flux (sometimes referred to as the “superficial velocity”) is proportional (as in the case of saturated flow) to the driving force, namely the hydraulic gradient, thus: v w = −K w ( θ ) where  dh dx  (2.3)  vw = water volume flux (m3 water/m2/s) Kw = hydraulic conductivity (m/s) x = distance (m)  Flow and Transport in Unsaturated Beds  22  The negative sign shows that flow is in the direction from higher head to lower head. Equation 2.3 is Darcy’s Law for water flow in an unsaturated porous medium. It was originally developed by Darcy (in 1856) for saturated porous media, where Kw was constant and equal to the saturated hydraulic conductivity [Hillel, 1998; Lal and Shukla, 2004]. It was later adapted for use in unsaturated porous media, where Kw is a varying function of the water content. The hydraulic conductivity depends on the properties of both the fluid and the porous medium [Zheng and Bennett, 2002], thus: K w (θ ) =  where  k (θ )ρw g µw  (2.4)  k = intrinsic permeability (m2) µw = water viscosity (kg/m/s)  The retention and distribution of water and solutes in a porous medium depend on the hydraulic properties of the medium. These properties culminate in the “soil water characteristic curve” and the “hydraulic conductivity function.” They are highly dependent on pore size distribution and hence on fragment and/or aggregate size distributions, the mode of ore placement, the degree of saturation of the ore heap, and solution application methods and rates (boundary conditions). Insight into the effects of varying some of these initial conditions (e.g., by adjusting the water content of the ore prior to stacking) and the boundary conditions (by adjusting the application methods and rates) and incorporating them into a flexible mathematical model will help heap operators to effect changes in the flow regimes and flow path distributions in the heap to improve leach efficiency [Orr, 2002]. These properties are discussed below.  2.2.1.1  The Soil Water Retention Curve  The soil water characteristic curve (SWCC) describes the relationship between the matric potential and the volumetric water content in a porous medium [Freeze and Cherry, 1979; Fetter, 1993; Hillel, 1998]. If the water in a saturated soil is  Flow and Transport in Unsaturated Beds  23  put under a slight suction, initially no air inflow occurs until the suction is increased to a certain critical value at which the largest surface pore starts draining and its water content is displaced by air [Freeze and Cherry, 1979; Fetter, 1993; Hillel, 1998].  This critical suction is referred to in soil science  literature as the air entry suction (AES) [Freeze and Cherry, 1979; Hillel, 1998]. AES is generally very small in coarsely textured and well-aggregated materials, and is much larger in dense, poorly aggregated and finely textured soils, because coarsely textured soils have a higher number of larger pores than finely textured soils. As suction increases, the relatively large pores drain first, followed by those of intermediate size and finally by small pores. At a very high suction, only the water that is adsorbed onto the soil surfaces remains, and this is referred to as residual water content [Freeze and Cherry, 1979; Hillel, 1998]. This residual water can only be removed by evaporation. The SWCC may be obtained in two ways; namely, by desorption and by sorption.  In desorption, the sample is  initially saturated and increasing suction (in a step-wise manner) is gradually applied to dry the soil while taking successive measurements of wetness versus suction.  In sorption, an initially dried soil sample is gradually wetted while  reducing the suction incrementally [Hillel, 1998]. The amount of water remaining at equilibrium at any given suction is a function of the sizes and volumes of the pores and the amount of water adsorbed onto the particles.  Hence, it is a  function of the matric potential [Hillel, 1998]. This function is usually measured experimentally in the laboratory using Tempe cells and graphically represented by a water content vs. matric potential curve called the soil-moisture retention, or characteristic, curve [Arya and Paris, 1981; Hillel, 1998; Arya et al., 1999]. The field SWCC determination consists of measurements of soil water pressure by tensiometers installed at depths of interest and of water content by gravimetric sampling or neutron or gamma attenuation techniques [Arya and Paris, 1981]. The resultant curves (water content vs. matric potential) obtained from wetting and the drying SWCC tests are usually not identical.  The equilibrium soil  wetness at a given suction is greater in desorption (drying) than in sorption  Flow and Transport in Unsaturated Beds  24  (wetting). This phenomenon is called hysteresis [Hillel, 1998]. Hysteresis in SWCC is attributed to four causes; namely, non-uniformity of pore spaces, contact angle differences between an advancing meniscus and a receding one, encapsulation of air in dead-end pores, and differential changes in the soil due to swelling, shrinkage and/or aging phenomena [Hillel, 1998].  2.2.1.2  Permeability of Unsaturated Porous Medium  The hydraulic conductivity function describes the relationship between the unsaturated hydraulic conductivity and the water content [Freeze and Cherry, 1979; Fetter, 1993, Hillel, 1998].  In an unsaturated porous medium, the  presence of two phases of fluid (i.e., water and air) in the voids of the medium tends to impose opposite effects on the magnitude of hydraulic conductivity with respect to each phase. In the case of the water phase, the presence of air in the pores of the medium decreases the magnitude of the coefficient of permeability as the degree of saturation decreases. As water content or degree of saturation and matric suction are related, the coefficient of permeability is a function of water content (or degree of saturation), and is therefore also a function of the medium matric suction. The coefficient of permeability decreases as the matric suction increases in magnitude. The hydraulic conductivity of a porous medium is at its maximum at 100% saturation. The unsaturated hydraulic conductivity is a non-linear function of both degree of saturation and matric suction. This relationship is illustrated in Figure 2.1 for finely and coarsely textured materials [Lal and Shukla, 2004]. The same trends were found in the work of Newman et al. [1997] and O’Kane et al. [1999]. As shown in Figure 2.1, at saturated conditions, the coarsely textured material has a higher hydraulic conductivity than the finely textured material due the higher number of large pores in the coarse material.  However, with a little  suction, the large pores drain their water, and the hydraulic conductivity of the coarsely textured material decreases due to a decrease in the cross-sectional area available for water flow, increased tortuosity and viscosity of the system  Flow and Transport in Unsaturated Beds  25  [Hillel, 1998; O’Kane et al., 1999].  The hydraulic conductivity of the finely  textured material will remain near its saturated permeability as the suction increases, due to the higher number of small pores in the finely textured material. In other words, the finely textured material has greater water retention ability than the coarse material. As the suction increases further, the fine material will also start draining its water. The rate at which hydraulic conductivity decreases for any unsaturated porous medium is a function of the gradation of the material, and tends to decrease more rapidly for uniform materials than for well-graded materials [Hillel, 1998; O’Kane et al., 1999].  Figure 2.1  Typical relationship between permeability and degree of saturation for finely textured and coarsely textured materials [after Lal and Shukla, 2004].  Another important observation in Figure 2.1 is that even though the saturated hydraulic conductivity of the coarsely textured material is higher than that of the finely textured material, the conductivities of both the coarse and fine materials decrease steeply with increasing suction.  The implication of this trend in  unsaturated bed hydrology is that transport processes taking place in wet conditions are faster than in dry conditions in the same porous medium [Hillel, 1998].  Flow and Transport in Unsaturated Beds  26  The permeability-saturation functions of the coarsely textured and finely textured materials cross at a specific degree of saturation, as shown in Figure 2.1. During the construction of a heap, segregation of coarse and fine materials may result which can affect the entire hydrology of the heap. Solution application rates that are greater than the saturated permeability of the finely textured material will lead to preferential flow through the coarse material. The high solution flux will create a pressure profile leading to saturation conditions where the effective permeability of the coarse material is greater than that of the fine material [O’Kane et al., 1999]. For an irrigation rate that is lower than the saturated hydraulic conductivity of the fine material, preferential flow in the finely textured material will result. The low application rate creates a pressure profile where the coarse material drains while the fine material retains water and its effective permeability is higher than that of the coarse material at a given suction. O’Kane et al. [1999] obtained from experimental work some interesting results which illustrate the dependence of preferential flow on segregation due to particle size differences in an unsaturated porous bed. They filled a laboratory column with inert material intentionally segregated as shown in Figure 2.2. The surface of the column material was uniformly irrigated with water using a peristaltic pump. At steady state, the amounts of water collected from below the fine and coarse fractions were compared to illustrate the preferential flow due to the flux applied to the surface. When the applied flux was slightly higher than the saturated permeability of the finely textured material, approximately 5% and 95% of the water was collected from the fine and coarse compartments, respectively, as shown in Figure 2.2. On the other hand, when the flux applied at the surface was less than the saturated permeability of the finely textured material, approximately 68% and 32% of the water was collected from the fine and coarse compartments, respectively.  Flow and Transport in Unsaturated Beds  27  Figure 2.2  Preferential flow induced by irrigation rate and particle size [after O’Kane et al., 1999]  The direct measurements in the laboratory and in the field of both SWCC and permeability functions are tedious, time-consuming and expensive [Arya and Paris, 1981; Hwang and Powers, 2003], and the results are variable, prone to error, and may only be applicable over a narrow range of saturation [Arya et al., 1999]. As a result, efforts have been and are being made to develop indirect techniques of estimating these functions.  Two main approaches have been  adopted. The first is to relate easily measurable soil properties (including particle size distribution, particle density, bulk density, and structure- and texture-related  Flow and Transport in Unsaturated Beds  28  parameters) to the water contents at specified pressures using pedo-transfer function (PTF) techniques [Looney and Falta, 2000; Arya and Paris, 1981; Arya et al., 1999a,b; Nimmo, 1997; Zhuang et al., 2001]. The second involves the estimation of the functions from power-law curves with constants that must be determined empirically. Several such formulae have been proposed including those of Brooks and Corey [1964], Mualem [1976] and van Genuchten [1981]. Although both the PTF and power-law based formulations are used to predict the SWCC and permeability functions, both lack the physical basis to account for the effects of texture and packing characteristics (structure) of the porous medium. Furthermore, each only represents the functions over a certain limited range, and fails to trace the shape of the curves over certain ranges, especially near saturation [Hillel, 1998]. Mbonimpa et al. [2002] cautioned that these predictive models should not replace testing, but rather they should be additional tools for analysis and decision-making during the preliminary stages of a project. The understanding of this complex flow behaviour under unsaturated conditions is required for meaningful operation strategies and modeling of flow of solution and transport of solutes in unsaturated porous media.  The causes and  implications of preferential flow during heap leaching operations will ultimately flow from this understanding.  2.2.1.3  Preferential Flows  As discussed earlier, the primary purposes of irrigating the heap are to transport solutes for leaching to the various parts of the heap, to transport the products of reaction out of the heap and to maintain heap temperatures within the limits of microbial activity in the case of heap bioleaching.  However, one deleterious  aspect of heap irrigation under unsaturated conditions is the preferential flow of solution [Orr, 2002; Dixon and Petersen, 2003].  This is the situation where  regions of the heap are by-passed by the flowing liquid. Researchers of tricklebeds [e.g., Funk et al., 1990] and metallurgists using heap and dump leaching [Howard,1968; Jacobsen, 1972; Roman et al., 1974; Woolery et al., 1978; Murr  Flow and Transport in Unsaturated Beds  29  et al., 1982; Eriksson and Destouni, 1997] have long been aware of the challenges posed by preferential flow of solution.  Several recent field  investigations confirm that solution flow in heaps and dumps tends to concentrate in preferred pathways and thus bypass much of the ore bed, in extreme cases leading to channelling [Eriksson and Destouni, 1997; Orr, 2002; Orr and Vesselinov, 2002; Petersen and Dixon, 2003]. Preferential flow is initiated, promoted and influenced by different heap structures, by pre-treatment of the ore, by the composition of the solution (lixiviant) and the solution application rate and schedule [Orr, 2002; Orr and Vesselinov, 2002]. For example, the practice of irrigating leaching solution over the surface of the heap from drip emitters, on coarse grids of 50 to 100 cm spacing (with flow from each dripper corresponding to the application rate to be spread over the entire surrounding area of roughly 0.25 to 1 m2) results in solution influx from point sources at relatively high local flux. For the unsaturated flow that is desired in heaps, such point sources may lead to the formation and maintenance of preferential flow channels [Dixon and Petersen, 2003]. Observations at copper heap leaching operations in Chile lend support to the occurrence of this phenomenon, and simple experiments undertaken by Dixon and Petersen [2003] also seem to confirm this. There is ample evidence of channelling of the solution between clusters and/or individual particles [Schlitt and Nicolai, 1987; Eriksson and Destouni, 1997; Orr, 2002; Orr and Vesselinov, 2002; Dixon and Petersen, 2003].  Preferential flow of solution during heap  leaching leads to unusually long times of extraction, less than expected extraction of the metal or mineral of interest, and/or the development of hot spots, which can lead to destruction of microbes in heap bioleaching. During heap decommissioning studies, solution channelling has been pointed out as one of the factors responsible for the delay in meeting solute concentration limits in the effluent [Catalan and Li, 1998]. Preferential flow in porous media can be attributed to one of three generic causes; namely, short-circuiting flow, finger flow and funnel flow [Freeze and  Flow and Transport in Unsaturated Beds  30  Cherry, 1979]. A comprehensive coverage of the mechanisms of formation of preferential flows has been provided by Freeze and Cherry [1979].  Howard  [1968] attributed non-uniform solution flow to the solution encountering impermeable and/or semi-permeable layers in the dump or heap which form due to the inclusion of clay, compaction due to the weight of loading trucks, or impermeable layers of iron precipitates. Despite preferential flow (channelling), the heap is nonetheless fully wetted [Dixon and Petersen, 2003], as moisture is drawn by capillary action into the entire pore space of the heap, albeit to a point that is generally not enough to induce the solution to flow. This implies that stagnant pore solution would exist in the bed. Reagents and dissolved products must therefore migrate into particle pores or into clusters (through these stagnant films of solution) by diffusion, which, being a slow process, retards the delivery of reagents to the reaction front and the removal of dissolved products from the reaction front through the stagnant films surrounding the particles or clusters to reach the flowing channels and be removed from the heap.  2.2.2  Solute Transport  Solutes in the heap may either be added (for example, sulphuric acid to copper heaps or cyanide to precious metal heaps) or generated by desired or undesired reactions within the heap. Undesired reactions include consumption of cyanide by non-precious metals and consumption of sulphuric acid by gangue silicates in copper oxide or sulphide leaching. The desired reactions comprise dissolution of the desired metals and generation of acid in sulphide leaching. Other reactions such as loss of water due to evaporation, precipitation, adsorption and other equilibrium reactions also affect reagent, heat and product balances and have varying effects on the hydraulic properties of the heap as mentioned in Chapter 1. Solutes in solutions infiltrating porous media move with the flowing water as well as within the medium due to concentration gradients. Understanding the transport of solutes in heap leaching is important to optimizing the use of lixiviant  Flow and Transport in Unsaturated Beds  31  and hence the irrigation system and the transport of dissolved metals out of the heap. The transport of solutes in a porous medium involves two main mechanisms; namely, convection or advection, and hydrodynamic dispersion.  Each is  discussed briefly in the following paragraphs.  2.2.2.1  Advection  Advection is the process whereby solutes are transported by the bulk mass of the flowing fluid [Freeze and Cherry, 1979]. In the absence of other processes, the solution and the solute will move at the same rate and in the same direction with a piston displacement [Schwartz and Zhang, 2003; Lal and Shukla, 2004]. The driving force is the hydraulic gradient [Spitz and Moreno, 1996]. The mass flow of water in the porous medium carries with it an advective flux of solutes proportional to their concentration [Zheng and Bennett, 2002; Schwartz and Zhang, 2003] as given by Equation 2.5: na = v w c = −K w c where  dh dx  (2.5)  na = advective molar flux (kmol/m2/s) vw = water volume flux (m3 water/m2/s) c = solute concentration (kmol/m3 water)  In a highly permeable medium, advection is the most important transport mechanism by which solution travels to the discharge point, and depends on the hydraulic gradient imposed by the boundary conditions [Freeze and Cherry, 1979; Spitz and Moreno, 1996].  As mentioned above, the solution usually  follows the path of least resistance. Unknown preferred pathways, especially at the large scale, can lead to errors during transport predictions. Advective flow is more complex when the density (and/or viscosity) of the solution varies with solute concentration [Freeze and Cherry, 1979].  Flow and Transport in Unsaturated Beds  32  2.2.2.2  Hydrodynamic Dispersion  Hydrodynamic dispersion is the process that causes a zone of mixing to develop between a fluid of one composition which is adjacent to, or being displaced by, a fluid of different composition [Schwartz and Zhang, 2003].  Hydrodynamic  dispersion, if present in flow through a porous medium, will cause the mass of solute to spread beyond the region that can be reached by advection alone. This spreading phenomenon causes dilution of the solute [Freeze and Cherry, 1979]. Hydrodynamic dispersion occurs by two processes; namely, molecular diffusion and mechanical dispersion.  2.2.2.3  Molecular Diffusion  Aqueous diffusion is a spontaneous process whereby dissolved mass is transported from a higher potential to a lower potential by random molecular motion, and so it does not depend on bulk flow of solution. The driving energy responsible for diffusion of ions, atoms, or molecules is due to gradients in chemical potential or (for dilute solution) concentration [Lim et al., 1998]. As such, diffusion is irreversible [Spitz and Moreno, 1996] and tends to decrease concentration gradients and restore homogeneity.  The movement of solutes  occurs in the opposite direction of the concentration gradient, and may be expressed according to Fick’s law [Zheng and Bennett, 2002]: ndiff = −Ddiff (θ )  dc dx  where  ndiff = diffusive molar flux (kmol/m2/s)  and  Ddiff = effective coefficient of molecular diffusion (m2/s)  (2.6)  Diffusion coefficients in porous media are generally smaller than those in pure solution because collisions with the solid matrix due to the tortuous flow paths in the porous medium hinder the diffusion process [Freeze and Cherry, 1979; Zheng and Bennett, 2002]. The tortuous pore passages in the porous medium render the actual diffusion length longer than the direct path length.  In  Flow and Transport in Unsaturated Beds  33  unsaturated porous media, such as heaps, the effective diffusion coefficient depends on the degree of saturation of the medium and decreases with decrease in water content due to increase in the diffusion path length [Lim et al., 1998]. Other factors contributing to this decrease in diffusion with decrease in wetness include increased viscosity of the liquid adjacent to soil particles, increased ionic interaction along small pores and water films, and interconnectivity between the water-filled pores and the water films [Hillel, 1998, Lim et al., 1998]. Diffusion is a slow process and it is often considered negligible compared to advection at high flow velocities and in media with high hydraulic conductivities.  2.2.2.4  Mechanical Dispersion  Mechanical dispersion is a microscopic process that occurs because of velocity variations in a porous medium compared to the average pore water velocity, and it is caused by three mechanisms. The first mechanism is the result of variations in the travel velocities of molecules at different points across the pores due to the drag exerted on the fluid by the roughness of the pore surfaces. According to Poiseuille’s Law, the laminar flow velocity varies with radial distance across a capillary tube, thus [Hillel, 1998]:   r2  u(r ) = 2v 1 − 2   R  where  (2.7)  u = velocity (m/s) v = volume flux or superficial velocity or average velocity (m/s) r = radial distance from the centre of tube (m) R = tube radius (m)  At the tube wall, the velocity u = 0. At the center of the tube, the velocity u is at its maximum and equal to twice the average velocity v. The velocity of the solute molecule carried by the advective (convective) stream therefore depends on its position within the pore passage.  Flow and Transport in Unsaturated Beds  34  The second mechanism is the result of differences in pore sizes, shapes and pore surface areas along which the solution flows, leading to different fluid bulk velocities. Another variant of Poiseuille’s Law gives the volumetric flow through a capillary tube as a function of the pressure head gradient, thus [Hillel, 1998]: Qw = πR 2v w =  πR 4 ρw g dh 8 µw dx  (2.8)  In any given cross-sectional area, the pore diameters may vary from micrometres to centimetres, and the volumetric flow rate is proportional to radius to the fourth power.  This results in widely varying bulk liquid velocities, thus causing  mechanical dispersion. The third mechanism is the result of fluctuations in the flow paths of a fluid element compared to the mean flow direction caused by tortuosity, branching and inter-fingering of pore channels. The pore passages in porous media are tortuous such that the actual path length of diffusion is greater than the apparent straight-line distance. In an unsaturated porous medium, the tortuosity increases as the moisture content decreases. In heaps, fluid particles follow tortuous flow paths around ore particles, and therefore the velocities vary in direction from one point to another. Also, the magnitude of the velocity varies, and is greater in the large pores than in the small pores. For isolated pores or ‘dead-end’ pores the velocity is even zero [Hillel, 1998; Dixon and Petersen, 2003]. Solute flux by mechanical dispersion is also given according to Fick’s Law, thus [Spitz and Moreno, 1996; Schwartz and Zhang, 2003]: ndisp = −Ddisp (θ )  dc dx  where  ndisp = molar flux by mechanical dispersion (kmol/m2/s)  and  Ddisp = coefficient of mechanical dispersion (m2/s)  Flow and Transport in Unsaturated Beds  (2.9)  35  Macroscopically, dispersion is similar to diffusion; however, unlike diffusion, it occurs only as a result of flow [Freeze and Cherry, 1979; Lal and Shukla, 2004]. As a result of this macroscopic similarity, the two processes are additive and the sum of mass fluxes due to molecular diffusion and mechanical dispersion yield the net mass flux due to hydrodynamic dispersion, thus: nd = ndiff + ndisp = −Ddiff (θ )  dc dc dc − Ddisp (θ ) = − D( θ ) dx dx dx  where  nd = molar flux by hydrodynamic dispersion (kmol/m2/s)  and  D = coefficient of hydrodynamic dispersion (m2/s)  (2.10)  Finally, the net molar flux is the sum of advective and dispersive fluxes, thus: n = na + nd = v w c − D  2.2.2.5  dc dx  (2.11)  Dispersion and Scale  The spreading of solute in the direction of bulk flow is known as longitudinal dispersion (LD) while spreading in the directions normal to the bulk flow is called transverse dispersion (TD) [Freeze and Cherry, 1979; Zheng and Bennett, 2002]. LD is normally much stronger than TD [Freeze and Cherry, 1979].  One  characteristic feature of the dispersive process is that it causes spreading of the solute in directions normal to the flow path as well as along the flow path. The amount of mechanical dispersion is a function of the average linear velocity in the direction of flow [Fetter, 1993; Schwartz and Zhang, 2003]. Thus, the coefficients of hydrodynamic dispersion in the longitudinal and transverse directions are often taken as linear functions of water volume flux, with coefficients which are functions of volumetric water content, thus: DL = α L (θ ) v w + εD (θ ) D0  (2.12)  DT = αT (θ ) v w + εD (θ ) D0  (2.13)  Flow and Transport in Unsaturated Beds  36  where  DL = longitudinal coefficient of hydrodynamic dispersion (m2/s) DT = transverse coefficient of hydrodynamic dispersion (m2/s) αL = longitudinal dispersivity (m) αT = transverse dispersivity (m) εD = effective diffusivity factor (–)  and  D0 = molecular diffusivity (m2/s)  Thus, if a solute is transported along the flow path, it spreads in all directions in the horizontal plane [Freeze and Cherry, 1979]. The total mass of solute would be conserved in the domain of flow but the mass occupies an increasing volume of the porous medium. Even though the medium may be isotropic (directionally uniform) with respect to textural properties and hydraulic conductivity, mechanical dispersion is directionally dependent. Dispersion is thus anisotropic, being stronger in the direction of flow (longitudinal dispersion) than in directions normal to flow (transverse dispersion) [Freeze and Cherry, 1979; Zheng and Bennett, 2002]. This implies that 1D transport models may only be useful in the interpretation of laboratory column data, but not in analyzing field data. At high solution flows where mechanical dispersion is the dominant mechanism of dispersion, the migrating solute cloud is ellipsoidal. In the same light, at low solution flow rates, where molecular diffusion is the dominant mechanism of dispersion; the migrating solute cloud is circular or spherical [Freeze and Cherry, 1979]. As noted above, microscopic heterogeneity and tortuosity of a porous medium are responsible for dispersion.  These factors determine the dispersivity in  laboratory scale experiments, resulting in longitudinal dispersivities in the range of 0.1 to 1 cm [Zheng and Bennett, 2002]. However, dispersivities two to four orders of magnitude larger are required to account for observed dispersion in field scale experiments [Zheng and Bennett, 2002].  This disparity arises  because dispersion at the field scale is caused by macroscopic heterogeneities rather than pore-scale heterogeneity [Zheng and Bennett, 2002].  Thus, the  macroscopic heterogeneities mask the microscopic heterogeneities, rendering  Flow and Transport in Unsaturated Beds  37  the effect of the latter minor in comparison. These macroscopic heterogeneities may be characterized by values of hydraulic conductivity and porosity [Zheng and Bennett, 2002].  Though the scale-dependence of dispersion has been  documented for work on both saturated and unsaturated media, most common solute transport models are based on governing equations with constant dispersion coefficients. This may be so because researchers ignore its variability with scale, and/or their data fail to capture this information.  If a constant  dispersion coefficient is applied to simulate data at multiple locations for a multilocation system, model predictions will under- or overestimate the experimental data [Pang and Hunt, 2001]. Notable studies of scale-dependent dispersion includes those of Aral and Liao [1996], Huang et al. [1996], Logan [1996], Hunt [1998], and Pang and Hunt [2001]. Aral and Liao [1996] developed a 2D advection equation with a timeHuang et al. [1996] presented analytical  dependent dispersion coefficient.  solutions for scale-dependent dispersion, assuming that dispersivity increases linearly with distance until some distance after which dispersivity reaches an asymptotic value. Logan [1996] provided an analytical solution for 1D transport models incorporating rate-limited sorption and first order decay under timevarying boundary conditions, assuming an exponentially increasing coefficient of dispersion. Hunt [1998] developed 1D, 2D and 3D analytical solutions of a scaledependent dispersion equation at unsteady state with an instantaneous source and at steady state with a continuous source. Pang and Hunt [2001] presented a simpler model assuming that dispersivity increases linearly with distance, and that the coefficient of bulk diffusion is negligible. The linear [Aral and Liao, 1996; Hunt, 1998; Pang and Hunt, 2001] and exponential [Logan, 1996] scaledependencies of dispersivities in 1D (assuming constant volumetric water content) are given, respectively, thus: α = ax  Flow and Transport in Unsaturated Beds  (2.14)  38    bx  α = aL 1 − exp −   L   where  (2.15)  a, b = empirical constants (–) x = scale length (m) L = typical observation length (m)  2.2.2.6  Sorption  When a porous medium is contacted with a solution containing dissolved matter, it frequently happens that certain solutes are removed from solution and immobilized within or onto the solid matrix by electrostatic or chemical forces [Zheng and Bennett, 2002].  This process is referred to as sorption.  The  opposite process, whereby solute particles are detached from the solid matrix into solution, is called desorption. Sorption comprises both adsorption (onto solid surfaces) and absorption (into particle pores). Sorption of solute onto the solid phase during solute displacement in the heap decreases the concentration of the solute in solution. When a solution containing a particular solute is contacted with a porous solid matrix, the concentration in solution will fall as concentration on the solid phase rises, until equilibrium is established. A plot of the equilibrium concentration in solution against the equilibrium concentration on the solid is called an adsorption isotherm [Zheng and Bennett, 2002], and is often represented by a power-law equation, thus: q = K ac a  (2.16)  ∂q = K aac a −1 ∂c  (2.17)  with the following gradient:  where  q = solid phase (adsorbed) concentration (kmol/kg solid)  Flow and Transport in Unsaturated Beds  39  Ka = adsorption equilibrium constant (variable) a = exponent (–)  The equilibrium constant and exponent must be determined for each species in each porous medium. Equation 2.16 is called the Freundlich isotherm.  For  certain species, particularly at low concentrations, a may be taken as equal to one, thus giving a linear isotherm. The nature of the Freundlich isotherm assumes that the solid matrix has infinite sorption capacity, whereby the sorbed concentration increases indefinitely with solute concentration. In practice this is not the case, and there is an upper limit to adsorption.  A pseudo-linear sorption isotherm which accounts for the  maximum sorption capacity of the solid matrix, called the Langmuir isotherm, is given thus [Zheng and Bennett, 2002]: q = qmax  K ac 1 + K ac  (2.18)  with the following gradient: ∂q Ka = qmax ∂c (1 + K ac )2  where  (2.19)  qmax = maximum sorption capacity (kmol/kg solid)  If sorption occurs during solute transport in a heap then the rate of mass accumulating on the solid matrix by sorption is equal to the rate at which the interstitial solution in the matrix loses its solute.  These sorption isotherm  expressions are empirical, and the parameters can only be obtained by fitting them to experimental data.  2.3  Review of Current Methods of Modeling Heap Hydrology  In a recent review of heap and dump models, Dixon [2003] mentioned that while most models accounted for wetting and advection in heap hydrology and  Flow and Transport in Unsaturated Beds  40  transport modeling, only HeapSim addressed the issues of variable saturation and diffusion explicitly. The ANSTO conceptual model combined the Richards equation with the formalism of van Genuchten and Mualem [van Genuchten, 1981; Mualem, 1976] (discussed in Chapter 3 below) to predict saturation as a function of depth and temperature.  2.3.1  Solution Flow in Heaps  As discussed previously, lateral variation in levels of saturation due to uneven solution distribution at the heap surface under widely spaced drip emitters has been largely overlooked in heap wetting studies [Dixon, 2003].  Under point  sources of irrigation, the only mechanism for solution ingress into the zones between the drip emitters is by capillary action [Dixon, 2003]. It is also known from unsaturated flow studies in porous media that capillary forces tend to decrease quite dramatically as the relative wetting-phase permeability increases with an increasing degree of saturation. Hence, not much lateral flow of solution will be expected away from the channels following the initial heap wetting. This has also been confirmed experimentally and mathematically [Dixon, 2003; Petersen and Dixon, 2003]. Dixon [2003] established that although the entire bed attains some level of saturation, the actual downward flow of solution is restricted to within just a few centimetres of the axis of flow from under the drip emitter. This effect is accounted for in HeapSim, segregating the liquid hold-up into flowing and stagnant portions, each with its own uniform liquid mass fraction. The separation distance between adjacent flow channels is taken as a function of drip emitter spacing. Thus, the true lateral saturation profile should be predicted from formal hydrological modeling. So far all heap models (except HeapSim) assume uniform application of solution over the entire surface, and consequently, uniform saturation and downward velocity across any horizontal plane of the heap, which is inadequate in light of the experimental and mathematical evidence [Dixon, 2003; Petersen and Dixon, 2003].  Flow and Transport in Unsaturated Beds  41  2.3.2  Solute Transport in Heaps  Most of the models reviewed by Dixon [2003] accounted for solute transport by advection only. However, there is clear evidence of diffusion and dispersion of solutes during heap leaching [Bouffard and Dixon, 2001; Dixon, 2003; Petersen and Dixon, 2003]. For example, in sulphide heaps, where sulphide is ultimately oxidized to sulphate, many leaching reactions require additional acid for ferrous oxidation. This is especially true for chalcocite where the first stage of leaching generates neither sulphate nor elemental sulphur [Dixon, 2003; Petersen and Dixon, 2003].  The effect of pore diffusion on the overall kinetics of metal  extraction depends on the length of the diffusion path, and this can have a significant impact on the extraction rate in heaps with poor solution distribution [Dixon and Petersen, 2003]. However, simulations undertaken by Bouffard and Dixon [2003], show that the pore distance did not significantly affect pyrite oxidation during the heap leaching of a refractory gold ore.  The acid  concentration gradient across the radial distance is not large enough to limit the reaction. This is because acid is generated locally during the bio-oxidation of pyrite, as opposed to chalcocite bioleaching where no acid is generated locally and acid must diffuse from the advection channels to the reaction fronts to effect reaction [Dixon and Petersen, 2003]. As mentioned in Section 2.3.2, dispersive solution transport occurs in both the directions parallel (longitudinal) to and normal (transverse) to the flow axis. According to a criterion established by Mears [1971], axial dispersion effects are insignificant for particle beds of great depth and low levels of convection – qualities shared by most heaps.  This is the justification for neglecting axial  dispersion in some heap leach models (e.g., that of Sanchez-Chacon and Lapidus [1997]). However, the wide variability of pore sizes in heap material, effects of tortuosity, variations in granular solid surface characteristics, and solute concentration gradients across the interfaces of mobile and stagnant solutions could result in wide variations in the velocity of solute in the heap. Thus, the validity of adopting the Mears criterion for heaps must be reviewed to reflect the  Flow and Transport in Unsaturated Beds  42  current state of insight into heap hydrology. Transverse dispersion can be very important, especially in those situations where solution flow channels are separated by large pockets of stagnant solution. With this new understanding, Petersen and Dixon [2003] concluded that the rate-limiting step in chalcocite heap leaching is the macro-diffusion of acid. HeapSim accounts for this explicitly [Dixon, 2003].  The original version of HeapSim used the ‘Turner Structure’  concept where solution is assumed to flow downward through the heap, with diffusion of solutes into and out of side branches which may be defined as having linear, cylindrical or spherical geometry [Dixon, 2003]. Bouffard and Dixon [2001] investigated three solute transport models, which are essentially variants of the ‘Turner Structure’ concept, and found all of them to be adequate for modeling column hydrology from laboratory tracer test data. However, these models are empirical, and their applicability to full-scale heap data remains untested [Dixon, 2003].  2.4  Conclusions and Focus of this Study  Most heap models assume that the solution flows vertically downward through a heap as a front (in ideal plug flow) [Sanchez-Chacon and Lapidus, 1997; Bouffard and Dixon, 2001; Dixon et al., 1993; Dixon and Hendrix, 1993a,b] at constant velocity given by the product of the superficial velocity (typically between 5 and 15 L/m2/h) and the degree of saturation (which is also assumed to be constant). The works of Casas et al. [1993, 1998] and others (e.g., Decker and Tyler [1996]; Montero et al. [1994]; Leahy et al. [2005a,b, 2006]) identified the need for 3D configuration and resolution of water, heat and solute transport, but dwelt on uniform saturation and axial dispersion only. Others [Bennett et al., 2003a,b; McBride et al., 2005; Cross et al., 2006a,b] have neglected solute dispersion completely. As such, the hydrological parts of these models have been handled at best in a 1D configuration. The problems of channelling and the recent findings from experiments and mathematical modeling of variable degrees of saturation and spreading of applied solution during heap leaching have prompted this study from a hydrological point of view in at least a 2D  Flow and Transport in Unsaturated Beds  43  configuration. As hydrology is considered the backbone of any heap leach model [Bouffard and Dixon, 2001; Dixon, 2003], it is essential to develop a comprehensive hydrodynamic model based on sound hydrology in unsaturated porous media. Darcy’s Law (combined with suitable constitutive relationships such as the van Genuchten and Mualem equations) provides an accurate relationship for determining the velocity by accounting for the pressure gradient associated with hydrostatic head (the sum of gravitational and capillary forces), the viscosity of the solution and the local permeability. The formulation of a 2D axisymmetric model for heap leaching based on Darcy’s Law, and its solution strategies, are given in the next chapter.  Flow and Transport in Unsaturated Beds  44  Chapter 3 MODEL DEVELOPMENT In Chapter 2, the mechanisms that impact on solution flow and solute transport were reviewed. In this chapter, the principles and conventions used in hydrology and soil physics relevant to modeling of water flow and solute transport in unsaturated porous media have been employed to develop a 2D axisymmetric model and its solution strategies. Formulation of the 2D axisymmetric model of water flow in heap leaching operations and the numerical techniques of solving the model to obtain conservative and stable solutions are presented in Section 3.1.  Formulation and solution strategies for the associated solute transport  model are given in Section 3.2.  3.1  The Water Flow Model  The movement of water and solutes in porous media is governed by equations (usually referred to as governing equations) based on the principle of the conservation of mass. Thus, for mass conservation of matter in a certain fixed volume, the mass input is equal to the sum of the mass output, plus the mass accumulated within the volume.  Also, Darcy’s Law provides an accurate  relationship for determining the velocity of water through a porous medium by accounting for the pressure gradient associated with hydrostatic head, the viscosity of the water, and the local permeability of the porous medium. For effective design and operation of drip systems for the irrigation of heaps, there is the need to predict water movement through the ore bed. Modeling of water flow and solute transport begins with quantifying the fluid flow through the medium. Let us assume a heap surface being irrigated by a set of drip emitters spaced at regular intervals 2R, and each discharging leach solution at an equal rate of U kg/m2/h. Due to the symmetry of the drip emitter layout, the surface of the heap  Model Development  45  can be subdivided into identical, nearly cylindrical elements of diameter 2R and depth Z, with the drip emitter placed at the origin of the vertical (longitudinal) axis. The solution flow pattern for the heap can be conveniently described by analyzing the flow in this single representative, cylindrical volume. Thus because of the axial symmetry around the vertical axis, the infiltration process can be viewed as an axisymmetric flow as shown in Figure 3.1. 2R  2R Solution Pipe Drip emitters  Z  2R  Figure 3.1  A schematic of a heap under irrigation from drip emitters  Assuming the density of water to be virtually constant (it only varies ± 2% over the entire temperature range of stability), then water volume represents water mass and the 2D axisymmetric equation of water volume continuity (EOC) (its derivation is given in Appendix A) is written thus:  − ∇ ⋅ vw = −  where  ∂v w ,z 1 ∂(r v w ,r ) ∂θ Mw − = sw − ∂z r ∂r ∂t ρw  (3.1)  vw = water volume flux (a vector) (m3 water/m2/s) z = depth (m) r = radius (m) θ = volumetric water content (m3 water/m3) t = time (s) sw = net rate of water production (kmol/m3/s)  Model Development  46  ρw = water density (assumed to be 1000 kg/m3 water) Mw = water molecular mass (18.015 kg/kmol)  Expressing Darcy’s Law (Equation 2.3) in a 2D axisymmetric configuration, the water volume flux, assuming no pressure gradients within the air phase may be written in isotropic, unsaturated porous media thus:  vw =  k k rw ( ρ w g∇z + ∇pc ) = K w k rw (∇z + ∇hc ) µw  (3.2a)  k k rw  ∂p   ∂h   ρ w g + c  = Kw k rw 1 + c  µw  ∂z  ∂z    (3.2b)  k k rw ∂ pc ∂h = K w k rw c µw ∂r ∂r  (3.2c)  v w ,z =  v w ,r =  where  k = intrinsic permeability (m2) µw = water viscosity (kg/m/s) g = gravitational acceleration (9.81 m/s2) Kw = hydraulic conductivity = k ρw g µw (m/s) krw = relative water permeability (–) pc = capillary pressure (Pa) hc = capillary head = pc ρw g (m)  Two additional equations (called constitutive relations) required to solve Equations 3.1 and 3.2 are the ore/soil water retention curve and the permeabilitypressure head relations.  3.1.1  Constitutive Relations  Constitutive relations express interrelations and constraints of the physical processes, variables and parameters as functions of a set of primary variables selected to make the governing equations (e.g., Equations 3.1 and 3.2) solvable [Looney and Falta, 2002].  Model Development  47  One of the most popular functions for describing the ore/soil water retention curve (relating the water content and matric potential) is the Brooks and Corey (BC) equation [Brooks and Corey, 1964; Corey, 1977; Lal and Shukla, 2004]  (hc hc ,0 )− λ Se =  1   where  Se =  θ − θr θs − θr  or  hc > hc ,0 hc ≤ hc,0  (3.3)  θ = θr + (θs − θr )Se  and where λ = an exponent related to the pore size distribution (–) hc,0 = a capillary head scaling factor, often referred to within the context of the BC model as the “air entry” head (m) θs = saturated volumetric water content (also the effective macroporosity) (m3 water/m3) θr = residual (fully drained) volumetric water content (m3 water/m3) A disadvantage of the BC expression is the abrupt change of the Se curve at hc = hc,0, denoting the pressure at which the largest pore drains [Looney & Falta, 2002].  Some continuously differentiable (smooth) equations have been  proposed to improve the description of the soil moisture characteristic curve near saturation.  One such equation with attractive properties is the closed-form  equation of van Genuchten [van Genuchten, 1980]:  Se =  where  1 [ 1 + (hc hc ,0 )n ]m  (3.4)  n = an exponent related to the pore size distribution (–) m = an exponent with values over the range 0 < m < 1 (–), and often taken as m = 1 – 1/n  As can be seen from Equation 3.4, the van Genuchten equation is continuous in both Se and dSe/dhc at hc = hc,0.  Model Development  48  Using Mualem’s capillary tube bundle model [Mualem, 1976], van Genuchten derived the well-known van Genuchten-Mualem (VGM) formula which defines the relative hydraulic permeability, krw, as a function of capillary head, hc: k rw =  {1 − (hc hc,0 )n −1[1 + (hc hc,0 )n ]− m } 2  (3.5)  m  [1 + (hc hc,0 )n ] 2  Combining Equations (3.4) and (3.5) and simplifying by taking n = (1–m)–1 yields the most widely used van Genuchten-Mualem formula defining the relative water permeability, krw, as a function of water saturation, Se: S 2 (1 − y )2 S < 1 e = e where 1 Se ≥ 1 1  k rw  1  y = (1 − Sem )m  (3.6)  Equation 3.4 may be rearranged to express capillary head as a function of Se, thus: 1   y  mn hc = hc,0    Se   (3.7)  Thus, according to the VGM formalism, the relative water permeability and capillary head may be specified as unique functions of effective saturation, Se (the fraction of macro-porosity filled with flowing water) as expressed in Equations 3.6 and 3.7, respectively. Given these relationships, Darcy’s Law may be recast in terms of θ thus:  v w ,z = Uw θ − Dw  ∂θ ∂z  v w ,r = − Dw  ∂θ ∂r (3.8)  where  Model Development  Uw =  K w k rw θ  and  Dw = −K w k rw  dhc dθ  49  These equations represent the advection and dispersion of water volume. Uw has units of velocity (m/s) and Dw has units of diffusivity (m2/s) and takes the following form: K w hc,0 (1 − y )2 y mn − m Dw = 1 1 θs − θr mnSemn + 2 1  1  where  1  y = (1 − Sem )m  (3.9)  Also, each approaches zero asymptotically at zero effective saturation. This is important because it implies that both quantities may be defined as zero at effective saturations less than zero, thus:  Uw = 0  and  Dw = 0  @  θ ≤ θr  (3.10)  and that is important because it allows the problem of wetting a dry (or nearly dry) bed to be solved easily and accurately. This model degenerates at saturation since Dw → ∞ as θ → θs. However, this problem is easily avoided in practice by limiting Dw to some finite maximum value (for example, say, the value at Se = 0.99). In this way, the differential Dw·∂θ will function as the differential phreatic head ∂h upon reaching saturation. It also bears noting that, at low values of Se (as would be expected in heaps) and high values of n, the VGM model approaches the Brooks-Corey (BC) model asymptotically: 1  1  2  k rw ≈ Se2 (mSem )2 = m 2Sem  + 21  = m 2Seψ  −  1  −1  hc ≈ hc,0Se mn = hc ,0Se λ  Of course, this implies that the values of the hydraulic conductivity Kw assumed by the VGM and BC models will differ by a factor of m2 in order to fit the same set of data. Hence, it follows that neither value is necessarily the “true” value, and should therefore be viewed only as a fitting parameter at low saturations.  Model Development  50  3.1.2  Discretization of the Water Transport Model  The governing and constitutive equations (Equations 3.1, 3.2, 3.6 and 3.7) are highly non-linear and analytical solution techniques usually involve assumptions that tend to oversimplify the system and therefore, numerical solution procedures must be used. The entire heap may be divided into an array of identical cylindrical columns (excluding the sloped sides for simplicity). The height and radius of each cylinder are divided into annular control volumes of vertical thickness ∆z and radial thickness ∆r, respectively.  It is assumed that each control volume is  representative of the heap; thus, characteristics such as particle size distribution, grade, porosity, etc are the same as that of the heap.  Figure 3.2  i -1, j -1  i -1, j  i -1, j +1  i , j -1  i, j  i , j +1  i +1, j -1  i +1, j  i +1, j +1  Schematic diagram of a discretised 2D axisymmetric grid  Consider the system of control volumes and boundaries shown in Figure 3.2 [Patankar, 1980]. Under a point source (such as a dripper emitter), water or fresh leach solution enters only the upper boundary of the cylinder at i = j = 1, and spreads both outward and downward from there to the control volumes below and away from the centre. Thus the entire column of heaped material is divided into (I×J) control volumes, where I is the number of vertical divisions and  Model Development  51  J is the number of radial divisions. Patankar’s technique [Patankar, 1980] of defining vector quantities at the control volume boundaries and scalar quantities within the volumes has been adopted here. Across the (i,j)th control volume, approximating the time differential with a forward finite difference, the discrete equation of water volume continuity may be written thus:  v w( i ,−z1, j ) − v w( i ,,zj ) ( j − 1) v w( i ,,rj −1) − ( j ) v w( i ,,rj ) δw θ ( i , j ) Mw ( i , j ) sw = − + ∆z ( j − 21 ) ∆r ∆tw ρw  (3.11)  where the delta notation has the following meaning:  δw θ ( i , j ) = θt( i , j ) − θt(−i ,∆j )tw and where ∆tw is the timestep of the water volume transport problem. Volume fluxes may be written thus, defining arithmetic mean values of volume dispersivity at the nodal boundaries: v w( i ,,zj ) = Uw( i , j )θ ( i , j ) −  v w( i ,,rj ) = −  Dw( i , j ) + Dw( i +1, j ) θ ( i +1, j ) − θ ( i , j ) 2 ∆z  Dw( i , j ) + Dw( i , j +1) θ ( i , j +1) − θ ( i , j ) 2 ∆r  (3.12a)  (3.12b)  Defining the following rate coefficients (all with units of inverse time):  γ w( i , j ) =  Uw( i , j ) ∆z  βw( i,,zj ) =  Dw( i , j ) + Dw( i +1, j ) 2 (∆z )2  βw( i ,,rj ) =  Dw( i , j ) + Dw( i , j +1) 2 (∆r )2  then volume fluxes may also be rewritten thus:  Model Development  v w( i ,,zj ) = ( ∆z )[ γ w( i , j )θ ( i , j ) − βw( i,,zj ) (θ ( i +1, j ) − θ ( i , j ) )]  (3.13a)  v w( i ,,rj ) = −( ∆r )βw( i,,rj ) (θ ( i , j +1) − θ ( i , j ) )  (3.13b)  52  and the discrete 3D axisymmetric equation of water volume continuity, Equation 3.11, may be expanded thus:  δw θ ( i , j ) ∆t w    − jj −−11 βw( i,,rj −1)θ ( i , j −1)  − (γ w( i −1, j ) + βw( i,−z1, j ) )θ ( i −1, j )    2     +  + (γ w( i , j ) + βw( i,−z1, j ) + βw( i,,zj ) )θ ( i , j )  +  + ( jj −−11 βw( i,,rj −1) + j −j 1 βw( i ,,rj ) )θ ( i , j )  = ωw( i , j ) (3.14) 2 2     ( i , j ) ( i +1, j )   j i j ( , ) − βw ,z θ  − j − 1 βw ,r θ ( i , j +1)      2  where the source term is defined thus: ωw( i , j ) =  Mw ( i , j ) sw ρw  Along the heap surface, vertical water volume fluxes are specified directly, and along the heap base, a zero vertical moisture gradient may be assumed, while radial volume fluxes at the axis and outer radius are also taken as zero, leading to the following very succinct boundary conditions: ( 0, j )  (1, j ) w  ω  v M = w sw(1, j ) + w ,z ∆z ρw  (3.15) βw( 0,,zj ) = βw(I,,zj ) = βw( i,,r0 ) = βw( i,,rJ ) = 0 (Setting the β coefficients equal to zero, and thereby defining the diffusion coefficients as zero, at specific nodal boundaries effectively renders those boundaries impervious to diffusive flow.  This is the most convenient way to  define zero-flux boundary conditions when the control volume technique is employed.) For a point source at relatively low flow rates, the entire column flux uw may be fed to the control volume at the origin, and the upper boundary conditions are taken thus: v w( 0,,z1) = J 2uw  Model Development  and v w( 0,,zj ≠1) = 0  if  J 2uw ≤ K w(1,1)  (3.16)  53  At higher flow rates (or finer radial grid sizes), where the ponding criterion would be exceeded at the central volume surface, the applied flow must “spill over” to concentric nodes until the entire flow rate is accommodated, thus: j −1    1  2  J uw − ∑ (2n − 1) K w(1,n )  K w(1, j )  v w( 0,,zj ) = 0 n =1  2 j − 1     (3.17)  where the double-line notation denotes that the function in the centre is bounded by the values on the left (at the minimum) and right (at the maximum).  3.1.3  Solution of the Water Transport Model  The entire system of (I×J) balance equations may be expressed in matrix notation thus: 1 [δw θ ] + [ ∆w ,z ][θ ] + [ ∆w ,r ][θ ] = [ωw ] ∆t w  (3.18)  For an implicit solution in time which is accurate to O(∆t), the matrix equation may be recast thus:   1   + [ ∆w ,z ] + [ ∆w ,r ] [δw θ ] = [ωw ] − ([ ∆w ,z ] + [ ∆w ,r ])[θt − ∆t w ]  ∆t w   (3.19)  The ∆w,z and ∆w,r matrices on the LHS may be factored into separate terms:   1  1  1 + [ ∆w ,z ] + [ ∆w ,r ] = ∆tw  + [ ∆w ,z ]  + [ ∆w ,r ]  − ∆tw [ ∆w ,z ][ ∆w ,r ] ∆t w  ∆t w  ∆tw    1  1 + [ ∆w ,r ]  ≅ ∆tw  + [ ∆w ,z ]    ∆tw  ∆t w  (3.20)  Neglecting the final term only introduces O(∆t3) error, and thus does not affect a problem which is, at best, only O(∆t2) accurate.  Hence, we may write the  following:  Model Development  54   1  1  1  + [ ∆w ,z ]  + [ ∆w ,r ] [δw θ ] = {[ωw ] − ([ ∆w ,z ] + [ ∆w ,r ])[θt − ∆t w ]} ∆ tw  ∆t w  ∆tw   (3.21)  Finally, we can solve this system in two steps, using the solution from the first step to find the final solution in the second step, thus:    1 1  {[ωw ] − ([ ∆w ,z ] + [ ∆w ,r ])[θt − ∆tw ]} + [ ∆w ,z ] [δw θ ]⊗ = ∆t w   ∆t w  (3.22a)   1   + [ ∆w ,r ] [δw θ ] = [δw θ ]⊗  ∆t w   (3.22b)  It is important to note that the entire system of lines in either direction may be solved simultaneously by ordering the θ values such that either the ∆w,z matrix or the ∆w,r matrix becomes globally tri-diagonal. This is akin to stringing the lines together into one long, continuous line in both directions, and only requires reordering between the steps. In this way, both matrices can be defined using only a 2D array of 3×(I×J) elements, and the solution is very efficient. Finally, simulation results suggest that the solution order z → r (as shown above) is to be preferred over r → z. Due to the extreme non-linearity of the coefficients Uw and Dw, the problem is naturally very sensitive to time-step, and running r → z is significantly more so.  Also, this problem is easily formulated for O(∆t2)  accuracy in principle, but much higher time-step sensitivity overwhelms any gains in accuracy in practice.  3.1.4  Testing the Water Transport Model  Given the highly nonlinear nature of the constitutive relations in the water transport model, no analytical solution is available to compare with numerical simulations. Hence, in order to demonstrate the veracity of the water transport model code, flow from a point source was simulated at very high and very low capillary head in a 4 × 1 m column to test whether the model would produce the  Model Development  55  expected results. For very high capillary head one would expect the water to spread out uniformly across the column, while for very low capillary head one would expect the water to remain confined near the vertical axis. These results are confirmed in Figures 3.3a and 3.3b, respectively.  0.35  0.6  0.3 0.5  0.25  0.1  0.4  0.3  0.2  Effective saturation, Se  0.15  Effective saturation, Se  0.2  0.05 0.1  6  Figure 3.3  Radius, r  b)  31  31  21  21  a)  0 26  11  Depth, z  16  -0.05 S1 26  Depth, z  16  11  6  1  1  0  S1  Radius, r  Predicted column wetting patterns at extreme values of capillary head in a 4 × 1 m column: a) hc,0 = 1.0 m; b) hc,0 = 0.001 m (all other parameters as determined in Chapter 5).  Finally, running the model to steady state under both sets of parameters gave integrated vertical fluxes at every depth which were identical to the applied flow rate, thus confirming that the model is perfectly conservative, as expected.  3.2  The Solute Transport Model  As discussed in Chapter 2, solute transport in unsaturated flow depends on the combined phenomena of advection and hydrodynamic dispersion, both of which depend on the flow field. Since solute fluxes in porous media are controlled by physical transport, chemical interactions and biological processes, a general model of solute transport in an unsaturated flow field should comprise the equations for flow and solute displacement, with terms accounting for physical transport, chemical, or biological interactions of the solute with the flow  Model Development  56  environment.  Thus, the production and consumption of solute species are  accounted for in the transport equations as source and sink terms, respectively.  3.2.1  Formulation of the Solute Transport Model  For a 3D axisymmetric system, the equation of solute continuity for the transport of a single aqueous solute (derived in Appendix A) is given as follows: − ∇ ⋅ nA = −  where  ∂ nA,z 1 ∂(r nA,r ) ∂ = (θ c A ) − s A − ∂r ∂t ∂z r  (3.23)  nA = molar flux of solute A (a vector) (kmol/m2/s)  cA = molar concentration of solute A (kmol/m3 water) sA = net molar production rate of solute A (kmol/m3/s) Molar flux is given by Fick’s Law corrected for bulk flow, which, assuming dilute solution, may be written thus:  n A = v w c A − D A ⋅ ∇c A where  (3.24)  DA = dispersivity of solute A (a symmetric tensor) (m3 water/m/s)  According to the formalism of Bear [Bear, 1972], dispersivity takes different values in the direction of flow (L, longitudinal) and perpendicular to the direction of flow (T, transverse), each of which may be represented with linear functions of volume flux:  where  DA,L = α Lv w + εD DA,0  (3.25a)  DA,T = αT v w + εD DA,0  (3.25b)  vw = v w2 ,z + v w2 ,r = resultant water volume flux (m3 water/m2/s) αL = aqueous longitudinal dispersivity coefficient (function of vw) (m) αT = aqueous transverse dispersivity coefficient (function of vw) (m) εD = aqueous molecular diffusivity coefficient (function of θ) (–)  Model Development  57  DA,0 = molecular diffusivity of solute A in water (m2/s) Hence, in the L-T coordinate plane, the dispersivity tensor only has non-zero values along the diagonal, and Fick’s law may be written simply as follows: n A,L = v w c A − DA,L  n A,T = −DA,T  ∂c A ∂L  (3.26a)  ∂c A ∂T  (3.26b)  However, as volume may flow at any angle β to the vertical, in the z-r coordinate plane the dispersivity tensor has non-zero off-diagonal values, and Fick’s law must be written thus: n A,z = v w ,z c A − DA,z  ∂c A ∂c − Dzr A ∂r ∂z  (3.27a)  n A,r = v w ,r c A − DA,r  ∂c A ∂c − Drz A ∂z ∂r  (3.27b)  where DA,z = α zv w + εD DA,0  DA,r = α r v w + εD DA,0  Dzr = α zr v w = α rzv w = Drz  and where αz =  α Lv w2 ,z + αT v w2 ,r v w2  αr =  α Lv w2 ,r + αT v w2 ,z v w2  α zr =  (α L − αT ) v w ,zv w ,r = α rz v w2  It bears noting here that only the values of molecular diffusivity are related to specific solutes. Dispersivities are solely functions of the prevailing flow field. Also, molecular diffusion is usually an insignificant source of solute motion relative to hydrodynamic dispersion, and is therefore typically ignored in hydrological transport problems.  Model Development  58  3.2.2  Discretization of the Solute Transport Model  Across the (i,j)th control volume in Figure 3.2, approximating the time differential with a forward finite difference, the discrete 3D equation of solute continuity may be written thus: nA( i,−z1, j ) − nA( i,,zj ) ( j − 1) nA( i,,rj −1) − ( j ) nA( i,,rj ) δ A (θc A )( i , j ) = − s A( i , j ) + ∆z ( j − 21 ) ∆r ∆t A  (3.28)  where the delta notation has the following meaning:  δ A (θc A )( i , j ) = (θc A )(t i , j ) − (θc A )(t i−,∆j )t A and where ∆tA is the timestep of the solute transport problem. Molar fluxes may be written thus, defining arithmetic average central differences for the “shear” solute fluxes and adopting Patankar’s “Power Law” upwinding scheme to ensure that information related to advection is only propogated in the direction of flow [Patankar, 1980]:  n  (i , j ) A,z  DA( i,,zj ) + DA( i,+z1, j ) c A( i +1, j ) − c A( i , j ) = v c − φ(Pe ) 2 ∆z (i , j ) ( i +1, j ) ( i , j +1) ( i +1, j +1) ( i , j −1) D + Dzr cA + cA − cA − c A( i +1, j −1) − zr 2 4 ∆r  (3.29a)  DA( i,,rj ) + DA( i,,rj +1) c A( i , j +1) − c A( i , j ) 2 ∆r ( i +1, j ) ( i +1, j +1) ( i −1, j ) cA + cA − cA − c A( i −1, j +1) 4 ∆z  (3.29b)  (i , j ) (i, j ) w ,z A  (i , j ) A,z  n A( i,,rj ) = v w( i ,,rj )c A( i , j ) − φ(Pe(Ai ,,rj ) ) Drz( i , j ) + Drz( i , j +1) − 2  where the resultant water volume fluxes taken for calculating the dispersivities are defined within each node based on arithmetic averages thus: 2  v  Model Development  (i , j ) w   v ( i , j ) + v w( i ,+z1, j )   v w( i ,,rj ) + v w( i ,,rj +1)   +  =  w ,z    2 2      2  59  and where the upwinding function φ is defined as a function of the local boundary Péclet number to give the proper magnitude and direction of advection relative to dispersion thus: − Pe Pe < −10  (1 + 0.1 Pe )5 − Pe − 10 ≤ Pe < 0  φ(Pe ) =  5 0 ≤ Pe < 10  (1 − 0.1 Pe )  0 Pe ≥ 10 Defining the following rate coefficients (all with units of inverse time):  γ A( i,,zj ) =  v w( i ,,zj ) ∆z  γ A( i,,rj ) =  v w( i ,,rj ) ∆r  βA( i,,zj ) = φ(Pe(Ai ,,zj ) )  DA( i,,zj ) + DA( i,+z1, j ) 2 (∆z )2  β A( i,,rj ) = φ(Pe(Ai ,,rj ) )  βzr( i , j ) =  Dzr( i , j ) + Dzr( i +1, j ) 8(∆z )( ∆r )  βrz( i , j ) =  DA( i,,rj ) + DA( i,,rj +1) 2 (∆r )2  Drz( i , j ) + Drz( i , j +1) 8(∆r )( ∆z )  allows the Péclet numbers to be defined thus:  Pe(Ai ,,zj ) =  2 ( ∆z ) v w( i ,,zj ) γ A( i,,zj ) = DA( i,,zj ) + DA( i,+z1, j ) β A( i,,zj )  Pe(Ai ,,rj ) =  2 ( ∆r ) v w( i ,,rj ) γ A( i,,rj ) = DA( i,,rj ) + DA( i,,rj +1) βA( i,,rj )  With these, the solute fluxes may be rewritten thus: n A( i,,zj ) = ( ∆z )[ γ A( i,,zj )c A( i , j ) − β A( i,,zj ) (c A( i +1, j ) − c A( i , j ) ) − βzr( i , j ) (c A( i , j +1) + c A( i +1, j +1) − c A( i , j −1) − c A( i +1, j −1) )] n A( i,,rj ) = ( ∆r )[ γ A( i,,rj )c A( i , j ) − βA( i,,rj ) (c A( i , j +1) − c A( i , j ) ) − βrz( i , j ) (c A( i +1, j ) + c A( i +1, j +1) − c A( i −1, j ) − c A( i −1, j +1) )]  (3.30a)  (3.30b)  Combining gives the final equation for solute concentration within the (i,j)th control volume:  Model Development  60  δ A (θ c A ) ( i , j ) ( ∆t ) A   − jj−−11 (γ A( i,,rj −1) + β A( i,,rj −1) ) c A( i , j −1)   − (γ A( i,−z1, j ) + β A( i,−z1, j ) ) c A( i −1, j )   2     +  + ( β A( i,−z1, j ) + γ A( i,,zj ) + β A( i,,zj ) ) c A( i , j )  +  + [ jj−−11 β A( i,,rj −1) + j −j 1 (γ A( i,,rj ) + β A( i,,rj ) )] c A( i , j )  2 2     ( i , j ) ( i +1, j )   j β c − A, z A  − j − 1 β A( i,,rj ) c A( i , j +1)    2    = ω A( i , j )   β zr( i −1, j ) (c A( i −1, j +1) + c A( i , j +1) − c A( i −1, j −1) − c A( i , j −1) )     − β zr( i , j ) (c A( i , j +1) + c A( i +1, j +1) − c A( i , j −1) − c A( i +1, j −1) )  −  + j −1 β ( i , j −1) (c ( i +1, j −1) + c ( i +1, j ) − c ( i −1, j −1) − c ( i −1, j ) )    j − 1 rz A A A A 2    − j −j 1 β rz( i , j ) (c A( i +1, j ) + c A( i +1, j +1) − c A( i −1, j ) − c A( i −1, j +1) )  2    (3.30)  where the source term is defined thus: ωA( i , j ) = s A( i , j ) Along the heap surface, incoming water flow and solute concentration are specified directly, no flow crosses the axis or the perimeter, and no dispersion occurs across the boundaries. These considerations lead to the following boundary conditions:  (1, j ) A  ( i ,0 ) A,r  ω  v w( 0,,zj )  =s  (1, j ) A  +  =β  ( i ,J ) A,r  =β  ∆z  c A( in ) (3.31)  γ  ( 0, j ) A,z  =β  ( 0, j ) A,z  =β  (I, j ) A,z  =β  ( 0, j ) zr  =β  (I, j ) zr  =β  ( i ,0 ) rz  =β  ( i ,J ) rz  =0  Also, it is important to note that, although Dzr = Drz for any control volume (a result of the symmetry of the dispersion tensor), βzr ≠ βrz as these values are associated with the boundaries between control volumes, and as such represent the arithmetic averages of different Dzr or Drz values from adjacent control volumes.  Model Development  61  3.2.3  Solution of the Solute Transport Model  The entire system of (I×J) equations may be expressed in matrix notation thus:   δ A (θc A )   ∆t  + [ ∆ A,z ][c A ] + [ ∆ A,r ][c A ] = [ωA ] − ([ ∆ zr ] + [ ∆ rz ])[c A ]   A  (3.32)  Since the zr and rz coefficient matrices cannot be forced into tri-diagonal form, they must be consigned to the RHS of the equation, and will thus operate on old concentration values only.  The alternative to this would be to iterate to  convergence, although since these shear terms are expected to be very small, it should make little difference. Also, applying the discrete product rule to the delta term on the LHS: δ A (θc A ) = θt ⋅ δ A c A + δ Aθ ⋅ c A,t − ∆t A results in the following:   δ Aθ    θ   [δ A c A ] + [ ∆ A,z ][c A ] + [ ∆ A,r ][c A ] = [ω A ] −    + [ ∆ zr ] + [ ∆ rz ]  [c A,t − ∆t A ]  ∆t A    ∆t A    (3.33)  Finally, the source term may be differentiated thus:  [δ A ω A ] = [ω ′A ][δ A c A ]  where  ω ′A =  ∂ω A ∂s A = ∂c A ∂c A  and thus  [ω A ] = [ω A,t − ∆t A ] + [ω ′A ][δ A c A ]  giving the final equation for solute concentration in matrix form:  θ − ω ′A   [δ A c A ] + [ ∆ A,z ][c A ] + [ ∆ A,r ][c A ]  ∆t A   δ θ   = [ω A,t − ∆t A ] −   A  + [ ∆ zr ] + [ ∆ rz ]  [c A,t − ∆t A ]   ∆t A    Model Development  (3.34)  62  The source term derivative in the coefficient of the accumulation terms on the LHS ensures stability by eliminating the possibility of negative concentration values for reactant solutes at large solute timesteps. For an implicit solution in time which is accurate to O(∆t), this equation may be recast thus:    θ − ω ′A      ∆t  + [ ∆ A,z ] + [ ∆ A,r ]  [δ A c A ] A     δ θ   = [ω A,t − ∆t A ] −   A  + [ ∆ A,z ] + [ ∆ A,r ] + [ ∆ zr ] + [ ∆ rz ]  [c A,t − ∆t A ]   ∆t A    (3.35)  The ∆A,z and ∆A,r matrices on the LHS may be factored into separate terms:  θ − ω ′A    + [ ∆ A,z ] + [ ∆ A,r ]  ∆t A   θ − ω ′A  =   ∆t A   −1  −1    θ − ω ′A    θ − ω ′A    θ − ω ′A      ∆t  + [ ∆ A,z ]   ∆t  + [ ∆ A,r ]  −  ∆t  [ ∆ A,z ][ ∆ A,r ]} (3.36) A  A  A         θ − ω ′A  ≅   ∆t A   −1     θ − ω ′A    θ − ω ′A      [ ∆ ] [ ∆ ] + + A z , A , r     ∆t  ∆t A  A       As before, neglecting the final term only introduces O(∆t3) error, and thus does not affect a problem which is, at best, only O(∆t2) accurate. Hence, we may write the following:  θ − ω ′A     ∆t A   −1    θ − ω ′A    θ − ω ′A       [δ A c A ] [ ∆ ] [ ∆ ] + + A z A r , ,    ∆t  ∆t  A  A        δ θ   = [ω A,t − ∆t A ] −   A  + [ ∆ A,z ] + [ ∆ A,r ] + [ ∆ zr ] + [ ∆ rz ]  [c A,t − ∆t A ]   ∆t A    (3.37)  or  Model Development  63    θ − ω ′A    θ − ω ′A       [δ A c A ] + + [ ∆ ] [ ∆ ] A r A z , ,    ∆t  ∆t    A A        δ Aθ    θ − ω ′A   =  + [ ∆ A,z ] + [ ∆ A,r ] + [ ∆ zr ] + [ ∆ rz ]  [c A,t − ∆t A ]  [ω A,t − ∆t A ] −     ∆t A     ∆t A    (3.38)  Finally, as before, we can solve this system in two steps, using the solution from the first step to find the final solution in the second step, thus:   θ − ω ′A     [δ A c A ] ⊗ + [ ∆ ] A z ,  ∆t    A      δ Aθ    θ − ω ′A   =  + [ ∆ A,z ] + [ ∆ A,r ] + [ ∆ zr ] + [ ∆ rz ]  [c A,t − ∆t A ]  [ω A,t − ∆t A ] −     ∆t A     ∆t A     θ − ω ′A      ∆t  + [ ∆ A,r ]  [δ A c A ] = [δ A c A ] ⊗  A    (3.39)  (3.40)  Again, it is important to note that the entire system of lines in either direction may be solved simultaneously by ordering the cA values such that either the ∆A,z matrix or the ∆A,r matrix becomes globally tri-diagonal. This is akin to stringing the lines together into one long, continuous line in both directions, and only requires re-ordering between the steps.  In this way, both matrices can be  defined using only a 2D array of 3×(I×J) elements. Generalizing this problem to several solutes is also relatively straightforward. Assuming that the source term for any solute is potentially a function of all solutes, then differentiation of the source term vector gives a Jacobian matrix thus:  Model Development  64  ′ ′ ω12  δω1   ω11  δω   ω ′ ′  2  =  21 ω 22  M   M M    δω M  ω M′ 1 ω M′ 2  where  ′ ≡ ω mn  ∂ω m ∂s m = ∂c n ∂c n  L ω1′N   δc1  L ω 2′ N   δc 2  O M  M    ′  δc N  L ω MN  and thus  δω m ≡ ω m − ω m,t − ∆t A  Applying this to solute m gives the following:  θ   [δ A c m ] + [ ∆ m,z ][c m ] + [ ∆ m,r ][c m ]  ∆t A   δ θ   ′ ][δ A c n ] −   A  + [ ∆ zr ] + [ ∆ rz ]  [c m,t − ∆t A ] = [ω m,t − ∆t A ] + ∑ [ω mn  ∆t  n  A    (3.41)  or ′   θ − ω mm  [δ A c m ] + [ ∆ m,z ][c m ] + [ ∆ m,r ][c m ]  ∆t A   δ θ   ′ ][δ A c n ] −   A  + [ ∆ zr ] + [ ∆ rz ]  [c m,t − ∆t A ] = [ω m,t − ∆t A ] + ∑ [ω mn   n ≠m   ∆t A    (3.42)  For an implicit solution in time which is accurate to O(∆t), this equation may be recast thus:   θ − ω mm  ′    + [ ∆ m,z ] + [ ∆ m,r ] [δ A c m ]    ∆t A   ′ ][δ A c n ] = [ω m,t − ∆t A ] + ∑ [ω mn  (3.43)  n ≠m   δ θ   −   A  + [ ∆ A,z ] + [ ∆ A,r ] + [ ∆ zr ] + [ ∆ rz ] [c m,t − ∆t A ]   ∆t A   which, upon factorization, gives the following two-step solution:  Model Development  65    θ − ω mm  ′    [δ A c m ] ⊗ [ ∆ ] + m z ,    ∆ t A     ′ ][δ A c n ] [ω m,t − ∆t A ] + ∑ [ω mn    n ≠m ′   θ − ω mm =     δ θ    ∆t A  −   A  + [ ∆ A,z ] + [ ∆ A,r ] + [ ∆ zr ] + [ ∆ rz ]  [c m,t − ∆t ] A     ∆t A     (3.44)    θ − ω mm  ′    [ δ A c m ] = [δ A c m ] ⊗ + [ ∆ ] m , r    ∆ t  A    (3.45)  This set of equations may be solved iteratively, in either Gauss-Jacobi or GaussSeidel mode. In Gauss-Jacobi mode, [δcm] values for all solutes are determined assuming previous [δcn≠m] values on the RHS. Next, the system is solved again with the updated values of [δcn≠m] on the RHS, and this is repeated until no further changes in [δcm] values are obtained. Finally, the concentration values themselves are updated. In Gauss-Seidel mode, updated values of [δc1] are immediately applied towards the solution of [δc2], then updated values of both [δc1] and [δc2] are applied towards the solution of [δc3], and so on. Alternatively, the system may simply be solved once without iteration, using either mode. Finally, the off-diagonal (RHS) Jacobian factors may be eliminated from Equations 3.44 and 3.45 with no loss of stability and little loss of accuracy.  3.2.4  Testing the Solute Transport Model  In order to demonstrate the veracity of the solute transport model, transport of a single solute undergoing first-order chemical reaction in a 4 × 1 m column was simulated. Solute transport was superimposed over steady-state saturation and water velocity profiles calculated using the hydrological parameters determined in Chapter 5. For these simulations, the results of which are shown in Figure 3.4, dispersion parameters were varied from zero to 100 times those values determined in Chapter 5.  As expected, the reaction rate is virtually uniform  throughout the column with very high dispersion, while the solute is confined to the advective flow path of the water in the column with zero dispersion.  Model Development  66  1  0.45  0.9  0.4  0.25 0.2 0.15  0.7  Effective saturation, Se  0.3  0.6 0.5 0.4 0.3 0.2  0.05  0.1  21  11  16  0  Radius, r  31  S1  a)  Depth, z  26  31  21  26  Depth, z  0 16  11  6  6  1  1  0.1  Relative solute concentration, c/ c0  0.8  0.35  S1  b)  Radius, r  0.9  0.8  0.8  0.7 0.6 0.5 0.4 0.3 0.2  0.7 0.6 0.5 0.4 0.3 0.2  Relative solute concentration, c/ c0  1  0.9 Relative solute concentration, c/ c0  1  0.1  Figure 3.4  6  d)  21  16  11  Radius, r  0  31  S1  Depth, z  26  31  21  c)  26  Depth, z  0 16  11  6  1  1  0.1  S1  Radius, r  Predicted relative solute concentrations for a first-order chemical reaction (rate constant k = 1/day) at extreme values of dispersion parameters in a 4 × 1 m column: a) steady-state saturation profile used in the simulations (all hydrological parameters as determined in Chapter 5); b) Baseline dispersion parameters (as determined in Chapter 5); c) Dispersion parameters set to zero; d) Dispersion parameters set to 100 times the baseline.  Finally, simulation with no chemical reaction confirms that the solute transport model, like the water transport model, is perfectly conservative, as expected.  Model Development  67  Chapter 4 EXPERIMENTAL METHODOLOGIES 4.1  Introduction  The objectives of the laboratory test program are to determine the hydrological properties and characteristics that influence solution flow and solute transport during heap leaching, to measure these properties from experimental work and to determine the efficiency of extraction of copper from a typical heap leach ore varying some of the important parameters determined. While the experimental design, set-up and operations might not necessarily duplicate industrial practice, particularly with respect to ore agglomeration and stacking, and heap confinement, the procedures used were meant to obtain and model an ideal heap leach operation as a preliminary step. The guiding principle for the design of the laboratory experimental set-up was to maintain the fundamental challenges associated with modeling the full-scale heap leach facility while simultaneously creating experiments that would be affordable and would be completed in a practical time frame.  In designing the experimental program, four aspects  required some attention. Firstly, a standard ore sample crushed to a consistent particle size range and suitable as a heap leach feed was required to ensure the same physico-chemical environment for all tests. Secondly, careful attention was given to the geometries of the reactors to obtain the correct initial and boundary conditions for the modeling exercise. Thirdly, the tracer used to generate the breakthrough curves should be inert to the ore and its concentration in the bed should be uniquely measurable. Finally, the measurement procedures of the required parameters needed to be established and tested to ensure reproducibility. The criteria and procedures used for each aspect are discussed in the following section.  Experimental Methodologies  68  4.2  Materials, Equipment and Procedures  The criteria used to select the various pieces of equipment and experimental materials, and preparations for the various stages of the experimental work are discussed in the following subsections.  4.2.1  Ore  It was necessary to use real ore crushed to specifications of a heap leach operation, rather than, say, crushed silica or glass beads, because the physical properties (namely, the texture, structure, porosity and particle size distribution) and the hydraulic properties (namely, permeability and matric suction) of a real ore and silica or glass beads are very different. Furthermore, it was necessary because the results of this work would be applied to model normal heap leaching practice. Mantoverde Mine in Northern Chile provided about 12 tonnes of lowgrade copper oxide ore for this work. The ore was crushed for their heap leach operations. Its leaching behavior is well known and therefore it is an ideal ore for model validation purposes. Thus, using it for every stage of the experimental work makes the modeling work all the more credible. The copper in the ore occurs mainly as oxidized copper minerals, including brochantite, malachite and chrysocolla.  Small quantities of sulfide minerals are also present, including  chalcocite and, at depth, chalcopyrite. The principal gangue minerals are quartz, orthoclase, sericite, chlorite and specularite, while minor to trace quantities of limonite and calcite have been found in fractures and cavities [Zarate and Kelly, 1995; Zarate et al., 2003]. On receiving the sample, its particle distribution and chemical analysis were undertaken. The chemical composition and the particle size distribution curve of the as-received ore are given in Table 4.1 and Figure 4.1, respectively. Thereafter, the sample was agglomerated in 20 kg batches with water to 6.5% moisture content in a cement mixer and stored in plastic bags.  Experimental Methodologies  69  After loading the ore into the columns, prior to actual tracer concentration measurements (i.e. tests requiring measurement of bulk electrical conductivity of the material with TDR probes and/or four-electrode sensors), indigenous solutes were washed (under unsaturated conditions) out of the ore by irrigating with water.  Their presence would have interfered with the electrical conductivity  measurements. To further mask the effects of any leached minerals that could affect the bulk electrical conductivity of the material during the tracer tests, the background bulk electrical conductivity was elevated with a solution of 0.1 g/L KCl. Chemical composition of the ore (major metallic elements only)  Table 4.1  Element  Cu  Al  Fe  Ca  Mg  Mn  K  Na  Mass % 0.64 6.88 11.63 1.34 2.96 0.12 5.19 1.89  Cumulative Percent passing (%)  100  80  60  40  20  0 0.01  0.1  1  10  100  Particle Size (mm)  Figure 4.1  Particle size distribution curve of the as-received ore.  Experimental Methodologies  70  4.2.2  Reactor Geometries  This work is about developing a 3D axisymmetric model, so all the tracer tests were undertaken in cylindrical columns.  Though commercial heaps are not  confined, laboratory heap leach test work is normally conducted in cylindrical columns.  Two cylindrical columns, namely, the Small Column and the Large  Column, were used for all the tracer tests. The Small Column had an internal diameter of 100 mm and a height of 1200 mm. It was designed to be a 1D reactor to minimize the effects of heterogeneity of the ore, and contained four sets of four-electrode sensors and mini-tensiometers. The Large Column had an internal diameter of 1 m and a height of 4 m. It was designed to investigate the scale (vertical and radial) dependence of the hydraulic parameters from the tracer tests. There were four levels of five TDR probes each (making 20 probes total) in the Large Column to measure the resident bulk electrical conductivity and volumetric water content in the vertical and radial directions. Thus, for a given tracer test, a total of twenty breakthrough curves could be generated for the various positions, and each could be inverted to determine the scale dependence of the particular hydraulic parameters.  For those hydraulic  parameters that would be found to be scale dependent, any scale beyond the height and radius in this column would be estimated. The ore in the Large Column was also acid leached after washing out residual KCl from the tracer tests. This was to simulate leaching in 3D axisymmetric coordinates.  4.2.3  Tracer  To obtain breakthrough curves simulating only the solution flow and solute transport, an inert tracer was required.  Alkaline chlorides have been used  successfully by a number of investigators as inert tracers [Decker and Tyler, 1999; Inoue et al., 2000] to obtain breakthrough curves from local (resident) and global (flowing) effluent tests in soils and ores. Analytical grade KCl supplied by Fisher Scientific was used as the tracer over the concentration range of 0.1 to 2.5  Experimental Methodologies  71  g/L. TDR has an upper limit of about 10 dS/m for obtaining attenuation of the signal for the calculation of the concentration of the pore solution [Nichol, 2002; Or et al., 2003]. The electrical conductivity of a solution of 2.5 g/L KCl (the highest concentration used in this work) is 4.5 dS/m, well below the 10-dS/m upper limit of the TDR probes. The lower concentration value of 0.1 g/L KCl (equivalent to 0.2 dS/m) was set as mentioned above to elevate the background bulk electrical conductivity to mask the effects of leaching during the tracer tests. Running this background solution through the column also minimized any possible adsorption of the tracer by the ore, when solutions of high tracer concentration were used [Muñoz et al., 1997]. De-ionized water was used to prepare all tracer solutions except those used in the Large Column, where tap water was used due to the enormous volumes required.  The solutions for the tracer tests were prepared with potassium  chloride.  4.2.4  Measuring Devices and Data Acquisition Systems  As set out in the study objectives, it was necessary to measure the resident or local tracer concentration distribution in the vertical and radial directions. This cannot be obtained from measurements of tracer concentration in the column effluent alone. The non-destructive methods of measuring porous media salinity include buried porous electrical conductivity sensors, four-electrode sensors, time-domain reflectometry (TDR) probes and electromagnetic induction systems [Inoue et al., 2000].  These methods measure the bulk porous medium solute concentration  rather than the solution concentration of individual ions.  TDR probes  [Kachanoski et al., 1992; Decker and Tyler, 1999; Inoue et al., 2000; Nichol, 2002] and four-electrode sensors [Inoue et al., 2000] have been successfully used to measure resident tracer concentrations and water content in soils and ores. Due to size and installation limitations, the TDR probes were used only in the Large Column to measure the resident tracer concentrations, while four-  Experimental Methodologies  72  electrode sensors were used in the Small Column for the same purpose. Flowthrough conductivity probes measured the concentration of the tracer in the effluent from both columns. Mini-tensiometers attached to transducers were used to measure the pore water pressure changes during infiltration in the Small Column. Changes in the mass of ore due to changes in irrigation flux in the Small Column were measured by an Ohaus electronic balance. All the measuring devices were connected to some form of data acquisition system to allow continuous and automatic data reading and storage during the experiments. The specific system used for each device is discussed in the following sections. 4.2.4.1  TDR (Time-Domain Reflectometry) Probes  A detailed review of the background information for TDR is given in Appendix B. Only its application is discussed here. The material presented in Appendix B may be considered elementary in hydrology and soil sciences. It is presented, however, because hydrometallurgists are typically not familiar with TDR technology. The TDR probes used in this work were supplied by Campbell Scientific (Edmonton, AB).  Each probe consisted of three 316 stainless steel rods of  length 300 mm, with outside diameters of 4.76 mm and spacing between the rods of 22 mm.  The probes were connected via coaxial cables to three  multiplexers (model MXDS 50, Campbell Scientific).  The second and third  multiplexers were connected to the first multiplexer while the latter was connected via a 0.5-m long, 50 Ω coaxial cable to a cable tester (model 1502B, Tektronix, Beaverton, OR). WinTDR 2003 software by Utah State University [Or et al., 2003] controlled the settings of the TDR probes, recorded and analyzed the waveforms based on curve-fitting algorithms described by Or et al. [2003]. The TDR measurements of water content and the apparent medium electrical conductivity were made every 10 minutes and saved directly to a computer connected to the cable tester.  Experimental Methodologies  73  4.2.4.2  Four-Electrode Sensors  The four-electrode sensors used in this investigation consisted of four 316stainless rods of 3.2 mm outside diameter and 80 mm length. They were set in cast epoxy at both ends to ensure that they remained parallel during use in the coarse ore bed. The two inner and two outer rods were spaced at 16 and 32 mm, respectively.  Four four-electrode sensors were fabricated and used in the  Small Column. One major disadvantage of four-electrode sensors is that soilspecific calibration is required. When a current flows through the four-electrode sensor which is in contact with a porous medium, the ratio, I/V, of the electric current, I, flowing through the outer rods to the voltage difference, V, between the inner rods, is measured. The ratio, I/V, is inversely proportional to the electrical resistance of the medium or directly proportional to its electrical conductivity (EC).  The magnitude of the electric  current through the two outer rods was obtained from I = V1/Rf, where Rf is a known resistance inserted in the circuit, as shown in Figure 4.2. The ratio V1/V2 was automatically measured by a datalogger (model CR 1000, Campbell Scientific), which also supplied the required AC at 60 Hz. The proportionality constant between the output value and the bulk electrical conductivity (ECb) depends on the shape and length of the rods, inter-rod spacing, and construction of the sensor. It was determined for all four sensors from calibration experiments by measuring known EC values of solutions of potassium chloride at a known reference temperature. The four sensors were then calibrated for the heap leach material used in this work within the Small Column. The four four-electrode sensors were connected via a relay-multiplexer (model AM16/32, Campbell Scientific) to the datalogger.  LoggerNet v3.1.2 software  (Campbell Scientific) coupled with a custom-written program (in CR-Basic) controlled the settings of the datalogger for data generation and collection from the four four-electrode sensors, the four mini-tensiometer transducers, and the two flow-through conductivity probes at predetermined time intervals.  Experimental Methodologies  74  V1  CR1000 datalogger  Rf  Outer electrodes  V2  Inner electrodes  Figure 4.2 4.2.4.3  The four-electrode sensor circuit  Flow-through Conductivity Probes  Flow-through electrical conductivity probes (model 545 Multi-cell (Pt), Amber Science Inc; OR) were coupled to EC meters (model 19101-00, Cole-Parmer Instruments) which were calibrated before use with a standard solution of 692 ppm NaCl of conductivity 1.413 dS/m (Thermo Electron Corporation, MA) and connected to the datalogger for continuous data recording. 4.2.4.4  Tensiometers  A photograph of a mini-tensiometer with its transducer (model SKT850T, Earth Systems Solutions LLC, CA) is shown in Figure 4.3. The tensiometer comprises four main parts, namely, a ceramic tip with an air-entry pressure of 100 kPa, flexible clear tubing, pressure transducer and jet-filled cup (reservoir).  The  pressure transducer, which can measure capillary pressures between 0 and −100 kPa, was attached to the datalogger for automatic recording to the data acquisition system. The tensiometer was filled with de-aerated distilled water by sucking it through the ceramic tip using a vacuum pump. This ensured that the ceramic cup was fully saturated and that all air bubbles were removed. When good contact with the porous unsaturated medium was achieved, the negative pore-water pressure exerted tension on the ceramic tip and this tension was transferred through the water in the tensiometer shaft to the pressure transducer, which converted the tension to an electrical signal. This electrical signal was  Experimental Methodologies  75  read and stored by the datalogger. Good contact between the tensiometer tip and the porous medium is required for correct functioning of the instrument.  Figure 4.3 4.2.4.5  A photograph of a mini-tensiometer with transducer  Electronic Balance  An electronic balance (Model Ohaus 31, supplied by Ohaus Corporation, NJ) of sensitivity 0.002 kg was used to monitor the changes in mass of water retained by the ore in the Small Column during irrigation. A custom-written Excel macro program read and stored the variations in the mass of the column (recorded by the balance) directly to a computer.  4.3  Test Procedures  In order to achieve the study objectives, the experimental program was subdivided into two main groups, namely, the Tracer Tests and the Heap Leach Performance Tests. The details of the various experiments are presented in the following sections.  Experimental Methodologies  76  4.3.1  Tracer Tests  The Tracer Tests encompassed the displacement of potassium chloride as an inert miscible tracer through samples of agglomerates of the low-grade copper ore in the Small and Large Columns. Sensors (TDR and/or four-electrode) to measure water content and tracer concentration, and tensiometers to measure matric potential, were placed during column filling. While solute mobility within the actual heap leach would be influenced by physico-chemical factors (including dissolution, precipitation, adsorption and desorption) and physical processes (comprising advection and hydrodynamic dispersion), the tracer tests would only address the physical processes. The emphasis in this set of experiments was on flow of solution and transport of potassium chloride as the inert tracer.  4.3.2  Small Column Apparatus  Three sets of experiments were undertaken in the Small Column. The specific objectives of these experiments were: To determine moisture retention characteristics of the ore using a succession of steady state infiltration experiments. To obtain breakthrough curves from tracer tests for the estimation of 1D hydraulic and solute transport parameters. To generate a final drain-down curve to validate the hydraulic parameters. The experimental setup comprised an acrylic column with an internal diameter of 100 mm and a height of 1200 mm, filled with agglomerated low-grade copper ore and placed on a balance to monitor the mass changes due to changes in water content of the ore bed. The first 100 mm of the column was filled with coarse gravel to serve as a drainage layer. The ore was packed in the column above the drainage layer to a height of 1100 mm. The four four-electrode sensors and the four tensiometers were installed in pairs at heights of 200, 400, 600 and 800 mm above the drainage layer as the column was being filled with the ore. The  Experimental Methodologies  77  column was set in a metal frame for support and then loaded onto the balance. The column was closed at its upper end with a tightly fitted PVC cover. The cover had two ports, one at its centre for irrigation of the ore in the column, and the other for ‘breathing’. A photograph and a schematic of the Small Column in operation are shown in Figures 4.4 and 4.5, respectively.  Figure 4.4  A photograph of the Small Column  The pump was calibrated to deliver solution to the surface of the ore at certain specific fluxes.  Solution concentrations of 0.1, 1.0, 2.0, 2.5 g/L KCl were  prepared and used.  During the irrigation, variations in water content, matric  potential, bulk electrical conductivity and electrical conductivity of the effluent with time were measured using the balance, the four tensiometers, the four fourelectrode sensors and the flow-through conductivity probe, respectively. Data were recorded every 10 minutes by the data acquisition system. It was assumed that water retained on the gravel in the drainage layer, the sensors, and any part  Experimental Methodologies  78  of the column apparatus itself was negligible compared to that retained by the ore.  Pump  C0  4-electrode sensors  C1 Feed Solution  Column  Tensiometers and pressure transducers Balance  Drainage layer Datalogger (CR1000)  Flow-through Conductivity probe  EC Meter  Laptop 2  Laptop 1 Multiplexer AM16/32  Collection bucket Figure 4.5 4.3.2.1  Schematic diagram of the Small Column  Determination of Hydraulic Functions  The purpose of the first set of experiments in the Small Column was to determine the moisture retention characteristics of the heap leach ore.  This was required  for the water and solute models. It involved simultaneous measurement of matric pressure and mass of water retained during infiltration of 0.1 g/L KCl solution.  Experimental Methodologies  79  A solution of 0.1 g/L KCl was used to irrigate the ore bed in the Small Column at a flux of 30 L/m2/h (a relatively higher flux than those common in heap leaching practice). After establishing steady state drainage (by constant mass reading on the balance), the irrigation rate was decreased in steps to 15, 10, 7 and 3 kg/m2/h, respectively, attaining steady state drainage at each flux, to obtain drying retention data. The irrigation flux was then increased in steps to 7, 10, 15 and 30 L/m2/h, respectively, again attaining steady state drainage at each flux, to obtain wetting retention data.  The tensiometers and the balance measured  variations in matric potential and water content, respectively, with time. Only the measurements taken during steady state drainage were used to determine the ore-moisture retention characteristic function. 4.3.2.2  1D Tracer Tests  The purpose of the second set of experiments undertaken in the Small Column was to characterize the material hydraulically by steady state tests from tracer experiments. The procedure in this set of experiments involved establishing a steady-state flow at fluxes of 7, 10 and 15 L/m2/h over the ore bed in the Small Column with 0.1 g/L KCl solution. Thereafter the solution concentration was changed to 2.5 g/L, maintaining the same flux. The local bulk EC of the ore column was measured at four depths by the four- electrode sensors, while the global EC of the effluent was measured in the outflow tube at the bottom of the column. Breakthrough curves were generated from these data. 4.3.2.3  Column Drain-Down and Decommissioning  At the end of all the tests, the column was allowed to drain down completely, and the mass was measured continuously to give a final drain-down curve. Finally, the column was emptied and rinsed out. The wet ore and the slurry from rinsing were dried to determine the exact mass of dry ore in the column. The column and all the ancillary equipment and wet gravel of the drainage layer were  Experimental Methodologies  80  remounted on the balance to determine the mass of the apparatus without the ore.  4.3.3  The Large Column  The Large Column was 4 m in height and comprised four 1 × 1 m cylindrical sections plus a bottom section which served as the stand. All parts were made of HDPE. The individual sections were bolted together using flange fittings with rubber gaskets between the flanges to seal the joints. The space above the stand section was filled with two layers of coarse silica gravel (supplied by Target Products Ltd, Vancouver, BC), each of a height of 10 mm, but of different particle size distributions. The particle size distributions were 80% passing 25 mm and 12.5 mm for the bottom and the top layers, respectively. The silica gravel layer constituted a drainage layer for the column material during irrigation.  The  bottom section also had a perforated drain and spout for the effluent to exit the column. As each section was installed, it was filled with the ore to half its height (0.5 m). Then, five TDR probes were laid horizontally in a rectangular pattern on top of the ore bed, one across the center, and one to the east, one to the west, one to the north and one to the south at 0.2, 0.2, 0.3 and 0.3 m from the center, respectively, as shown Figures 4.6 and 4.7. Finally, it was filled to the top and the next section was installed. The top (fourth) section was covered with a flat round cover of the same material. It had nine ports connected by tubing to the solution delivery pumps. Eight of the nine ports were distributed every 45° at half the internal radius of the column and the ninth port was at the centre. The column effluent exited via the spout on the bottom section, through a flow-through conductivity probe and into a large plastic container.  Experimental Methodologies  81  Figure 4.6 4.3.3.1  A photograph of the Large Column, with TDR probe layout  Experimental Procedure  For all tracer tests, there were two tracer solutions of concentrations 0.1 and 2.5 g/L KCl, giving electrical conductivities of 0.2 and 4.5 dS/m, respectively.  The  ore bed was irrigated with solution of electrical conductivity 0.2 dS/m to mask background effects due to residual indigenous solutes. The irrigation pumps were calibrated to deliver solution at preset flow rates of 6, 9 and 12 L/m2/h. The tracer test had two profiles, namely, increasing bulk conductivity of the ore bed (salination) using the 4.5 dS/m solution pumped uniformly through all nine ports, and decreasing bulk conductivity of the ore-bed (desalination) using the 0.2 dS/m solution, combining the flows from all nine delivery tubes through the central port. The TDR probes measured the bulk electrical conductivity and volumetric water content, while the flow-through conductivity probe measured the EC of the effluent, at a preset time interval. The data were saved on their respective data acquisition systems.  Experimental Methodologies  82  Multiplexer MXDS 50  C0 Pump C1  A  Feed Solution  A  TDR cable Tester Tektronix 1502B  TDR probes  Computer  Column Drainage layer Datalogger (CR1000)  Section A-A (probe plan) 3N 3C  EC Meter  3W  3E 3S  Laptop 1 Multiplexer AM16/32  Flow-through Conductivity probe  Collection bucket Figure 4.7  Schematic diagram of the Large Column and instrumentation, with TDR probe layout  4.4  Performance Tests  The Heap Leach Performance Tests comprised leaching copper from low-grade oxidized copper ore with dilute sulfuric acid solution.  These tests were  undertaken in three mini-columns (described below) and in the Large Column following completion of all tracer tests. All leach solutions were prepared using analytical grade sulfuric acid. The details of the various experimental apparatus and procedures are given in the following sections.  Experimental Methodologies  83  4.4.1  Mini-Columns  The leach tests were undertaken in three C-PVC cylindrical columns, each of internal diameter 96 mm and height 500 mm. In each column, a perforated plate covered with a thin layer of glass wool to prevent loss of fine solids, rested on top of a spacer, which rested atop a bottom plate used to seal the column. Thus, the ore was loaded onto the glass wool in the column.  A 100 mm layer of  polyethylene balls (5 mm in diameter) and 5 mm layer of sponge covered the surface of the ore bed to help distribute the irrigating solution uniformly across the surface area of the column. The effective height of the bed of ore in all three columns was 270 mm.  The three columns were individually irrigated by  peristaltic pumps at the same flux of 10 L/m2/hr with sulphuric acid solutions of concentration 5, 10 and 20 g/L, respectively. The effluent from each column drained into a receiving container from which it was harvested periodically, weighed and sampled for analysis of pH and concentrations of free acid and copper. The copper concentration was determined using ICP while free acid was determined after complexing hydrolyzable metals with sodium oxalate.  At the  end of the test, the irrigation pumps were shut down and the columns were allowed to drain completely. The residue was removed, dried and sampled for analysis for residual copper.  4.4.2  Acid Leaching in the Large Column  The purpose of acid leaching the sample in the Large Column after the tracer tests were completed was to determine the extraction of copper from the sample (simulating 3D axisymmetric geometry) under a point source of irrigation. 4.4.2.1  Experimental Procedure  The sample in the column was rinsed with tap water for one week and allowed to drain to eliminate KCl from the tracer tests. Irrigation of water was resumed at 10 L/m2/h until the influent and effluent flow rates were equivalent. Thereafter, 10 g/L H2SO4 solution was pumped at 10 L/m2/h through the central discharge tube  Experimental Methodologies  84  to irrigate the surface of the sample, thus simulating a single drip emitter. The effluent was harvested every 24 hours. Solution samples were taken for the analysis of pH, free acid and ICP multi-element analyses. 4.4.2.2  Kinetics of Copper Leaching  Modeling leaching data from the Column leach experiments requires kinetic data of copper leaching from the ore. Only the oxidized copper fraction in the ore was expected to leach. As a result another test was undertaken to determine the leaching rate of copper from the ore, as described below. 4.4.2.3  Experimental Procedure  Three bottle-roll tests were undertaken with H2SO4 solutions at concentrations of 5, 10 and 20 g/L. Each bottle containing 400 g (dry mass) of agglomerated ore and 1900 kg of sulfuric acid solution was rolled for 7 days. Solution samples were taken at specific times. The solids in the samples were returned to the bottles after centrifuging and removing the supernatant. The solution samples were analyzed for pH, free acid and copper concentrations. At the end of the experiments, the slurries were filtered, and the residues were washed, dried, weighed, and analyzed for mass balance purposes. The head assay of the ore for acid-soluble and total copper was performed. The acid-soluble fraction was determined by leaching 200 g of the agglomerated ore in 1200 mL of 20 g/L H2SO4 for a period of five days in a bottle roll. The slurry was filtered and the residue was washed and dried. Both the original filtrate and wash solution were weighed, sampled and analyzed for copper.  The dried  residue was also weighed and analyzed for copper for mass balance purposes.  Experimental Methodologies  85  Chapter 5 RESULTS AND DISCUSSION – TRACER TESTS The theoretical formulation of the model for water flow and solute transport in heaps and the procedure for its numerical solution were presented in Chapter 3. Experiments designed to determine the parameters of the model were described in Chapter 4. The analyses, interpretation and discussion of the results of the experimental work, and the predictions of the hydraulic and transport parameters of the inert solute transport model, are presented in this chapter.  5.1  The Small Column  As mentioned in Chapter 4, three types of experiments were performed on lowgrade oxide copper ore from Mantoverde in the Small Column, namely: Steady-state infiltration tests for the determination of hydraulic functions Steady state tracer tests for the determination of solute transport parameters Drain-down  5.1.1  Determination of Hydraulic Parameters  The steady-state infiltration tests involved feeding a constant flow of solution to the Small Column for a period of about two days (until the mass of the column reached a constant value), measuring the capillary pressure using tensiometers installed at various levels (as mentioned in Section 4.3.1), and then repeating the procedure at a different solution flow rate. The entire test involved stepping both downward and upward through five flow rates as follows: 30, 15, 10, 7, 3, 7, 10, 15, and 30 L/m2/h. By noting the mass of the entire column apparatus on a digital balance, and subtracting the known tare weight and the mass of dry ore, the total amount of moisture in the column could be measured directly as a function of flow rate. The data for the drying and wetting cycles are shown in Tables 5.1 and 5.2, respectively. As expected, the moisture levels increased  Results and Discussion – Tracer Tests  86  with increasing flow rate. It was assumed that the densities of all solutions used (thus solutions of KCl (between 0.1 and 2.5 g/L KCl) and sulphuric acid (between 5 and 20 g/L H2SO4) used in the various tests) in this study were equal to 1000 kg/m3. Table 5.1  Moisture retention with irrigation flux (drying cycle)  Flux (L/m2/h) Mass water content (%) Volumetric water content (%)  30  15  10  7  3  12.38 19.58  11.92 18.86  11.73 18.56  11.56 18.29  11.24 17.79  The bulk density of the ore was measured as 1582 kg/m3. This value was used to convert all values of water content from mass to volumetric units. Table 5.2  Moisture retention with irrigation flux (wetting cycle)  Flux (L/m2/h) Mass water content (%) Volumetric water content (%)  3  7  10  15  30  11.24 17.79  11.62 18.38  11.76 18.60  11.91 18.83  12.40 19.61  Figure 5.1 shows the final drain-down curve for the Small Column. From this curve, the residual volumetric water content θr was estimated to be 9.2% by mass, which is equivalent to 14.55% by volume.  Also, the saturation water  content θs was estimated to be 24% by mass, which is equivalent to 37.97% by volume, based on mass measurement of a saturated column of ore. These values have been used in all of the calculations of effective saturation Se reported in this thesis.  Results and Discussion – Tracer Tests  87  0.185  Average Volumetric Water Content 3 3 (m water/m )  0.180 0.175 0.170 0.165 0.160 0.155 0.150 0.145 0.140 0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  Time (days)  Figure 5.1  Final drain-down curve for the Small Column  Plotting the 1D applied flux at steady state against the effective saturation values, Se, gives the hydraulic permeability relationship for the ore. The volumetric water content values in Tables 5.1 and 5.2 may be converted to effective saturation by the following formula: Se =  θ − θr θs − θr  (5.1)  The parameters in Equation 5.1 maintain their definitions from Section 3.2.1. As shown in Figures 5.2 and 5.3, these functions are virtually identical for both the drying (decreasing flow rate) and wetting (increasing flow rate) cycles, suggesting very little hysteresis in the hydraulic behaviour of this ore in the moisture range investigated. The data in both cases are well fitted by a power law function, with coefficients of determination (R2) of about 0.99, suggesting that the Brooks-Corey model is adequate to describe the hydraulic characteristics of  Results and Discussion – Tracer Tests  88  the Mantoverde ore, as would be expected at the low saturation levels typical of beds of coarse material. The Brooks-Corey model is defined as follows: v w = K w Seψ  (5.2)  hc = hc ,0 Se−1/ λ  (5.3)  where vw is the local water volume flux (m/s), Kw is the saturated hydraulic conductivity (m/s), hc is the local capillary head (m), hc,0 is the “air-entry head” and λ and ψ are material-specific constants. The air-entry head has the biggest impact in determining the degree by which water will spread out under a point source such as a drip emitter in a heap. Very low values of hc,0 associated with coarsely crushed solids ensure very little capillary head, and therefore, very little lateral spreading of water in heaps. 35  5.2315  Solution Flux = 99073Se 2  R = 0.9922  25  2  Solution flux (L/m /h)  30  20  15  10  5  0 12%  14%  16%  18%  20%  22%  24%  Effective saturation, Se (–)  Figure 5.2  Moisture retention curves with best power-law fits (drying)  Results and Discussion – Tracer Tests  89  35 Solution Flux = 101530Se5.2632 R2 = 0.9888  25  2  Solution flux (L/m /h)  30  20  15  10  5  0 12%  14%  16%  18%  20%  22%  24%  Effective saturation, Se (–)  Figure 5.3  Moisture retention curves with best power-law fits (wetting)  The capillary pressure data, samples of which are shown in Figure 5.4, while noisy, appear to be sufficient for determining the relationship between capillary head and effective saturation. Air entry pressure head hc,0 is usually obtained from soil water retention curve measurement undertaken in a Tempe cell as mentioned in Chapter 2 (Section 2.2.1.1). Several attempts to obtain the soil water characteristic curve of this heap leach sample failed due to bubbles forming continuously on the porous plate. The bubbles resulted from the wide particle size distribution of the ore sample. Tempe cells are good for uniformly fine materials; for example, the fine soils or clays used for soil liners.  Other methods like neutron or gamma  attenuation techniques were not tried because of safety concerns and licensing requirements. The noise might be due to the effects of ore stratification, structural or textural variability, and/or non-uniform distribution of water in the column of ore, any of which could lead to medium discontinuity in certain parts of the ore bed. These  Results and Discussion – Tracer Tests  90  effects could have been exacerbated by the fact that the porous tips of the tensiometers were about the same diameter as the largest particles in the bed. There may also be some hysteresis due to diurnal temperature variations detected in the trends of the raw data. A logical pattern is, however, hard to discern.  Nonetheless, power-law fits appear to be sufficient, thus further  reinforcing the fact that the Brooks-Corey model is sufficient to describe the hydraulic behavior of the Mantoverde ore. 1.70 Capillary pressure = 0.8242Se-0.3294 R2 = 0.5788  Capillary pressure (kPa)  1.60  1.50  1.40  1.30 Drying Wetting 1.20 12%  14%  16%  18%  20%  22%  24%  Effective saturation, Se (–)  Figure 5.4  Capillary head vs effective saturation with best power-law fits a) 800 mm depth  Results and Discussion – Tracer Tests  91  1.90 -0.2181  Capillary pressure = 1.1564Se 2  Capillary pressure (kPa)  R = 0.8411 1.80  1.70  1.60  Drying Wetting 1.50 12%  14%  16%  18%  20%  22%  24%  Effective saturation, Se (–)  Figure 5.4  b) 600 mm depth  1.50 Capilary pressure = 0.5531Se-0.4808 R2 = 0.6685  Capillary pressure (kPa)  1.40  1.30  1.20  1.10 Drying Wetting 1.00 12%  14%  16%  18%  20%  22%  24%  Effective saturation, Se (–)  Figure 5.4  c) 400 mm depth  Results and Discussion – Tracer Tests  92  1.80 -0.302  Capillary pressure = 0.9111Se 2  R = 0.3989  Capillary pressure (kPa)  1.70  1.60  1.50  1.40  1.30 Drying Wetting 1.20 12%  14%  16%  18%  20%  22%  24%  Effective saturation, Se (–)  Figure 5.4  d) 200 mm depth  The parameters of the well-known Brooks-Corey and van Genuchten-Mualem models as estimated from the degree of saturation and capillary pressure data are given in Tables 5.3 and 5.4. Average power-law fits are also plotted based on log-mean averages of both the capillary head fits and their derivatives as shown Figures 5.5 and 5.6, respectively. Either method gives nearly the same results, although the average based on the derivatives will be used for further calculations since it is only the derivative which is involved in the water flow model.  Results and Discussion – Tracer Tests  93  Table 5.3  Hydraulic conductivity parameters derived from the Small Column Drying Wetting log mean  Test 2  vw coefficient (L/m /h) vw exponent  99073 5.23  101530 5.26  100294 5.25  Brooks-Corey (BC) parameters Kw (m/s) ψ  0.0275 5.23  0.0282 5.26  0.0279 5.25  van Genuchten-Mualem (VGM) parameters Kw (m/s) 0.157 0.154 0.160 m 0.421 0.423 0.420  Table 5.4  Capillary pressure parameters derived from the Small Column  Probe depth pc coefficient (kPa) pc exponent hc,0 (m) BC exponent, λ VGM exponent, n  200 mm 400 mm 600 mm 800 mm 0.91 –0.30 0.093 3.31 7.86  0.55 –0.48 0.056 2.08 4.94  1.16 –0.22 0.118 4.59 10.88  log mean log mean (hc) (dhc/dSe)  0.82 –0.3294 0.084 3.03 7.21  0.83 –0.33 0.085 3.01 7.14  0.80 –0.33 0.082 3.01 7.14  Brooks-Corey (BC) parameters are given directly from the power-law fits of the permeability-saturation and capillary pressure data. The relationship between these parameters and the van Genuchten-Mualem (VGM) parameters actually used in the model are given as follows: ψ 1 m= −   2 4  −1  n=  λ m  K w( VGM ) =  Kw(BC ) m2  (5.4)  It is important to note that the effective values of hydraulic conductivity Kw are different for the two models, suggesting that Kw should be viewed as more of a fitting parameter than a physical property of the bed per se.  Results and Discussion – Tracer Tests  94  1.9 1.8  Capillary pressure (kPa)  1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 12%  200 mm 400 mm 600 mm 800 mm mean p mean dp  14%  16%  18%  20%  22%  24%  Effective saturation, Se (–)  Figure 5.5  Capillary pressure power-law fits based on average hydraulic parameters from the Small Column: capillary pressure  The trends obtained in Figure 5.5 are similar to those obtained by Inoue et al. [2002] for volumetric water contents in soils from 10 and 20%. This is expected as the lower the water content of a porous medium, the higher its suction (matric) potential, as discussed in Section 2.3.1.  Results and Discussion – Tracer Tests  95  -1.0 -1.5 -2.0  dpc/dSe (kPa)  -2.5 -3.0 -3.5 200 mm 400 mm 600 mm 800 mm mean p mean dp  -4.0 -4.5 -5.0 -5.5 -6.0 12%  14%  16%  18%  20%  22%  24%  Effective saturation, Se (–)  Figure 5.6  Capillary pressure power-law fits based on average hydraulic parameters from the Small Column: derivative of capillary pressure  5.1.2  Estimation of Transport Parameters  All the basic hydraulic parameters have been determined in the previous sections (as given in Tables 5.3 and 5.4). The parameters which remain to be determined for predicting solute transport (as discussed in Section 3.2.2) are the dispersivity coefficients, αL and αT, and the diffusivity coefficient, εD. In porous bed flow, solutes are transported by the mechanisms of advection (bulk fluid motion) and dispersion, including mechanical dispersion and molecular diffusion. These driving forces are expressed in Fick’s Law, as given in Equation 3.24, with dispersivity coefficients as defined in Equation 3.25. There is no satisfactory technique to determine dispersivity coefficients directly from experiment, and the common practice is to fit miscible displacement data to analytical or numerical solutions of solute transport models. The most common models used for the inverse estimation of this parameter are the Advection  Results and Discussion – Tracer Tests  96  Diffusion Equation (ADE) and the Dual Porosity Model (DPM). The former is usually used for systems which include only single domain transport media, and the latter is used for systems which include at least two domains for water flow and solute transport [Coats and Smith, 1964; Decker and Tyler, 1999; Forrer et al., 1999 and Nützmann et al., 2002]. From the inert tracer transport tests run in the Small Column (as given in Section 4.3.1) and the 1D version of the solute transport equation (Chapter 3), the 1D transport parameters in the heap material were estimated.  The strategy  employed here consisted of using the breakthrough curves (BTC) of the effluent solution to eliminate the local fluctuations of the measured data and to estimate representative dispersivity parameters of the ore. While solute transport in large diameter columns may require resolution in a 2D configuration, solute transport in small diameter columns may be adequately represented in a 1D configuration. Thus for flow in the vertical direction in the Small Column tests, only the coefficient of longitudinal dispersivity, αL, may be relevant.  Transverse  dispersivity, αT, in such cases, will be considered negligible. Hence from these tests, the longitudinal dispersivity as a function of water content was estimated. The longitudinal dispersivity, αL, for the Small Column data was determined indirectly by fitting observed effluent concentration profiles of the inert solute from solute transport studies, conducted at steady state, run at nominal flow rates of 7, 10, and 15 L/m2/h. These nominal rates were corrected by calculating the area between the tracer curves and the concentration axis, as shown in Figure 5.7. This integral bears the following relationship to the flow rate and the water content: ∞  Zθ  ∫ (1 − C ) dt = v 0  (5.5)  w  where Z is the effective depth of the ore bed in the column, C is the normalized concentration of the solute at time t, θ is the average volumetric water content,  Results and Discussion – Tracer Tests  97  calculated from the mass of water retained (noted on the balance) during the experimental run at a particular flux data, and vw is the actual volume flux.  1.0 Flux=16.2 L/m2/h 2  Normalized Concentration  Flux=7.04 L/m /h 0.8  Flux=10.34 L/m2/h  0.6  0.4  Data Model fit Data Model fit Data Model fit  0.2  0.0 0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  Time (days)  Figure 5.7  Small Column normalized effluent tracer data and model fits  From the curve of effluent concentration vs time, the area between the curve and the concentration axis was calculated using the trapezoid rule.  With known  values of Z and θ, vw was calculated directly from Equation 5.3. On this basis, given the three known values of volumetric water content from the balance, the three nominal flow rates were corrected to 7.04, 10.34, and 16.2 L/m2/h, respectively. Thus, slight drifts of the peristaltic pumps from the set-points for particular fluxes over time could introduce slight variations in the irrigation flux. The 1D version of the solute transport model developed in this investigation (Chapter 3) was run at various values of αL for the corrected irrigation flux to obtain a concentration vs time curve data. The corresponding area between the curve and the concentration axis was calculated with the trapezoid rule, as done for the data. For each run, the sum of square error (SSE) from the actual data was calculated. The SSE versus αL was plotted. Using a third order polynomial  Results and Discussion – Tracer Tests  98  empirical regression model to fit this plot, the value of αL at the minimum SSE was determined for each irrigation flux. The corrected fluxes and optimal values of the longitudinal dispersivity along with the goodness of fit are given in Table 5.5 and shown in Figure 5.7. Table 5.5  5.1.2.1  Variation of dispersivity with applied flux in the Small Column R2  θ (%) vw (L/m2/h)  αL (m)  18.33 18.51 18.90  0.1539 0.9998 0.1934 0.9999 0.2639 0.9997  7.04 10.34 16.20  Water content as function of volumetric flux  It is interesting to note that, over the water content range investigated, the volumetric water content increased with irrigation flux.  Under steady state  conditions of irrigating the heap from a point source in the 1D column, the water content would be uniform throughout the column. Thus, it is assumed that there is little or no variation in water content implying no gradient of water content with depth. With these assumptions, Equations 3.2b and 3.8 reduce to:  v w = K w k rw  (5.6)  Using the best-fit values of VGM parameters m and Kw (from Table 5.3), and θ (from Section 5.1.1), the predicted values of vw were calculated from the water content values as given in Table 5.5 and shown in Figure 5.8. The coefficient of determination (R2) of the model fit to the observed data is 0.99, further affirming that the model represents the data well.  Results and Discussion – Tracer Tests  99  18  2  Volumetric Flux (L/m /h)  16  14  12  10  8  Model Data Linear (Model)  6  4 0.183  0.184  0.185  0.186  0.187  0.188 3  0.189  0.19  3  Volumetric Water Content (m water/m )  Figure 5.8  Linear correlation of volumetric water content with normalizationcorrected volume flux vw  5.1.2.2  Dispersivity as a function of volume flux  Another useful observation, given in Figure 5.9, is a near-perfect linear relationship of dispersivity with the flow rate. Since these tests were 1D, it is presumed that the flow rate corresponds to the resultant velocity. Hence, the following function for αL will be adopted in the 2D axisymmetric solute transport model:  α L = 0.06929 + 0.01201v w (L/m2 /h) = 0.06929 + 43240 v w (m/s )  Results and Discussion – Tracer Tests  (5.7)  100  0.27  0.25  2  Longitudinal dispersivity (m /s)  Long. dispersivity = 0.012(Volume Flux) + 0.0693 2 R =1  0.23  0.21  0.19 data Linear (data) 0.17  0.15 6  8  10  12  14  16  18  2  Volume Flux (L/m /hr)  Figure 5.9  Linear correlation of longitudinal dispersivity with normalizationcorrected volume flux vw  The relationship between water content and dispersivity is the subject of controversy in the literature.  It has been shown by various mathematical  analyses, [Nützmann et al., 2002; Toride et al., 2003; and Javaux and Vanclooster, 2003; inter alia] that dispersivity decreases as volumetric water content increases.  However, the experimental data from these works could  hardly be fitted as such. For example, the data generated by Toride et al. [2003] have a wide water content range (about 0.08 ≤ θ ≤ 0.38), and various groups of datapoints could be fitted differently, especially Figures 7 and 8 in Toride et al., [2003]. Fitting them separately could lead to contradictory conclusions. In the lower water content region (0.08 ≤ θ ≤ 0.18), the dispersivity could be said to increase with increasing water content. In the present work, with the limited water content range relevant to safe and successful operation of a heap leach facility, this relationship suggests that an increase in water content increases the dispersivity (Table 5.5 and Figure 5.9), which corroborates the work of  Results and Discussion – Tracer Tests  101  Vanderborght et al. [1997] who observed an increase in dispersivity with increasing flow rate (flux). In heap leaching (and hence the adopted experimental protocol in this investigation), the bottom boundary condition does not allow high ranges of water content. Javaux and Vanclooster [2003] discussed the two contradicting impacts of flow rate on dispersivity. They attributed increased flow rate with increased dispersivity to pores in the transport domain becoming larger and thus increasing the variation of mixing time of the solute particles. The decrease of dispersivity with increase in flow rate they attributed to reduced tortuosity, making flow paths shorter resulting in narrow arrival time distribution.  Ore bed stratification in the  column and in the full scale heap leach pad could also mask the full effect of water content on dispersivity. Thus quantitative relationship of water content and dispersivity for the range of water contents relevant to heap leaching (incorporating ore-bed structural or textural variability, non-uniform moisture content and scale) certainly warrants further study from both the mathematical and experimental points of view. 5.1.2.3  Diffusivity  Molecular diffusion is often the dominant mechanism of solute transport when hydraulic gradients and/or hydraulic conductivity are low.  Diffusion may be  controlled by manipulating the diffusion coefficient and/or the pathways for diffusion. The diffusion coefficient has been known to depend on water content in unsaturated porous media [Lim et al., 1998]. The diffusivity coefficient, εD, is difficult to measure owing to the fact that molecular diffusion generally plays an insignificant role in 1D experiments, and is not even very important in 2D experiments. Hence, for completeness of the modeling effort, data were fitted from the literature [Lim et al., 1998]. Figure 5.10 shows the results of several studies, and the solid curve, fitted for the purpose of this work, represents the following very simple correlation:  Results and Discussion – Tracer Tests  102  ε D = 3 .3 θ 2  (5.8)  This curve passes through the solid symbols, which represent natural graded materials. The open symbols represent artificial materials such as glass beads, and are therefore not relevant to heap leaching.  In their work, Lim et al. [1998]  found a general increasing trend in the normalized diffusion coefficient with increasing saturation. This trend is expected and may be attributed mainly to decreasing diffusion path length (reduced tortuosity) associated with increasing water content.  Another factor that may contribute to this trend is a decrease in  the viscosity of the liquid next to the ore particles.  It bears noting that this  function has virtually no effect on either the 1D or the 2D simulations undertaken with the measured hydraulic values of the Mantoverde copper oxide ore. This will be discussed further in Chapter 7.  1  Diffusivity coefficient  0.8  0.6  0.4  0.2  0 0  0.1  0.2  0.3 3  0.4  0.5  3  Volumetric water content (m water/m )  Figure 5.10 Diffusivity function vs volumetric water content from several studies  [adapted from Lim et al., 1998]  Results and Discussion – Tracer Tests  103  Finally, the 1D solute transport model results are compared with the Small Column internal tracer responses measured by the four-electrode sensors at depths of 800, 600, 400 and 200 mm as shown in Figure 5.11.  Normalized Concentration  1.0 200 mm data 200 mm model 400 mm data 400 mm model 600 mm data 600 mm model 800 mm data 800 mm model 1000 mm data 1000 mm model  0.8  0.6  0.4  0.2  0.0 0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0  Time (days)  Figure 5.11 Small Column tracer curves at various depths and their  corresponding model fits (at hc,0 = 0.05 m, vw = 16.2 L/m2/h) The effluent concentration data (as measured by the flow-through conductivity probe) were used to estimate the transport parameters and hence the model fits this data best.  The degree of fitness gets worse the farther away the four-  electrode sensor is from the point of effluent data measurement. Obviously, there was a significant degree of diurnal variation in the experimental tracer curves from the four-electrode sensors. BTCs obtained from probes at depth 200 mm below the surface (the topmost) in the Small Column experienced more heterogeneous conditions because of the impact of solution dripping at the centre of the surface of the bed. Also, the tailing effect observed in the Small Column data could result from a fraction of pores with relatively smaller water velocities in comparison with the mean velocity in the dominant water flow  Results and Discussion – Tracer Tests  104  domain.  Such effects may be modeled using a two-domain transport model  such as the dual-porosity transport model [Decker and Tyler, 1999].  These  heterogeneities would affect the responses of the probes to various extents and contribute to deviations between the predicted and the observed concentration data. However, within this limitation, it would appear that the simulations are very close to the experimental results, which lends credence to the parameters determined from the tracer tests.  5.2  Estimation  of  Parameters  for  the  3D  Axisymmetric  Transport Model Concentration profiles obtained from the Large Column effluent tracer concentration data at various applied fluxes were fitted with the 3D axisymmetric solute transport model to estimate the remaining hydraulic parameters. Three nominal volume fluxes were applied to the column: 6, 9, and 12 L/m2/h. The inert tracer was distributed from all nine feed points at the top of the column, and was rinsed from the central feed point only, at the same overall applied flux. Tests irrigated from multiple points are denoted “m” for multiple feed points, and those irrigated by the single point in the centre are denoted “s” for single feed point. All 3D simulations carried out in support of the Large Column experiment employed a 16 × 32 grid of radial and axial control volumes, respectively. Normalized (i.e., 0 to 1) concentration outputs were generated for those nodes which corresponded most closely to the placements of the TDR probes. The probes were placed at depths of 0.5, 1.5, 2.5, and 3.5 m, and at radial distances from the centre, r = 0 (C), 2R/5 (E and W) and 3R/5 (N and S). Hence, the nodes (i,j) as stated in Table 5.6, were queried. In all of the figures that follow, experimental results from the TDR probes and flow-through conductivity probes are compared with model results obtained with the aid of the set of hydraulic and solute transport parameters determined from the Small Column experiments and presented in Tables 5.3 and 5.4.  Results and Discussion – Tracer Tests  105  Grid points queried for comparison of model and TDR results  Table 5.6  r = 0 (C) r = 2R/5 (E,W) r = 3R/5 (N,S)  Depth z = 0.5 m z = 1.5 m z = 2.5 m z = 3.5 m  (4,1) (12,1) (20,1) (28,1)  (4,6) (12,6) (20,6) (28,6)  (4,10) (12,10) (20,10) (28,10)  The normalized concentrations of the tracer as measured at the various levels and positions were plotted alongside the corresponding simulated concentrations obtained from the model run at an irrigation flux of 6 L/m2/h in the Large Column, as shown in Figures 5.12a through 5.12d and summarized in Table 5.7. While in some cases, the TDR probe concentration profiles lie on the simulated data, in other cases, there are deviations. The possible reasons for these deviations are discussed at the end of this section. 1.0 14 13  Normalized Concentration  0.8 12  4S 4N 4W 4E 4C  0.6  11  0.4 10  0.2  9 8 7  6 1  0.0 0  2  4  6  8  10  12  14  16  Time (days)  Figure 5.12 Simulated (at 16 radial nodes) and observed tracer profiles in the  Large Column at an irrigation flux of 6 L/m2/h a) 0.5 m depth (TDR probe level 4)  Results and Discussion – Tracer Tests  106  1.0 16  Normalized Concentration  0.8 15  3W 3S 3N 3C 3E  0.6 14  0.4  13  12  0.2  11 10 9  8 1  0.0 0  2  4  6  8  10  12  14  16  Time (days)  Figure 5.12 b) 1.5 m depth (TDR probe level 3) 1.0  Normalized Concentration  0.8 16  0.6  15  2W 2N 14  0.4 13 12  0.2  11 10 8 9 1  0.0 0  2  4  6  8  10  12  14  16  Time (days)  Figure 5.12 c) 2.5 m depth (TDR probe level 2)  Results and Discussion – Tracer Tests  107  1.0  Normalized Concentration  0.8  0.6 1C 1E 1S 1W  16 15  0.4 14 13 12  0.2  11 8  9  10  1  0.0 0  2  4  6  8  10  12  14  16  Time (days)  Figure 5.12 d) 3.5 m depth (TDR probe level 1)  The increased axial and radial distribution of water (degree of saturation as shown in Figures 5.13a to 5.13c) with increase in irrigation flux implies the dispersivity ratio increases with increase in flux, since the solute is in the liquid phase. The optimal values of αT/ αL when the tracer effluent data from the Large Column at the three fluxes were fitted to the model show increase in dispersivity ratio, αT/αL with increase in flux, as given in Table 5.8.  Results and Discussion – Tracer Tests  108  Table 5.7  Summary of concentration profiles of TDR probes Depth (m)  0.5  1.5  2.5  3.5  Probe ID  Expected Position  Observed Position  4C 4E 4N 4S 4W 3C 3E 3N 3S 3W 2C 2E 2N 2S 2W 1C 1E 1N 1S 1W  c_4(1) c_4(6) c_4(10) c_4(10) c_4(6) c_12(1) c_12(6) c_12(10) c_12(10) c_12(6) c_20(1) c_20(6) c_20(10) c_20(10) c_20(6) c_28(1) c_28(6) c_28(10) c_28(10) c_28(6)  c_4(1) c_4(5) c_4(8) c_4(8) c_4(7) c_12(7) c_12(4) c_12(9) c_12(9) c_12(11) NA NA c_20(10) NA c_20(12) c_28(7) c_28(7) NA c_28(12) c_28(12)  NA: For the TDR probes in these positions, the concentration data were dominated by noise and could therefore not be worked on further. Table 5.8  Optimal αT/αL ratios for various solution fluxes of the tracer data Flux (L/m2/h)  Optimal αT/αL  6 9 12  2 3 4  Results and Discussion – Tracer Tests  109  Column height  0.35-0.40 0.30-0.35 0.25-0.30 0.20-0.25 0.15-0.20 0.10-0.15 0.05-0.10 0.00-0.05  Column diameter  Figure 5.13 Se profile in the Large Column a) vw = 6 L/m2/h  Column height  0.40-0.45 0.35-0.40 0.30-0.35 0.25-0.30 0.20-0.25 0.15-0.20 0.10-0.15 0.05-0.10 0.00-0.05  Column diameter  Figure 5.13 b) vw = 9 L/m2/h  Results and Discussion – Tracer Tests  110  Column height  0.40-0.45 0.35-0.40 0.30-0.35 0.25-0.30 0.20-0.25 0.15-0.20 0.10-0.15 0.05-0.10 0.00-0.05  Column diameter  Figure 5.13 c) vw = 12 L/m2/h  This observation amounts to a linear relationship between flux and αT/αL in the Large Column. Since αL was found to have a linear relationship with volume flux, this would imply that αT has a quadratic relationship with volume flux. Thus it would appear that while αL has to do with momentum, αT has more to do with kinetic energy. In all subsequent simulations, this relationship was accounted for to correspond to the flux of irrigation. The effluent data from the Large Column fitted with the solute model using hc,0 = 0.082 m (the derivative log-mean value from the capillary data, given in Table 5.3) is shown in Figure 5.14. It shows that the tracer solution was reaching the bottom of the column well before expected, and that the column was behaving as if the value of hc,0 was much smaller (and, therefore, the degree of water spreading was much less) than that predicted from the Small Column tests. While this could indicate channeling, the smooth monotonic shapes of the experimental curves would seem to rule out solution bypassing, which typically results in tracer curves with an abrupt change of slope.  Results and Discussion – Tracer Tests  111  1.0 0.9 Data Model  Normalized Concentration  0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 12 L/m2/h  6 L/m2/h  9 L/m2/h  0.0 0  2  4  6  8  10  12  14  16  18  20  Time (days)  Figure 5.14 Model and experimental effluent tracer responses from the Large  Column at fluxes of 6, 9 and 12 L/m2/h assuming hc,0 = 0.082 m Here again, the best-fitting calibrated value of hc,0 was obtained using an iterative process of adjusting hc,0 by trial and error until the simulated dependent variable (effluent solute concentrations) matched the experimental data run at a flux of 6 L/m2/h. The optimal value of hc,0 obtained from this exercise was 0.05 m. It was validated for the effluent data obtained from the runs at 9 and 12 L/m2/h. The data for all three solution fluxes and their corresponding simulations at hc,0 = 0.05 m are in Figure 5.15.  It is clear from this figure that these fits better  represent the observed data than the models run at hc,0 = 0.082 m. Thus, the optimal value of hc,0 was decreased from 0.082 to 0.05 m for use in all of the solute transport model fits except where stated otherwise.  The initial value  (0.082 m) was estimated from Small Column tensiometer (capillary pressure) data which experienced a high degree of scatter as shown Figure 5.4. The scatter in the data probably relates to the coarse nature of the ore, with maximum particle size around 15 mm which is comparable to the 12 mm diameter of the  Results and Discussion – Tracer Tests  112  porous tips of the tensiometers themselves. Furthermore, the very low value of hc,0 is not surprising. This is the result of the coarse particle size of the heap leach material, and simply emphasizes the fact that coarse-sized porous media are generally characterized by low suction potentials. This is further discussed in Chapters 6 and 7. 1.0 0.9 Data Model  Normalized Concentration  0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 12 L/m2/h  6 L/m2/h  9 L/m2/h  0.0 0  2  4  6  8  10  12  14  16  18  20  Time (days)  Figure 5.15 Model and experimental effluent tracer responses from the Large  Column at fluxes of 6, 9 and 12 L/m2/h assuming h0 = 0.05 m Predicted axial velocity (volume flux) profiles using the hydraulic and transport parameters from the Small Column experiments (given in Tables 5.3 and 5.4) are shown for all three tracer experiments in Figure 5.16.  Also, as mentioned  previously, velocities are well distributed evenly from the single point source given the specific set of hydraulic parameters taken (especially hc,0). These plots show that the higher the flux, the wider the volume covered (in both axial and radial directions) by the flow under the emitter.  Results and Discussion – Tracer Tests  113  2.0E-05 1.8E-05 1.6E-05  1.2E-05 1.0E-05 8.0E-06 6.0E-06  z-normal volume flux, vz  1.4E-05  4.0E-06  S31  S29  S27  S25  S23  S21  Radius  S19  S17  S15  S13  S9  S11  S7  S5  S3  25  a)  31 S1  0.0E+00 19  13  7  1  2.0E-06 Depth  2.0E-05 1.8E-05 1.6E-05  1.2E-05 1.0E-05 8.0E-06 6.0E-06  z-normal volume flux, vz  1.4E-05  4.0E-06  S31  S29  S27  S25  S23  S21  Radius  S19  S17  S15  S13  S11  S9  S7  S5  S3  25  b)  31 S1  0.0E+00 19  13  7  1  2.0E-06 Depth  2.0E-05 1.8E-05 1.6E-05  1.2E-05 1.0E-05 8.0E-06 6.0E-06  z-normal volume flux, vz  1.4E-05  4.0E-06  S31  S29  S27  S25  S23  S21  Radius  S19  S17  S15  S13  S11  S9  S7  S5  S3  25  13  19  0.0E+00 31 S1  c)  7  1  2.0E-06 Depth  Figure 5.16 Axial velocity profiles for the single point irrigation in the Large  Column tests predicted using the parameters determined from the Small Column tests a) vw = 6 L/m2/h, b) vw = 9 L/m2/h, c) vw = 12 L/m2/h Results and Discussion – Tracer Tests  114  The key value of the TDR data would seem to be in pointing up the very different tracer responses at different locations within the column, particularly from probes which were placed to give similar responses. Several general observations present themselves based on the TDR data: The magnitude of, and variation in, the degree of wetting with column depth (measured by the TDR probes) is much larger than predicted, and may be a product of inaccurate probe readings.  Topp’s equation (Equation B3 in  Appendix B) used for volumetric water content determination might also contribute to the differences between the predicted values and the data. Topp’s equation was originally developed empirically for relatively fine agricultural soils, and may not bear the same coefficients when applied to heap leach material as shown by Decker and Tyler [1999].  At depths of 0.5  and 3.5 m, the relative wetting values of the various probes are more or less as predicted, while the other two levels (i.e. depths 1.5 and 2.5 m) give more random values. It is clear from Figures 5.14 and 5.15 that, within the range of volume flux investigated in this work, the higher the irrigation flux, the faster the solute exited the column in the effluent. In heap leaching, hydraulic factors, kinetics of leaching and availability of solution will be taken into account to determine an economic irrigation flux. Based on the placement of the probes and the axisymmetric irrigation pattern, it was expected that the E and W probes at each depth would give the same response, as would the N and S probes. Obviously, this was not the case. From the data, it would appear that the flow axis shifted laterally down the column towards E and slightly towards N. However, given the reasonably good correspondence between the measured and predicted tracer profiles with column radius, it would appear that the flow only shifted some 10 cm. Notwithstanding the poor conductivity resolution from some of the probes (and the complete lack of data from the N probe at 3.5 m), it is apparent that the responses of some of the probes at some locations are simulated fairly  Results and Discussion – Tracer Tests  115  well by the model using the parameters from the Small Column. Also, the model appears to have captured at least the broad outlines of the tracer response at the different depth levels. However, in general, it is fair to say that the data generated by the TDR probes in this study were not of a sufficiently high quality for model validation purposes. The relatively large size (length and width) of the TDR probes might also overshadow the required fine grid coordinates, thus yielding average values over a large area spanning several coordinates in the horizontal plane.  Another possible  source of error could be that bed shifting during irrigation may have resulted in at least one of the three rods not remaining parallel to the others. In this case, a smaller size TDR probe could be used. However, there is a limitation on the minimum length of the probes to avoid signal attenuation [Nichol, 2002]. Thus four-electrode sensors may be a better choice for measuring solute concentrations over relatively small areas. Given all this, it must be noted that only the flow-through conductivity probe in the effluent operated reliably enough to give simulation-ready results (as shown in Figure 5.15).  Results and Discussion – Tracer Tests  116  Chapter 6 RESULTS AND DISCUSSION – PERFORMANCE TESTS 6.1  Introduction  A comprehensive heap leach model must account for the two main components of the process, namely, reaction processes at the solution-solid interface and transport of dissolved species in the infiltrating solution within and between the solid particles comprising the heap. The analyses of the data generated from the inert tracer tests conducted on the Mantoverde copper ore held in cylindrical columns of various sizes, and hydrological model fitting of the data, were discussed in Chapter 5. Also, the optimal values of parameters that could not be determined directly by experiments were estimated and discussed in that chapter.  In the present  chapter, the data generated when the same ore was leached with dilute sulphuric acid in what are termed ‘Performance Tests’ are analyzed and discussed, and fitting of the data by the reactive solute transport model is undertaken. The solutes in this case were H+ from the sulphuric acid (in the infiltrating solution) and Cu2+ from the reaction of copper minerals in the ore with the acid solution. The reactive solute transport model is essentially the same as the inert solute transport model but with terms accounting for chemical sources (whether injected or generated as a result of the reaction) and rates of reaction (rates of consumption of H+ and production of Cu2+). The data from the Performance Tests were used to determine the reaction parameters.  6.1.1  The Concept of Leaching Large Particles  The Performance tests comprised irrigating columns of the ore continuously from a point source (located at the centre), with dilute sulfuric acid solution of known concentration for a specified period of time. The acid reacted with the copper minerals in the ore particles and the metal ions (copper) were transported by the  Results and Discussion – Performance Tests  117  solution through the bed of ore to the outlet of the column. During the reaction, the acid protons reacted with the copper minerals according to Reactions R2.3 to R2.9 (in Section 2.2.2). Conceptually, the protons react initially with the copper minerals located at the particle surface and thereafter progress inward towards the core of the particle, leaving a layer of copper-depleted material surrounding the core of unreacted copper ore, as shown schematically in Figure 6.1. Protons must diffuse from the bulk of solution through this copper-depleted layer to the reaction front. Thus, for large particles, diffusion through this layer could be the rate-controlling step for the dissolution reaction. In the same context, Cu2+ forms at the reaction front and must diffuse from the reaction zone through the copperdepleted zone to the bulk of solution surrounding the ore particle.  Un-reacted core  Reacted rim Cu  2+  Cu ore particle +  H  Liquid phase (around particle)  Figure 6.1  Schematic of leaching copper from a round particle of copper ore  The acid reacts simultaneously with the copper minerals and other mineral phases, generally referred to as gangue. The acid and copper ions are thus considered to move from their original positions by the processes of advection and hydrodynamic dispersion within the steady state flow field of solution. In order to model the reaction and transport of copper from the ore due to heap leaching, the model developed in Chapter 3 has to be modified to account for continuous acid consumption and copper ion generation with time and in space. Results and Discussion – Performance Tests  118  In other words, as opposed to the inert solute transport considered earlier in modeling the results of the tracer tests, the solutes in the Performance tests (H+ and Cu2+) are not inert. Also, certain reaction parameters are required, thus necessitating various batch tests on the material to estimate them.  6.2  Reactive Solute Transport Model  The 2D axisymmetric model used to fit the data generated in the Performance Tests is general and comprises three material balances, namely, acid and dissolved copper balances in the leach solution, and unreacted solid phase (residual) copper in the ore.  6.2.1  Dissolution and Transport of Copper  In general, the leaching of copper from the ore involved an increase of the concentration of copper in the solution and a decrease of copper from the solid phase. The leaching of copper and transport of dissolved copper through the heap of ore are governed by several well-defined processes. These include dissolution of copper, which is a chemical reaction occurring between the lixiviant and the mineral particles of the ore, causing copper in the solid phase to be transferred into solution as cations, and sorption (including adsorption and precipitation) of copper (into the solid phase) occurring while the dissolved copper in solution migrated through the heap of ore. Usually, the adsorption phenomenon is represented by equilibrium isotherms such as those discussed in Equations 2.16 to 2.19. In the use of the model for Cu2+ transport in this work, it is assumed that Cu2+ is an inert solute, thus dissolved Cu2+ neither adsorbs nor reacts to form secondary solid phases. This assumption is valid for solution pH less than 4, because the sorption capacity of Cu2+ is then low and the solubility of secondary copperbearing minerals is high [Eriksson and Destouni, 1997].  However, this  assumption for acid heap leaching of a copper ore requires some caution. Neither the acid (in the influent solution) nor the Cu2+ (generated by the reaction)  Results and Discussion – Performance Tests  119  could be considered inert in the system. During the transport of Cu2+ in solution in the heap, any spreading of the solution outside the continuous channels of solution flow may result in Cu2+ encountering basic matter in the ore, such as limestone (CaCO3) or dolomite ((Ca,Mg)CO3). Acid in solution spreading to such areas could be totally neutralized, maintaining a largely basic environment and thus dropping its load of copper ions. In such a situation, Cu2+ would react to form secondary Cu precipitates such as basic copper sulphates/carbonates and simple copper hydroxides and oxides, as given in Equations (R6.1) to (R6.4), and thus form a solid-phase copper rim around the main solution channels. 3 CuSO4 + 2 Ca(OH)2 → Cu3(SO4)(OH)4 + 2 CaSO4  (R6.1)  4 CuSO4 + 3 Ca(OH)2 → Cu4(SO4)(OH)6 + 3 CaSO4  (R6.2)  CuSO4 + Ca(OH)2 → Cu(OH)2 + CaSO4  (R6.3)  2 CuSO4 + H2O + 2 CaCO3 → Cu2(OH)2CO3 + 2 CaSO4 + CO2  (R6.4)  The basic copper sulphates formed in (R6.1) and (R6.2) are antlerite and brochantite, respectively.  The pH of the environment determines the stability of  a particular phase. It is worth noting that calcium sulphate (gypsum) formed from these reactions can also impact on the hydrology of the bed with time and cause permeability problems.  It was further assumed that, apart from the minerals  (copper and the gangue) identified as acid consumers, the acid was not consumed by any other phase. The precipitated basic copper phases would generally re-dissolve once the pH of the solution reaching them fell below 4 [Eriksson and Destouni, 1997; Martens et al., 2003].  Thus, to minimize the loss of Cu2+ from solution (or its delay in  showing up in the effluent), the pH of the effluent should be maintained below 4. Two possible mechanisms may be responsible for Cu loss from solution. One is hydrological  immobilization  and  the  other  is  chemical  immobilization.  Hydrological immobilization is caused by the slow movement of solution in large  Results and Discussion – Performance Tests  120  parts of the heap because most of the irrigation solution is flowing through preferential channels.  Chemical immobilization is caused by precipitation of  dissolved products, forming secondary solid phases, or adsorption, or ion exchange on, for example, ferric oxyhydroxides [Eriksson and Destouni, 1997]. Hydrological immobilization and chemical immobilization may be distinguished by running a bottle-roll test on the ore.  If the loss persists in the bottle-roll  experiment, then it is likely to be chemical immobilization.  Cu ions are also  known to adsorb on aluminum or manganese hydroxides. This adsorption may be pH dependent. The solute transport model developed in this investigation does not incorporate the effects of pH and adsorption-desorption on the transport of Cu2+, at this stage.  Hence, no tests to determine the equilibrium isotherm for adsorption of  copper ion on the ore or its precipitation from solution were undertaken. A further study is required for comprehension of these mechanisms and their eventual inclusion in the reactive solute transport model.  6.2.2  Kinetics of Leaching  The kinetics of copper leaching and gangue acid consumption were represented by the following empirical power-law expressions, respectively:  where  φ sCu(a) = −sCu(s) = ρb kCuCCu(s) CHη2SO4 (a)  (6.1)  φ sH2SO4 (a) = − ρb kGACCH2SO4 (a) − 0.75 ρb kCuCCu(s) CHη2SO 4 (a)  (6.2)  sCu(a) = net rate of dissolved copper production, kmol/m3/s sCu(s) = net rate of solid copper production, kmol/m3/s  sH2SO4(a) = net rate of sulfuric acid production, kmol/m3/s ρb = ore bulk density, kg ore/m3 kCu = copper leaching rate constant kGAC = gangue acid consumption rate constant CCu(s) = concentration of unleached copper in the ore, kmol/kg ore  Results and Discussion – Performance Tests  121  CH2SO4(a) = concentration of acid, kmol/m3 water φ, η = empirical constants Empirical rate equation 6.1 states that the leaching rate of copper is proportional to the concentration of solid copper remaining to be leached to a certain constant power φ and to the concentration of acid to a certain power η. Braun et al. [1974] proposed a similar but simpler leaching model to represent the reaction between copper and sulphuric acid, assuming second order reaction kinetics of copper leaching, thus fixing the values of both φ and η to unity. Montero et al. [1994] and Muñoz et al. [1997] adopted the model of Braun et al. [1974], but both took the model a step further by accounting for copper adsorption, which they both represented with the Freundlich equilibrium isotherm. In the case of gangue acid consumption, a first-order acid consumption rate is adopted here, thus fixing the values of φ and η to zero and unity, respectively. The reason for the difference in the mathematical expressions used for copper leaching and gangue acid consumption is discussed below. Assuming the pH of the environment was always kept below 4, where there would be little or no loss of copper from solution due to chemical immobilization, and the change in the concentration of unleached copper in the solid phase at any time t after leaching can be obtained from mass balance as follows: dCCu(s) dt  =  sCu(s) ρb  (6.3)  It is important to remember that, in the solution of these equations, the following assumptions were made: Density, viscosity and temperature gradients of the fluid did not affect the distribution of solute velocities. The porous medium was incompressible, homogeneous and isotropic with respect to coefficients of longitudinal and transverse dispersivity.  Results and Discussion – Performance Tests  122  The influence of immobile water in the heap was negligible. The process was undertaken under a steady-state hydraulic regime. Thus, depending on the water flux irrigated onto the surface of the ore in the column and the hydraulic properties of the ore, the moisture content, θ, was determined. The consumption of sulfuric acid from the leach solution occurred due to dissolution of both copper and gangue matter while the solution is being transported through the heap.  6.2.3  Numerical Solution  The heap was discretized as a uniform grid of 2D axisymmetric control volumes having nodes and boundaries as defined by Patankar [1980]. The properties characterizing the medium and the initial and boundary conditions were defined for each cell. The scalar and vector properties were defined within the nodes and across the nodal boundaries, respectively, as detailed by Patankar [1980]. The model comprising the flow equation and the solute transport equations were solved using the strategies discussed in Chapter 3.  6.2.4  Calibration of the Model  Calibration of the model comprised determining the values of the parameters characterizing it.  The relatively large number of parameters in the model  required that some of them should be estimated separately. The hydraulic (Kw, m, n, hc,0) and dispersive (αT, αL) parameters were estimated from the results of the inert tracer tests undertaken in cylindrical columns of ore. The chemical reaction parameters comprising the copper and gangue leaching kinetics were estimated from the results of acid leaching the ore in bottle-roll tests (Section 4.3.2).  Results and Discussion – Performance Tests  123  6.2.5  Estimation of Hydraulic and Dispersive Parameters  The values of unsaturated hydraulic conductivity Kw, VGM parameters m and n, and the ratio of transverse to longitudinal dispersivities, αT/αL were determined from the results of tracer tests conducted in both the Small and Large Columns, as detailed in Chapter 5. It was assumed that these values would remain largely unchanged as a result of leaching with dilute sulfuric acid solution.  6.2.6 6.2.6.1  Estimation of Chemical Parameters Estimation of Leach Kinetic Parameters  The data used to estimate the parameters of the kinetic expressions (Equations 6.1 and 6.2) were obtained from leach tests undertaken by bottle-rolling specified masses of ore (Section 4.3.2). The elements that contributed significantly to sulfuric acid consumption during the experiments are copper, aluminum, calcium, iron, magnesium, manganese, potassium and sodium.  All elements leached  from the ore except copper in this case were referred to as gangue elements. The relative conversions of copper and acid from the three experiments with initial acid concentrations of 5, 10 and 20 g/L sulfuric acid are shown in Figures 6.2, 6.3 and 6.4, respectively. The conversion data generally show more rapid extraction around the initial moments of leaching, demonstrating high initial conversion rates followed by relatively low conversion rates. For heap leaching ores, this trend may be due to two main reasons. The first reason has to do with elements/minerals in the fine fractions of the ore leaching faster than those in the coarse fractions by virtue of the high surface areas of the fine particles. The second reason may be that the elements/minerals are associated with different mineral phases whose leaching rates can be classified as fast and slow. The rates of dissolution of copper and the constituents of the gangue phases and their evolution with time have to be quantified to estimate the parameters in the kinetic models.  Results and Discussion – Performance Tests  124  The bottle roll data were modeled using Equation 6.1.  Integrating using the  trapezoid rule, and noting that, for the bottle-roll tests, the concentration of acid was not constant but decreased monotonically (as shown in Figures 6.2 through 6.4), the kinetic parameters of leaching copper from the ore were determined. The values of kCu, φ and η in Equation 6.1 were determined by minimizing the sum of square error between the predicted and observed extraction values of copper for all three data sets, using the Solver routine in Microsoft Excel, assuming an ore bulk density value of 1,582 kg/m3 as previously determined. At a total sum of square error of all three data sets of 0.00696, the optimized values of kCu, φ and η for copper dissolution from this ore in acid solution of concentration between 5 and 20 g/L sulphuric acid are given Table 6.1. As can be seen from Figures 6.2 through 6.4, the copper extraction data are well predicted by the model. Table 6.1  Kinetic copper oxide leaching parameters from the ore by bottle roll-tests Parameter  Value  kCu (1/s) kGAC (m3/kg/s) φ η  1.53 × 1012 1.24 × 10–8 4.89 0.363  Results and Discussion – Performance Tests  125  0.05  100% 90%  0.04  70% 0.03  60% 50% Copper Data  40%  0.02  Copper leach model Acid Data  30%  Acid Concentration (mol/L)  Extraction of Copper (%)  80%  0.01  20% 10% 0% 0  20  40  60  80  100  120  140  160  0.00 180  Time (hr)  Copper extraction from bottle-roll test at 5 g/L initial acid  100%  0.10  90%  0.09  80%  0.08  70%  0.07  60%  0.06 Copper Data Copper Model Acid Data  50% 40%  0.05 0.04  30%  0.03  20%  0.02  10%  0.01  0% 0  20  40  60  80  100  120  140  160  Acid Concentration (mol/L)  Extraction of Copper (%)  Figure 6.2  0.00 180  Time (hr)  Figure 6.3  Copper extraction from bottle-roll test at 10 g/L initial acid  Results and Discussion – Performance Tests  126  Extraction of Copper (%)  0.20  90%  0.18  80%  0.16  70%  0.14  60%  0.12  Copper Data Copper Model Acid Data  50%  0.10  40%  0.08  30%  0.06  20%  0.04  10%  0.02  0% 0  20  40  60  80  100  120  140  160  Acid Concentration (mol/L)  100%  0.00 180  Time (hr)  Figure 6.4  Copper extraction from bottle-roll test at 20 g/L initial acid  Due to the low conversion of acid-consuming gangue minerals in this ore (aluminum, calcium, iron, magnesium, manganese, potassium and sodium) the acid consumption by all the gangue elements was lumped into a single rate equation; specifically, Equation 6.2. The decrease in copper extracted in the 5 g/L initial acid bottle-roll test is noted in Figure 6.2. This has been ascribed to precipitation of copper from solution due to depletion of acid from the bottle after 96 hours of leaching. Thus, the results suggest that the key to preventing metal loss due to chemical immobilization is to maintain acidic conditions in the effluent. The copper and gangue extraction data from the bottle-roll tests have shown some interesting results that are worthy of the attention of hydrometallurgists. The first one has to do with differentiating between the acid consumed by copper and that consumed by the gangue phase elements. This was accomplished by determining the stoichiometric sum of the acid consumed by all the major gangue elements and the remaining acid taken to be that consumed due to leaching of  Results and Discussion – Performance Tests  127  copper. The average acid concentration was plotted against this calculated rate of acid consumption by all gangue elements, and the result is shown in Figure 6.5.  Each curve in Figure 6.5 shows a rapid initial rate of acid consumption at  high acid concentration (the steep portion of each curve). This corresponds to the leaching of copper oxides.  Once the copper oxides are dissolved, the  diminished rate of acid consumption in each test falls onto a single straight line. This line represents first-order acid consumption by all of the gangue elements lumped together, and the slope of this line represents the first-order rate constant. The best-fit slope is 0.0094/h (with a degree of fit represented by R2 = 0.993), giving a rate constant (rate of acid consumption per unit mass of ore due to the gangue) of kGAC = 1.24 × 10–8 m3 liquid/kg ore/s during bottle roll leaching.  Gangue Acid consumption rate (mmol/L/hr)  12 Gangue Acid Consumption rate = 9.3912(Acid concentration) + 0.1161 R2 = 0.9927 10  8 5 g/L Acid 10 g/L Acid  6  20 g/L Acid Linear (1st order Model )  4  2  0 0.00  0.02  0.04  0.06  0.08  0.10  0.12  0.14  0.16  0.18  0.20  Acid concentration (mol/L)  Figure 6.5  Calculated rates of acid consumption during bottle roll leach tests  The rate of acid consumption by copper oxides can therefore be estimated as the difference between the points shown in Figure 6.5 and the corresponding point on the straight line at the same acid concentration. These estimated values are plotted against the average reaction rate of copper oxides in Figure 6.6. The  Results and Discussion – Performance Tests  128  high correlation coefficients (R2 > 0.996) confirm the goodness-of-fit for the three curves by straight lines of slopes 0.7660, 0.8363 and 0.8148, for the 5, 10 and 20 g/L initial acid concentrations, respectively. These slope values represent the moles of acid consumed per mole of copper leached.  12  Acid consumption rate (mol/s)  Acid consumption rate = 0.766(Copper reaction rate) - 0.0098 R2 = 0.9996 for 5 g/L Acid 10  Acid consumption rate = 0.8363(Copper reaction rate) - 0.2462 2 R = 0.9958 for 10 g/L Acid  8  Acid consumption rate = 0.8148(Copper reaction rate) - 0.6545 2 R = 0.9999 for 20 g/L Acid  6  4 5 g/L Acid 10 g/L Acid 20 g/L Acid  2  0 0  2  4  6  8  10  12  Copper reaction rate (mol/s)  Figure 6.6  Rate of acid consumption by copper oxides  This has an important implication about the mineralogy of the major copper minerals that were leached from the ore by dilute sulfuric acid. Assuming all copper is in the form of brochantite, CuSO4·3Cu(OH)2, this ratio should be 0.75. Since other oxide minerals of copper are present in the ore [Zarate et al., 1995, 2003], each will contribute to this ratio in accordance with its stoichiometric ratio, as well as its relative abundance. The reactions of the common oxides of copper with dilute sulphuric acid were given by Equations R2.3 to R2.9 in Section 2.2.2. The calculated ratios of moles acid consumed per mole of copper leached of the various oxides are given in Table 6.2. Attempts to characterize samples of the ore mineralogically using X-ray Diffraction and Reitveld analysis to determine the relative proportions of the copper containing minerals failed to yield useful  Results and Discussion – Performance Tests  129  information due to the low copper concentration in the ore. Values consistently near 0.75 obtained for all three tests suggest that brochantite is the dominant copper oxide mineral in this ore. Table 6.2  Stoichiometric sulphuric acid required to leach copper oxides  Copper Oxide Mineral Name  Formula  Reaction No.  Acid/Copper Mole Ratio  Malachite Azurite Cuprite Tenorite Chryosocolla Brochantite  CuCO3.Cu(OH)2 2Cu(CO3).Cu(OH)2 Cu2O CuO CuSiO3.2H2O CuSO4.3Cu(OH)2  R2.3 R2.4 R2.5 R2.7 R2.8 R2.9  1 1 0.5 1 1 0.75  6.2.6.2  Estimation of Leach Kinetic Parameters in Small and Large Columns of Ore  The kinetic parameters (kCu, φ and η for copper leaching and kGAC for gangue leaching) obtained for Equations 6.1 and 6.2 would then be used in the reactive transport model to simulate the reaction between the acid and the ore, and the transport of residual acid and dissolved copper through the bed and into the effluent solution. The values of kCu and kGAC were obtained from the data from the bottle-roll tests were further modified to be consistent for input sA in Equation 3.23. For copper, the following equation was used: φ 4.89 sCu(a) = ρb kCuCCu(s) CHη2SO 4 (a) = 1582 × 1.53 × 1012 CCu(s) CH02.363 SO 4 (a)  (6.4)  By stoichiometry, 0.75 mole of acid is consumed with every mole of copper leached from brochantite.  This rate, combined with the molar rate of acid  consumed by gangue leaching, gives the following expression for acid:  (  φ sH2SO4 (a) = − ρb kGACCH2SO4 (a) + 0.75k CuCCu(s) CHη2SO4 (a)  (  )  4.89 = −1582 1.24 × 10 −8 CH2SO4 (a) + 0.75 × 1.53 × 1012 CCu(s) CH02.363 SO 4 (a)  Results and Discussion – Performance Tests  )  (6.5)  130  Due to the nature of heaps as unsaturated fixed beds, the rate and extent of leaching should be smaller than the corresponding values from agitated systems with the same ore particle size distribution, such as bottle-rolling tests.  It is  important to note that the leach kinetic models used here are for the microscopic level [Braun et al., 1974, Montero et al., 1994]. However, they are being used for the macroscopic level at which measurable, differentiable and continuous variables can be determined and boundary conditions can be stated and solved. The macroscopic effects of the microscopic factors are retained in the form of effective coefficients or parameters whose numerical values can be determined experimentally [Montero et al., 1994]. The kinetic parameters estimated from the bottle-roll tests were thus adjusted to fit the data of the acid leach tests from the 10-cm mini-columns and the Large Column. Thus again, the values of the leach kinetic constants (kCu and kGAC) were adjusted until the predicted solute concentration profiles matched the observed concentration profiles of acid and copper in the effluent from these column tests. As expected, the kinetic parameters must be reduced significantly to fit the data.  6.3  Results and Discussion of the Reactive Solute Transport Model  The best-fit kinetic parameters for both copper and gangue dissolution as obtained for the acid and copper concentrations in the effluent from the 10-cm mini-columns (irrigated with solutions of concentrations 5, 10 and 20 g/L sulphuric acid) are given in Table 6.3 and the fits to the observed data are given in Figures 6.7 and 6.8.  The best-fit kinetic parameters for both copper and  gangue dissolution as obtained for the acid and copper concentrations in the effluent from the Large Column (irrigated with solutions of acid concentration of 10 g/L sulphuric acid) are also given in Table 6.3 and the fits to the observed data are given in Figure 6.9. It should be noted that during the sampling of the effluent for acid and copper analyses, the entire effluent received in the  Results and Discussion – Performance Tests  131  containers during a particular time interval was removed and weighed, thoroughly mixed and sampled for analysis. A new container of known mass replaced the one removed to receive a new batch of effluent. This was taken into account in calculating the solute (acid and copper) concentrations in the simulated effluent. Table 6.3  Optimal values of kinetic parameters for copper and gangue leaching from the10-cm mini-columns and the Large Column  Acid Concentration (g/L H2SO4)  10-cm mini-columns kCu  kGAC  Large Column kCu  kGAC  kCu(BR)/10 kGAC(BR) /3.5 5 kGAC(BR)/3.0 kCu(BR)/10 kGAC(BR) /2.5 kCu(BR)/10 10 kCu(BR)/10 kGAC(BR)/5.0 20 kCu(BR) and kGAC(BR) are the best-fitting values of kCu and kGAC from the bottle-roll acid leach tests. While the best-fitting value of kCu is kCu(BR)/10 in the 10-cm mini-columns and the Large Column, the best-fitting value of kGAC varies between kGAC(BR)/5 and kGAC(BR)/3 for the 10-cm mini-columns and is kGAC(BR)/2.5 for the Large Column. The reason for the higher value of kGAC in the Large Column than in the minicolumns may have to do with a higher degree of acid dispersion in the former than in the latter. It is reasonable that the values of the kinetic parameters are higher in the bottle-roll data than in the column data since the former system is an agitated reactor and the latter is a fixed bed unsaturated reactor irrigated from a drip emitter. For the model fitting of the 10-cm mini-columns leach data, the acid data fitted better than the copper data.  In particular, simulations in the first two days  suggest that more copper was exiting the columns than was observed. This trend was also seen in the results of the Large Column acid leach data as shown in Figure 6.9. This is the problem that will characterize the transport of reactive solutes in fixed beds such as heaps. The copper leached from the upper layer might precipitate (as mentioned above) or adsorb on the ore in the lower levels  Results and Discussion – Performance Tests  132  due to depletion of acid from the moving solution. (The loss of copper from solution might be exacerbated by the presence of basic matter in the ore, such as limestone or dolomite). As can be seen from Figure 6.8 the model followed the data closely on days 2, 1.5 and 1, for the 5, 10 and 20 g/L acid irrigations, respectively.  In Figure 6.9, the model followed the Large Column copper  extraction closely only after day 8. The water content profile in the Large Column during the acid leach is shown in Figures 6.10a to 6.10c for periods of 3.5, 7.5 and 15.5 days’ of leaching. As expected, these profiles are essentially identical, since the water content profile reaches steady state fairly quickly.  The  corresponding copper and acid concentration profiles are shown in Figures 6.11 and 6.12, respectively. The water content shows a conical shape leaving the ‘shoulders’ with the lowest moisture content and with water spreading from the central axis to the side of the column. (The bed of ore was initially irrigated with water at the same flux used to irrigate the acid solution prior to irrigation of the acid solution, so there is no difference between water content profiles for 3.5, 7.5 and 15.5 days).  The copper profiles follow the water profiles closely as the  dissolved copper is expected to be inert in the heap. If this inert characteristic is removed, then the copper profile will be expected to be similar to that of the acid (Figure 6.12). The acid spreading (Figure 6.12a to 6.12c) is limited to the vertical zone under the drip emitter, because of its reaction with gangue components in the ore. Due to insufficient acid reaching the peripheral areas of the column, dissolved copper in solution in these areas would precipitate to form secondary copper phases. This precipitation of dissolved copper from solution is the cause of delay in the copper exiting the column in the effluent as seen in Figures 6.8 and 6.9.  The precipitated (or adsorbed) copper would re-dissolve once acid  arrived at these peripheral areas (generally below pH 4) with prolonged irrigation of acid solution. For the model to capture the delay of copper exiting the column (heap) in the early days of irrigation, the mechanisms and factors responsible for the precipitation and/or adsorption of copper ions from solution should be incorporated into the model.  Results and Discussion – Performance Tests  133  0.16  Acid Concentration (mol/L)  0.14 0.12 5 g/l Acid Data 10 g/l Acid Data 20 g/l Acid Data Model  0.10 0.08 0.06 0.04 0.02 0.00 0  2  4  6  8  10  12  Time (days)  Figure 6.7  Effluent acid concentrations from the 10-cm mini-columns  Throughout the 15 days of running the acid leach test in the Large Column, the effluent samples contained virtually no acid. The simulations that fit the observed copper profile show extremely low acid concentrations (as shown in Figure 6.9) which would be below the analytical detection limit of the titration method used for sulfuric acid analysis. On the whole, the reactive model, which is an extension of the inert tracer transport model, represented the acid and copper concentrations in the effluent from the 10-cm mini-columns run with 5, 10 and 20 g/L and the Large Column run with 10 g/L acid.  Results and Discussion – Performance Tests  134  0.10  Copper Concentration (mol/L)  0.09 0.08 0.07 5 g/l Acid Data 10 g/l Acid Data 20 g/l Acid Data Model  0.06 0.05 0.04 0.03 0.02 0.01 0.00 0  2  4  6  8  10  12  Time (days)  Effluent copper concentrations from the 10-cm mini-columns  Figure 6.8 0.06  4.5E-04 Cu Data  0.05  4.0E-04 3.5E-04  0.04  3.0E-04 2.5E-04  0.03 2.0E-04 0.02  1.5E-04 1.0E-04  0.01  Acid Concentration in Effluent (mol/L)  Cu Concentration in Effluent (mol/L)  Cu Model Acid Model  5.0E-05 0.00  0.0E+00 0  2  4  6  8  10  12  14  16  Time (days)  Figure 6.9  Effluent copper and acid concentrations from the Large Column  Results and Discussion – Performance Tests  135  Column height  0.24-0.25 0.23-0.24 0.22-0.23 0.21-0.22 0.2-0.21 0.19-0.2 0.18-0.19 0.17-0.18 0.16-0.17 0.15-0.16 0.14-0.15  Column diameter  Figure 6.10 Predicted water content profile in the Large Column a) after 3.5 days of acid leaching  Column height  0.24-0.25 0.23-0.24 0.22-0.23 0.21-0.22 0.2-0.21 0.19-0.2 0.18-0.19 0.17-0.18 0.16-0.17 0.15-0.16 0.14-0.15  Column diameter  Figure 6.10 b) after 7.5 days of acid leaching  Results and Discussion – Performance Tests  136  Column height  0.24-0.25 0.23-0.24 0.22-0.23 0.21-0.22 0.2-0.21 0.19-0.2 0.18-0.19 0.17-0.18 0.16-0.17 0.15-0.16 0.14-0.15  Column diameter  Column height  Figure 6.10 c) after 15.5 days of acid leaching  0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 6.11 Predicted aqueous Cu concentration profile in the Large Column a) after 3.5 days of acid leaching (Cu concentration in mol/L)  Results and Discussion – Performance Tests  137  Column height  0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Column height  Figure 6.11 b) after 7.5 days of acid leaching (Cu concentration in mol/L)  0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 6.11 c) after 15.5 days of acid leaching (Cu concentration in mol/L)  Results and Discussion – Performance Tests  138  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 6.12 Predicted acid concentration profile in the Large Column a) after 3.5 days of acid leaching (acid concentration in mol/L)  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 6.12 b) after 7.5 days of acid leaching (acid concentration in mol/L)  Results and Discussion – Performance Tests  139  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 6.12 c) after 15.5 days of acid leaching (acid concentration in mol/L)  Results and Discussion – Performance Tests  140  Chapter 7 SENSITIVITY OF COPPER EXTRACTION TO HEAP AND SOLUTE PARAMETERS 7.1  Introduction  The experimental results described in Chapters 5 and 6 are limited to data generated from laboratory scale tests and the estimation of parameters from the solute transport model. The information obtained helped to identify the main hydraulic, solute transport and chemical reaction parameters of the model developed for heap leaching of copper from the Mantoverde copper oxide ore. In this chapter, the sensitivity of the copper extraction from this ore to variations in the parameters of the model is undertaken. The aim of the sensitivity analysis was to estimate the rate of change in the output of the reactive model (copper extraction) with respect to changes in the values of some of the model inputs (parameters). Such knowledge would be useful for: evaluating the range of applicability of the model, determining parameters that can be manipulated at low cost to improve process performance, and understanding the behavior of the system being modeled. The choice of a sensitivity analysis method depends to a great extent on the sensitivity measure employed, the desired accuracy in the estimates of the sensitivity measure, and the computational cost involved. The method used in this work comprised running the original model for a set of input/parameter combinations and estimating the sensitivity/uncertainty of the model output. This approach is often used to evaluate the robustness of a model, by testing whether  Sensitivity of Copper Extraction to Heap and Solute Parameters  141  the model response changes significantly with changes in model parameters and/or the structural formulation of the model. The primary advantage of this approach is that it accommodates both qualitative and quantitative information regarding variation in the model. However, the sensitivity information obtained depends to a great extent on the choice of the sample points, especially when only a small number of simulations can be performed. The impacts on copper extraction of such factors as emitter spacing, irrigation rate, heap height, ore permeability, dispersion parameters, particle reaction kinetics and acid concentration have been estimated by simulation of copper extraction from the Large Column, using the parameter values in Table 7.1 as the baseline. The fraction of copper extracted from the heap, F, was calculated using the formula:  F=  vw  t  CCu(s),0 ρb Z ∫0  CCu(a)dt  (7.1)  The integral in Equation 7.1 was computed using the trapezoid rule on the timecopper concentration data generated from the model simulations utilizing the optimized parameters.  Sensitivity of Copper Extraction to Heap and Solute Parameters  142  Table 7.1  Parameter values used as the baseline case for simulations in the Large Column  Parameter Height of heap, Z Radius (emitter half-spacing), R Bulk density of heap, ρbulk Hydraulic Conductivity, Kw Saturated Water Content, θs Residual Water Content, θr Van Genuchten parameter, m Van Genuchten parameter, n Air entry capillary head, hc,0 Molecular diffusion coefficient, D0 ε coefficient ε exponent α ratio, αT/αL α slope α intercept Irrigation flux, vw Copper concentration in feed solution Acid concentration in feed solution Copper leaching rate constant, kCu Gangue leaching rate constant, kGAC  Table 7.2  Baseline value  Units  4 0.5 1582 0.1570 0.3797 0.1455 0.4213 7.137 0.05 1.86 × 10–9 3.0 0.5 1/30 43240 0.06929 10 0 0.1 1.53 × 1011 4.96 × 10–9  m m kg/m3 m/s m3/m3 m3/m3  m m2/s  L/m2/h kmol/m3 kmol/m3 /s 3 m /kg/s  Ratio of quantity of copper extracted from heaps of varying emitter spacing with time  Time (days)  0.25  Emitter half-spacing (m) 0.50  4 15  1.00 1.00  1.80 0.98  Sensitivity of Copper Extraction to Heap and Solute Parameters  1.00 2.40 0.64  143  7.2 7.2.1  Hydraulic Properties Emitter Spacing  For short emitter spacing (small R), a lower volume of solution contacts the surface of the heap per unit time, spreading to a comparatively wider area/volume and thus effecting more copper leaching and extraction. For large emitter spacing, the high solution flux to the surface would cause the solution to exit the heap quickly without much spreading. Thus larger R would show initial higher copper extraction but would be overtaken by the copper extraction from the smaller emitter spacing as shown in Figure 7.1. Also shown in Table 7.2 above is the change in the ratios of quantity of copper extracted with time and emitter spacing. Furthermore, as seen from the acid profiles in Figures 7.2a to 7.2c, the acid spreads relatively more in the narrower column than the wider columns.  The fact that the acid zone has reached the bottom of the wider  column supports the earlier assertion that initial rate of extraction of copper in the wider column is higher than from the narrower column. Once the copper in the path of the acid in the wider column is extracted, any further acid passing through would only be consumed partly by gangue with the remaining acid exiting the column without spreading to wider zones to effect copper extraction. Though the acid has not broken through to the bottom of the narrower columns, the solution spreading has caused more copper extraction from these columns than in the wider ones. A danger of wider emitter spacing is the high volume of solution hitting the surface of the heap could result in formation and maintenance of preferential flow (channeling) of solution in the heap, as discussed in Chapter 2.  Sensitivity of Copper Extraction to Heap and Solute Parameters  144  0.20 F (R = 0.25 m)  0.18  F (R = 0.50 m) F (R = 1.00 m)  Fraction Copper extracted  0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0  2  4  6  8  10  12  14  16  Time (days)  Simulation of varying emitter spacing on copper extraction  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column height  Figure 7.1  Column diameter  Figure 7.2  Acid profiles for copper extraction with varying emitter spacing a) R = 0.25 m (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  145  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  b) R = 0.50 m (acid concentration in mol/L)  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column height  Figure 7.2  Column diameter  Figure 7.2  c) R = 1.00 m (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  146  7.2.2  Heap Height  The initial rate and extent of copper extraction increased with decrease in bed height (all other parameters remaining unchanged) as illustrated in Figure 7.3. As shown in the corresponding acid profiles in Figures 7.4a to 7.4c, the acid has broken through the shorter heap with high concentration exiting the column while the bottoms of the taller heaps were not reached by acid to effect copper extraction. The extent of copper extraction from the 4m and 6m-high heaps surpassed the 2m-high heap after days 6 and 8, respectively, because of wider acid coverage in the 4 and 6 m-high columns and hence more copper-acid contact. This corroborates the practice of leaching heaps to near completion in lifts and then loading another lift onto the already leached one before proceeding with leaching or removing the already leached heap and replacing it with a similar short heap. However, process economics is the main factor in the selection of heap height. The higher capital and operating costs associated with operating shorter heaps usually outweigh the advantage of faster metal recoveries. 0.18 Z=2m Z=4m Z=6m  0.16  Fraction Copper extracted  0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0  2  4  6  8  10  12  14  16  Time (days)  Figure 7.3  Simulation of varying heap height on copper extraction  Sensitivity of Copper Extraction to Heap and Solute Parameters  147  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.4  Acid profiles for copper extraction with varying heap height a) Z = 2 m (acid concentration in mol/L)  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.4  b) Z = 4 m (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  148  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.4  7.2.3  c) Z = 6 m (acid concentration in mol/L)  Irrigation Flux  In this section, the effect of varying the solution flow rate on the extraction of copper from the ore in the Big Column was simulated for solution fluxes of 5, 10 and 15 L/m2/h. Simulation results are shown in Figure 7.5.  High solution flow  rates resulted in more irregular flow patterns as compared to low flow rates. Consequently, the extent of longitudinal and transverse dispersion is higher under high flow rates, leading to higher copper extraction. The higher solution flux has caused more acid dispersion to the reaction sites to effect the leaching of copper and enough acid to maintain the dissolved copper in solution to the effluent from the heap. The acid profiles (Figures 7.6a to 7.6c) confirm these findings, showing increased axial and radial acid spreading with increase in irrigation flux.  Sensitivity of Copper Extraction to Heap and Solute Parameters  149  0.18 Vw = 5 L/m2/h Vw = 10 L/m2/h Vw = 15 L/m2/h  0.16  Fraction Copper extracted  0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0  2  4  6  8  10  12  14  16  Time (days)  Simulation of varying irrigation flux on copper extraction  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column height  Figure 7.5  Column diameter  Figure 7.6  Acid profiles for copper extraction with varying irrigation flux a) vw = 5 L/m2/h (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  150  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  b) vw = 10 L/m2/h (acid concentration in mol/L)  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column height  Figure 7.6  Column diameter  Figure 7.6  c) vw = 15 L/m2/h (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  151  7.2.4  Hydraulic Conductivity  The permeability of the heap is an essential process parameter that controls the irrigation and aeration (in the case of heap bioleaching) of the bed. During ore stacking, segregation may take place (particularly in the absence of agglomeration), resulting in coarse particles at the bottom and finer particles at the top. This situation leads to an increase in permeability with depth. A counter effect is when fines at the top are mobilized by the irrigation solution and plug the void spaces between the larger particles towards the bottom of the heap, resulting in a decrease in permeability with depth. Acid leaching also generally opens up the bed while leaching under alkaline conditions can result in plugging by fine precipitates. The implication is that the permeability of the ore bed will vary in time and space during the life of the leach pad. The effects of fine particles can be reduced by agglomeration. Wetting and agglomeration of the ore prior to stacking is known to increase the average bed permeability by 1 to 2 orders of magnitude for ore crushed to between 3/8” and 1” [Bartlett, 1996]. The extraction of copper from this ore is not sensitive to varying Kw between 0.1570 m/s ± 50%. This is expected as the absolute permeability and wetting phase viscosity, which are the factors upon which hydraulic conductivity depends, are without effect on the steady-state saturation profile if the heap permeability is uniform, as pointed out by Dixon [2003], and hence variations of hydraulic conductivity have no effect on spreading of moisture and copper extraction from the heap as shown in Figure 7.7. Furthermore, the acid profiles (Figures 7.8a to 7.8c), show that there is basically no difference in the acid distribution due to changes in the value of hydraulic conductivity.  Sensitivity of Copper Extraction to Heap and Solute Parameters  152  0.18 Kw = 0.0785 m/s  0.16  Kw = 0.1570 m/s Kw = 0.3140 m/s  Fraction Copper extracted  0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0  2  4  6  8  10  12  14  16  Time (days)  Simulation of varying heap permeability on copper extraction  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column height  Figure 7.7  Column diameter  Figure 7.8  Acid  profiles  for  copper  extraction  with  varying  hydraulic  conductivity a) Kw = 0.0785 m/s (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  153  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  b) Kw = 0.1570 m/s (acid concentration in mol/L)  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column height  Figure 7.8  Column diameter  Figure 7.8  c) Kw = 0.3140 m/s (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  154  7.2.5  Air Entry Capillary Head  Decreasing the air entry head hc,0 (all other parameters remaining constant) increases the rate and extent of copper extraction from the heap, as shown in Figure 7.9. This is a hydrological parameter that affects the transport of solutes within the heap. The coarser the particle size of the ore, the lower hc,0, and the faster the solute transport. This is shown in the acid profiles (Figures 7.10a to 7.10c), with acid breaking through the heap with a lower hc,0 value with lower radial spreading than the heap with a higher hc,0 value. The wider acid spreading will eventually cause higher copper extraction from heaps with higher hc,0 values. This was also discussed in Chapter 5 for inert solute transport, where a lower value of hc,0 was required for earlier breakthrough of irrigation solution (to fit the inert tracer curves in Figures 5.13 and 5.14). 0.18 hc,0 = 0.02 m 0.16  hc,0 = 0.05 m hc,0 = 0.08 m  Fraction Copper extracted  0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0  2  4  6  8  10  12  14  16  Time (days)  Figure 7.9  Simulation of varying air entry head on copper extraction  Sensitivity of Copper Extraction to Heap and Solute Parameters  155  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.10 Acid profiles for copper extraction with varying air-entry head a) hc,0 = 0.02 m (acid concentration in mol/L)  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.10 b) hc,0 = 0.05 m (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  156  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.10 c) hc,0 = 0.08 m (acid concentration in mol/L)  7.3 7.3.1  Transport Properties Dispersion Coefficients  Increasing either the longitudinal dispersivity αL or the transverse-to-longitudinal dispersivity ratio αT/αL, causes the acid to disperse more widely in the heap and hence increases the extent of copper extraction (as shown in Figures 7.11 and 7.12). From these figures, the extraction of copper from the ore is more sensitive to αL than to αT/αL. From the acid profiles (Figures 7.13a to 7.13c) for αL (αT maintained constant), there is no change in the radial acid spreading, but a substantial increase in axial spreading with increase in αL. In the case of αT/αL, there is a small decrease in axial spreading and a marginal increase in the radial spreading of acid with increase in αT/αL as shown in Figures 7.14a to 7.14c. Thus, copper extraction from the heap is more sensitive to changes in αL than to αT/αL, because of the low αT/αL values used for fitting the data.  Sensitivity of Copper Extraction to Heap and Solute Parameters  157  0.20 alpha slope = 21620 m  0.18  alpha slope = 43240 m alpha slope = 86480 m  Fraction Copper extracted  0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0  2  4  6  8  10  12  14  16  Time (days)  Figure 7.11 Simulation of varying the slope of αL vs vw on copper extraction 0.20 alpha ratio = 1.6667%  0.18  alpha ratio = 3.3333% alpha ratio = 6.6667%  Fraction Copper extracted  0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0  2  4  6  8  10  12  14  16  Time (days)  Figure 7.12 Simulation of varying the transverse-to-longitudinal dispersivity ratio  αL/αT on copper extraction  Sensitivity of Copper Extraction to Heap and Solute Parameters  158  Column height  0.07-0.08 0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.13 Acid profiles for copper extraction with varying longitudinal  dispersivity a) slope of αL vs vw = 21620 m (acid concentration in mol/L)  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.13 b) slope of αL vs vw = 43240 m (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  159  Column height  0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Column height  Figure 7.13 c) slope of αL vs vw = 86480 m (acid concentration in mol/L)  `  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.14 Acid profiles for copper extraction with varying transverse-to-  longitudinal dispersivity ratio a) αT/αL = 1.6667% (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  160  Column height  `  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.14 b) αT/αL = 3.3333% (acid concentration in mol/L)  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.14 c) αT/αL = 6.6667% (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  161  7.3.2  Diffusion Coefficient  Variations in the molecular diffusion coefficient (with all other parameters remaining constant) do not have a significant effect on copper extraction from the ore. The effect of D0 on copper extraction should be similar to that of αL and αT/αL (Figure 7.11and 7.12), but to a lesser extent. Thus, a higher value of D0 should cause more acid dispersion in the heap and hence higher copper extraction. But its effect is negligibly small, as shown in Figure 7.15. It was pointed out in Chapter 5 that D0 is not really a critical parameter in solute transport modeling. This justifies the assertion not to undertake any extensive experimental work to estimate its value in this investigation, but to estimate it from the literature. It also justifies selecting a common value of D0 for all solutes, which renders the model code much more efficient by eliminating the need to evaluate transport parameters for each solute individually. 0.20 D0 = 0  0.18  D0 = 1.86E–09 m2/s  Fraction Copper extracted  0.16  D0 = 3.37E–09 m2/s  0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0  2  4  6  8  10  12  14  16  Time (days)  Figure 7.15 Simulation of varying molecular diffusivity coefficient on copper  extraction  Sensitivity of Copper Extraction to Heap and Solute Parameters  162  7.4  Chemical Reaction Parameters  7.4.1  Particle Reaction Rate Constants  Increasing the gangue dissolution rate constant kGAC decreases the rate of copper extraction as shown in Figure 7.16, as it decreases the quantity of acid available for copper leaching and its subsequent maintenance in solution. The acid profiles in Figures 7.17a to 7.17c confirm this, with decreasing acid levels with depth of the heap with increase in kGAC. 0.25 kGAC = 2.48E–09 m3/kg/s kGAC = 4.96E–09 m3/kg/s  Fraction Copper extracted  0.20  kGAC = 9.92E–09 m3/kg/s  0.15  0.10  0.05  0.00 0  2  4  6  8  10  12  14  16  Time (days)  Figure 7.16 Simulation of varying gangue dissolution rate constant on copper  extraction Although the value of kCu used is always higher than the value of kGAC, the copper extraction is more sensitive to kGAC than to kCu, because of low copper proportion in the ore and its faster extraction rate than the gangue. The latter consumes acid ad infinitum after all copper has been extracted from a zone.  Sensitivity of Copper Extraction to Heap and Solute Parameters  163  Copper leaching and gangue acid consumption are parallel reactions with the latter being much slower than the former. The higher the ratio kCu/kGAC, the more selective the leach for copper and hence the lower the capital and operating costs of a solution purification plant (e.g., a solvent extraction facility) downstream of leaching.  Column height  0.07-0.08 0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.17 Acid profiles for copper extraction with varying gangue acid  consumption rate constant a) kGAC = 2.48 × 10–9 m3/kg/s (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  164  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Column height  Figure 7.17 b) kGAC = 4.96 × 10–9 m3/kg/s (acid concentration in mol/L)  0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.17 c) kGAC = 9.92 × 10–9 m3/kg/s (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  165  On the other hand, as expected, increasing the copper leaching rate constant kCu increases the rate of copper extraction as shown in Figure 7.18. Increasing the copper reaction rate constant decreases the acid concentration with depth in the heap as shown in Figures 7.19a to 7.19c. 0.20 kCu = 7.64E+10 /s  0.18  kCu = 1.53E+11 /s kCu = 3.05E+11 /s  Fraction Copper extracted  0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0  2  4  6  8  10  12  14  16  Time (days)  Figure 7.18 Simulation of varying the copper leaching rate constant on copper  extraction  Sensitivity of Copper Extraction to Heap and Solute Parameters  166  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.19 Acid profiles for copper extraction with varying copper leaching rate  constant a) kCu = 7.64 × 1010 /s (acid concentration in mol/L)  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.19 b) kCu = 1.53 × 1011 /s (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  167  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.19 c) kCu = 3.05 × 1011 /s (acid concentration in mol/L)  7.4.2  Acid Concentration  Increasing the concentration of acid in the irrigation solution increases the rate of copper extraction as shown in Figure 7.20, as more acid is available for leaching and maintaining copper in solution. The acid profiles shown in Figures 7.21a to 7.21c confirm that the there is a marginal increase in radial acid concentration and a substantial increase in the axial acid concentration with depth, resulting in increase copper extraction from the heap. The acid concentration in the irrigation solution will depend on the relative rates of metal leaching and gangue dissolution, and hence on the degree of selectivity desired. There are operating cost implications for use of excess acid in leaching and for downstream treatment processes to purify the desired metal from the pregnant solution as mentioned above.  Sensitivity of Copper Extraction to Heap and Solute Parameters  168  0.25 0.05 M Acid 0.10 M Acid 0.20 M Acid  Fraction Copper extracted  0.20  0.15  0.10  0.05  0.00 0  2  4  6  8  10  12  14  16  Time (days)  Figure 7.20 Simulation of effects of varying acid concentration on copper  extraction  Sensitivity of Copper Extraction to Heap and Solute Parameters  169  Column height  0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.21 Acid profiles for copper extraction with varying acid concentration a) CH2SO4 = 0.05 mol/L (acid concentration in mol/L)  Column height  0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.21 b) CH2SO4 = 0.10 mol/L (acid concentration in mol/L)  Sensitivity of Copper Extraction to Heap and Solute Parameters  170  Column height  0.13-0.14 0.12-0.13 0.11-0.12 0.1-0.11 0.09-0.1 0.08-0.09 0.07-0.08 0.06-0.07 0.05-0.06 0.04-0.05 0.03-0.04 0.02-0.03 0.01-0.02 0-0.01  Column diameter  Figure 7.21 c) CH2SO4 = 0.20 mol/L (acid concentration in mol/L)  7.5  Practical Applications of Sensitivity Analysis for Heap Leaching  While the rate and extent of copper extraction from a heap are functions of all the parameters discussed above, emitter spacing, heap height, solution irrigation flux and lixiviant (acid) concentration are the only ones that can be practically controlled to enhance the process. Dispersivity parameters can be influenced by one or more of the above-mentioned parameters; particularly, the irrigation flux. The chemical reaction parameters are functions of the reactivity of the minerals and the gangue with the lixiviant and the various prevailing conditions. Thus the extraction of copper from the heap may be enhanced by: reducing the emitter spacing reducing the height of the bed  Sensitivity of Copper Extraction to Heap and Solute Parameters  171  increasing the solution irrigation flux to the extent that does not cause solution channeling and/or threaten the structural integrity of the heap increasing the concentration of acid to the point that will not compromise selectivity of copper leaching over gangue dissolution.  Sensitivity of Copper Extraction to Heap and Solute Parameters  172  Chapter 8 CONCLUSIONS AND RECOMMENDATIONS 8.1  Summary and Conclusions  A 2D axisymmetric model was developed to describe the leaching of metals and the transport of water and solutes during heap leaching operations. This model was solved using control volume discretization and approximate factorization. The model encompasses various physico-chemical processes that occur during heap leaching. These comprise (for copper leaching from simple oxide ores) transport of acid in the irrigated solution through the heap, chemical reaction between the acid and leachable mineral phases (metal oxides and gangue) in the ore, and transport of the dissolved phases (metal ions) and residual acid through the bed of ore to the effluent.  The leaching process was modeled by an  empirical rate law while solute transport was characterized by parameters for advection and dispersion in the ore bed under unsaturated conditions. These parameters were identified, analyzed and estimated from experimental data. The model was calibrated by means of data generated from specially designed laboratory experiments. The first set of experiments comprised inert tracer tests on heaps of Mantoverde copper ore in cylindrical columns of different sizes. For these tests, miscible displacement experiments were undertaken to measure and/or estimate the unsaturated hydraulic properties (namely, the water retention – matric pressure head relationship and hydraulic conductivity of the heap) and inert solute transport characteristics (consisting of advective and dispersive parameters) from breakthrough curves (BTCs) generated from resident and effluent solute concentrations. The comparison of parameters at several depths to assess the solute transport model is an obvious advantage of in situ measurement of BTCs.  The second set of experiments consisted of batch  agitated leach experiments to estimate kinetic parameters of metal (copper) and gangue leaching in acid solutions. Conclusions and Recommendations  The third (final) set of experiments was 173  designed to evaluate the capability of the model to predict reactive solute transport and comprised leaching the same ore packed in columns of different sizes and irrigated by dilute sulphuric acid solution. The parameters estimated from the first two sets of experiments were used successfully to model the leaching, and transport of residual acid and dissolved copper in the effluents from these columns. Different techniques were employed to determine the parameters, including direct estimation and iterative adjustment of model parameters by trial and error until  the  simulated  dependent  variable  (for  example,  effluent  solute  concentrations or resident solute concentrations as measured by a probe in the ore bed) matched the corresponding experimental data. Furthermore, sensitivity of copper extraction from the ore by heap leaching was simulated to determine parameters that can be manipulated economically to enhance the process and to shed more light on the applicability of the model. The experimental work and the subsequent simulations show that hydraulic properties influence the water retention and distribution and hence solute transport characteristics of the heap. Opportunities thus exist for improving the extraction of metals from heap leaching by manipulating the system design and operating strategies which affect the hydraulic characteristics of the ore bed. The Brooks-Corey model adequately described the hydraulic characteristics of the Mantoverde ore. Also as expected, the lower the water content of the ore bed, the higher its suction (matric) potential and within the flux range investigated in this work, the higher the irrigation flux, the faster the solute exited the column in the effluent. It was observed that within the limited water content range relevant to safe and successful operation of a heap leach facility, dispersivity increased with an increase in irrigation flux and hence water content.  High solution flow rates  resulted in more irregular flow patterns as compared to low flow rates.  Conclusions and Recommendations  174  Consequently, the degree of longitudinal and transverse dispersion is higher under high flow rates. A sensitivity analysis showed that faster copper extraction results from higher solution flux. The empirical rate model used for copper and gangue dissolution represented the data well. The predicted leach kinetic parameters from the agitated leached system were readily modified for the fixed–bed leach data, confirming the adequacy of the empirical approach. Due to the fixed bed nature of heaps and unsaturated irrigation, the rate and extent of leaching were smaller than the corresponding values from the agitated system (bottle-rolling) with the same ore particle size distribution. Thus, as expected, the kinetic parameters of dissolution from the heap (fixed-bed) were 10% of the values estimated for copper, and 20 to 40% for gangue, in the agitated leach system. The sensitivity analysis showed that copper extraction from heap leaching the ore is not particularly sensitive to hydraulic conductivity and solute molecular diffusion coefficient. It is slightly sensitive to air entry head, the magnitude of longitudinal dispersivity, the ratio of transverse to longitudinal dispersivity, and the kinetics of copper leaching, but is highly sensitive to drip emitter spacing, heap height, irrigation flux, acid concentration, and the kinetics of gangue dissolution. On the whole the modeling approach captured the trends of the water distribution and transport of both inert and reactive solutes in heaps, clearly demonstrating transport of solutes to be the rate limiting step in the process. Close agreement between observed data and predicted results lends credence to the model and its solution methodology. The model can be used to simulate the impact of ore hydraulic properties and irrigation fluxes on water flow and solute transport during the extraction of metals from a heap.  Since solute transport is a  hydrological issue depending on the hydraulic parameters of the heap, hydrology of the system should be incorporated in heap leaching models.  There is,  therefore, a clear need to incorporate the ore properties in the design and  Conclusions and Recommendations  175  operating strategies of heaps for efficient lixiviant distribution, leaching and recovery of metals. It can also be used for optimal design and various process enhancement strategies for efficient operation, as well as economic evaluations and potential environmental impact assessments of heap leaching operations. Apart from its use in heap leach facility design and operation, the model developed in this work can also be incorporated into heap bioleaching models such as HeapSim [Dixon and Petersen, 2003].  It may also be used for  modeling: reagent rinsing from heaps prior to de-commissioning [Decker and Tyler, 1996] soil washing by flushing with acid to remove heavy metals from contaminated soils [Brusseau, 1996] and contaminated agricultural soils in-pit tailings management facilities (TMF) for disposal of uranium tailings [West et al, 2003]  8.2  Recommendations  This investigation is an initial attempt to incorporate ore hydraulic characteristics into the modeling of typical heap leaching operations. Further improvements to this approach will be required on more targeted fieldwork to complement progress with the modeling and scenario analyses. Particular areas requiring further investigation and possible extension of the model include: the effects of chemical, biological and physical changes on solute transport during heap leaching, including the effects of formation of such phases as sulfur and jarosite (and associated hydrophobicity), changes in pH, neutralization of acid, precipitation of copper hydroxides and degradation of ores during the life of the heap the effects of changes in viscosity and density of solutions recycled from solvent extraction plants for heap leaching  Conclusions and Recommendations  176  the effects of solution flow, air flow, temperature distribution and heat transport capacity of solution and/or air for optimal microbe activity during heap bioleaching the effects of two-phase flow, by investigating the disturbance caused by oxygen to the flow field and the effects (if any) of dissolution of oxygen gas on the rates of chemical reactions the influence of possible stagnant or immobile solution (or a fraction of the heap with lower solution velocity than the mean solution velocity in the dominant flow domains) in the heap and its relationship with the mean flow velocity of the solution through the heap the influence of particle size and textural structure of ores on water flow, solute transport and kinetics of mineral extraction the effect of the ionic strength on reactive solute transport the effects of evaporation of water and rainfall on the extraction of metals from heaps Given the complexities of the physical and chemical reactions and their impact on the hydrological and chemical environment of heaps, modeling cannot fully replace laboratory and pilot plant experimentation for optimal design and operation of heaps.  The modeling effort should instead supplement the  experimental work, in order to save time and money, and to provide a means for the more intelligent interpretation of the results of laboratory and pilot plant trials.  Conclusions and Recommendations  177  References Abeluik, R., Identification of unsaturated solute transport parameters. Groundwater Management: Quantity and Quality (Proceedings of the Benidorm Symposium, October 1989), IAHS Publ. No. 188, 1989, pp 261–270. Amey, E.B. Gold, Chapter 33 in Minerals Year Book, Volume 1, US Department of Interior, US Geological Survey, Metals and Minerals, 1999, pp 33.1–33.14 Aral, M.M. and Liao, B. Analytical solutions for two-dimensional transport equation with time dependent dispersion coefficient. J. Hydrological Engng. 1 (1996) 20–32. Arya, L. M., and Paris, J. F. A physicoempirical model to predict the soil moisture characteristic form particle-size distribution and bulk density data. Soil Sci. Soc. Am. J. 45 (1981) 1023–1030. Arya, L. M., Leij, F.J., Shouse, P.J. and van Genuchten, M.T. Relationship between the hydraulic conductivity function and particle-size distribution. Soil Sci. Soc. Am. J. 63 (1999) 1063–1070. Arya, L. M., Leij, F.J., van Genuchten, M.T. and Shouse, P.J. Scaling parameter to predict soil water characteristic from particle-size distribution data, Soil Sci. Soc. Am. J. 63 (1999) 510–519. Bartlett, R.W. Metal extraction form ores by heap leaching. Metall. Mater. Trans., 28B (1997) 529–545. Batarseh, K.I. and Stiller, A.H. Modelling the role of bacteria in leaching of lowgrade ores. AIChE J. 40 (1994) 1741–1756. Bear, J. Dynamics of Fluids in Porous Media, Elsevier, NY, 1972. Bennett, C.R., Cross, M., Croft, T.N., Uhrie, J.L., Green, C.R. and Gebhardt, J.E. A comprehensive copper stockpile leach model: background and model formulation. In Hydrometallurgy 2003, Proceedings of the 5th International Symposium Honoring Professor Ian M. Ritchie, Volume 1: Leaching and Solution Purification (C.A. Young, A.M. Alfantazi, C.G. Anderson, D.B. Dreisinger, B. Harris and A. James, eds.), TMS, Warrendale, PA, 2003, pp 315–328. Bhakta, P and Arthur, B. Heap biooxidation and gold recovery at newmont mining: first year results, J. Metals 54 (2002) 31–34. Bhakta, P. Ammonium thiosulfate heap leaching. In Hydrometallurgy 2003, Proceedings of the 5th International Symposium Honoring Professor Ian M. Ritchie, Volume 1: Leaching and Solution Purification (C.A. Young, A.M.  References  178  Alfantazi, C.G. Anderson, D.B. Dreisinger, B. Harris and A. James, eds.), TMS, Warrendale, PA, 2003, pp 259–267. Bouffard S.C. and Dixon, D.G. Mathematical modeling of pyretic refractory gold ore heap biooxidation: model development and Isothermal column simulations. In Hydrometallurgy 2003, Proceedings of the 5th International Symposium Honoring Professor Ian M. Ritchie, Volume 1: Leaching and Solution Purification (C.A. Young, A.M. Alfantazi, C.G. Anderson, D.B. Dreisinger, B. Harris and A. James, eds.), TMS, Warrendale, PA, 2003, pp 275–288. Bouffard, S.C and Dixon D.G. Investigative study into the hydrodynamics of heap leaching processes, Metall. Mater. Trans. 32B (2001) 763–776. Bouffard, S.C., Understanding the Heap Biooxidation of Sulfidic Refractory Gold Ores. Ph.D. Thesis, The University of British Columbia, 2003. Braun, R.L., Lewis, A.E. and Wadsworth, M.E. In-place leaching of primary sulfide ores: laboratory leaching data and kinetic model. Metall. Trans. 5 (1974) 1717–1726. Brierley, C.L. Bioheap leaching: commercial status, operating considerations and future prospects. In Biomine ’99 and Water Management in Metallurgical Operations, Conference Proceedings, Australian Mineral Foundation, August, 1999. Brooks, R.H. and Corey, A.T. Hydraulic properties of porous media. Hydrology Paper No. 3, Colorado State University, Fort Collins, 1964. Brusseau, M.L. Evaluation of simple methods for estimating contaminant removal by flushing. Ground Water 34 (1996) 19–22. Buyuktas, D. and Wallender, W.W. Numerical simulation of water flow and solute transport to tile drains. J. Irrigation Drainage Engng. 128 (2002) 49–56. Catalan, L.J.J. and Li, M.G. Rinsing copper oxide heaps with sulfuric acid: a pilot plant study. In Tailings and Mine Waste ’98, Balkerma, Rotterdam, 1998, pp 881–893. Coats, K.H. and Smith, B.D. Dead-end pore volume and dispersion in porous media. J. Soc. Petrol. Eng. 4 (1964) 73–84. Corey, A. Mechanics of Heterogeneous Fluids in Porous Media, Water Resources Publications, Fort Collins, CO, 1977. Cross, M., Bennett, C.R., Croft, T.N., McBride, D. and Gebhardt, J.E. Computational modeling of reactive multi-phase flows in porous media: applications to metals extraction and environmental recovery processes. Miner. Engng. 19 (2006) 1098–1108.  References  179  Cross, M., Bennett, C.R., McBride, D., Taylor, D.A. and Gebhardt, J.E. Computational modelling of heap leach processes. In Sohn International Symposium Advanced Processing of Metals and Materials, Volume 3, Thermo and Physicochemical Principles (F. Kongoli and R.G. Reddy, eds.), TMS, Warrendale, PA, 2006, pp 323–337. Dai, X., Chee, K.C., Jeffrey, M. and Breur, P. A comparison of cyanide and thiosulfate leaching for the recovery of gold from a copper containing ore. In Hydrometallurgy 2003, Proceedings of the 5th International Symposium Honoring Professor Ian M. Ritchie, Volume 1: Leaching and Solution Purification (C.A. Young, A.M. Alfantazi, C.G. Anderson, D.B. Dreisinger, B. Harris and A. James, eds.), TMS, Warrendale, PA, 2003, pp 123–136. Dalton, F.N., Henkelrath, W.N., Rawlins, D.S. and Rhoades, J.D. Time domain reflectometry. Simultaneous measurement of soil water and electrical conductivities with a single probe. Science 224 (1984) 989–990. Decker, D.L. and Tyler, S.W. Evaluation of flow and solute transport parameters for heap leach recovery materials. J. Environ. Qual. 28 (1999) 543–555. Dixon, D.G. Heap leach modelling – the current state of the art. In Hydrometallurgy 2003, Proceedings of the 5th International Symposium Honoring Professor Ian M. Ritchie, Volume 1: Leaching and Solution Purification (C.A. Young, A.M. Alfantazi, C.G. Anderson, D.B. Dreisinger, B. Harris and A. James, eds.), TMS, Warrendale, PA, 2003, pp 289–314. Dixon, D.G. Personal communication, 2008. Dixon, D.G. and Petersen, J. Comprehensive modelling study of chalcocite and heap bioleaching. In Copper 2003, Volume VI-2, Hydrometallurgy of Copper (P.A. Riveros, D.G. Dixon, D.B. Dreisinger and J. Menacho, eds.), TMS, Warrendale, PA, pp 493–517. Dixon, D.G. and Hendrix, J.L. A general model for leaching of one or more solid reactants from porous ore particles. Metall. Trans. 24B (1993) 157–169. Dixon, D.G. and Hendrix, J.L. A mathematical model for leaching of one or more solid reactants from porous ore pellets. Metall. Trans. 24B (1993) 1087–1102. Dixon, D.G., Dix, R.B. and Comba, P.G. A mathematical model for rinsing of reagents from spent heaps. In Extraction and Processing for the Treatment and Minimization of Wastes 1994 (J.P. Hager, B.J. Hansen, W.P. Imrie, J.F. Pusateri and V. Ramachandran, eds.), TMS, Warrendale, PA, pp 701–713. Eriksson, N. and Destouni, G. Combined effects of dissolution kinetics, secondary mineral precipitation, preferential flow on copper leaching from mining waste rock. Water Resources Res. 33 (1997) 471–483.  References  180  Fetter, C.W. Contaminant Hydrogeology, 1993, Macmillan, NY. Forrer, I. and Kasteel, R. Longitudinal and lateral dispersion in an unsaturated field soil. Water Resources Res. 35 (1999) 3049–3060. Fredlund, M.D. The Role of Unsaturated Soil Property Functions in the Practice of Unsaturated Soil Mechanics, Ph.D. Thesis, University of Saskatchewan, 1999. Freeze, R.A., and Cherry, J.A. Groundwater, Prentice Hall, Englewood Cliffs, NJ, 1979. Funk, G.A., Harold, M.P. and Ng, K.M. A novel model for reaction in trickle beds with flow maldistribution. Ind. Engng. Chem. Res. 29 (1990) 738–748. Gardner, W.R., Hillel, D. and Benyamini, Y. Post irrigation movement of soil water: I. Redistribution. Water Resources Res. 6 (1970) 851–861. Giese, K., Tiemann, R. Determination of complex permittivity from thin sample time domain reflectometry: improved analysis of the step response waveform. Adv. Mol. Relax. Processes 7 (1975) 45–59. Gupta, C.K. and Mukherjee, T.K. Hydrometallurgy in Extraction Processes, Volume 1, CRC Press, Boca Raton, FL, 1990. Harrington, J.G., Bartlett, R.W. and Prisbrey, K.A. In Hydrometallurgy Fundamentals, Technology and Innovation (J.B. Hiskey and G. Warren, eds.), SME, Littleton, CO, 1993, pp 691–708. Harrington, J.G., Prisbrey, K.A. and Bartlett, R.W. In Biohydrometallurgical Technologies (A.E.. Torma, J.E. Wey and V.L. Lakshmanan, eds.), TMS, 1993, pp 521–530. Heimovaara, T. J., Focke, A. G., Bouten, W., and Verstraten, J. M. Assessing temporal variations in soil water composition with time domain reflectometry. Soil Sci. Soc. Am. J. 59 (1995) 689–698. Hillel, D. Environmental Soil Physics, Academic Press, 1998. Hook, W.R., Livingstone, N.J., Sun, Z.J. and Hook, P.B. Remote diode shorting improves measurement of soil water by time domain reflectometry. Soil Sci. Soc. Am. J. 56 (1992) 1384–1391. Howard, E.V. Chino uses radiationlogging for studying dump leaching processes. Mining Engng. (1968) 70–74. Huang, K. van Genuchten, M.T. and Zhang, R. Exact solutions for onedimensional transport with asymptotic scale-dependent dispersion. Applied Math. Modell. 20 (1996) 298–308.  References  181  Hunt, B. Contaminant source solutions with scale-dependent dispersivities. J. Hydrol. Engng. 3 (1998) 268–275. Hwang, S.I. and Powers, S.E. Using particle-size distribution models to estimate soil hydraulic properties. Soil Sci. Soc. Am. J. 67 (2003) 1103–1112. Inoue, M., Simunek, J., Shiozawa, S. and Hopmans, J.W. Simultaneous estimation of soil hydraulic and solute transport parameters from transient infiltration experiments. Adv. Water Resources 23 (2000) 677–688. Jacobsen R.H. Optimization of dump and heap leaching – theoretical and practical studies. Preprint T-IV-a-3, Joint Meeting, MMIJ-AIME, Tokyo, 1972. Kachanoski, R.G., Pringle, E. and Ward, A. Field measurement of solute travel times using time domain reflectometry. Soil Sci. Soc. Am. J. 56 (1992) 47–52. Kelly, S.F., Selker, J.S. and Green, J.L. Using short soil moisture probes with bandwidth time domain reflectometry instruments. Soil Sci. Soc. Am. J. 59 (1995) 97–102. Lal, R. and Shukla, M.K. Principles of Soil Physics, CRC Press, Boca Raton, FL, 2004. Leahy, M.J., Davidson, M.R. and Schwarz, M.P. A model for heap bioleaching of chalcocite with heat balance: bacterial temperature dependence. Minerals Engng. 18 (2005) 1239–1252. Leahy, M.J., Davidson, M.R. and Schwarz, M.P., A model for heap bioleaching of chalcocite with heat balance: mesophiles and moderate thermophiles. Hydrometallurgy 85 (2006) 24–41. Leahy, M.J., Davidson, M.R. and Schwarz, M.P. A two-dimensional CFD model for heap bioleaching of chalcocite. ANZIAM, 46 (2005) C439–C457. Li, M.G. Decommissioning of sulphuric acid-leached heaps by rinsing. In Tailings and Mine Waste ’95, Balkerma, Rotterdam, 1995, pp 295–304. Lim, P.C., Barbour, S.L. and Fredlund, D.G. The influence of degree of saturation on the coefficient of aqueous diffusion. Can. Geotech. J. 35 (1998) 811–827. Logan, J.D. Solute transport in porous media with scale-dependent dispersion and periodic boundary conditions. J. Hydrology 184 (1996) 261–276. Looney, B.B. and Falta, R.W. (eds.) Vadose Zone Science and Technology Solutions, Volume II, Battelle, 2000.  References  182  Madsen B.W. and Wadsworth, M.E. A mixed kinetics dump leaching model of ore containing a variety of copper sulfide minerals. USBM Report of Investigations No. 8547, 1981. Mallants, D., Vanclooster, M., Toride, N., Vanderborght, J., van Genuchten, M.Th. and Feyen, J. Comparison of three methods to calibrate TDR for monitoring solute movement in undisturbed soil. Soil Sci. Soc. Am. J. 60 (1996) 747–754. Marchbank, A. R., Kirby, E., and Roberto, F. The application of biologically mediated leaching methods – a comprehensive review. Can. Metall. Q. 96 (2003) 143–147. Martens, W., Frost, R.L., Kloprogge, J.T. and Williams, P.A. Roman spectroscopic study of the basic copper sulphates-implications for copper corrosion and “bronze disease”. J. Raman Spectroscopy 34 (2003) 145–151. Mboimpa, M, Aubertin, M. Chapuis, R.P. and Bussiere, B. Practical pedotransfer function for estimating the saturated hydraulic conductivity. Geotech. Geolog. Engng. 20 (2002) 235–259. McBride, D., Cross, M., Croft, T.N., Bennett, C.R.and Gebhardt, J.E. Modelling variably saturated flow in porous media for heap leach analysis. In Computational Analysis in Hydrometallurgy, (D.G. Dixon and M.C. Dry, eds.), 2005, pp 45–59. Mears, D.E., The role of axial dispersion in trickle-flow laboratory reactors, Chem. Engng. Sci. 26 (1971) 1361–1366. Misra, S. and Parker, J.C. Parameter estimation for coupled unsaturated flow and transport, Water Resources Res. 25 (1989) 385–396. Montero, J.P., Munoz, J.F., Abeliuk, R and Vauclin, M. A solute transport model for the acid leaching of copper in soil columns. Soil Sci. Soc. Am. J. 58 (1994) 678–686. Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Res. 12 (1976) 513–522. Muir, D.M. and Aylmore, M.G. Thiosulphate as an alternative to cyanide for gold processing – issues and impediments. Miner. Proc. Extr. Metall. 113 (2004) C2– C12. Munoz, J.F., Rengifo, P. and Vauclin, M. Acid leaching of copper in saturated porous material: parameter identification and experimental validation of a twodimensional transport model. J. Contaminant Hydrology 27 (1997) 1–24.  References  183  Murr, L.E. Schlitt, W.J. and Cathles, L.M. Experimental observations of solution flow in the leaching of copper-bearing wastes. In Interfacing Technologies in Solution Mining (W.J. Schlitt, ed.), SME, NY, 1982, pp 271–290. Nadler, A., Dasberg, S. and Lapid, I. Water content and electrical conductivity determination of layered soil profiles using time domain reflectometry. Soil Sci. Soc. Am. J. 55 (1991) 938–943. Newman, L.L., Herasymuik, G.M., Barbour, S.L., Fredlund, D.G. and Smith, T. The hydrogeology of waste rock dumps and a mechanism of unsaturated preferential flow. In Proceedings of the 4th International Conference on Acid Rock Drainage, Vancouver, BC, 1997, pp 553–565. Nichol, C. Transient Flow and Transport in Unsaturated Heterogeneous Media: Field Experiments in Mine Waste Rock. Ph.D. Thesis, The University of British Columbia, 2002. Nimmo, J. R. Modelling structural influences on soil water retention. Soil Sci. Soc. Am. J. 61 (1997) 712–719. Nützmann, G., Maciejewski, S. and Joswig, K. Estimation of water saturation dependence of dispersion in unsaturated porous media: experimental and modelling analysis. Adv. Water Resources 25 (2002) 565–576. O’Kane, M., Barbour, S.L., and Huang, M.D. A framework for improving the ability to understand and predict performance of heap leach piles. In Copper 99, Volume IV, Hydrometallurgy of Copper (S.K. Young, D.B Dresinger, R.P. Hackl, and D.G Dixon, eds.), TMS, Warrendale, PA, 1999, pp 409–419. Or, D., Jones, S.B., VanShaar, J.R., Humphries, S. and Koberstein, L. WinTDR Soil Analysis Software Users Guide, Utah State University Soil Physics Group Version 6.0, 2003. Orr, S. Enhanced heap leaching, part 1: Insights. Mining Engng. (Sept. 2002) 49–56. Orr, S. and Vesselinov, V. Enhanced heap leaching, part 2: Applications. Mining Engng. (Oct. 2002) 33–38. Pacchepsky, Y.A. and Rawls, W.J. Soil structure and pedotransfer functions. Eur. J. Soil Sci. 54 (2003) 443–451. Pacchepsky, Y.A., Timlin, D. and Varallyay. Artificial neural networks to estimate soil water retention from easily measurable data. Soil Sci. Soc. Am. J. 60 (1996) 727–733. Pang, L. and Hunt, B. Solution and verification of scale-dependent dispersion model. J. Contaminant Hydrology 53 (2001) 21–39.  References  184  Pantelis, G., Ritchie, A.I.M. and Stepanyants, Y.A. A conceptual model for the description of oxidation and transport processes in sulphidic waste rocks dumps. Applied Math. Modell. 26 (2002) 751–770. Parker, J.C. and van Genuchten, M.Th. Flux-averaged and volume-averaged concentrations in continuum approaches to solute transport. Water Resources Res. 20 (1984) 866–872. Patankar, S. Numerical Heat Transfer and Fluid Flow, Taylor & Francis, 1980. Petersen, J. and Dixon, D.G. TeckCominco Z2K Heap Leach Simulation Package, Reference Manual, May/June 2002. Peterson, J., and Dixon, D.G. The dynamics of chalcocite heap bioleaching. In Hydrometallurgy 2003, Proceedings of the 5th International Symposium Honoring Professor Ian M. Ritchie, Volume 1: Leaching and Solution Purification (C.A. Young, A.M. Alfantazi, C.G. Anderson, D.B. Dreisinger, B. Harris and A. James, eds.), TMS, Warrendale, PA, 2003, pp 351–364. Qi, M; Lorenz, M. and Vogelpohl, A. Mathematical solution of the twodimensional dispersion model. Chem. Eng. Technol. 25 (2002) 693–697. Roman, R. J. and Bhappu, R.B. Heap hydrology and the decommissioning of a leach heap. Chapter 64 in Hydrometallurgy Fundamentals, Technology and Innovation (J.B. Hiskey and G. Warren, eds.), SME, Littleton, CO, 1993. Roman, R.J., Brenner, B.R, and Becker, G.W. Diffusion model for heap leaching and its application to scale-up. Trans. SME-AIME 256 (1974) 247–251. Roth, K., Schlin, R., Fluhler, H. and Attinger, W. Calibration of time domain reflectometry for water content measurement using a composite dielectric approach. Water Resources Res. 26 (1990) 2267–2273. Rudolph, D.L., Kachanoski, R.G., Celia, M.A., LeBlanc, D.R. and Stevens, J.H. Infiltration and solute transport experiments in unsaturated sand and gravel, Cape Cod, Masschusetts: experimental design and overview. Water Resources Res. 32 (1996) 519–532. Sanchez-Chacon, A. E. and Lapidus, G. T. Model for heap leaching of gold ores by cyanidation. Hydrometallurgy 44 (1997) 1–20. Schlitt, W.J. and Nicolai, L.F. Non-vertical solution flow and its applications in heap and dump leaching. Miner. Metall. Proc. (Feb. 1987) 1–7. Schwartz, F. W. and Zhang, H. Fundamentals of Groundwater, 2003.  References  185  Songrong, Y., Jiyuan, X., Guanzhou, Q., and Yuehua, H. Research and application of bioleaching and biooxidation technologies in China. Minerals Engng. 15 (2002) 361–363. Spitz, K. and Moreno, J. A Practical Guide to Groundwater and Solute Transport Modeling, J. Wiley and Sons, NY, 1996. Stollenwerk, K.G and Kipp, K.L. ACS Symposium Series: Chemical Modeling of Aqueous Systems II, ACS, Washington, DC, 1990, pp 241–257. Straits Resources Limited-Heap Leach SX-EW http://www.straits.com.au/technology/sx_ew/sxewhl.htm  Technology,  Topp, G.C. Electromagnetic determination of soil water content: measurements in coaxial transmission lines. Water Resources Res. 16 (1980) 574–582. Topp, G.C., Yanuka, M., Zebchuk, W.D. and Zeglin, S. Determination of electrical conductivity using time domain reflectometry. Soil and water experiments in coaxial lines. Water Resources Res. 24 (1988) 945–952. Vallochi, A.J. Use of temporal moment analysis to study reactive solute transport in aggregated porous media, Geoderma 46 (1990) 233–247. Van Genuchten, M.Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (1980) 892–898. Vanclooster, M., Gonzalez, C., Vanderborght, J., Mallants, D., Diels, J., An indirect calibration procedure for using TDR in solute transport studies. In Symposium and Workshop on Time Domain Reflectometry in Environmental, Infrustructure and Mining Aapplications, Northwestern University, Evanston, IL, September 7–9, 1994, pp 215–226. Vanclooster, M., Mallants, D. and Feyen, D.J. Determining the local-scale of solute transport parameters using time domain reflectometry (TDR). J. Hydrology 148 (1993) 93–107. Vanclooster, M., Mallants, D., Vanderborght, J., Diels, J., Van Orshoven, J., and Feyen, J. Monitoring solute transport in a multi-layered sandy lysimeter using time domain reflectometry. Soil Sci. Soc. Am. J. 59 (1995) 337–344. Veras, C.A. and Videla, N.M. A new method for heap leaching testing. In Copper 2003, Volume VI-2, Hydrometallurgy of Copper (P.A. Riveros, D.G. Dixon, D.B. Dreisinger and J. Menacho, eds.), TMS, Warrendale, PA, 2003, pp 147–159. Visser, W.C. Progress in the knowledge about the effect of soil moisture content on plant production. Tech. Bull. No. 5, Inst. Land Water Management, Wageningen, The Netherlands.  References  186  Vogel, T., Gerke, H.H., Zhang, R. and van Genuchten, M.Th. Modeling flow and transport in two-dimensional dual permeability system with spatially variable hydraulic properties. J. Hydrology 238 (2000) 78–89. Ward, A.L., Kachanoski, R.G., and Elrick. D.E. Laboratory measurements of solute transport using time domain reflectometry. Soil Sci. Soc. Am. J. 58 (1994) 1031–1039. Weast, R.C. CRC Handbook of Chemistry and Physics, 57th ed., CRC Press, Boca Raton, FL, 1976. Woolery, R.G., Hansen, D.J., Weber, J.A. and Ramachandran, S. Heap leaching of uranium: a case history. Mining Engng. (March 1978) 285–290. Wosten, J.H.M., Pacchepsky, Y.A. and Rawls, W.J. Pedotransfer functions: bridging the gap between available basic soil data and missing soil hydraulic characteristics. J. Hydrology 251 (2001) 123–150. Yazdani, J., Barbour, L., and Wilson, W. Soil water characteristic curve for mine waste rock containing coarse material. Proceedings of CSCE Annual Conference, London, Ontario, 2000, pp 198–202. Zarate, G. and Kelly, R. The metallurgy of the Mantoverde project. Hydrometallurgy 39 (1995) 307–319. Zarate, G., Trincado, L., Troncoso, U. and Vargas, C. (2003) Process improvements at Mantoverde heap leach-SX-EW plant. In Copper 2003, Volume VI-2, Hydrometallurgy of Copper (P.A. Riveros, D.G. Dixon, D.B. Dreisinger and J. Menacho, eds.), TMS, Warrendale, PA, 2003, pp 811–819. Zegelin, S.J. and White, I. Improved field probes for soil water content and electrical conductivity measurement using time domain reflectometry. Water Resources Res. 25 (1989) 2367–2376. Zheng, C. and Bennett, G.D. Applied Contaminant Transport Modeling, 2nd ed., J. Wiley and Sons, NY, 2002. Zhuang, J., Jin, Y. and Miyazaki, T. Estimating water retention characteristics from soil particle-size distribution using a non-similar media concept. Soil Sci. 166 (2001) 308–321.  References  187  Appendix A DERIVING THE EQUATIONS OF CONTINUITY The schematic diagram of the cylindrical differential element used to formulate the equations of continuity of water flow and solute transport in a 2D axisymmetric configuration is shown in Figure A1.  ∆r r  z  ∆z  z + ∆z  Figure A.1  Schematic of a cylindrical differential element  Appendix A – Deriving the Equations of Continuity  188  A.1  Deriving the Equation of Continuity for Water Flow  The mass rate of water flow into the upper z-normal surface at z: 2πr∆r ρw v w ,z z The mass rate of water flow out of the lower z-normal surface at z + ∆z: 2πr∆r ρw v w ,z z + ∆z The mass rate of water flow into the inner r-normal surface at r: 2πr∆z ρw v w ,r  r  The mass rate of water flow out of the outer r-normal surface at r + ∆r: 2πr∆z ρw v w ,r  r + ∆r  The mass rate of accumulation of water in the element, assuming constant mass density: 2πr∆r∆z ρw  ∂θ ∂t  The net mass rate of production of water in the element by chemical reaction:  2πr∆r∆z Mw sw By conservation of mass of water in the element:  mass rate of water   net mass rate of   mass rate of   mass rate of    −   =   −   accumulati on water production water flow in water flow out         Inserting the mass rates into this equation gives:  Appendix A – Deriving the Equations of Continuity  189    ∂θ 2πr∆r∆z  ρw − Mw sw  ∂t    (  )  (  = −2πr∆r ρw v w ,z z + ∆z − ρw v w ,z z − 2π∆z rρw v w ,r  r + ∆r  − rρw v w ,r  r  )  Dividing by the volume of the element and the mass density gives:  v w ,z z + ∆z − v w ,z z 1 rv w ,r ∂θ Mw sw = − − − ∆z r ∂t ρw  r + ∆r  − rv w ,r  r  ∆r  Finally, taking the limit as the volume of the element approaches zero and recognizing the definition of the first derivative gives the differential equation of continuity:  ∂v ∂θ Mw 1 ∂(rv w ,r ) sw = − w ,z − − = −∇ ⋅ vw r ∂t ρw ∂z ∂r  A.2  Deriving the Equation of Continuity for Solute Transport  The molar rate of solute transport into the upper z-normal surface at z: 2πr∆r n A,z z The molar rate of solute transport out of the lower z-normal surface at z + ∆z: 2πr∆r nA,z z + ∆z The molar rate of solute transport into the inner r-normal surface at r: 2πr∆z nA,r  r  The molar rate of solute transport out of the outer r-normal surface at r + ∆r: 2πr∆z nA,r  Appendix A – Deriving the Equations of Continuity  r + ∆r  190  The molar rate of accumulation of solute in the element: 2πr∆r∆z  ∂(θc A ) ∂t  The net molar rate of production of solute in the element:  2πr∆r∆z sA By conservation of moles of solute in the element:  molar rate of solute   net molar rate of   −     accumulati on   solute production  molar rate of  molar rate of     −   =   solute transport in   solute transport out  Inserting the molar rates into this equation gives:  (  )  (    ∂(θc A ) 2πr∆r∆z  − s A  = −2πr∆r nA,z z + ∆z − nA,z z − 2π∆z rnA,r   ∂t  r + ∆r  − rnA,r  r  )  Dividing by the volume of the element gives:  nA,z z + ∆z − nA,z z 1 rnA,r ∂(θ c A ) − sA = − − ∆z r ∂t  r + ∆r  − rnA,r  r  ∆r  Finally, taking the limit as the volume of the element approaches zero and recognizing the definition of the first derivative gives the differential equation of continuity: ∂n ∂(θc A ) 1 ∂(rnA,r ) − s A = − A, z − = −∇ ⋅ nA ∂t ∂z r ∂r  Appendix A – Deriving the Equations of Continuity  191  Appendix B TIME DOMAIN REFLECTOMETRY B.1  Overview  In order to characterize solute transport, its travel-time density functions at different depths in the heap must be measured. This may be accomplished by measuring the resident and/or flux-averaged concentration of a conservative tracer at different depths of an ore column [Freeze and Cherry, 1979].  The  movement of the tracer through unsaturated field soils was traditionally monitored by soil coring or vacuum solution sampling techniques [Kachanoski et al., 1992; Vanclooster et al., 1995]. Soil sampling gives resident concentration and mass distribution of the tracer at depth. Simple salts such as NaCl, KCl, CaCl2 and NaBr are typically used as inert tracers, and their concentrations are monitored by measurements of electrical conductivity by soil coring and/or vacuum solution sampling. This method is laborious, time-consuming, expensive and destructive, and only one distribution of solute can be obtained at a particular time and location. Other traditional methods use radioisotopes such as 3H, 29  Br and  131  I,  36  Cl which can be monitored using radioactivity detectors such as the  Geiger–Muller detector [Freeze and Cherry, 1979; Kachanoski et al., 1992; Vanclooster et al., 1995] but may be limited in their application due to cost and the radiation hazards involved. In the late 1970’s, Time Domain Reflectometry (TDR) was developed and has been established over the years to give high-quality data for water content and bulk electrical conductivity in a porous medium.  The method measures the  velocity of an electromagnetic pulse propagating along a probe placed in contact with the soil. The pulse velocity is related to the dielectric properties of the soil and hence the water content in the soil [Kachanoski et al., 1992; Vanclooster et al., 1993, 1995].  Appendix B – Time Domain Reflectometry  192  The attenuation of the electromagnetic wave, which is propagated along TDR rods in the soil, is related to the bulk electrical conductivity of the soil. Based on this assumption and the knowledge of the relationship between the soil bulk electrical conductivity and solute resident concentration, solute transport can be monitored by TDR [Freeze and Cherry, 1979; Kachanoski et al., 1992; Vanclooster et al., 1995]. The TDR method is scale-independent and leaves the surrounding soil relatively undisturbed, allowing continuous and automated sampling at different depths in the soil. As such TDR may be adopted to yield the data required for the evaluation of different solute transport models in heterogeneous soils. Another advantage of TDR is the elimination of the cost and time to analyze hundreds of samples for the tracer, as is the case with soil coring and vacuum solution sampling techniques.  Kachanoski et al. [1992],  Vanclooster et al. [1993, 1995], Mallants et al. [1996], Decker and Tyler [1999] and Nichol [2002] have all used TDR probes in laboratory columns.  B.2  A Brief History of the Evolution of the TDR Technique  TDR was initially introduced by Topp et al. [1980] to measure volumetric water content in soils. In 1984, Dalton et al. [1984] showed that TDR probes used for water content measurement could also measure electrical conductivity in the same volume of soil. Zegelin and White [1989] and Nadler et al. [1991] have demonstrated that using triple-wire probes, without a balancing transformer, yielded more accurate conductivity measurements than two-wired probes with a balancing transformer. Nadler et al. [1991] also gave an improved technique for interpreting TDR waveforms for estimating electric conductivity, including a calibration procedure to determine the cell constant of TDR probes. In 1992, Kachanoski et al. [1992] obtained electrical conductivity data from TDR to determine solute travel times in a column of soil.  Heimovaara et al. [1995]  proposed a modification to account explicitly for the impedance of the sample by taking the difference of the total impedance and the impedance due to the combined series resistance in the cable, connectors and the cable tester.  Appendix B – Time Domain Reflectometry  193  B.3  TDR Theory  When an electromagnetic wave is propagated along a TDR probe inserted in a soil sample, the propagation velocity is determined by the time response of the system to a pulse generated by the cable tester, and the soil dielectric constant [Or et al., 2003]: vp =  2Lp tp  2  t c c  εb =   =  p   2L  v   p  p where  (B1)  2  (B2)  vp = propagation velocity (m/s) Lp = probe length (m) tp = pulse reflection time (s) c = speed of light in a vacuum (3 x 108 m/s)  and  εb = soil dielectric constant (–)  From Equation 2, the dielectric constant of a medium can be defined as the square of the ratio of propagation velocity in vacuum to that in the medium. The soil dielectric constant depends on the water content of the soil [Vanclooster, 1993]. The dielectric constants of the components of an unsaturated soil (water, air and soil) are given in Table B.1. The large disparity of the values of the dielectric constants makes the method relatively insensitive to soil composition and texture, and thus, a good method for soil water measurement [Or et al., 2003].  B.3.1 Volumetric Water Content Two basic approaches have been used to establish the relationship between the soil dielectric constant and the volumetric water content. The first approach, by Topp et al. [1980], is empirical, and involves fitting a third order polynomial to the  Appendix B – Time Domain Reflectometry  194  observed data of the dielectric constants and water content for a number of soils (mainly with water contents limited to the range θ < 0.5 in mineral soils) [Or et al., 2003]. As such, it fails to adequately describe the relationship at water contents above 0.5 and for soils high in organic matter [Or et al., 2003]. The empirical relationship derived by Topp et al. [1980] is given as follows: θ = − 5.3 × 10 −2 + 2.92 × 10 −2 εb − 5.5 × 10 −4 εb2 + 4.3 × 10 −6 εb3 Table B.1  (B3)  Dielectric constants of soil materials (20–25ºC) [Or et al., 2003] Material  εb  Water 80.4 – 78.5 Air 1.0 Fused Quartz (SiO2) 3.78 Sandy Soil (dry) 2.55 Loamy Soil (dry) 2.51 The second formulation (called the mixing model approach), which is physically based on the phases of the components of the unsaturated porous medium, was proposed by Roth et al. [1990]. It involves the dielectric constants of the threephase system (i.e. soil, water and air) and the porosity of the soil [Or et al., 2003]: θ=  where  εbβ − (1 − η )εsβ − ηεaβ εwβ − εaβ  (B4)  εs = dielectric constant of solid phase εw = dielectric constant of aqueous phase εa = dielectric constant of gaseous phase εb = bulk dielectric constant of the three-phase system η = porosity of the soil β = empirical parameter related to soil geometry (–1 ≤ β ≤ 1), usually taken as 0.5 for an isotropic medium [Or et al., 2003]  Appendix B – Time Domain Reflectometry  195  Thus the porosity of the soil must be known or estimated when using this approach. In the range of water contents of less than 0.5, the mixed model and the empirical model produce similar calibration curves [Or et al., 2003].  B.3.2 Electrical Conductivity As surface waves (known as transverse electromagnetic, or TEM, waves) propagate along TDR probes buried in a bed of soil, the signal energy is attenuated in proportion to the bulk electrical conductivity along the travel path [Vanclooster; 1993; Or et al., 2003]. This proportional reduction in signal voltage is the basis for the simultaneous measurement of the bulk soil electrical conductivity and volumetric water content using TDR. It is also assumed that the bulk electrical conductivity of the soil is proportional to the impedance load of the ongoing signal and the signal load occurring after multiple reflections, thus [Vanclooster; 1993; Or et al., 2003]: EC = α  where  Z0 K = ZT ZT  (B5)  Z0 = impedance load before the wave enters the soil (characteristic probe impedance) (Ω) ZT = ‘lumped’ impedance load of the soil-probe system after multiple reflections (or the total resistance of the cable tester, coaxial cable, and TDR probe inserted in a sample) (Ω) α = calibration constant  and  K = probe constant given by [Vanclooster, 1993; Or et al., 2003]:  K=  Z0 120πLp  Heimovaara et al. [1995] modified Equation 5, which was used by Dalton et al. [1984], Topp et al. [1988], and Nadler et al. [1991], as follows [Vanclooster, 1994; Heimovaara et al., 1995]:  Appendix B – Time Domain Reflectometry  196  EC =  where  K ZT − Zc  (B6)  Zc = impedance due the combined series resistance in the cable, connectors and cable tester (typically 50 Ω)  Zc increases with increasing cable length. A direct calibration of the TDR probe in different solutions of known EC could be used to determine Zc and K by simple linear regression.  This modification is necessary, especially for use of long  cables in automated TDR systems [Heimovaara et al., 1995]. A typical TDR trace (reflectogram) resulting from a step voltage pulse propagated along a probe inserted into a dielectric medium is given in Figure 1.  The  reflectogram, reproduced by means of a high-frequency oscilloscope, is analyzed to estimate the dielectric properties of the material along the probe [Vanclooster, 1993].  Figure B.1  A schematic diagram of TDR wave showing points of measurement for water content and load resistance [Ward et al., 1994]  Appendix B – Time Domain Reflectometry  197  The impedance ZT is the total opposition to the flow of electrical energy in the transmission line. It comprises the opposition to alternating current (AC), called reactance, and opposition to direct current (DC), called resistance, RT [Ward et al., 1994]. At zero or low frequencies, ZT is only dependent on RT. Since the lowest frequency of the TDR is associated with the longest travel time, ZT should be measured at the longest travel time [Ward et al., 1994]. Therefore, the total resistance ZT at the end of the transmission line (typically a coaxial cable), which is equal to the impedance at low frequencies, can be obtained at very long times, and is given thus [Heimovaara et al., 1995; Or et al., 2003]: ZT = Zc  where  (1 + ρ∞ ) (1 − ρ∞ )  (B7)  ρ∞ = voltage reflection coefficient at long times (i.e., the point of the waveform where the multiple reflection have ceased and the waveform has reached a stable level) [Heimovaara et al., 1995; Or et al., 2003].  The reflection coefficient of this stable level is defined as follows [Heimovaara et al., 1995; Or et al., 2003]: ρ∞ =  where  V∞ − VA VA  (B8)  VA = amplitude of the incident voltage wave generated by the cable tester V∞ = final voltage amplitude in the transmission line after all multiple reflections have ceased  These are illustrated in Figure B.1. Two main approaches exist for determining EC from TDR measurements, namely, the Cell Constant Approach and the Thin-Section Approach.  Appendix B – Time Domain Reflectometry  198  In the Cell Constant Approach, ρ∞ and ZT can be determined by Equations 7 and 8, respectively [Heimovaara et al., 1995; Or et al., 2003].  Equation 5 or its  modified version, Equation 6, can then be used to calculate EC. The Thin-Section Approach was originally proposed by Giese and Tiemann [1975], and has been shown to be effective in quantifying the bulk EC of the soil [Or et al., 2003]:  cε  Z  2V  EC =  0  0  A − 1 L  Z   p  c  V∞ where  (B9)  ε0 = permittivity of free space (8.9×10–12 F/m)  The unique ability of TDR to simultaneously measure θ and EC in the same soil volume has led to its use in monitoring and quantifying the behavior of salts and ionic solutes in the soil environment. The electrical conductivity of the bulk soil depends on the electric conductivity of the soil solution and contributions of the ionic species on the solid particle exchange sites, and the tortuosity of the flow path of the electrical current [Or et al., 2003]. The latter attribute is a function of the liquid-solid-air phase geometry, and thus changes in response to θ and alteration of the particle or pore arrangement. Under constant water content (static or steady flow) conditions, electrical conductivity of the soil water is strongly correlated to the concentration of ionic species in the soil solution. These features of TDR have been exploited by many soil scientists (Kachanoski et al. [1992] and Decker and Tyler [1999], inter alia) to estimate the transport of ionic solutes (serving as inert tracers) through columns of soils or ores. The time history of soil bulk electrical conductivity, EC(t), as estimated from the time history of impedance of TDR waveforms at infinity, Z(t), is approximately linearly related to electrical conductivity of the soil water, ECw(t). The latter, in turn, is linearly related to resident concentration, Cr(t), of a particular solute or an inert tracer in the soil solution [Ward et al., 1994]:  Appendix B – Time Domain Reflectometry  199  EC (t ) =  K = a ⋅ ECw (t ) + b = ACr (t ) + B Z ( t ) − Zc  (B10)  If the soil initially has a low electrical conductivity with solute concentration Cri (corresponding to an attenuation of 1/Zi) and a solution of high EC and concentration Cr0 (corresponding to an attenuation of 1/Z0) of an inert tracer is added to the soil to displace the initial soil water, then the relative concentration of the solute added to the soil is given as follows [Kachanoski et al., 1992]: 1 1 − Cr (t ) − Cri EC (t ) − ECi Z (t ) − Zc Zi − Zc = = 1 1 Cr 0 − Cri EC0 − ECi − Z0 − Zc Z i − Zc  (B11)  If Cri is sufficiently small and can be considered negligible then Equation 11 reduces to Equation 12 [Vanclooster et al., 1993, 1994]: 1 1 − Cr (t ) EC (t ) − ECi Z (t ) − Zc Zi − Zc = = 1 1 Cr 0 EC0 − ECi − Z 0 − Zc Z i − Z c  (B12)  The empirically estimated value of Cr(t)/Cr0 from the impedance at infinity can be used to estimate solute transport parameters, provided that 1/Z0 is known. The value of 1/Z0 corresponds to the attenuation at infinity along the horizontal axis of the TDR trace, when the soil water in the TDR region is completely saturated with solution of concentration Cr0. The value of 1/Z0 reflects the reduction of electrical conductivity in the porous medium compared to the electrical conductivity of water, ECw, an effect induced by the tortuosity of the porous medium [Ward et al., 1994].  B.4  Orientation of Probe Installation  TDR probes can be installed vertically (in the direction of solution flow) or horizontally (in a plane normal to the direction of flow). There are merits and  Appendix B – Time Domain Reflectometry  200  demerits of both orientations, and the preference for one over the other will depend on the experimental setup and the type of data required.  B.4.1 Vertical Installation When TDR probes are installed vertically, the indirect calibration reduces to equating the specific mass of tracer applied, the difference in EC before and after the solute application [Ward et al., 1994]. Using this approach, the calibration coefficients drop out if simple transport models are considered. Kachanoski et al. [1992] have used vertically-installed TDR probes to estimate the mass flux leaving the nearly cylindrical volume around the TDR rods during steady-state field and column breakthrough experiments. The limitations of vertically installed probes include the assumption that all solute applied will immediately be detected by the TDR probe, ignoring any possible horizontal bypass of solute in a transport experiment [Ward et al., 1994]. Furthermore, analysis of the trace for moisture content in multi-layered soils or deep in a homogeneous soil are difficult and at risk of error, due to multiple reflections or incomplete attenuation of the electromagnetic wave [Ward et al., 1994]. To overcome these problems and to obtain solute resident concentration and moisture content data representative for each measurement depth, horizontally installed TDR probes are used, although the indirect calibration of the horizontal installation is more involved.  B.4.2 Horizontal Installation When the probes are installed horizontally, one can use a step- (continuous-) type boundary condition at the soil surface, during solute displacement experiments.  Thus, at infinite time, all solute-free water will be completely  displaced and the calibration constant 1/Z0 [Ward et al., 1994] can be read directly from the TDR screen. The limitations of this approach include the huge amounts of tracer solution required to displace all soil water, and physical nonequilibrium effects which may increase the time required to distribute the solute uniformly in the soil, and which make it practically impossible to continue application of the solute, especially in field experiments [Ward et al., 1994]. Appendix B – Time Domain Reflectometry  201  The pulse boundary condition is an alternative that can be applied with horizontally installed TDR probes.  More accurate estimates of transport  parameters have been obtained by this approach [Vallochi, 1990].  In this  approach, the calibration constant is obtained by equating the time-integrated TDR signal to the total solute mass applied at the soil surface [Vanclooster et al., 1993, Ward et al., 1994].  The limitations in adopting this approach include  assumption of complete solute mass recovery, and no differences between solute resident and flux concentrations. These restrictions make it difficult to assess solute transport parameters, especially when the soil is structured and contains micro-scale heterogeneity [Vanclooster et al., 1994]. Vanclooster et al. [1993, 1995], Mallants et al. [1996], Decker and Tyler [1999] and Nichol [2002] used horizontally-installed TDR rods in soil or ore columns.  B.5  Limitations of TDR  Limitations of the TDR method include the relatively high cost of the equipment, potential limited applicability under highly saline conditions due to excessive signal attenuation, and the fact that soils with high bound water and high organic matter contents require soil-specific calibrations [Kachanoski et al., 1992; Owino et al., 2000]. High sample electrical conductivity (obtained with solutions of high salt concentrations) reduces the quality of the TDR waveform through loss of signal amplitude, which leads to eventual failure of conventional TDR methods. Owino et al. [2000] observed that soil solutions with electrical conductivities higher than 9.9 dS/m represented the practical limit for TDR. Two main methods are used to surmount this problem, namely, waveform differencing by remote diode shorting, and covering the probe with a highly resistant coating [Hook et al., 1992; Kelly et al., 1995; Nichol, 2002].  Appendix B – Time Domain Reflectometry  202  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0067786/manifest

Comment

Related Items