Published on

Hydro-mechanical-chemical coupled modeling of copper heap leaching with iCP



  • Name
    Orlando Silva


Several coupled Hydro-Mechanical-Chemical (HMC) processes are involved during the construction and operation of a heap leach pad. Due to the complexity of the problem, common existing heap leaching models solve the different phenomena in an uncoupled or partially coupled way. Also, often modeling of leaching of copper oxide ores assumes that mineral dissolution/precipitation is described by the diffusion-controlled shrinking core model and the system is composed by a single copper-containing mineral and a single mineral representing the gangue. Although these simplifications help to solve the problem, they could result in poor predictions of the composition of the pregnant solution (PLS) and consequently, an over- or under- estimation of the copper recovery and acid consumption. Further, improving the design of heap leach pads requires performing an integrated analysis of the HMC phenomena. Mechanics stability and mineral recovery efficiency cannot be estimated accurately without considering multiple coupling of the complex chemistry, soil mechanics and hydrodynamics. Here we present a coupled HMC model that accounts for the kinetics of several mineral species and the evolution of a muti-lift heap. Numerical simulations were performed with Comsol-PhreeqC (iCP), an interface that allows solving coupled unsaturated flow, soil mechanics and reactive transport using the Sequential Non Iterative Approach (SNIA).

Numerical implementation in iCP

The HMC governing equations are solved with the numerical interface Comsol-PhreeqC (iCP) (Nardi et al., 2013). iCP is a Java interface that couples the general purpose finite element framework Comsol Multiphysics ( and the geochemical simulator PHREEQC (Parkhurst and Appelo, 1999; In brief, the tool is based on implementing an operator splitting technique for solving the hydro-mechanical problem (a coupled or uncoupled system of partial/ordinary equations and algebraic differential equations) on Comsol and the chemical system of algebraic and differential equations on PHREEQC.

Simulation of a multi-lift heap leach pad

The present HMC model was applied to simulate the construction and operation of a synthetic copper heap leach pad consisting of five ore layers settled on a terrain with a slope of 2%. Each layer is simulated sequentially, mimicking the sequence of construction and operation of a real heap leach pad. We assumed that the construction of each layer is performed during 10 d, followed by an irrigation period of 90 d at a constant irrigation rate of 7.5 l/h-m2. The stacked ore consists of an artificial mineral mixture made from Chrysocolla (20% w/w), Chalcopyrite (10% w/w), Pyrite (8% w/w), Albite (30% w/w) and K-Feldspar (32% w/w), whose chemical evolution is assumed to be governed by kinetics reactions. Equilibrium precipitation-dissolution reactions of Jarosite-K, Jurbanite, Basaluminite, FeOH3(a) and Schwertmannite were also taken into account.

The figure below shows the evolution of concentrations of selected dissolved species at the output of the drainage system. The evolution of pH shows the arrival of the acidic solution to the liner system in a time of around 60 d (top left). Kinetic dissolution of Chrysocolla and Chalcopyrite are triggered by the infiltration of the acidic solution and, as a result, an increase of Cu concentrations is observed in the outlet flow (top right). The chemical evolution of impurities like Si and Al (bottom) clearly shows that under this low pH conditions kinetic dissolution of alumina-silicates is also favoured. This process is highly dependent on the pH and, over a long period of time, porosity can be significantly modified. As expected, the composition of the PLS responds to the irrigation and drainage rates of the system (shown in grey colour).

Evolution of pH (a) and concentrations of Cu (b), Si (c) and Al (d) in the pregnant leach solution.

Chrysocolla and oxide minerals leach rapidly under heap leach conditions. Chalcopyrite hardly reacts at all in heap leaching. The video below shows that the Chalcopyrite concentration (middle) practically does not change, while the Chrysocolla (top) is almost completely dissolved at the end of the leaching of each layer. The latter has a major impact on the porosity of the ore bed (bottom).

Evolution of Chrysocolla (top) and Chalcopyrite (mid) concentrations and porosity due to chemical reactions (bottom).

The copper recovery and acid consumption curves (see the figure below) directly reflect the different dissolution rates of the mineral species considered in the model, specially the Chrysocolla (fast dissolution rate) and the Chalcopyrite (slow dissolution rate). The first layer provides a copper recovery of about 30%, while the copper recovered from the leaching of subsequent layers is about 8%. At the end of the operation (500 d), approximately the 60% of the copper has been recovered. The copper recovering increases due to the leaching of residual mineral from previous layers, depending on the thickness of the loaded layers and the irrigation times. The competition between Chrysocolla and Chalcopyrite is also evidenced by the acid sulphuric consumption. The acid consumption during the leaching of a new layer is progressively lower than the acid consumption of previous leaching stages. This is because the mass of mineral of the new layer is lower and the impact of acid consumption from previous leached layers is minor.

Predicted copper recovery (top) and acid consumption (bottom).

Changes in porosity and permeability due to mineral dissolution/precipitation and irrigation time can affect significantly the hydromechanical behavior of the heap leach pad. First, the hydrodynamics behavior of the system is conditioned by the time that takes the irrigation solution to saturate the heap as it flows to the drainage system. After a short period of time, a quasi-steady state is achieved for every layer. Secondly, the mechanical effects are traduced in an increase in the vertical stress developed mainly on the rock foundation as the layers are progressively loaded. This vertical stress increase also reaches the top of the third layer, where the vertical displacement reaches a maximum. The evolution of liquid saturation, Darcy velocity field and vertical stress is shown in the next video.

Evolution of liquid saturation and Darcy velocity field (top) and vertical stress (bottom) in the leach pad.

There are zones of the ore where the porosity increases due to the mineral dissolution and other zones where the consolidation of the layers leads to a porosity decrease. In the present application case, chemical predominates over hydromechanical effects. The figure below (top) shows the combined hydromechanical and chemical effects on the global porosity at an observation point (white dot in the previous video). The consolidation of the lifts related to the porosity decrease is illustrated in the figure below (bottom) in terms of the vertical displacement at the cross section A-A´ shown in the previous video. The maximum vertical displacement is obtained at mid height, which seems to be in agreement with field observations.

Evolution of porosity at the centre of the first lift (top) and vertical displacement (bottom) at the cross section A-A’.


The results show that changes in porosity and permeability due to mineral dissolution/precipitation and irrigation time affect significantly the hydromechanical behavior of a heap leach pad, and vice versa. A strong advantage of the present model is that it can treat multiple chemical species simultaneously, providing solutions that are best for interpreting all the chemical interactions within a leaching system. The stability of the pad and the ore recovery efficiency can be predicted with more accuracy because the model considers the interaction between variable saturated flow, mechanical deformation and chemical reactions during the construction and operation of the heap.


I would like to express my gratitude and appreciation to Albert Nardi and Mayu Tincopa who laid the foundations for this research. I also want to acknowledge the partners of the iMaGe consortium ( for supporting the iCP development.


Nardi, A., Idiart, A., Trinchero, P., de Vries, L. M. and Molinero, J. (2013). ‘Interface Comsol-PHREEQC (iCP), an efficient numerical framework for the solution of coupled multiphysics and geochemistry’, Computers & Geosciences, submitted.

Parkhurst, D. L. and Appelo, C. A. J. (1999). ‘User's guide to PHREEQC (version 2)-A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations’, U.S. Geological Survey Water-Resources Investigations, vol Report 99–4259, 312 p.