A Posteriori Analysis of Analytical Models for Heap Leaching Using Uncertainty and Global Sensitivity Analyses

The heap leaching of minerals is one of the more commonly used processes in the mining industry. This process has been modeled by many authors. However, the validation, verification, and implementation of these models are difficult since there is uncertainty about the operating conditions and the leaching model parameters. This work uses the uncertainty quantification, based on uncertainty and sensitivity analysis, for studying the model strength against uncertainties in heap leaching. The uncertainty analysis (UA) is used to quantify the effect of the magnitude of the uncertainties of the input variables on the recovery of heap leaching. Global sensitivity analysis (GSA) is used to study the nature of connections between the recovery and input variables of the leaching model. In addition, GSA facilitates the detection of whether a leaching model is over-parameterized. The information obtained allows studying some applications of the kinetic model. The Mellado et al. kinetic model is used as an example. The UA results indicate that the kinetic model can estimate the recovery behavior considering the full range of uncertainties of input variables. The GSA indicates that the kinetic model is over-parameterized on the uncertainties range considered; this conclusion contradicts the results when the local sensitivity analysis is used. However, the model shows a good correlation between the results of GSA and the kinetic behavior of heap leaching. In addition, the kinetic model presents versatility because it allows the determination of operating regions for heap leaching.


Introduction
The hydrometallurgical process of heap leaching is one of the most commonly used processes in the mining industry.Heap leaching allows the processing of low-grade metallic minerals (for example, less than 1% copper) [1,2]; nonmetallic minerals [3,4], and potentially yttrium and heavy rare earth elements [5].In fact, heap leaching is often the preferred method for extracting metal from low-grade ore deposits as it provides a low capital cost compared to other methods [6].This low cost is due to the fact that energy-intensive comminution is not required.However, this is contrasted with often slow and inefficient recovery [7,8].In addition, small changes in metal extraction efficiency drive considerable changes in operating cash flow margins and affect project viability [9].The performance of heap leaching depends on many input variables (operational and design), i.e., its optimization is complex.A thorough review on heap leaching has been published [10], and one important conclusion is that the long-term success of heap leaching as a technology choice requires a concerted high-level cross-disciplinary engineering approach from research and development to design to construction to operation and to closure.
Heap leaching is a mineral processing technology whereby piles of crushed or agglomerate materials are leached with various chemical solutions that extract valuable minerals [11].These chemical solutions are a weak sulfuric acid solution for copper, diluted alkaline cyanide solution for gold, and water for nitrate minerals.The valuable-containing minerals are irrigated with the chemical solution which dissolves the valuable metal from the mineral, and the pregnant leach solution (PLS) passes down through the ore pile and is recovered at the base of the heap.The valuable material is extracted from the PLS using different technologies, and the chemical solution is recycled back onto the leach pile.The most common methods for recovery of valuable minerals are solvent extraction and the electro winning processes for copper, crystallization operations for nitrates, and either activated carbon or precipitation with zinc for gold.
In the literature, there is a range of scientific research on modeling and optimizing the leaching process; some examples will be briefly presented.De Andrade Lima [12] reported a model developed to simulate the transient evolution of the dissolved chemical species in the heap and column isothermal leaching processes.Bouffard and Dixon [13] published a study about heap leaching planning considering an optimization approach.Bouffard and Dixon [14] compared the rates of pore diffusion and cyanide gold dissolution in coarse, porous gold oxide ore particles.Cross et al. [15] described a computational model for the analysis of multiphase flows in reactive porous media targeted at metal recovery through stockpile leaching and environmental recovery processes.Sheikhzadeh et al. [16] numerically and experimentally studied an unsaturated flow of liquid in a bed of uniform and spherical ore particles.Galvez et al. [17] performed a computational evaluation of a phenomenological model to optimize the copper leaching.They developed an unsteady and two-dimensional model based on the mass conservation equations of the liquid phase in the bed and the particles.Mellado and Cisternas [18] presented an analytical-numerical scheme for solving the heap leaching problem, but the model still contained mathematical complexities.Bouffard and Dixon [19] reviewed heap leaching of refractory gold ores and analyzed the effect of operating parameters (e.g., flow rate, particle size, height and aeration flow).Bennett et al. [20] and McBride et al. [21] presented a comprehensive copper heap leach model based on computational fluid dynamics together with its parameterization and validation against laboratory column test data and plant data.The modeling of heap leaching also includes soft computing such as an artificial neural network.Hoseinian et al. [22] applied an artificial neural network and a genetic algorithm to predict copper recoveries and determine the optimized conditions of column leaching of a copper ore.
A good model development should include model verification, validation and uncertainty quantification [23].Model verification is used to ensure that the model is behaving properly.For example, the model can be compared with other models or with known analytical solutions.Model validation is the comparison with experimental data.Uncertainty quantification studies the effect of the uncertainties on the model.Uncertainty and sensitivity analyses can be used to provide a general overview of the effect of the uncertainties.
Several aforementioned studies helped to understand the behavior of leaching systems.However, the model is usually validated based on experimental data, and verification and uncertainty quantification are not considered in the model development.For example, Ordoñez et al. [24,25] performed an experimental validation of a phenomenological model for the leaching of caliche with the aim of assessing their level of mathematical adjustment.However, Liddell [26] indicates that a simple adjustment or comparison with experimental values is not enough to validate a leaching model, especially empirical and hybrid models.Specifically, Liddell [26] reported that a data set of titanium leaching from ilmenite (FeTiO 3 ) could be represented by four shrinking core models (spherical particles under reaction control, spherical particles under layer diffusion control, cylindrical particles under reaction control and cylindrical particles under layer diffusion control).
Mellado et al. [27] developed a model that combined aspects of phenomenology with empirical expression.Their results are comparable to the ones obtained from more complex models.However, the Mellado et al. model considered that the recovery at infinite time is equal to one, which is not observed in practice.Therefore, an expression for the recovery at infinite time was added to the previous model [28].The Mellado et al. model has been validated in several ways.First, it was used to model the recovery of copper minerals from northern Chile in pilot heaps measuring 3, 6 and 9 m in height, obtaining a determination coefficient of 99.74 [28].In other words, the Mellado et al. model was able to explain 99.74% of the variability in the recovery.Ding et al. [29] applied the Mellado et al. equations to model the heap leaching of uranium ore.The kinetic model gave the fitting precision of more than 95% and prediction precision of more than 93%.The Mellado et al. model was also verified with the results obtained from a simulation of heaps following the Dixon and Hendrix [30,31] phenomenological model, solved by using an analytic-numerical procedure as proposed in the literature [18].Different values of particle size and superficial velocity were used to adjust the Mellado et al. model, showing the fit to be very good (0.28% deviation).Then, the model was used for prediction (both interpolation and extrapolation) with very similar results (0.29% deviation).All these studies have been performed on one variable at a time (or one-factor-at-a-time, OAT).OAT has been demonstrated to be inefficient for the study of model sensitivity [32].Uncertainty quantification has not been applied to the Mellado et al. model.Saltelli et al. [33] indicated that sensitivity analysis had been claimed to be an essential ingredient in all phases of model building (model identification, model calibration, and model corroboration).Sensitivity analysis allows studying the nature of the connections between the output variables and input variables.Additionally, sensitivity analysis allows identifying if a model is over-parameterized.There are two types of sensitivity analysis: Local sensitivity analysis and global sensitivity analysis (GSA) [34,35].Local sensitivity analysis has several limitations such as linearity and normality assumptions, as well as local variations.GSA considers the entire range of input variations and does not assume linearity and normality of the input variables.
Sensitivity analysis can be another tool for model development that can be used in addition to verification and experimental validation.Mellado et al. [36] applied a local sensitivity analysis to completely evaluate an analytical heap leaching model, as previously described in the literature [27,28].Later, these results were used for the optimization of the irrigation flow rate [37].McBride et al. [21] performed a local sensitivity analysis to validate the computational fluid dynamics model proposed.Local sensitivity analysis had been applied in numerous studies of heap leaching [4,24,31,38].The most popular local sensitivity analysis practice is that of "one-factor-at-a-time".This consists of analyzing the effect of varying one model input variable at a time while keeping all others fixed.Although the shortcomings of "one-factor-at-a-time" are known from the statistical literature, it is still widely used.In fact, Salteli [32] introduced geometric proof of the inefficiency of "one-factor-at-a-time", with the purpose of providing the modeling community with a convincing and possibly definitive argument against "one-factor-at-a-time".The inadequacy of "one-factor-at-a-time" is not limited to sensitivity analysis (SA), but to uncertainty analysis (UA) as well.The elementary statistics of the model output, for instance, its maximum, or mode, can be misinterpreted via "one-factor-at-a-time".
This work aims to apply UA and GSA as tools for uncertainty quantification of heap leach models.Specifically, the goals are identification of the key inputs that influence the uncertainty of model recovery, a better understanding of the effects of inputs variables on the output variables, and accomplishing model simplification based on detection of whether the model is over-parameterized.
The paper is organized as follows.Section 2 describes the models and methodologies used.In Section 2.1, the kinetic model studied in this paper is described.In Section 2.2 the uncertainty analysis (UA) is presented.In Section 2.3, a brief description of the GSA is presented, and Section 2.4 presents the regionalizations of variables.In Section 3, the results are given together with the discussion, including UA (Section 3.1), GSA (Section 3.2), and regionalization (Section 3.3).Finally, in Section 4, some conclusions are presented.

Mellado Model
Mellado et al. [27] considered that the behavior of the heap could be modeled using a system of first-order equations such as: where y is a dynamic quantity, such as the concentration or the recovery (R t ), k τ is the kinetic constant and n τ is the order of the reaction.The subscript τ represents a time scale that depends on the phenomenon to model.To solve Equation ( 1), an initial condition is required.Then, Mellado et al. introduced a delay, i.e., a time ω where the R t begins to change (R t (ω) = 0).The general solution to this problem is known (see Mellado et al. [27] for the general solution), for n τ = 1 the solution is: Dixon and Hendrix [30,31] considered that the leaching phenomenon occurs at different scales of size and time, and that different phenomena participate in the leaching process.It is possible to represent these phenomena using Equation ( 2) with expressions associated with the particle and height of heap leaching.At the particle level, they considered that the phenomenon was dominated by the diffusion of the dissolved products from the mineral particles to the bulk of the leaching solution.Then, they defined a dimensionless diffusion time [30,31]: where D Ae is the effective diffusivity of the solute within the particle pores (cm 3 /cm•s), ε o is the porosity of the particle, t is the time (s), r the particle radius (cm).At the bulk level, Dixon and Hendrix considered that the heap is represented by the porous bulk formed of particles through which the leaching solution flows at a constant rate.Then, they defined the dimensionless time of the bulk [30,31] where u s is a superficial velocity of lixiviant flow through the bed (cm 3 /cm•s), ε b is the volumetric fraction of the bulk solution in the bed, and Z is the height of the heap (cm).Then, Equation ( 2) can be written for both dimensionless times as [27,28]: where ω represents the delay in the scale of τ.Following the phenomenological model of Dixon and Hendrix [30,31], it is possible to demonstrate that [27] Minerals 2018, 8, 44 5 of 16 For the inclusion of both scales in the processing of the heap, it was assumed that the total recovery is the sum of the recoveries on both scales, that is R = R τ + R θ .In addition, it was considered that that each recovery has an asymptotic behavior with increasing time, this is Therefore, the total recovery is: where λ is defined as the kinetic weight factor.In addition, Mellado et al. [28] developed a mathematical expression for R ∞ : This new expression shows that the R ∞ is influenced by the radius of the particle (r), height of the heap leaching (Z) and three coefficients of mathematical adjustment (α, β and γ).
Finally, the heap leaching model is given by [28]:

Uncertainty Analysis
The UA corresponds to determining the uncertainty in the output variables as a result of uncertainty in the input variables.The UA is usually done using the probability theory where the uncertainty is represented by probability distribution functions (PDF), but there are other available theories such as the fuzzy theory.In the probability theory, the process is performed in four steps.First, the type of PDF and the magnitude of uncertainty is determined for each input variable; in other words, the input uncertainty is characterized.Then, for each input variable, a sample is generated from the PDF.This process is usually performed using random sampling, but other methods such as Latin hypercube sampling are available when small samples are required.Thirdly, values of the output variables are determined for each element of the sample.If the model is costly, this stage may be demanding regarding computational calculations.Finally, the results are analyzed using graphics, descriptive statistics, and statistical tests to characterize the behavior of the output variables.
As commented before, the uncertainties of input variables can be represented through a PDF.When the input variables have epistemic uncertainties, the uncertainty can be represented by a uniform distribution.Design and operating variables present this type of uncertainty.When input variables have stochastic uncertainties, the normal distribution is usually used to represent this type of uncertainty.

Sensitivity Analysis
The SA can be defined as the study of how the uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input [39].There are two types of SA: Local sensitivity analysis and GSA.The first quantifies the rate of change of the model output due to small variations in the uncertainties of the input variables.One of the disadvantages of local sensitivity analysis is that the choice of the assessment point could sharply influence the results, especially when the sampling space of the input variables is affected by a considerable uncertainty.However, GSA considers the full range of uncertainties of the input variables represented by PDF.
GSA techniques have been widely used in various engineering areas and are of great importance in determining the most significant variables in a model.The general objectives of GSA are [40,41] Minerals 2018, 8, 44 6 of 16 (a) the identification of the significant and insignificant input variables and the possible reduction of the dimensions of an optimization problem; (b) the improvement in the understanding of the model behavior (highlighting interactions among factors and finding combinations of factors that result in high or low values for the model output); (c) mapping the output behavior in function of the inputs by focusing on a specific domain of inputs if necessary; and (d) calibrating some model inputs using some available information.
There are several methods to perform the GSA, but among these methods, those based on variance decomposition stand out, due to their versatility and efficiency [42].Therefore, in this study, we will use such methods, specifically those based on the Sobol method [43].In variance-based methods, the variance V of a model Y = f (X 1 , X 2 , . . . ,X k ) is fully decomposed into terms depending on the k uncertain input variables and their interactions.The number of interactions increases exponentially with the number of input variables, which is why all interactions are rarely determined.More often practitioners are satisfied with computing just k first order effects (S i ) and k total effects (ST i ), the latter describing synthetic interactions among input factors.The first order sensitivity index measures only the main effect contribution of each input variable on the output variance.It does not take into account the interactions among variables.Two factors are said to interact if their total effect on the output is not equal to the sum of their first order effects.A model without interactions is said to be additive.The first order indices sum to one in an additive model with orthogonal inputs.For non-additive models, information from all the interactions as well as the first order effect is searched for.
The first-order sensitivity index (S i ) is important when the objective is to determine the most important input uncertainties.The total sensitivity index (ST i ) is important when the objective is to reduce the uncertainty in the output model.If the first-order sensitivity index of the ith input factor is negligible, the uncertainty in X i does not affect the uncertainty in the output model Y.Therefore, X i is non-influential or unimportant.This does not determine any information about input interactions or high-order sensitivity indices.If the total sensitivity index (ST i ) is also small, then apart from being unimportant, X i does not interact with other factors (high-order effects of X i are negligible).The implication of small values of S i and ST i is that the uncertainty in X i has no effect on the uncertainty in Y. Thus, in subsequent analyses, X i can be fixed to its nominal value (mean or median) and further research, measurement, analysis and data gathering can be directed at other factors.Conversely, regardless of the magnitude of ST i , a large value of the first-order sensitivity index S i implies that X i is influential.The arithmetic difference between ST i and S i indicates the magnitude of the interactions between X i and other factors.

Regionalization of Input Variables
The regionalization is based on dividing the response of the model into two subsets according to some criteria [44].Here, one set brings together the values of input variables that provide the desired outcome (group B) and another set brings together the undesired outcome (B), i.e., for each input variable X i of model where R is the total number of simulations and m is the total number of input variables of the model.To apply the regionalization, one should first identify the effect of the influential input variables on the responses of the model.The goal of the regionalization is to identify operational regions of the process, as described in Lucay et al. [45].

Results and Discussion
The results are arranged in three subsections.UA results, which address the UA at 30 and 100 days of leaching.GSA results, which were obtained using the Sobol-Jansen method [32,46], included in the package "Sensitivity" [47] of R [48].Finally, the third part considered the application of regionalization analysis.
A good distribution function to represent epistemic uncertainties is uniform PDF.Thus, uniform PDF with the following magnitudes were used: 3%, 5%, 3%, 5%, 3%, 5%, 5%, 5% and 5% for the kinetic weight factor, volumetric fraction of the solution of the bulk, height of heap, delay, kinetic constants, effective diffusivity, porosity of mineral, radius of particle, and superficial velocity, respectively.The UA was analyzed twice: after 30 and 100 days of leaching.The results are shown in Figures 1 and 2.
Minerals 2018, 8, x FOR PEER REVIEW 7 of 16 A good distribution function to represent epistemic uncertainties is uniform PDF.Thus, uniform PDF with the following magnitudes were used: 3%, 5%, 3%, 5%, 3%, 5%, 5%, 5% and 5% for the kinetic weight factor, volumetric fraction of the solution of the bulk, height of heap, delay, kinetic constants, effective diffusivity, porosity of mineral, radius of particle, and superficial velocity, respectively.The UA was analyzed twice: after 30 and 100 days of leaching.The results are shown in Figures 1 and 2.  Figures 1 and 2 show that the leaching time affects the leaching kinetic uncertainty.The histograms show that the uncertainty in the input variables produces a greater uncertainty in the recovery after 30 days compared to the uncertainty after 100 days.In fact, the difference between the highest value and the lowest recovery value after 30 days is approximately double that observed after 100 days.The test QQ-plot indicates that when the leaching time is equal to 30 days, the recovery presents a normal PDF with a standard deviation of 1.25 and a mean of 55.5% of recovery.The test QQ-plot indicates that when the leaching time is equal to 100 days, the recovery moves away from normal behavior.These results show that the kinetic model can estimate the recovery of heap leaching considering the full range of uncertainties of input variables.The change of recovery A good distribution function to represent epistemic uncertainties is uniform PDF.Thus, uniform PDF with the following magnitudes were used: 3%, 5%, 3%, 5%, 3%, 5%, 5%, 5% and 5% for the kinetic weight factor, volumetric fraction of the solution of the bulk, height of heap, delay, kinetic constants, effective diffusivity, porosity of mineral, radius of particle, and superficial velocity, respectively.The UA was analyzed twice: after 30 and 100 days of leaching.The results are shown in Figures 1 and 2.  Figures 1 and 2 show that the leaching time affects the leaching kinetic uncertainty.The histograms show that the uncertainty in the input variables produces a greater uncertainty in the recovery after 30 days compared to the uncertainty after 100 days.In fact, the difference between the highest value and the lowest recovery value after 30 days is approximately double that observed after 100 days.The test QQ-plot indicates that when the leaching time is equal to 30 days, the recovery presents a normal PDF with a standard deviation of 1.25 and a mean of 55.5% of recovery.The test QQ-plot indicates that when the leaching time is equal to 100 days, the recovery moves away from normal behavior.These results show that the kinetic model can estimate the recovery of heap leaching considering the full range of uncertainties of input variables.The change of recovery behaviour over time is due to recovery tending to present asymptotic behaviour.Figures 1 and 2 show that the leaching time affects the leaching kinetic uncertainty.The histograms show that the uncertainty in the input variables produces a greater uncertainty in the recovery after 30 days compared to the uncertainty after 100 days.In fact, the difference between the highest value and the lowest recovery value after 30 days is approximately double that observed after 100 days.The test QQ-plot indicates that when the leaching time is equal to 30 days, the recovery presents a normal PDF with a standard deviation of 1.25 and a mean of 55.5% of recovery.The test QQ-plot indicates that when the leaching time is equal to 100 days, the recovery moves away from normal behavior.These results show that the kinetic model can estimate the recovery of heap leaching considering the full range of uncertainties of input variables.The change of recovery behaviour over time is due to recovery tending to present asymptotic behaviour.

Global Sensitivity Analysis
The realization of GSA considered six cases.The first case applies GSA to the kinetic model using the nominal values and the uncertainties used in UA, with leaching time between one and 100 days.In cases two to six, the GSA was realized by modifying the nominal value of one input variable and maintaining the others in their nominal values.The same uncertainties were used in all cases.The second case considered a height of the heap of 300 cm and 900 cm.The third case considered a particle radius of 1.5 cm and 3.5 cm.The fourth case considered a superficial velocity of 9.0 cm 3 /cm 2 •d and 57.0 cm 3 /cm 2 •d.The fifth case considered an effective diffusivity 0.06 cm 3 /cm 2 •d and 0.112 cm 3 /cm 2 •d.The sixth case considered a porosity of the particle of 0.01 and 0.06.
Figure 3 shows the results of the first case.It can be observed that u s and ε b are influential on recovery throughout the leaching time while the influence of Z on recovery increases over time.The influence of λ on recovery increases until approximately day 18 and then decreases over time.The variable k_θ presents a moderate influence on recovery, while input variables r and ω present a low influence on recovery.The remaining input variables do not influence the recovery.These GSA results allow a reduction of the number of input variables of the kinetic model.For example, if non-significant variables are set at their nominal values, Equation ( 11) can be rewritten as: Evidently, these results depend on the uncertainty range of the input variables of the kinetic model.It is important to indicate that these results contradict the results when local SA is applied, which concludes that the model of Equation (11) does not have extra variables [36].
Minerals 2018, 8, x FOR PEER REVIEW 8 of 16 9.0 cm 3 /cm 2 •d and 57.0 cm 3 /cm 2 •d.The fifth case considered an effective diffusivity 0.06 cm 3 /cm 2 •d and 0.112 cm 3 /cm 2 •d.The sixth case considered a porosity of the particle of 0.01 and 0.06. Figure 3 shows the results of the first case.It can be observed that and are influential on recovery throughout the leaching time while the influence of on recovery increases over time.The influence of on recovery increases until approximately day 18 and then decreases over time.The variable k_θ presents a moderate influence on recovery, while input variables and present a low influence on recovery.The remaining input variables do not influence the recovery.These GSA results allow a reduction of the number of input variables of the kinetic model.For example, if non-significant variables are set at their nominal values, Equation ( 11) can be rewritten as: Evidently, these results depend on the uncertainty range of the input variables of the kinetic model.It is important to indicate that these results contradict the results when local SA is applied, which concludes that the model of Equation (11) does not have extra variables [36].Figures 4 and 5 show the results for the second case; this is for a heap height of 300 and 900 cm.These figures show that when the heap height increases, the of the heap height decreases.These results agree with the kinetic behavior of heap leaching shown in Figure 6.This figure shows that when the height of heap decreases, the recovery is more sensitive to changes in the heap height.Figures 4 and 5 show the results for the second case; this is for a heap height of 300 and 900 cm.These figures show that when the heap height increases, the ST i of the heap height decreases.These results agree with the kinetic behavior of heap leaching shown in Figure 6.This figure shows that when the height of heap decreases, the recovery is more sensitive to changes in the heap height.
Figures 4 and 5 show the results for the second case; this is for a heap height of 300 and 900 cm.These figures show that when the heap height increases, the of the heap height decreases.These results agree with the kinetic behavior of heap leaching shown in Figure 6.This figure shows that when the height of heap decreases, the recovery is more sensitive to changes in the heap height.Figures 3, 7 and 8 show the index when the particle radius is 2.5, 1.5 and 3.5 cm.The effect of is small and only at times less than 50 days, which is consistent with what was observed in the recovery (see Figure 6).The effect of increases as its value increases and it produces a change in the index of the other input variables.This occurs because the uncertainties were considered as    Figures 3, 7 and 8 show the index when the particle radius is 2.5, 1.5 and 3.5 cm.The effect of is small and only at times less than 50 days, which is consistent with what was observed in the recovery (see Figure 6).The effect of increases as its value increases and it produces a change in the index of the other input variables.This occurs because the uncertainties were considered as Figures 3, 7 and 8 show the ST i index when the particle radius is 2.5, 1.5 and 3.5 cm.The effect of r is small and only at times less than 50 days, which is consistent with what was observed in the recovery (see Figure 6).The effect of r increases as its value increases and it produces a change in the ST i index of the other input variables.This occurs because the uncertainties were considered as a percentage of the nominal value.Then, at a higher nominal value of a variable, greater uncertainty was considered.This shows that the level of uncertainty plays an important role and that its value must be carefully determined.In all cases analyzed, it is observed that the number of input variables of the kinetic model can be reduced.These reductions depend on the uncertainties range of the input variables.This result is different to the one obtained when local SA is applied [36].Therefore, it is recommended that GSA must be used and not local SA.In all cases analyzed, it is observed that the number of input variables of the kinetic model can be reduced.These reductions depend on the uncertainties range of the input variables.This result is different to the one obtained when local SA is applied [36].Therefore, it is recommended that GSA must be used and not local SA.In all cases analyzed, it is observed that the number of input variables of the kinetic model can be reduced.These reductions depend on the uncertainties range of the input variables.This result is different to the one obtained when local SA is applied [36].Therefore, it is recommended that GSA must be used and not local SA.

Regionalization
The information obtained with GSA can be used to study some applications of the model.For example, to know the operating conditions necessary to obtain a recovery higher than 73.5% (obtained when the leaching time is approximately equal to 35 days, and heap height is 300 cm).The regionalization of the input variables can give some orientations.The results of GSA at day 35 indicate that the regionalization must be done only for the height of heap, superficial velocity and volumetric fraction of the bulk solution in the bed; the other input variables can be fixed in their current operational values (nominal values).This reduction of the heap leaching model helps to obtain operating regions more defined and easier to interpret.The groups B (desired outcome) and B (undesired outcome) are defined as follows.
where E is the space formed by the uncertainties of three influential input variables, while R(x) is the recovery of heap leaching in the point x.The regions B and B are shown in Figure 15.
For example, Figure 15 (left) shows the values of the height of heap, superficial velocity, and volumetric fraction of the bulk solution in the bed when the desired outcome is obtained (recovery higher than 73.5%).In fact, if the heap leaching operating conditions are Z = 300 cm, u s = 27 cm 3 /cm 2 •s and b = 3.0% (region B), the recovery calculated with the Mellado et al. model is 72.7%, but if Z = 300 cm, u s = 29 cm 3 /cm 2 •s and b = 2.9% (region B), the recovery calculated with the Mellado et al. model is 74.1%.The regionalization analysis can be complemented with planes like the ones shown in Figure 16.
Figure 16 shows the plane Z = 300 cm including both regions (B ∪ B).This figure allows us to clearly operate the heap leaching when the leaching time is equal to 35 days.For example, when u s = 28 cm 3 /cm•d, b begins to increase.Figure 16 shows that the recovery of heap leaching decreases.Therefore, the regionalization of influential input variables allows us to control the uncertainty of recovery of heap leaching.This application indicates that the kinetic model presents versatility.
Minerals 2018, 8, x FOR PEER REVIEW 13 of 16 For example, Figure 15 (left) shows the values of the height of heap, superficial velocity, and volumetric fraction of the bulk solution in the bed when the desired outcome is obtained (recovery higher than 73.5%).In fact, if the heap leaching operating conditions are = 300 cm, = 27 cm /cm • s and = 3.0% (region ), the recovery calculated with the Mellado et al.
model is 72.7%, but if = 300 cm, = 29 / • s and = 2.9% (region ), the recovery calculated with the Mellado et al. model is 74.1%.The regionalization analysis can be complemented with planes like the ones shown in Figure 16.
Figure 16 shows the plane = 300 cm including both regions ( ∪ ).This figure allows us to clearly operate the heap leaching when the leaching time is equal to 35 days.For example, when = 28 cm /cm • , begins to increase.Figure 16 shows that the recovery of heap leaching decreases.Therefore, the regionalization of influential input variables allows us to control the uncertainty of recovery of heap leaching.This application indicates that the kinetic model presents versatility.

Conclusions
This work shows that UA and GSA are useful tools to test the model robustness against uncertainties in leaching models.The kinetic model proposed by Mellado et al. [27,28] was used as an example.The UA indicates that the kinetic model can estimate the recovery behavior considering the full range of uncertainties of input variables.Different types of PDF for recovery were observed after 35 and 100 days of operation.In addition, the uncertainty in the recovery at day 100 was much less than the uncertainty at day 35.The GSA indicates that the kinetic model is over-parameterized on the uncertainties range considered; this conclusion contradicts the results when local SA is used.The level of uncertainty affects the results, which is why the range of uncertainty must be determined with some degree of certainty.Additionally, the kinetic model presents versatility because it allows determining operating regions for controlling the uncertainty in the recovery of the leaching process.The regionalization of the most significant variables (higher values of ST i ) can be a tool to determine operating conditions or understand the relationship between variables that produce undesired or desired effects.

Figure 1 .
Figure 1.Uncertainty analysis at 30 days of leaching.

Figure 2 .
Figure 2. Uncertainty analysis at 100 days of leaching.

Figure 1 .
Figure 1.Uncertainty analysis at 30 days of leaching.

Figure 1 .
Figure 1.Uncertainty analysis at 30 days of leaching.

Figure 2 .
Figure 2. Uncertainty analysis at 100 days of leaching.

Figure 2 .
Figure 2. Uncertainty analysis at 100 days of leaching.

Figure 3 .
Figure 3. Global sensitivity analysis (GSA) of the model with a heap height equal to 600 cm.

Figure 3 .
Figure 3. Global sensitivity analysis (GSA) of the model with a heap height equal to 600 cm.

Figure 4 .
Figure 4. GSA of the model with a heap height equal to 300 cm.Figure 4. GSA of the model with a heap height equal to 300 cm.

Figure 4 . 16 Figure 5 .
Figure 4. GSA of the model with a heap height equal to 300 cm.Figure 4. GSA of the model with a heap height equal to 300 cm.Minerals 2018, 8, x FOR PEER REVIEW 9 of 16

Figures 3 -
Figures 3-5 also show that the relative weight of the significant variables changes if the nominal values of the height of the heap change.To complement the GSA results, Figure6shows the recovery versus time for all the cases studied.

Figure 5 .
Figure 5. GSA of the model at a heap height of 900 cm.

Figures 3 - 16 Figure 5 .
Figures 3-5 also show that the relative weight of the significant variables changes if the nominal values of the height of the heap change.To complement the GSA results, Figure6shows the recovery versus time for all the cases studied.

Figures 3 -
Figures 3-5 also show that the relative weight of the significant variables changes if the nominal values of the height of the heap change.To complement the GSA results, Figure6shows the recovery versus time for all the cases studied.

Figure 8 .
Figure 8. GSA of the model at radius 3.5 cm.

Figures 9 -
Figures 9-14 show the results of the other cases.Figures 3, 9 and 10 show the index when the is 9.0, 28.0 and 57.0 cm 3 /cm 2 •s.Figures 3, 11 and 12 show the index when the is 0.0864, 0.06 and 0.112 cm 3 /cm 2 •d.Figures 3, 13 and 14 show the index when the is 0.03, 0.01 and 0.06.The effect of the nominal value and the uncertainty is important in the results observed.

Figure 8 .
Figure 8. GSA of the model at radius 3.5 cm.

Figures 9 -
Figures 9-14 show the results of the other cases.Figures 3, 9 and 10 show the index when the is 9.0, 28.0 and 57.0 cm 3 /cm 2 •s.Figures 3, 11 and 12 show the index when the is 0.0864, 0.06 and 0.112 cm 3 /cm 2 •d.Figures 3, 13 and 14 show the index when the is 0.03, 0.01 and 0.06.The effect of the nominal value and the uncertainty is important in the results observed.

Figure 8 .
Figure 8. GSA of the at radius 3.5 cm.

Figures 9 -
Figures 9-14 show the results of the other cases.Figures 3, 9 and 10 show the ST i index when the u s is 9.0, 28.0 and 57.0 cm 3 /cm 2 •s.Figures 3, 11 and 12 show the ST i index when the D Ae is 0.0864, 0.06 and 0.112 cm 3 /cm 2 •d.Figures 3, 13 and 14 show the ST i index when the ε o is 0.03, 0.01 and 0.06.The effect of the nominal value and the uncertainty is important in the results observed.In all cases analyzed, it is observed that the number of input variables of the kinetic model can be reduced.These reductions depend on the uncertainties range of the input variables.This result is different to the one obtained when local SA is applied[36].Therefore, it is recommended that GSA must be used and not local SA.

Figures 9 -
Figures 9-14 show the results of the other cases.Figures 3, 9 and 10 show the index when the is 9.0, 28.0 and 57.0 cm 3 /cm 2 •s.Figures 3, 11 and 12 show the index when the is 0.0864, 0.06 and 0.112 cm 3 /cm 2 •d.Figures 3, 13 and 14 show the index when the is 0.03, 0.01 and 0.06.The effect of the nominal value and the uncertainty is important in the results

Figure 11 .
Figure 11.GSA of the model at D Ae 0.06 cm 3 /cm 2 •d.

Figure 14 .
Figure 14.GSA of the model at 0.06.

Figure 14 .
Figure 14.GSA of the model at 0.06.

Figure 14 .
Figure 14.GSA of the model at ε o 0.06.

Figure 15 .
Figure 15.Regionalization of the influential input variables Z, u s and b , B (left) and B (right).