Modelling particle scale leach kinetics based on X-ray computed micro-tomography images

Article history: Received 17 July 2015 Received in revised form 11 January 2016 Accepted 11 February 2016 Available online 13 February 2016 The apparent leach kinetics for an ore particle within a heap leaching system depend on the chemical conditions in the fluids around the particle, themass transport within the particle and the reaction kinetics at the surface of each mineral grain. The apparent rate kinetics thus depend upon the distribution of the mineral grains, in terms of both size and position, within the individual ore particles, as well as the evolution of this distribution. Traditionally this behaviour has been modelled using simplified relationships such as the shrinking core model. In this paper a method for simulating this evolution and the resultant kinetics based directly on 3D XMT images of the internal structure of the particles is presented. The model includes mass transport through the gangue matrix, surface reaction kinetics and the dissolution and subsequent evolution of the individual mineral grains within the ore particle. Different minerals and mineral associations will result in different surface reaction kinetics. One of the key inputs into thismodel is thus the distribution of the surface rate kinetics. Amethod for experimentally determining this distribution is presented. The simulation results are compared to the evolution of real particles as they undergo leaching as measured using a time sequence of 3D XMT images of a leaching column. It was found that these simulations are able to accurately predict both the overall leaching trends, as well as the leaching behaviour of mineral grains in classes based on their size and distance to the particle surface. The leaching behaviour did not follow that of a simple shrinking core approximation, with the actual spatial and size distribution of the grains, aswell as the distribution of their surface rate kinetics, all impacting the apparent leach kinetics. For the copper ore particles used in this work the best fit to the experiments was achieved at an intermediate value of the dimensionless group that characterises the relative importance of surface kinetics to diffusion indicating that both need to be considered for accurate modelling. This paper thus demonstrates that using 3D XMT to provide both structural and kinetic data and incorporating this information into a particle scale simulator provides an improved basis for predicting particle scale leach performance. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Heap leaching is an important hydrometallurgical technology, which extracts valuable metals such as gold, silver and copper from low grade ores. The increase in global metal demand and the decline in the high grade ore resources are making heap leaching an increasingly attractive alternative to conventional processing techniques such as flotation followed by smelting. There are, though, limitations in the use of heap leaching such as long processing times and relatively low extraction efficiencies, especially for primary sulphides (McKinney et al., 2007;Wills and Napier-Munn, 2006).
The chemical conditions and concentrations in the fluids around an ore particle, the mass (and heat) transport within the ore particle and surface reaction kinetics (often catalysed by biological processes in sulphide leaching) are the three main factors directly affecting leaching behaviour. The conditions around the individual particles will be strongly influenced by the inter-particle hydrodynamics and mass transport. The hydrodynamics and mass transport within the fluid surrounding the particles have been studied at both heap and column scale experimentally and via simulation (Mostaghimi et al., 2014;Petersen and Dixon, 2007;van Hille et al., 2010). The surface reaction rates have typically been studied using finely milled particles in stirred tanks (Córdoba et al., 2008;Koleini et al., 2011). However, there have been few studies coupling the processes occurring within the individual ore particles to the interaction with external conditions, internal mass transport processes and surface reaction kinetics. It is the interaction of all these processes which controls the apparent leach kinetics that are observed.
Understanding these apparent kinetics and how they depend upon the chemical conditions and the extent of extraction is an important consideration in the modelling and simulation of heap leaching. Current heap leaching simulators typically make use of simplified models to describe this relationship, often based on shrinking core models and variants of this model (Bartlett, 1998;Dhawan et al., 2012;Dixon and Petersen, 2003;Ogbonna et al., 2006). These models make a number of implicit and explicit assumptions about the leaching behaviour at the particle scale. The assumptions made by some or all of these models include that the particles are spherical and all the same size, the mineral grains are uniformly distributed throughout the particle, there is uniform porosity and that there is a sharp interface between the leached and unleached portions of the particle. These models need to be comparatively simple and easy to solve in order to be implemented within a heap scale simulator as they need to be used to predict the apparent kinetics at every point in the simulation domain at each time step. It is thus important to be able to understand the extent to which the various assumptions made in these models are valid and to what extent they impact the predictions.
The main aim of this paper is thus to develop a model that relaxes some of these assumptions in order to produce better predictions of the particle scale leaching kinetics of these systems. The purpose of this simulator is not to model an entire column or heap, but rather to model the performance of individual particles within a column. While beyond the scope of this particular paper, the ultimate purpose of this work will be to ascertain how the spatial and size distribution of the grains within the ore influence its apparent leach kinetics. This could then be used to develop apparent leach kinetic models that do a better job of incorporating mineralogical variability than the current particle scale models do and/or to understand how the fitting parameters in existing models depend upon the mineralogy and its texture.
Before the variability of the ore can be incorporated into the modelling a method for measuring this variability and its impact on leach performance is needed. The technique used in this paper is X-ray micro computed tomography (XMT or micro-CT), a non-destructive method for imaging the internal structure of materials, which has been extensively used and is now being applied to the study of geological (Blunt et al., 2013;Cnudde and Boone, 2013;Fonseca et al., 2013;Hall et al., 2010;Ketcham and Carlson, 2001) and minerals processing (Burlion et al., 2006;Ghorbani et al., 2011;Lin and Miller, 2005;Miller et al., 2003) systems.
Within an XMT system the sample is placed in an X-ray beam, and rotated through 360°while a series of micron level 2D projection images are collected by a detector. By using mathematical principles of tomography, these 2D images record the variation of X-ray attenuation within the sample and are reconstructed to produce a 3D image where each voxel (smallest volume element, equivalent to a 3D pixel) represents the X-ray attenuation at each single point (Landis and Keane, 2010;Lin et al., 2015).
Within the context of heap leaching XMT has been used to track the extent of leaching within ore particles (Safari et al., 2009), as well as morphological parameters such as the size and shape of particles and the distribution of the mineral grains within them, especially their relationship to features such as internal porosity and fractures (Kodali et al., 2011;Miller et al., 2003).
There has, though, been very limited work on the use of XMT to form the basis of modelling of heap leach behaviour. Ghorbani et al. (2013) modelled the dissolution of mineral grains within ore particles and compared their results to that obtained from XMT. In their work they measured the extent of leaching as a function of the distance to the grain surface by dividing the ore particles into a number of distance bands and compared their results to a shrinking core model and a one dimensional (spherical coordinates) reaction diffusion model solved numerically. They also compared the average extraction to a modified and unmodified version of the variable rate constant model. The work in this paper looks to extend this approach to not only using the XMT data to validate different particle scale leach models, but to also directly use the XMT data in the prediction of the performance. This means that the actual shapes of the particles and distributions of the mineral grains can be used in the predictions, rather than having to make assumptions such as spherical particles and uniformly distributed grains. The disadvantage of this approach is that it requires that that the mass transport within the ore particle be solved in 3D, which is computationally very expensive.
The first part of this paper describes the experimental methodology, including the image analysis applied to the XMT data. This is then followed by a description of the mathematical model and the computational approach used to solve it, while the final section of the paper compares the experimental and simulation results.

Methodology: experimental and 3D image analysis
The ore sample used in the experimental work came from Kennecott Copper's Bingham Canyon mine. The main sulphide minerals within the ore were chalcopyrite and pyrite, with small amounts of secondary sulphides also present (See Table 1). The leaching tests were carried out over a period of 200 days in a small scale leaching column with a 28 mm diameter and a 190 mm height. The leaching solution was 0.1 M H 2 SO 4 with 5 g/L Fe 3+ in the form of ferric sulphate. As ferric sulphate was added in the feed, the leaching did not rely on the biological leaching of the pyrite to produce ferric ions. The column leaching was carried out in a temperature controlled incubator at 60°C.
Note that in this paper when we refer to mineral grains it is the mineral grains of interest that we are referring to. In the specific experimental system used in this paper these minerals are the metal sulphides, though the methodology could be applied to any leach system. It should also be noted, though, that one of the shortcomings of XMT when applied to the copper leaching system is that the attenuation of Chalcopyrite and Pyrite are very similar to one another and that the recoveries presented in this work are the total amounts of sulphide leached rather than the copper extraction (see Fig. 1). The difference in the leach rate of these different minerals are accounted for by the fact that a distribution of surface leach rates is experimentally determined (see Section 4.2).
The reason for using these small columns in this work is because it is the effect of particle scale kinetics, and thus intra-particle phenomena, that are of interest. These columns are small enough that the particles within them can be scanned at high resolution without the need to unpack the particles from them, which would disturb the leach process. The small size also means that the wetting and reagent conditions experienced by each of the particles is likely to be similar. This is not the case in large columns or real heaps where channelling/preferential flow and other mass transfer effects can result in different particles experiencing quite different leach conditions. The leaching column was scanned at a sequence of time points (Day 0 dry scan and the end of leaching Days 5,11,23,43,53,83,118,136,168) using a Nikon Metris Custom Bay with a 1 mm aluminium filter to reduce the effect of beam hardening, an energy level of 89 kV, a 0.708 s exposure time and 2000 projections. The detector size is 2000 × 2000 pixels, which gave a linear voxel size of approximately 17 μm for the field of view used. XMT can easily be used to measure average extent of the sulphide leaching, while ICP analysis of the leachate can be used to measure average copper extraction. Fig. 1 shows the copper extraction versus the extent of sulphide leaching, which shows that initially the copper extraction is much quicker than then the average extent of sulphide leaching, which is probably due to the much quicker leach kinetics of the secondary sulphides and copper oxides, though once these are leached, the copper and sulphide recoveries approach one another. This will be because both chalcopyrite and pyrite are quite slow leaching and because eventually the leaching rate will become more diffusion dominated as mineral grains near the particle surface are depleted (as diffusion becomes more important, differences in the mineral surface kinetics become less important). Average leach behaviour, though, doesn't tell us much about the mechanisms at work and, especially, the relative importance of mass transport to surface kinetics, a crucial parameter if we wish to predict the apparent leach kinetics and its evolution. Fortunately, while chemical analysis of the leachate can only give us average extractions for the whole column, the XMT data can give us the location and time resolved leach rates for every mineral grain within the system, a vastly larger amount of information.
The conversion for the mineral grains can be calculated by the volume difference between each scanning time point and the initial grain volume in the Day 0 dry scan. Thresholding was applied to differentiate between the rock and the air phase, including internal porosity, and between the gangue and sulphide phases. All the mineral grains within the field of view were tracked individually over the course of the experiment. The segmentation of the individual ore particles within the column is achieved by using the Otsu algorithm to set the appropriate threshold between the air and the rock, with a marker base watershed algorithm used to distinguish the different ore particles. A tracking algorithm was then used to identify the particles across different time points as they move and rotate within the bed as the leaching progresses. The mineral grain segmentation uses the Maximum Entropy algorithm (Kapur et al., 1985). Since the translation and rotation of each ore particle relative to their location in the initial scan is known, the expected translation for every mineral grain in the original image can be calculated, which allows them to be identify in subsequent scans. The methodology for applying the segmentation and tracking these grains has been described in more detail in a previous paper (Lin et al., 2015). The algorithm labelled the features within the image, with the external air being given an integer label of 0, the gangue phase a label of 1, internal porosity a label of 2 and each of the sulphide grains a unique identifier starting from label 3 (Fig. 2). It should be noted that in this work, only the pores connected to ore surface were considered. Fig. 3 shows a sequence of XMT images that show the extent to which the grains are leached. Note that these figures are a virtual slice through the 3D volume of a single ore particle.
Since the XMT image analysis algorithm gives the location, initial size and the extent of leaching of every mineral grain within the particles in the small scale columns, it can be used to obtain the grain scale kinetics. In Section 4.2 a method for using this data to measure the variability in the surface kinetics, which is a key input into our modelling, is presented. This data is also used for validation of the modelling as the simulations not only predict the overall average leach behaviour of a particle, but also the mineral grain scale behaviour, which can be compared to the measured grain scale behaviour and how it varies with grain size and position within the particle. The actual 3D XMT images of ore particles together with the mineral grains within them are also used directly as the initial condition for the simulations, thus allowing the effect of particle shape and mineral grain distribution to be implicitly accounted for.

Mathematical model
The mathematical model takes the form of a set of partial differential equations that describe the motion of the reagents and which are coupled to a model for the mineral grain dissolution. These must then be solved numerically as the system geometry and its evolution is too complex for analytical solutions. This section describes the derivation of the governing equations and the methods used for numerically solving them.
As mentioned in the Introduction, the starting point for the assumptions made in this modelling are those used in the standard shrinking core model. The aim of this paper is to improve the predictions by relaxing certain of these assumptions by making use of XMT data, but with some of the assumptions maintained. In particular, the assumption Fig. 1. Average copper extraction as obtained from ICP measurements vs sulphide leached as obtained from XMT. Note that these are the average recoveries for the entire column rather than the particle and grain scale behaviour presented later (Unlike XMT analysis, chemical analysis of the leachate can only give average recovery). of spherical particles with uniformly distributed mineral grains is relaxed by incorporating the actual particle shapes and mineral grain distributions. The main assumptions that are maintained are that the mass transport is at quasi-steady state, that the surface kinetics are linear and that there is a uniform diffusivity. The last of these assumptions is maintained mainly because there is no data available as to the variability of this diffusivity at this sub-particle scale.

Reagent motion
Assuming that the reagent motion is predominantly diffusive, its flux, F, can be described by: where C is the reagent concentration and D is the effective diffusivity. A continuity equation is required to solve for the concentration within the ore particles as a function of time, t: In this work, this equation is simplified by making the approximation that the behaviour is quasi-steady state. In other words it is assumed that the accumulation of reagent species within the ore particles has an unimportant impact on the flux of the reagent through the particle, which implies that the diffusion timescale is shorter than the reaction timescale: The quasi-steady state assumption does not imply that that the flux remains constant with time (in the simulations presented below both the magnitude and the spatial distribution of the flux vary quite markedly with time). What this assumption does imply is that the current flux (and concentration) distribution depends only upon the current state of the mineral grains within the particle, as well as the external reagent concentration, but does not directly depend upon the flux history as the accumulation terms are ignored. Within the bulk of the ore particles the quasi-steady state assumption means that the concentration is governed by the following equation: The diffusivity will vary with the properties of the gangue through which it is moving, particularly as a function of the micro-scale porosity, which will typically be below the resolution at which the samples in this work are scanned. While the simulation framework that has been developed allows different regions to be assigned different diffusivities based, for instance, on differences in gangue mineralogy, without data on the variability of the diffusivity it is assumed constant in this paper. It is also possible that dissolution of gangue and/or deposition of species such as jarosite could alter the apparent diffusivity with time, but as there was no measurable change in the porosity in these experiments, there was no basis on which to include this effect, though again the simulator could incorporate these effects if the data on which to base it were available.
Chemical reactions occur at the surface of the metal sulphide. In general these reactions will be dependent on the concentration of both reactant and product species. For the purposes of this paper it will be assumed that there is a single limiting reagent and that the surface reaction kinetics are first order w.r.t. this reagent: where n̂is the outward unit normal for the sulphide grain boundary ∂ MS and k react is the first order surface rate constant for the reaction (note that unlike a volumetric first order rate constant, which has units of inverse time, the first order surface rate constant has units of length per unit time). Since the flux due to the reaction must match the diffusive flux into the surface: At the surface of ore particle there will be mass transport between the particle and the bulk fluid. The effect of boundary layers in the fluid and other related phenomena can be modelled using an external mass transfer coefficient, k ext : Similarly to the reaction flux, the mass transfer into the particle must match the flux within the particle: As the rates will all scale with the external concentration, C ext , it is useful to non-dimensionalise the concentrations using this quantity: The following set of equations can therefore be used to model the quasi-steady state reagent concentration within the ore particle: Boundary conditions: Surface of ore particle : The assumptions that the mass transport is quasi-steady state and that the surface reactions are linear are both potentially problematic, but, as was mentioned above, were chosen in this initial model development because they are both made in the most shrinking core type models. Relaxing both of these assumptions is the subject of ongoing work.
The set of equations is solved in 3D using the voxelised grid obtained from the XMT imaging. It is solved using a finite volume scheme for the discretisation, with the fluxes calculated on a staggered grid (a similar method is employed in Neethling and Cilliers (2003)). As the voxels are cubic, the approximation is very similar to that which would be obtained by simple finite differencing, though the finite volume scheme allows the complex shaped boundaries of both the ore particle and the mineral grains within it to be more naturally accommodated.
The problem is solved by matrix inversion as the assumption of first order kinetics means that the discretised equations are linear. As the XMT images are very large (10 8 to 10 9 voxels) and the numerical solution is obtained at the same resolution, the solver was implemented in parallel to both improve the computational times and to distribute the very large memory requirements over a number of nodes. The code was written in C++ and made use of MPI and the PETSc library (Balay et al., 2014) for the parallel implementation and matrix inversion. Fig. 4 shows a cross-section through a 3D image of the nondimensionalised reagent concentration. This simulation contained 300 million voxels and took about 48 min to solve on 50 cores. For subsequent concentration profiles as the grain distribution evolves the solution times are much quicker as the concentration profile from the previous time step can be used as quite an accurate initial guess for the next time step.
As it is being assumed that the surface kinetics are linear in a single rate limiting species, the exact identity of this species is unimportant. Even whether this rate limiting species is a reactant or a product does not impact the mathematical form. In these particular experiments, though, the limiting reagent is likely to be the ferric ions. In future work we will be relaxing the assumption of linear surface kinetics and including multiple reactant and product species. At this point the identity of these species and their relative concentrations in the surrounding fluid will impact the results.

Mineral Grain dissolution
The model can be used to calculate both the reagent concentration and its flux at each point in the system. In order to use this information to predict the dissolution process, the fluxes must be coupled to a model for the evolution of the mineral grains. If it is assumed that the ratio of the volume of the mineral grain leached to the reagent consumed is κ, then the rate of change of a grain's volume, V, can be expressed as follows: This equation on its own will only give the change of volume of the grain, but not the change in shape. In order to estimate both the change in shape and the change in volume, the equation must be applied at the voxel rather than the grain level. To avoid the dissolution of mineral within a voxel being binary, an additional scalar field S is introduced, which is the fraction of metal sulphide remaining within a voxel. Initially all the voxels that are identified as being within sulphide grains are assigned a value of S = 1, while S = 0 in all other voxels. The value of S in voxel i is then evolved by summing the fluxes over the surface of the voxel: where Δx is the edge length of the cubic voxels. Eq. (14) is the voxel level equivalent of Eq. (13), where the flux over a voxel face is assumed constant. As the fluxes are zero at all boundaries between sulphide voxels, the value of S will only change in voxels at the boundary of the mineral grains. When the value of S ≤ 0, the index for that voxel is changed from that of the mineral grain to that of the gangue, indicating that sulphide within that voxel is completely leached. To accurately capture the dissolution processes, the time step must be such that the fastest dissolving voxel takes at least one time step to disappear: Fig. 4. Reagent concentration within a single ore particle based on initial grain distribution. All reagent concentrations are non-dimensionalised using the external concentration. The non-symmetrical concentration profile is caused by both the irregularity of the particle shape and the non-uniformity of the mineral grain distribution. Note that this is a 2D slice through a 3D simulation and therefore much of the non-uniformity is out of plane.
As the disappearance of mineral voxels will impact the fluxes within the system, the concentration (and hence the fluxes) needs to be recalculated at every time step during which a metal sulphide voxel becomes completely leached. Between such steps, however, S i can be updated without needing to recalculate the concentration field due to the quasi-steady state assumption.

Ore particle case study
In addition to non-dimensionalising the concentration using the external reagent concentration, the simulation parameters also form dimensionless groups, which allows us to reduce the number of independent parameters that need to be examined. The most important of these dimensionless groups is a surface reaction Damköhler number, Da, which indicates the relative importance of the diffusive mass transport and the surface reaction rate: where l is a length scale associated with the particle, taken to be the equivalent spherical diameter of the particle. While the actual values of the flux will depend upon the specific values of the external concentration, diffusivity and rate constant, the shape of the internal concentration profile depends only on this dimensionless group and the current geometry of the system (assuming that the external mass transport coefficient is very large relative to the surface rate constant, an assumption made in this study). The leaching becomes diffusion limited when Da ≫ 1 and reaction limited when Da ≪ 1.
The simulations are carried out using a dimensionless time, t ⁎ , which is related to the actual time by the following relationship:

Diffusion limited and reaction limited cases
Fig . 5 shows the reagent concentration on a slice through the 3D simulation of an ore particle at four different time steps for a simulation with a high value of Da (Da = 100). The completely diffusion limited case (Da → ∞) (Santiago, 2001) cannot be directly simulated as it would require either a zero diffusion coefficient or an infinite reaction rate, but this simulation approaches this limit. In the diffusion limited case there is a distinct profile in the reagent concentration, with high reagent concentrations near the edges and low concentrations in the middle of the ore particle. This results in faster leaching of the mineral grains which are closer to the edge of the particle and slower leaching of those in the middle. Fig. 6a-c shows the 3D distribution of the mineral grains and their evolution with time. Given the 3D nature of these images, it is hard to see the change in size of the grains and their position relative to the surface. To address this, the average dissolution of grains has been plotted (Fig. 7a) at a particular distance from the surface as function of both that distance (x axis) and time (colour).
It is clear from this plot that the dissolution progresses from the outside inwards for the diffusion limited case. This inward progression of the leaching is similar to that predicted by the shrinking core model, though even when the system is very diffusion limited, as in this case, there is still not a distinct separation between the leached and unleached portions of the particle, with larger mineral grains that are Fig. 5. Reagent concentration for an ore particle at four different time steps with Da = 100.0 (diffusion limited). The white contour shows the outline of the ore particle. All reagent concentrations are non-dimensionalised using the external concentration.
close to the outer edge of the particle lasting longer than some smaller grains that are further in, with the shape of the ore particle and the distribution of the mineral grains further complicating the dissolution pattern. It is this variability in leach rate with mineral grain size that causes the variability in average extraction with distance seen in Fig. 7a.
At much lower values of Da (Da = 0.001) the reaction limited case is approached (Fig. 8), which is characterised by reagent concentrations throughout the particle that are virtually the same as those in the surrounding fluid (the zero concentrations in Fig. 8 are in the interior of mineral grains as all reactions occur at the surface of the grains). The  slight decrease in concentration with distance is due to the fact that, while Da is small, it is not zero. As the concentration of reagent is virtually constant in this case, the rate of dissolution of the individual grains is essentially independent of their position within the ore particle. The extraction rates will depend on the size of the mineral grains, though, as well as any variability in the surface reaction rates (assumed constant in these simulations, though this assumption is relaxed in Section 4.2). This can be seen in Fig. 6d-f and more clearly in Fig. 7b. Fig. 9 shows the experimentally obtained leaching profile based on XMT measurements through time. It is clear from the figure that the leaching profile does not follow the shrinking core behaviour as there is a gradual change in the average extent of dissolution with distance from the surface rather than a distinct separation between leached and unleached regions. The fact that leaching occurs in the centre of the particle (recovery increasing) at all times, albeit slower than towards the outside of the particle, indicates that in this system both diffusion and surface reaction rates play a role in the apparent kinetics. In a system with an intermediate value for δ, the observed rate will initially be dictated by the surface reaction kinetics as the grains near the surface, which will have minimal diffusion resistance to their leaching, will contribute most to the leaching rate. As these grains are depleted the grains deeper in the particle contribute more to the apparent leach rate and the effect of the diffusion resistance becomes stronger.
While in the idealised simulations above the surface kinetics rate constant for all the grains was assumed to be the same, the situation in the real particle is that the surface rate constants will vary from grain to grain due to differences in the mineralogy of the grains and their associations. A method for estimating the distribution of these rate constants is thus required before accurate simulations of the ore behaviour can be carried out.

Reaction kinetics distribution
In the leaching experiments the dissolution of every grain in the column is tracked as a function of time. This provides a lot of data as there are hundreds of thousands of grains in each column. Due to the large amount of available data, it is possible to divide the mineral grains into quite narrow categories based on their initial size and their distance from the nearest surface of the particle while still having quite a large amount of data in each category. If it is assumed that the size classes are narrow enough that the effect of particle size on the recovery within the category is small and that the distance from the surface can be used as a proxy for the chemical conditions experienced by the grain, then the remaining variability in the leach rate is mainly due to variability in the surface kinetics. Note that an assumption that the chemical conditions experienced by a grain within a category be constant with time over the time interval being considered is not required, only that the grains in a category experience a similar set of chemical conditions. The reason why distance from the surface can act as a reasonable proxy for the mass transfer resistance is that the mass transfer rate is proportional to the concentration gradient, which will be a strong function of the distance between the mineral surface and the particle surface. While it is a reasonable proxy it is not a perfect one as the distribution of neighbouring mineral grains will influence the concentration Fig. 8. Reagent concentration for an ore particle at four different time steps with Da = 0.001 (reaction limited). The white contour shows the outline of the ore particle. All reagent concentrations are non-dimensionalised using the external concentration. Fig. 9. Recovery against distance for an example ore particle based on image measurement. gradient being experienced by a particular grain as they will also be consuming reagent. The variability in the diffusivity of the reagent in the gangue species surrounding the particles will also impact the validity of the use of distance from the surface as a proxy for the mass transfer resistance.
At time t, the leaching rate of grain i of size V i (t) and in initial size category V i (0) when exposed to a concentration C can be expressed as follows, again assuming linear surface kinetics: where k i is the surface rate constant for grain i, A i is the surface area of grain i at time t, V i is the volume of grain i remaining at time t, and a is the proportionality in the relationship between volume and surface area. From Eq. (19), the reaction kinetics can be derived as: Integrating this equation gives: The ratio of the surface rate constant for grain i to the average surface rate constant for grains in the same size and distance to the surface category, k mean , can thus be calculated from the experimental data: where V i ð0Þ 1=3 −V i ðtÞ 1=3 is the value of this function at the average leach rate of the grains, given by: where N is the total number of grains in the category. Fig. 10 shows the cumulative distribution function of the rate constants for different size and distance categories. While the average rate constant varies markedly with position in the ore particle and with the size of the mineral grains, the distribution of rate constants relative to the mean is essentially the same in all categories. What this means is that, while the fastest leaching grains have leach rates orders of magnitude larger than the slowest leaching ones, what is independent of the grain distribution is the relative variability in this rate. That is, the proportion of grains that leach, for instance, twice as fast as the average rate for their particular size and distance category is the same for all categories (from Fig. 10 it can be seen that just over 10% of grains will have a leach rate at least twice as fast as the average for their size and distance class).
The slight scatter in the plots at high values of the rate constant is due to the fact that there is a finite range of rate constants that can be measured, since some of the fastest leaching grains disappear over the finite time interval (23 days) used in these measurements in the categories with the highest leach rates (small grains near the surface). It can also be observed that there is a small portion of grains with negative values of their relative rate constant (~5% in volume). This is mainly due to measurement error as there is some uncertainty in the measured volumes and thus a few grains with a rate around zero will appear to have a negative rate.
We can therefore use this average distribution as the basis of a Monte-Carlo type simulation in which the surface rate constant for each grain is assigned randomly based on the distribution, with the simulation repeated to ensure that the results are not strongly dependent on the rate constants assigned to specific grains (while there are a very large number of grains in the system, a lot of the grain volume is contained in a comparatively small number of larger grains). The cumulative distribution function can be converted to a probability density function and needs to be parameterised so that it can be easily input into the simulations. As can be seen from Fig. 11, a gamma distribution Fig. 10. Cumulative distribution function of the rate constants for different size and distance categories, with the average distribution function for all the plotted data. fits the data with positive reaction rates quite well. The blue area shows the portion with negative values (~5% in volume), which is mainly caused by the measurement error for non-leachable grains. 5% of the grains are thus assigned a rate constant of zero to account for the nonleaching portion. It can also be observed that there are peaks in the plot for the grains with higher reaction kinetics. These peaks are probably caused by grains composed of faster leaching minerals such as oxides and secondary sulphides (mainly covellite in this material).
While the analysis above ascribes this variability to surface kinetic variability, in reality it will also include the effect of variability in the local mass transport resistance in the gangue surrounding the grains. De-convolution of these two effects is hard without being able to, for instance, measure differences in the mineralogy and in the porosity surrounding the individual grains, though it would be interesting to achieve in order to better understand the relative importance of the different micro-scale mechanisms involved in heap leaching. From a purely simulation/prediction perspective, though, de-convoluting these effects is unimportant as this methodology is implicitly allowing the influence of both these mechanisms to be included in the simulation results.

Simulation results after applying gamma distributed kinetics
The apparent leach kinetics of the ore particles depend on a number of factors such as particle size, porosity and mineralogy. As Fig. 9 indicates that the actual behaviour is neither completely reaction nor diffusion limited, Da was tested at intermediate values, with a value of about 1 producing the closest fit. Fig. 12 shows the leaching profile using Da = 1 with the plot showing the average leaching as a function of distance at a number of time points. Comparing with Fig. 9, it can be observed that although the recovery against distance profile is not exactly the same as the profile obtained from the experiment, the general shape and trend of the profile are similar. The small scale variability in both Figs. 9 and 12 are mainly due to the leach kinetics of specific larger grains, which in the simulations are assigned randomly from the experimentally obtained distribution.
One of the main objective of this method is to simulate the leaching behaviour at longer times based on the initial response. The distribution of rate constants is found experimentally based on the measured leaching over the first 23 days using Eq. (23). The scaling factor for the time scale is fitted based on the first 53 days of leaching according to Eq. (17). Fig. 13 shows that the comparison between experimental measurements and simulated recovery is excellent. The test of this method is thus not how well it fits the initial leach response, but rather how well it fits at later times as this data was not considered in the calibration. Many models could be made to fit this data, but from our experience most have difficulty predicting the long term response based on the initial leach performance. One of the keys to the excellent performance of this approach is the use of the XMT data to measure the variability in the leach rate. Using this approach with only a single surface reaction rate leads to a serious over-estimation of the long term leach performance, meaning that it is not sufficient to know the size and spatial distribution of the grains and that surface kinetic variability also needs to be considered. Fig. 14 shows the parity between the experimental and simulated average recovery at the different measurement time points, with each of the sub-figures being for a different grain size and distance category. It shows that not only is the simulator able to predict the average behaviour, it is also able to predict the behaviour of specific categories of mineral grains.

Conclusions
In this paper a simulation methodology for predicting the grain scale behaviour within individual particles was introduced. This model made some of the same assumptions as most shrinking core, such as that the mass transport was pseudo-steady state and that the surface kinetics were linear. Where it differed markedly was that it made use of an XMT image as its simulation grid, which meant that the effects of the actual particle shape and the spatial and size distribution of the mineral grains were captured.
Another key feature of this model is that it included the effect of the variability of the surface rate kinetics in its predictions. This is important as it has a large impact on the long term evolution of the apparent kinetics. Fast leaching grains disappear first and thus the kinetics slow down with time. The variability in the particle kinetics thus impacts the evolution of the leach rate, with greater variability causing more of a slowdown relative to the initial rate. This leach variability comes from not only variability in grain size and position, but also variability in the surface kinetics. As the simulations inherently account for the effect of grain size and position, the effect of surface kinetic variability needs to be separately included. A methodology was developed by means of which the XMT data over the initial leach period could be used to decouple the effect of particle size and distance to the surface from the surface kinetic variability. It was found that approximately 5% of the grains had a zero rate constant, with the variability in the leach kinetics for the remaining grains following approximately a gamma distribution.
The balance between mass transfer and surface kinetics was characterised by the Damköhler dimensionless group, Da, with a value of around unity giving the best results. The simulation was carried out using an XMT images of the ore particles as its simulation grid. This required that the implementation be massively parallel given the high resolution and thus large number of voxels involved. Fig. 12. Recovery against distance for an example ore particle based on simulation (Da = 1). Fig. 13. Comparison between experimental XMT based measurement and simulation using gamma distributed kinetics. Note that the recovery is sulphide dissolution in the particles that were simulated rather than the overall recovery for the column.
The simulation was able to accurately predict not only the average leach behaviour, but also the behaviour of mineral grains within specific size and distance to the surface categories. It was found that the leaching behaviour, especially at the grain scale, does not follow that predicted by the shrinking core models and that the leach behaviour is most accurately simulated when the diffusion and reaction models are coupled to the realistic input mineral grain distributions, thus including the effects of size, spatial and surface kinetic distributions.
The purpose of this modelling is not to replace existing simpler particle scale leach models, but to act as a means of testing the predictions and assumptions of those models and to help to develop improved models. While there is still a very long way to go, the ultimate aim of this type of approach is to significantly reduce the current need for very time consuming column experiments by supplementing and/or replacing them by 3D scanning and detailed simulations in which mineralogical and textural effects can be accurately accounted for.
The assumptions of pseudo-steady state mass transport and first order surface kinetics are being relaxed in ongoing work by the authors and a study of the effect of these assumption will be presented in future papers.