Heap Leaching: Modelling and Forecasting Using CFD Technology

Heap leach operations typically employ some form of modelling and forecasting tools to predict cash flow margins and project viability. However, these vary from simple spreadsheets to phenomenological models, with more complex models not commonly employed as they require the greatest amount of time and effort. Yet, accurate production modelling and forecasting are essential for managing production and potentially critical for successful operation of a complex heap, time and effort spent in setting up modelling tools initially may increase profitability in the long term. A brief overview of various modelling approaches is presented, but this paper focuses on the capabilities of a computational fluid dynamics (CFD) model. Advances in computational capability allow for complex CFD models, coupled with leach kinetic models, to be applied to complex ore bodies. In this paper a comprehensive hydrodynamic CFD model is described and applied to chalcopyrite dissolution under heap operating conditions. The model is parameterized against experimental data and validated against a range of experimental leach tests under different thermal conditions. A three-dimensional ‘virtual’ heap, under fluctuating meteorological conditions, is simulated. Continuous and intermittent irrigation is investigated, showing copper recovery per unit volume of applied leach solution to be slightly increased for pulse irrigation.


Introduction
World ore deposits of copper are becoming lower in head grade with copper sulphides, particularly chalcopyrite being the predominant mineralogy [1]. Copper is typically recovered from these low-grade ore by means of heap leaching. Heap leaching involves trickling a leach solution though a bed of ore to extract target metals. The bed of ore (or heap) is built, often with crushed ore, and engineered to achieve percolation of the leaching solution to all regions. The leach solution interacts with the solid, where the lixiviant dissolves the target minerals from the solid into the solution, and the dissolved species are transported out of the heap for subsequent recovery. Bartlett [2] provides an overview of the process and Petersen [3] gives a review of the heap leach technology and emerging and potential future applications. Ghorbani et al. [4] reviews the current state, challenges and future of heap leach technology. The method of heap leaching is a low-cost technique that can often achieve high recoveries from oxide ores, typically greater than 70% depending upon the kinetics of reaction. However, the world copper resources are mainly chalcopyrite, and achieving high recovery from these copper sulphide deposits has been challenging due to the complexity of the leach kinetics leading to low leaching efficiency.
Over the past decades, there has been increasing interest in the modelling of the heap leach process. Typically, most heap leach operations will employ some form of modelling and forecasting to predict cash flow margins and project viability. However, there is considerable variation in modelling and forecasting tools from simple spreadsheets to advanced phenomenological models. The variety of different modelling approaches varies in complexity, range of input data, accuracy and time to setup and utilise. The more complex phenomenological models are not commonly employed forecasting tools, as they require the greatest amount of time and effort. Yet, accurate production modelling and forecasting are essential for managing production and potentially critical for successful operation of a complex heap leach site. Some heap leach operations operate for years to realise their eventual expected metal recovery, and the time and effort spent in setting up modelling tools initially may indeed contribute to profitability in the long term. Heap leaching of low grade complex copper sulphide ore bodies must contend with slow leach kinetics at ambient conditions, particularly for chalcopyrite and other copper-containing complex sulphide minerals, thus the control of leach reaction chemistry becomes critical in increasing the efficiency of recovery. The complex mineralogy of copper sulphide ore deposits, which often includes gangue chemistry and complex leach behaviour with many competing and interacting reactions, is challenging to predict and capture in any numerical model. In addition to the complex hydrothermal chemical interactions, the heap has strong structural heterogeneity, with localised pockets of high and low permeability regions, which often lead to preferential flow channels. Although it is acknowledged that the solution flow characteristics have a significant effect on the leach kinetics of the heap, this is often ignored in many modelling tools. Capturing the dynamics of a large heap leach operation, with width and depth of hundreds of meters and time scales of months and years is very challenging due to the complex interactions between the ore particles, fluid, thermal, bacterial and oxygen activity. Indeed, modelling the complexity of flow through the ore matrix is not a simple task due to variability in macro-pores, channelling of solution and preferential flow effects [5]. These effects are generally ignored in models but in reality, they can give rise to solute transport problems and limit reaction kinetics [6][7][8] Techniques based on soil physics to represent the preferential flow, such as dual continuum, dual or random porosity (or permeability), may be included in models.
To enable useful predictions of the heap leach process many different simulation approaches have been developed from simple spreadsheets to complex phenomenological models. There are considerable differences in the various heap leach models, and the modelling method employed to predict the heap leach operation very much depends on functionality requirements, together with cost, accuracy and time constraints. Many modelling approaches utilize laboratory data with some form of scale up to the commercial operation [8]. Marsden and Botz [9] discuss the advantages and disadvantages of different modelling approaches and classify the various models into functionality levels; a summary is given in Table 1. Marsden and Botz [9] detail the benefits of the various levels of modelling and conclude that, for general purpose production forecasts, models of level 3 and 4 provide the heap leach operator with reasonable estimates. However [9] also concluded that 'models of level 5 and 6 are particularly useful for providing detailed heap design, diagnostic and optimisation analysis's and for running what-if scenarios to investigate likely heap performance responses to input variable changes'. Indeed, providing a true representative 'virtual heap' which dynamically responds to operational and environmental changes requires all the complex interacting physics of the gas-solid-liquid-thermal heap matrix to be accounted for, thus arguably making the additional time and cost of a phenomenological model, setup and calibration, profitable for any long-term heap operation.
The phenomenological numerical modelling of the heap and bioleach process has been the focus of many investigations over recent years. Fluid flow models to capture the hydrodynamic behaviour of the heap have been presented by [10][11][12][13][14][15][16] amongst others. Computational fluid dynamics (CFD) technologies have enabled more complex multi-phase transport in the modelling of the heap leaching process [17][18][19][20][21][22] with much of this work focusing on capturing the kinetics in one-dimensional columns or two-dimensional slices. Mostaghimi et al. [23,24] applied CFD technology to capturing flow in a three-dimensional heap. McBride et al. [25][26][27] coupled flow-solid-gas and chemical interactions within a CFD model and applied the model to a three-dimensional heterogeneous heap with variable permeability, saturated conditions and solution traveling through preferential pathways. Local variation in permeability and torturous liquid flow paths through the micro-macro pore network has been captured by utilizing a dual pore-system [28] and by incorporating local voidage heterogeneities into the particle size distribution of the ore [29]. The diffusion of the dissolved mineral through the intra-particle network is often modelled by the shrinking core method, where the distribution of mineral within the ore particles is assumed to be homogeneous per particle size and the reaction takes place within a thin reaction zone that proceeds inwards. There has been much discussion on the validity of the shrinking core method [30], and although a simple shrinking core approximation does not account for the spatial and size distribution of the mineral within a particle, it has been shown to be a very useful tool to analyse commercial heap data [31]. Employing some of the same assumptions as the shrinking core model (such as pseudo-steady state and linear surface kinetics) the particle shape, spatial and size distribution of the mineral grains have been captured by X-ray tomography images to improve particle characterization and capture of the leach kinetics at the particle scale [30][31][32]. Others employ a diffusion-reaction equation [10,13]. The particle scale micro-model is often coupled to a fluid macro-model. The shrinking core model is normally coupled to the macro model by an explicit formulation as the implicit formulation is computationally expensive [17,19,25]. More recently Ferrier et al. [33] replaced the diffusion-reaction equation with a product of explicit functions, which significantly reduced computational time whilst retaining accuracy. In the micro-scale model, the mineral dissolution rate is dependent upon the current concentrations in the local fluid, where the macro models (fluid matrix, air matrix and thermal conditions) provide the inputs to the particle dissolution models. Copper sulphide ores have complex minerology and complex leaching behaviour, including oxidation, bacteria effects, gangue and precipitation. Watling [34] gives a review of the bioleaching of sulphide minerals, focusing on the characteristics of chalcopyrite leaching. Chalcopyrite leaching has been investigated by many [35][36][37][38][39][40] among others. Empirical rate laws for the dissolution of chalcopyrite were proposed by Kimball et al. [41] under a number of different conditions, oxidant concentrations and ambient temperatures. Li et al. [42] provided a review of the effect of oxidants. The slow reaction kinetics of chalcopyrite has been investigated by [39,40,[43][44][45][46][47][48]. Typically, the slow reaction rate is due to the formation of surface product layers and passivation and in general an increase in temperature will increase the rate of reaction. If the leach kinetics is limited by particle diffusion, the effect of temperature will be weak, but if the leach kinetics is chemically limited, the temperature effect will be significant [49]. Another important factor in achieving good metal recovery is the level of microbial activity in acid ferric sulphate leaching [50,51]. Chalcopyrite bioleaching has lattice energy much higher than secondary copper sulphide ores, with reaction rates strongly related to temperature. Reasonable rates of extraction are typically only achieved with temperatures above 50 • C. The bacterial oxidation of ferrous ions from pyrite to ferric ions generates heat exothermically within the heap creating high internal temperatures. The manipulation of aeration and irrigation rates can enable higher temperatures to be retained and allow heat to be optimized [52]. Panda et al. [53] carried out a pilot study of bioleaching of low grade chalcopyrite over 383 days and found that the frequency of solution transfer and acid concentration affected copper extraction but aeration and inoculum addition were less significant. Weather conditions such as heavy rain also affected the heap performance. The study also found that rest periods in the initial wetting period had a profound effect on the heap performance but only marginally thereafter. Continuous irrigation is typically applied in heap leaching but intermittent or pulse irrigation has been shown to improve copper recovery [38]. The method of combining irrigation and resting pulses is considered to decrease channelling effects due to solution draining from pores during the rest cycle, which allows the pore capillary forces to draw the next pulse of liquid into the mass of the ore. In a study of an engineered heap of low-grade chalcopyrite ore, lower irrigation rate or pulse leaching was recommended to enhance recovery [54]. When non-uniform percolation of leach solution through the heap occurs, local pockets of high metal inventory can develop. Rucker et al. [55] successfully investigated employing subsurface leaching, where solution is injected into dry regions, to enhance the performance of the heap.
Different strategies for heap bioleaching of chalcopyrite are studied by [56]. Panda et al. [57] provide a review on chalcopyrite leaching with particular emphasis on heap bioleaching. A brief introduction to the modelling of bioleach processes is given by Petersen [58] and Gebhardt et al. [59] give an overview of the leach kinetics of chalcopyrite dissolution in heap conditions using a CFD model. This contribution presents the coupling of the macro fluid-gas-thermal flow CFD algorithms with the micro-leach kinetics of a chalcopyrite leach. A shrinking core kinetic model for chalcopyrite is presented and implemented within the comprehensive multi-physics CFD framework, and the model utilizes temperature dependency based on the electrical conductivity of sulphur. The model is parametrized against experimental leach data and validated against a number of column leach tests at various temperatures. The calibrated model is then applied to a three-dimensional stockpile of similar material under fluctuating meteorological data, such as precipitation, evaporation and notably temperature. Chalcopyrite reaction, iron oxidation and internal temperatures are monitored for continuous and intermittent application rates.

Computational Model
The chalcopyrite kinetic model has been incorporated within an in-house multi-physics CFD framework. The model structure provides a three-dimensional finite volume unstructured mesh modular framework [60]. The governing equations for the macro variables, i.e., gas, fluid and heat, are discretised over a mesh of hexahedral elements representing the heap or solution domain, and implicitly solved at the centre of each finite volume, i.e., the mesh cell. The micro kinetic models are explicitly solved at a local level (mesh cell) as shown in Figure 1. The micro models representing local reaction kinetics are coupled with the macro models through sink and source terms as shown in Figure 2. Each mesh element contains three phases, a static solid phase of the ore body, a liquid phase of solution flowing predominately downwards, and a gas phase (of primarily air) flowing predominantly upwards. Within each mesh element, the chemical reactions are accounted for, both for each reaction sequence and also within each particle size fraction simultaneously. A schematic of the reaction stages occurring within each element for a copper sulphide leach is shown in Figure 2.
are discretised over a mesh of hexahedral elements representing the heap or solution domain, and implicitly solved at the centre of each finite volume, i.e., the mesh cell. The micro kinetic models are explicitly solved at a local level (mesh cell) as shown in Figure 1. The micro models representing local reaction kinetics are coupled with the macro models through sink and source terms as shown in Figure 2. Each mesh element contains three phases, a static solid phase of the ore body, a liquid phase of solution flowing predominately downwards, and a gas phase (of primarily air) flowing predominantly upwards. Within each mesh element, the chemical reactions are accounted for, both for each reaction sequence and also within each particle size fraction simultaneously. A schematic of the reaction stages occurring within each element for a copper sulphide leach is shown in Figure 2.   The chalcopyrite kinetic model has been incorporated within an in-house multi-physics CFD framework. The model structure provides a three-dimensional finite volume unstructured mesh modular framework [60]. The governing equations for the macro variables, i.e., gas, fluid and heat, are discretised over a mesh of hexahedral elements representing the heap or solution domain, and implicitly solved at the centre of each finite volume, i.e., the mesh cell. The micro kinetic models are explicitly solved at a local level (mesh cell) as shown in Figure 1. The micro models representing local reaction kinetics are coupled with the macro models through sink and source terms as shown in Figure 2. Each mesh element contains three phases, a static solid phase of the ore body, a liquid phase of solution flowing predominately downwards, and a gas phase (of primarily air) flowing predominantly upwards. Within each mesh element, the chemical reactions are accounted for, both for each reaction sequence and also within each particle size fraction simultaneously. A schematic of the reaction stages occurring within each element for a copper sulphide leach is shown in Figure 2.

Liquid-Gas-Thermal Transport
The computational algorithm for the solution of the variably saturated flow through the porous matrix is based on the mixed form of the classical Richards' equation. A detailed description of the algorithm is provided by [61]. In this formulation, the liquid-gas interaction is one-way, as it is assumed that the gas flow does not significantly affect the fluid but the fluid flow can influence the movement of the gas. Thermal gradients can influence both the liquid and gas flow, with buoyancy forces having a large impact on gas flow. Thermal effects on surface tension and liquid viscosity are accounted for within the pressure-moisture-conductivity relationships. Liquid and gas solute is transported throughout the domain.
The governing equation solved for variably saturated liquid flow is the mixed form of the Richards' Equation (1). It is written in terms of two unknown variables, moisture content (θ) and the pressure head (h), where K is the hydraulic conductivity.
z is the direction of gravity, t is time and S h is a sink/source term usually representing the boundary conditions or internal preferential flow paths that are represented as sink terms [28]. Equation (1) requires that pressure-moisture-conductivity relationships be defined to complete the model. The unsaturated soil hydraulic properties θ(h) and K(h) are in general highly nonlinear functions of the pressure head, can be temperature (T) dependent, and are represented in this work by the amended van Genuchten relationship as described by [16]. The saturated hydraulic conductivity is dependent upon the intrinsic permeability of the ore and the fluid properties that are in turn dependent upon temperature. The resulting equations are non-linear and require an iterative solution procedure. The governing equation solved for the solute transport of multiple solutes in the ore matrix is, where C i is the concentration of species i in the solution phase, q is the Darcy flux (q = −K(h)∇H, where H is the total hydraulic head, H = h + z), jk are (x,y,z) directions and S i the production or consumption of species i. The dispersion coefficient, D i,jk, is dependent upon the velocity components and longitudinal and transverse dispersivities. Assuming local thermal equilibrium, a single heat balance equation is solved for the fluid-solid-air matrix, Equation (3). The heat balance is predominately driven by the heat of reaction, boundary conditions, liquid and gas flow, evaporation and condensation, which are incorporated into the heat source term, where c p is the specific heat capacity, λ p is the thermal conductivity and ρ p is the matrix density, which are defined as in Equations (4) to (6), where the fluid, solid and air are represented by the subscript f, s, a respectively. The volumetrically averaged thermal conductivity is based on the inverse of the conductivities of the phases.
where m j is the mass of phase j, m p is the total mass and c j is the specific heat of phase j.
where λ j is the thermal conductivity of phase j.
where VF j is the volume fraction of phase j. The equation solved for the solution of gas flow assumes incompressible gas with a Boussinesq source term. The approximation is valid as local temperature differences are small and therefore changes in density are very small.
where k in is the intrinsic permeability of the porous media, k g (S) is the unsaturated permeability of the gas at liquid saturation, ε g is the volume fraction of the gas phase, µ g is the gas viscosity, and S g is the gas source term accounting for natural convection due to thermal gradients, volume changes due to saturation and other sources, such as boundary conditions.

Solid-Liquid Dissolution
A key aspect of the heap leach process is the dissolution of metals from the minerals of interest. In this work, the rate of dissolution is calculated using a shrinking core model to combine the effects of diffusion through the rock matrix and the kinetic reaction rates. The shrinking core equation is solved for each mineral in each characteristic particle size at each point in time. In this way, the ore material can be characterized by defining a particle size distribution, with a single heap grade or each particle size class having its own initial head grade. The equations solved are given below and [16] gives details of the shrinking core algorithm and coupling. The equation used to calculate the rate of dissolution of a given mineral is solved explicitly for each particle size at the start of each time step.
where r o is the initial particle radius, r m is the current mineral radius, A m comes from the kinetic rate equation for the current mineral, c o is the concentration of available reactant at particle surface, R is the gas constant, T is the temperature in Kelvin, D eff is the effective rock diffusion coefficient, ε p is the rock voidage, ρ ore is the ore density, M i is the molecular weight of the mineral and x i is the mass fraction of the mineral.
The value of A m comes from the general expression for the kinetic rate equations, such as those produced by [62], and the general form for which is where β is the fraction of mineral reacted and A,B are functions of the individual kinetic rate equation. The fraction of a particle reacted, β i for mineral i, is evaluated as, where the superscript 0 indicates the value at the previous time step. If there is insufficient reactant available to meet the demands of all the reactions, the reactant is shared out proportionally to all reactions. The particle fraction reacted, β i , and characteristic radius, (r m ) i , are then scaled back accordingly. The production and consumption of species i enters the transport Equation (2) as a source term. This approach allows minerals to react at different rates in an individual particle size, and although it requires the assumption that the minerals can be treated relatively independently, this is not unreasonable given the low concentration of reactive minerals present. This approach allows each mineral in each particle size fraction to be modelled using a single characteristic radius indicating how much has reacted. It can easily be related to experimental analysis of ore which is commonly given as mineral content by size classification, making validation easier. It also easily allows for different minerals to dominate the reactions at different stages of the leach cycle. This in turn allows the model to deal with multiple particle sizes and multiple minerals over large meshes without excessive memory usage in the overall CFD model framework. The mineral rate constants, A and B are unrelated to particle size and therefore allow the model to scale from small to large particle size distributions (e.g., from experimental column to 3D heap lifts).

Gas-Liquid Oxygen Mass Transfer
Oxygen transfer between the liquid and gas phase is driven by temperature according to Henry's law. Due to the large solution time steps and large surface area between the liquid and gas phases, it can be assumed that the transport rate is effectively instantaneous and an equilibrium state is obtained. Oxygen is partitioned between the liquid and gas phases by solving; where O g is the mass fraction of oxygen (O 2 ) in the gas phase, O f is the molar concentration of oxygen in the fluid phase, M O 2 is the molar weight of oxygen, M N 2 is the molar weight of nitrogen and K H is the Henry's law constant at a given temperature Henry's law constant, K H , is a function of temperature based on dissolved oxygen solubility against temperature data under atmospheric conditions. The partial pressure of oxygen can also be Minerals 2018, 8, 9 8 of 20 dependent on the level of salts in solution. An equilibrium state is solved using a Newton-Raphson scheme at the end of each time step.

Ferrous Oxidation
The ferric regeneration process in sulphide mineral bioleaching depends upon ferrous oxidising to ferric catalysed by bacterial action. Different classes of bacteria exist which tolerate highly acidic environments and high temperature. The model accounts for different classes of bacteria that are active in different temperature ranges. Bacteria populations are considered at the micro level depending upon the local conditions that contribute to growth and inhibition. The evaluated bacteria population numbers are used as modifiers on the rate equations for ferrous and/or sulphur oxidation.
An equilibrium relationship is assumed to determine the ferric ion concentration based on the concentrations of ferrous ions, dissolved oxygen and free acid, and an equilibrium constant, K c . The relationship is given by, where K c = 5.56 × 10 7 [63]. The ratio of ferric ions to ferrous is therefore given by An equilibrium relationship is solved using a Newton-Raphson scheme after all other chemical reactions have been completed in a time step.

Liquid-Solid Precipitation
Acid levels play an important role in bioleaching. The ferrous oxidation process requires acid along with oxygen, as do some other sulphide mineral leach reactions. Ferric iron salts are soluble at low pH, while at higher pH levels ferric iron salts can precipitate. The precipitation process creates acid as a by-product thus helping to buffer the acid levels in the heap. The ferric precipitation reactions are incorporated into the model by monitoring the local concentrations of iron salts and pH. At pH <2, the precipitation of insoluble hydronium jarosite is allowed. At pH >2, soluble ferric hydroxide is precipitated. The model assumes there is virtually no ferric in solution at pH >4 [64].
The hydronium jarosite precipitation reaction can be written as follows: Considering the stoichiometry for H 3 O + -jarosite, one mole of hydronium jarosite formed will remove 3 moles of ferric and 2 moles of sulphate and will add back 5 moles of H + (assuming that the sulfuric acid product dissociates to H + and HSO 4 − ). The relationship controlling hydronium jarosite precipitation can be written from the solubility product for the solid species on the right-hand side of the equation above.
The pH level at which H 3 O + -jarosite is precipitated is < 2. The KJar equilibrium value is calculated to give the same solubility as the hydroxide model at the pH at which jarosite starts to form.
The activity of ferric ion is described to be pH dependent [63]: with a heat of reaction ∆H = −19.83 Kcal/mol. Ferric hydroxide is generated by relating the maximum ferric in solution to pH using:

Leach Kinetics
The principle mineral reactions in the ore modelled in this paper are chalcopyrite (CuFeS 2 ) and pyrite (FeS 2 ). The shrinking core model is applied to each particle size distribution.

Pyrite Rate Kinetics
The pyrite reaction is, This is governed by the kinetic rate equation which is adapted from [64], where β is the fraction of pyrite reacted, r p is the pyrite grain radius (m), ζ is the particle shape function, T is the temperature (K) and R A is a rate constant based on known data.

Chalcopyrite Rate Kinetics
Dissolution of chalcopyrite is considered to occur through the following reaction: Munoz [65] reported that the mechanism by which chalcopyrite dissolves is by the continuous removal of iron and copper from the lattice with iron being removed first and then followed by removal of copper. The difficulty of leaching chalcopyrite under ambient conditions is generally considered to be surface passivation either from reaction products or adsorbed surface species layers, such as jarosite formation. The sulphur-rich layer presents a diffusion barrier to the continued dissolution reaction for copper. The initial kinetics of the reaction are uncertain and the sulphur or sulphur-rich layer is considered to form quickly. In general, the initial dissolution kinetics are considered to be quite fast, and then slow as the sulphur-rich layer forms. The rate of dissolution has been reported to be relatively insensitive to ferric concentration [66]. At ambient temperatures, the dissolution of chalcopyrite is limited and generally very minimal, and [64] suggests that the initial kinetic rate can be ignored under these conditions because the sulphur or sulphur-rich product layer forms quickly. Kaplun et al. [67] used kinetic rate and shrinking core analysis to investigate the activation energies associated with chalcopyrite leaching at moderate temperature. In this formulation, it is assumed that the rate-determining step in the reaction kinetics is the transport of electrons through the elemental sulphur or sulphur-rich layer.
A kinetic rate equation was developed from small scale leach tests for chalcopyrite by [66]. The fraction of chalcopyrite reacted was modelled using a shrinking core and was derived as shown in Equation (21).
The above rate expression and several other variations of the Arrhenius-type were a poor fit for experimental data at higher temperatures. Munoz et al. [66] further developed the kinetic leach rate equation for chalcopyrite by describing the rate-limiting step for the leaching of chalcopyrite as the electron transfer across the metal-deficient or sulphur-enhanced layer. The equation developed is: The molar density and electronic charge of chalcopyrite can both be considered constants for the type of simulations considered in this study. The parameters that are not constant are the free energy of the system and electrical conductivity of the sulphur layer. The free energy of the system can be defined as: where ∆G 0 R is the free energy of reaction as defined in Equation (24) ∆G 0 When Equation (24) is substituted into Equation (23) and all constants removed, the temperatures can be combined as in Equation (25), Fe 3+ (25) What remains for consideration is the electrical conductivity of the sulphur layer. The electrical conductivity of the sulphur layer is temperature dependent, and [68] showed a strong exponential relationship between temperature and electrical conductivity of elemental sulphur. Combining the predicted electrical conductivity relationship with Equations (25) and (22) gives Equation (26), which is more temperature sensitive. (1 − β)

Column Leach Tests
Experimental column leach tests were performed on 1.5 m (5ft) columns consisting of 41.64 kg of chalcopyrite containing ore of approximately 0.26% copper and 1% pyrite. The gas phase and liquid phase were given the temperature dependent properties of air and water respectively. The solid phase was assumed to have a constant density of 2950 kg·m 3 and thermal conductivity of 2.1 W·m −1 ·K −1 and specific heat capacity of 1172 J·kg −1 ·K −1 . A number of columns were leached, each column contained similar material and was held at a constant temperature, 25 • C, 35 • C, 50 • C and 75 • C. The particle size distribution, percent weight and percent copper concentration are given

Model Validation
Before the model can be used to predict recovery, the model must be calibrated to the reactivity characteristics of the ore. The material specific reaction rate modifier is used to characterise an ore in terms of reactivity. This rate modifier can then be applied to predict how this ore type will perform in other leach scenarios. The reaction rate modifier is mineral dependent and particle size independent.
One-dimensional flow was assumed for the column test and wall effects were assumed negligible. A mesh sensitivity study was undertaken and found the recovery varied only by approximately 2% between employing a very coarse mesh of 6 elements and a finer mesh of 60 elements. To enable a large number of simulations to be carried with little computational effort the test column was represented by 10 hexahedral elements. The ore was assumed to be initially dry, a fixed flux boundary conditions was applied to the top surface with a symmetry condition applied to the walls. A free drainage boundary condition was applied to the bottom surface. An adaptive time step scheme was employed. Further details on the solution algorithm can be found in [61]. The column test simulations take very little computational effort and complete in minutes.
The model was calibrated against experimental data for a base temperature of 35 • C. The reaction rate modifier was calibrated against copper recovery and particle diffusivity was calibrated against residual data. The copper, iron, and acid concentrations in the recovered solution compare well with the measured data, see Figure 3. In general, the predicted pH is slightly lower than the measured but the iron and copper concentrations are a good match apart from day 50 to 100 where much higher content is seen experimentally. As can be seen in Figure 4, the residual copper in chalcopyrite per particle size distribution is also in reasonable agreement with experimental analysis.
The model was then employed to predict recovery of column leach tests held at different temperatures, 25 • C, 50 • C and 75 • C. Figure 5 shows the copper predicted recovery against the experimental recorded recovery. The predicted recovery compares well with recorded data except for the case where the column is held at 50 • C. In this case, the predicted recovery is a reasonable match over the first 150 days but fails to predict the rate of increase in copper recovery during the latter part of the leach. At the end of the leach period, the predicted recover is 7% lower than both experimental leach tests. However, given the good agreement in recovery at the high and low temperatures the model can give a reasonable prediction of the effect of temperature changes on copper recovery. This is particularly important when simulating large scale heaps as local temperature variations within the heap will be significant.   The model was then employed to predict recovery of column leach tests held at different temperatures, 25 °C, 50 °C and 75 °C. Figure 5 shows the copper predicted recovery against the experimental recorded recovery. The predicted recovery compares well with recorded data except for the case where the column is held at 50 °C. In this case, the predicted recovery is a reasonable match over the first 150 days but fails to predict the rate of increase in copper recovery during the latter part of the leach. At the end of the leach period, the predicted recover is 7% lower than both experimental leach tests. However, given the good agreement in recovery at the high and low temperatures the model can give a reasonable prediction of the effect of temperature changes on copper recovery. This is particularly important when simulating large scale heaps as local temperature variations within the heap will be significant.

Three-Dimensional 'Virtual' Chalcopyrite Stockpile
A three-dimensional heap, 54 m by 24 m, comprised of three lifts each 9-m high of the ore employed in the column leach test, Table 2. The particle size of the ore was set to the experimental leach test ore distribution size. It is recognised that this small size particle distribution is not representative of typical leach pad particle sizes and in reality, this size distribution would likely cause permeability issues with compaction and migration of fines occurring. However, for simulation purposes and to enable calibrated data to be employed in the modelling work, no compaction was assumed and the intrinsic permeability of the ore was set using a random distribution using a uniform probability function with a maximum intrinsic permeability of 9.73 × 10 −10 m 2 and a minimum of 4.87 × 10 −13 m 2 giving an average hydraulic permeability of 5 × 10 −3 m/s at 20 • C which is typical of a well-draining material. Bouffard [69] gives a review of ore types, permeability and agglomeration practices in heap leaching.
The heap geometry was meshed employing 279,936 hexahedral elements. Each element was approximately 0.5 m × 0.5 m × 0.5 m. This is a very coarse mesh but will allow a number of scenarios to be investigated employing the virtual heap on a standard pc within a reasonable time period. A fixed flux boundary condition was applied to the active leaching surfaces, meteorological data was applied transiently to all external surfaces and the bottom surface was set to a free drainage condition. Further details on the adaptive time step scheme and solution algorithm are given in [61]. The two-year leach schedule took approximately 2.5 days to run on a Dell Intel(R) Core(TM), i7-2770, CPU @ 3.4 GHz (Dell Technologies, Round Rock, TX, USA).
A two-year leach period was simulated, with the second and third lifts being stacked after 90 days and 180 days, respectively. Initially, a simulation study was conducted on copper recovery assuming the same conditions but different external temperatures to illustrate the importance of capturing the thermal conditions within the heap. The annual meteorological data employed in the simulations is shown in Figure 6 where two external temperature ranges are considered. An intermittent leaching schedule is then investigated employing thermal reaction kinetics and the recovery per unit volume of solution applied is compared.

Thermal Effect on Recovery
Initially, the ore added to the heap was set to a temperature of 20 °C and simulations were performed assuming the low and high external temperatures shown in Figure 6. In both scenarios, isothermal reactions were assumed. As shown in Figure 7, the copper recovery was relatively high in both scenarios but slightly higher in the case with higher external temperatures. Initially both heaps are at 20 °C and reaction kinetics is similar producing local temperature increases within the reaction zone. The internal thermal kinetics dominate the copper recovery with external temperatures having only a small effect over time. Figure 8 shows the thermal contours within the heap at 17 months; regions of high temperature within the heap develop in areas where the leaching agent is reacting with pyrite and chalcopyrite minerals. In both simulations, the temperature in the reaction zones of the heap reached in excess of 45 °C.

Thermal Effect on Recovery
Initially, the ore added to the heap was set to a temperature of 20 • C and simulations were performed assuming the low and high external temperatures shown in Figure 6. In both scenarios, isothermal reactions were assumed. As shown in Figure 7, the copper recovery was relatively high in both scenarios but slightly higher in the case with higher external temperatures. Initially both heaps are at 20 • C and reaction kinetics is similar producing local temperature increases within the reaction zone. The internal thermal kinetics dominate the copper recovery with external temperatures having only a small effect over time. Figure 8 shows the thermal contours within the heap at 17 months; regions of high temperature within the heap develop in areas where the leaching agent is reacting with pyrite and chalcopyrite minerals. In both simulations, the temperature in the reaction zones of the heap reached in excess of 45 • C.

Continuous or Intermittent Leaching
The heap leach simulation assuming low external temperatures was repeated using an intermittent leach schedule of one day on followed by one day off. The first simulation assumed that the solution application rate and composition remained the same. The results were surprising as in the first year of the leach, grams of copper recovered per litre of leach solution applied remained approximately the same and copper recovery over the two-year leach period reduced to 33%. It is generally accepted that the inclusion of resting periods allows particles to remain wetted with leaching agent as solution drains down, allowing oxygenation and diffusion of reacting species to the bulk solution, thus improving leaching efficiency. However, looking at Figure 9 which shows the composition of the recovered pregnant solution, the pH levels remain high for most of the leach period and a much higher ratio of ferrous to ferric in the recovered solution is also seen. This indicates that the copper recovery is limited by the availability of acid and is not necessarily due to the reduced volume of liquid.

Continuous or Intermittent Leaching
The heap leach simulation assuming low external temperatures was repeated using an intermittent leach schedule of one day on followed by one day off. The first simulation assumed that the solution application rate and composition remained the same. The results were surprising as in the first year of the leach, grams of copper recovered per litre of leach solution applied remained approximately the same and copper recovery over the two-year leach period reduced to 33%. It is generally accepted that the inclusion of resting periods allows particles to remain wetted with leaching agent as solution drains down, allowing oxygenation and diffusion of reacting species to the bulk solution, thus improving leaching efficiency. However, looking at Figure 9 which shows the composition of the recovered pregnant solution, the pH levels remain high for most of the leach period and a much higher ratio of ferrous to ferric in the recovered solution is also seen. This indicates that the copper recovery is limited by the availability of acid and is not necessarily due to the reduced volume of liquid. The final simulation assumed the acid applied to the heap remained unchanged, i.e., for an intermittent leach where the volume of solution is halved the H 2 SO 4 in the leach solution was increased to 1.54 gpl. As can be seen in Figure 10, the grams of copper recovered per litre of solution applied almost doubled in the first year and significantly higher throughout the leach period. Increasing the acid in solution increased the copper recovery to 42%. Although the overall copper recovery is reduced, the grams of copper recovered per litre of solution applied is increased which may be advantageous when solution management is critical.  These simulations illustrate the advantage of employing a CFD phenomenological model as a diagnostic tool as well as forecasting. This study was undertaken employing a relatively coarse mesh as the aim was to investigate a range of scenarios in order to compare the effect of different operational conditions. For this study, a coarse mesh is probably sufficient for comparison purposes however, it should be noted that given the complexity of the interacting physics the recovery may be dependent on the mesh resolution. For the column test a 2% difference in recovery was seen when increasing the mesh resolution. If employing the model for forecasting purposes a much finer mesh would be recommended. The other critical aspect of the validity of the model results when employed for forecasting purposes is the characterization of the ore and parameterization of the These simulations illustrate the advantage of employing a CFD phenomenological model as a diagnostic tool as well as forecasting. This study was undertaken employing a relatively coarse mesh as the aim was to investigate a range of scenarios in order to compare the effect of different operational conditions. For this study, a coarse mesh is probably sufficient for comparison purposes however, it should be noted that given the complexity of the interacting physics the recovery may be dependent on the mesh resolution. For the column test a 2% difference in recovery was seen when increasing the mesh resolution. If employing the model for forecasting purposes a much finer mesh would be recommended. The other critical aspect of the validity of the model results when employed for forecasting purposes is the characterization of the ore and parameterization of the model. Although a computational model can simulate days, months of leaching in minutes, if the model has not been adequately calibrated to experimental data the results may be unsound. As duplicate experimental leach tests often show some fluctuations in recovery it is sensible to ensure a representative sample has been employed for calibration purposes. However, given the complexity of the interacting physics there is always likely to be some measure of error in the predicted results.
A key advantage of the model is sensitivity analysis of commercial heaps to understand which process parameters affect the recovery and to what extent. Such a sensitivity study on characterized and calibrated data can give insight into the effect of operational process changes and, therefore, provide a useful tool for 'what if' scenarios in heap operations. Changes in operational parameters may have unexpected knock on effects due to the interacting physics. The detailed spatial inventory data, as well as solution composition data, allows an analysis of all aspects of the heap and provides insight into the cause of the effect of any operational change.

Conclusions
Heap leaching requires balancing operational parameters to optimise a complex suite of interacting gas-solid-liquid-thermal physics that cannot be captured with simplistic models. Copper sulphide ores especially have complex mineralogy and complex leaching behaviour, including oxidation, bacteria effects, gangue and precipitation, where thermal conditions are critical. CFD phenomenological models, which capture the complex physics and respond to changing operational and environmental conditions, thus enable detailed heap design, diagnostic and optimisation analysis. This paper has described a comprehensive CFD model that enables capture of all these complex interactions, thus allowing operational scenarios to be investigated. The CFD model requires a large amount of data and calibration to enable accurate forecasting, but arguably, the time and effort are justified by the potential to provide valuable foresight in any long-term operation.
Author Contributions: The comprehensive CFD heap leach modelling tools, described in this paper, has been developed over a decade by Mark Cross and his 'multi-physics' team at Swansea University. Nick Croft and Diane McBride developed models and simulation tools for the analysis of the solution and gas flow through reactive porous media and the coupling of micro reaction models with macro transport models within the non-equilibrium finite volume framework. This research has been pursued collaboratively with James Gebhardt, who provided valuable knowledge on the complex process chemistry and extensive studies with collaborating industrial partners for the detailed validation and verification of the copper sulphide heap leach model. Diane McBride has recently extended the flow model to include freezing, thermal and channelling effects and wrote this paper.