**Multivariate block simulations of a lateritic nickel deposit and post-processing of a representative subset**

**J. Deraisme ^{I}; O. Bertoli^{I}; P. Epinoux^{II}**

^{I}Geovariances (Pty) Ltd, France

^{II}Centre Minier de Tiébaghi, Nouvelle, Calédonie

**SYNOPSIS**

Société le Nickel (SLN) exploits the Dome lateritic nickel orebody at its Tiébaghi operations in New Caledonia. The site geology has an obvious bearing on recovery and metallurgical performance, and the controls are extremely complex at all scales: lithology, oxidation-reduction conditions, mineralogy, and multivariate geochemistry. In that context, establishing an adapted recoverable resource estimation method that is efficient and transparent proves an interesting challenge. The method must address the key notions of: ]]>
> Support effect: the exploration data-set can only warrant the estimation of large panels 20×20×3 m^{3}, which are much larger than the selective mining units (SMUs) of 5×5×3 m^{3}

> Information effect: the selection at production stage uses estimates based on information much denser than that available for mine planning.

The problem is rendered more complex in this type of deposit by the fact that the additive variables used for any estimation process (metal accumulation and ore tonnage) are not the ones used to establish the selection at mining stage *(i.e.* the nickel grade of SMUs, which is the ratio of accumulation on tonnage). Eventually, the selection criteria involve not only Ni grades but also other constituents such as Al_{2}O_{3}, Fe_{2}O_{3}, MgO, and SiO_{2}.

The solution presented in this paper is the construction of a platform of SMU multivariate simulations in the saprolitic horizon, together with efficient post-processing aimed at producing the multivariate recoverable resource estimates. The efficiency constraints imposed by the number of blocks to be simulated (1.5 million) lead to resorting to direct block simulations. The underlying multivariate discrete Gaussian model (DGM), the validity of which is tested beforehand, is put to good use to mimic the selection process at the mining stage, by offering the ability to simulate in each block a composite value located at random within the block.

Finally, the paper presents the application of a scenario reduction algorithm to pick a representative subset of a few simulations to help appraise the risk attached to the downstream phases (reserve optimization, mine sequencing) of the project. The implementation presented here is for the first generation of the Scenario Reduction plugin built in Isatis software, where measuring the distance between the initial set of scenarios and the reduced set is based on Ni quantities only. Research is under way to adapt the measure to the true multivariate nature of the problem.

**Keywords:** direct block simulations, discrete Gaussian model, multivariate, scenario reduction.

**Introduction**

Société le Nickel (SLN) exploits the Dôme lateritic nickel orebody at its Tiébaghi operations located 385 km north of Nouméa on the western coastline of the Northern Province of New Caledonia (see Figure 1). Tiébaghi is SLN's largest active mining complex and relies on a fairly advanced UTM *(usine de traitementminerai:* ore treatment plant) to increase the nickel grade of marginal ore in order to augment production of traditional ore ('tradi') for shipment to the Doniambo pyrometallurgical processing plant in Nouméa. The geochemistry of the ore also impacts the nickel smelting and recovery processes at the Doniambo plant.

Owing to the confidential nature of the UTM design, all economic and technical parameters relevant to the ore resource characterisation process have been altered or disguised prior to publication of this paper.

The deposit geology significantly affects the recovery and metallurgical performance at the UTM, with controls that are extremely complex at all scales: lithology, oxidation-reduction conditions, mineralogy, and multivariate geochemistry.

It is within that context that SLN decided in 2011 to improve their recoverable resource estimation method to ensure robust models that would be produced capable of handling:

*> The traditional support effectÂ—the* exploration data-set can only warrant the estimation of large panels 20×20× 3 m^{3}, which are much larger than the selective mining units (SMUs) of 5×5× 3 m^{3}

*> The information effect*Â—the selection at production stage uses estimates based on information much denser than that available for mine planning

*> The problem of non-additivity of the selection variableÂ—the* additive variables used for estimation (metal accumulation and ore tonnage) are not the ones used to establish the selection at the mining stage (*i.e.* the SMU nickel grade, which is the ratio of accumulation on tonnage)

*>*The multivariate nature of the ore control selection criteria, which involve not only Ni but also other constituents such as Al

_{2}O

_{3}, Fe

_{2}O

_{3}, MgO, and SiO

_{2}.

The solution presented in this paper is to construct a platform of SMU multivariate simulations in the saprolitic horizon (which represents 90% of the recoverable resource at Dome) and construct an efficient post-processing method aimed at producing the multivariate recoverable resource estimates used as inputs to the pit optimization and mine sequencing phases.

**Regional and mine geology at Tiébaghi**

Nickel laterites represent nearly 60% of the world's Ni resources. Tiébaghi mine is located in the north of New Caledonia in the South Pacific. This island is known to contain about 25% of the world's nickel resources, developed from lateritic weathering profiles.

About one-third of New Caledonia's is covered by ultramafic rocks, mainly peridotites, emplaced as an ophiolite nappe by obduction of the oceanic mantle over the New Caledonia - Norfolk continental ridge during the late Eocene orogenic event (Cluzel *et al.,* 2001). A lateritic profile overlies the ultramafic rocks (see Figure 2). The distribution of ultramafic rocks (a large southeastern massif and smaller ones along the northwestern coast, with various elevations and a stepped lateritic surface) is the result of the combined effects of isostatic post-orogenic uplift, tectonics, and erosion.

]]> The Tiébaghi massif is a klippe of peridotite overlying Cretaceous basalts and Eocene-Paleocene sedimentary rocks. The massif forms a gently west-dipping plateau culminating at 584 m. The study area, called Dome, is located in the central-western part of the plateau and displays a complete and thick (30-80 m) lateritic weathering mantle.

Four main horizons make up the lateritic regoliths (see Figure 3), from the base to the top: rocky saprolites (isovolumic weathering), saprolites (isovolumic weathering), limonite (oxide of iron, closed porosity), and ferruginous/aluminous duricrust (Eggleton, 2001). Nickel is concentrated in secondary oxides and hydrous phylosilicates. Mineralization is controlled by the interplay of lithology, tectonics, climate, and geomorphology, and developed in several phases (Pelletier, 2003).

Mining operations at Tiébaghiaim are conducted to selectively extract three types of ore:

>* *The Ni-Co rich limonites, which are stockpiled for subsequent processing

>* *The intrinsically nickel-rich saprolites and saprolites that cannot be treated in the UTM because of low enrichment capacity or mineralogy preventing the solid-liquid separation ('Tradi' ore)

>* *The saprolite beneficiated by wet treatment in the UTM. These saprolites have undergone a sufficient degree of weathering to enable economic recovery and removal of gangue mineral by the UTM (UTM ore).

>* *The lateritic leaching processes induce intense chemical dissolution of the peridotite, as magnesia and part of silica are progressively leached while Al, Fe, Cr, Co, and Ni accumulate in the weathering mantles (Tardy, 1997; Nahon *et al.*, 1992; Anand and Paine, 2002). The ratio between iron-alumina and magnesia therefore indicates the intensity of weathering

>* *In addition to a high variability of peridotites with different compositions (plagioclase lherzolites, spinel dunite, harzburgite), the Tiébaghi massif exhibits a high degree of serpentinization, a hydrothermal phenomenon probably related to felsic and gabbro intrusions. These rocks introduce a significant amount of aluminum into the weathering profile. The presence of secondary aluminous phylosilicates such as kaolinite in the saprolites also affects the beneficiation characteristics at UTM

>* *Secondary silica may be present (chalcedony and quartz), resulting in the dilution of *in situ* nickel. Such silicification is due to the weathering of peridotite in which the low alumina content has restricted the formation of clays. Removal of silica by the UTM is an essential beneficiation step, increasing the recovered grade

>* *Structural phenomena and alteration have also led to the formation of a karst system, with alteration by groundwater having induced the formation of reduced horizons (gley) with specific mineralogical and chemical signatures (presence of carbonate and sulphides). These minerals affect the value and potential treatment of saprolites in the UTM.

Mineral processing and mineralogical studies have thus helped define the geochemical classification criteria of the nickel ore for Dome:

>* *The ratio (Fe_{2}O_{3} + Al_{2}O_{3}/MgO (IA_{2}) characterizes the degree of weathering, which affects the processing capacity and mass efficiency of the UTM

>* *Anomalous concentrations of alumina (Al_{2}O_{3}), iron (Fe_{2}O_{3}), and silica (SiO_{2}) are used to understand the mineralogy and indirectly predict nickel output from the UTM and subsequently the pyrometallurgical processing performance of the Doniambo plant.

It is therefore essential to estimate the nickel grade, together with all the secondary variables involved in the above calculations. These grades being non-additive, the variables to consider in any estimation process are Ts/m^{3} (proportion of altered material χ density) and accumulated variables i.e. Ni*Ts/m^{3}, Al_{2}O_{3}*Ts/m^{3}, Fe_{2}O_{3}*Ts/m^{3}, MgO*Ts/m^{3}, and SiO_{2}*Ts/m^{3}. The raw variables are then deduced by the ratios of the estimated accumulations on the estimated Ts/m^{3}. The grades are non-additive because they apply to the proportion of altered material for each composite, and even if these composites are of constant length (3 m), the proportion affected by alteration is variable. Consequently, the support over which these grades are measured becomes variable (Journel and Huijbregts, 1978).

The objective of the procedure developed for Tiébaghi is to estimate the recoverable proportions of 5×5×3 m^{3} SMUs within each 20×20×3 m^{3} panel when a cut-off is applied on Ni, and to attach to that proportion the grade of all five relevant elements being studied. An original procedure, based on direct block simulations and which is capable of handling the impacts of information and support effects as well as the non-additivity of the input variables, is presented below.

**Direct block simulation**

Block simulation may be performed either by averaging of point simulations or by a direct method. Direct block simulations are obtained by first generating block values that reproduce a predefined block variogram model, and then conditioning the block values directly from the point information. The obvious advantage of direct simulation is the gain in CPU performance and disk space. Owing to the efficiency constraints imposed by the number of blocks to be simulated at Tiébaghi (1.5 million), direct block simulation is preferable to re-blocking of point simulations.

The direct implementation can be performed in the framework of the discrete Gaussian model (DGM), where the deposit is partitioned into blocks of size *v* (SMUs), and each sample point (*x*) is considered random within each block, allowing a statistical characterization of the link between point and block information to be utilized for nonlinear estimation and simulations (Deraisme *et al.,* 2008).

In DGM, the variogram model attached to the block Gaussian equivalent values *Y _{v}* may be deduced from the variogram of point Gaussian values

*Y*regularized on blocks (Emery and Ortiz, 2005). The conditioning step is then achieved by simple kriging (SK) in the DGM framework, where

_{x}*Y*is obtained by linear regression from

_{x}*Y*(the correlation coefficient being the change of support coefficient

_{v}*r*that is linked to the physical loss of variability when changing the support of the distributions from point to block). This property of DGM can also be used to obtain a point simulated value within each simulated block value. The simulated points may be viewed as samples of the simulated reality, which can be used to mimic the grade control process by re-estimating the SMU values and comparing the resulting estimates to the simulated values. The only loss of optimality in that process concerns the location of the simulated point value, which is considered random within the block. The extension of the above framework to the multivariate setting is presented in Deraisme

*et al.*(2008).

One must note that the conditioning by SK imposes a very strict adherence to the stationarity hypothesis for the different grade distributions. The stationarity assumptions are well respected horizontally for all distributions, but obvious departures are observed along the vertical direction (see Figure 4), in keeping with the differential leaching of material in lateritic deposits.

This prompted tests for whether the conditioning step in the simulation algorithm (see the following section on Production and validation of multivariate simulations for the Dome orebody) would benefit from using ordinary kriging (OK), which is less sensitive to departures from ideal stationarity conditions, instead of the theoretically sound choice of SK. The tests concluded that simulations obtained via the two conditioning strategies were highly correlated realization-per-realization (0.95), and that SK-based realizations offer a better reproduction of the input variogram model as well as a better reproduction of the input statistical profile. SK was thus vindicated as the algorithm of choice for conditioning the simulations.

As for testing the applicability of the DGM, checking ratios of grade indicator variograms can theoretically test the applicability of the diffusive model (Rivoirard, 1994; Chiles *et al.,* 1999). Tests carried out on metric passes in the saprolitic horizon show that the ratios of cross-variograms to variograms of indicators increase with distance, confirming the applicability of the DGM. Moreover, the traditional test involving the ratio of the square root of the variogram to the first-order variogram (madogram) for Gaussian Ni grades further supports the plausibility of the bi-Gaussianity assumptions at play in the model (see Figure 5).

]]>

**Production and validation of multivariate simulations for the Dôme orebody**

Several tests were first carried out on bivariate simulations of Ts/m^{3} and Ni*Ts/m^{3} to help calibrate some key simulation parameters (number of bands in the turning band algorithm used for producing the non-conditional simulations, definition of the neighbourhood search parameters), and define the overall strategy to simulate the complete variable set. The number of simulations to be produced was a particular concern from a pragmatic viewpoint, as the implications of producing 50, 100, or 200 simulations are certainly not trivial from a performance viewpoint.

Initial tests were performed in the central area (see Figure 6) of the deposit by producing 100 simulations and testing how post-processing statistics of the first lot of 50 realizations compared to the statistics of the second lot of 50 simulations.

One of the comparisons performed is illustrated in Table I, which shows the average estimation error incurred in estimating the portion of the deposit above the economic cutoff for the different lot of simulations sampled at different drill spacings. The fact that the first lot of 50 realizations strictly reproduces the statistics for the second lot of 50 suggests that for the global characterization of the resource at Dome, working with 50 realizations is acceptable.

]]>

The initial step for the simulation of the six studied variables is to implement their histogram modelling via Gaussian anamorphosis. The histograms reconstructed by the anamorphosis functions offer a satisfactory fit of the weighted experimental histograms. The weighting of the distribution is mandatory as the data-set mixes different sampling grids (58×26 m^{2} for exploration, 20×20 m^{2} centred for mine planning, and 5×5 m^{2} for grade control) obtained via a cell declustering technique using a window of 20×20 m rotated at 45 degrees in keeping with the orientation of drilling (corresponding to step 4 in Figure 7).

Experimental variography is then performed on the Gaussian transforms in the two main directions of continuity in the horizontal plane (azimuth 135° and azimuth 45°) and along the vertical direction. A model is fitted using the linear model of coregionalization (LMC). The LMC is a generalization of the nested variogram model to the multivariate case. In this model, all simple and cross-variograms are linear combinations of the same elementary components. This can be interpreted by a decomposition of the variables themselves into linear combinations of independent random functions. The quality of the fitting obtained for the main variables (Ts/m^{3} and Ni*Ts/m^{3}) is shown in Figure 8.

^{3}support, the DGM change of support parameters are determined (see Deraisme

*et al.,*2008).

Fifty direct 5×5×3 m^{3} block co-simulations of Ts/m^{3}, Ni*Ts/m^{3} Al_{2}O_{3}*Ts/m^{3}, Fe_{2}O_{3}*Ts/m^{3}, MgO*Ts/m^{3}, and SiO_{2}*Ts/m^{3} are then produced within the saprolites. The number of turning bands used for the non-conditional simulations is set at 1000 to avoid potential streaking artefacts. The conditioning by SK is performed within a (300×200×12 m^{3}) neighbourhood rotated at azimuth 135° using eight sectors and seven 3 m composites per sector (total of 56 composites per neighbourhood).

To allow the incorporation of information effect incurred during grade control from a 10 m centred (10 m C) pre-exploitation ('pre-ex') drilling grid, a kriging of the variables of interest from the pre-ex sampling grid extracted from each realization is done following the production of the direct block simulation values and associated 3 m point values. The ranges of the neighbourhood utilized for that kriging are adapted (40×30×12 m^{3}) to reflect the density of information available at pre-ex stage.

A post-processing (utilizing all the realizations from the simulation platform) is then implemented to estimate the tonnage and metal quantities recoverable for each variable within 20×20×3 m^{3} panels. The same post-processing will be applied ultimately to each realization retained by scenario reduction.

In summary, the estimation of multivariate recoverable resources at 20×20×3 m^{3} panel level is performed for the series of economic cut-offs in the following manner:

1. For each 5×5×x3 m^{3} SMU simulated, 3 m composite values are extracted at random

2. A selection corresponding to the 10 m centered pre-ex sampling grid (10 m C) is made on the simulated points

3. 5×5×3 m^{3} SMU values are re-estimated from the 10 m C sampling set

4. Then, for each cut-off grade:

]]> - For each realization of the simulation platform, select the SMUs such that the estimated values kriged at step 3 are above the cut-off grade- For each 20×20×3 m

^{3}, calculate tonnage and metal contents for all variables and all 5×5×3 m^{3}smu making it past the cut-off- Looping on all realizations, calculate the recoverable estimates as the average of the values obtained at step b for all the realizations

5. Repeat for all cut-offs.

**Scenario reduction**

The above procedure allowed production of a robust estimate of the recoverable resources for the six variables of interest, which can be fed into the reserve calculation process. We present here the application of a scenario reduction algorithm that selects a representative subset of a few simulations out of the original platform to help appraise the risk attached to the downstream phases (reserve optimization, mine sequencing) of the project.

The implementation presented here is for the first generation of the scenario reduction plug-in built in Isatis^{®} software (see Bleines *et al.,* 2012), where measurement of the distance between the initial set of scenarios and the reduced set is based on Ni quantities only. Research is under way to adapt the measure to the true multivariate nature of the problem.

The central idea to the algorithm used is to characterize the difference between any two simulations by integrating, over all the panels, the difference in recoverable metal quantities between the two simulations for all the cut-offs. Once a dissimilarity matrix has been established for the 50 simulations in the platform, a combinatorial based on k-mean clustering is used to select the k (here, k=5) simulations that best capture the space of uncertainty as characterized by the 50 realizations. A detailed presentation of the underlying concepts is given by Armstrong *et al.* (2010), who base the minimization of the distance between the selected subset of simulations and the remainder set on a random sampling procedure of the combinatorial to be treated. The current selection procedure is based on a genetic sampling algorithm, which is presented in a companion paper in this volume. The parameters of the genetic sampling are optimized for the prior selection of three simulations out of 50 for which it is possible to perform an exhaustive sampling of all possibilities.

The set of simulations selected (29, 49, 2, 31, 32) comes with associated probabilities (76%, 18%, 2%, 2%, 2%). As can be seen, realization 29 can be viewed as a representative case, while realization 49 represents more marginal situations, and realizations 2, 31, and 32 can be considered as representative of the extremes.

]]> Post-processing of each one of these five realizations is then possible, and will help qualify the answer given by the resource estimate.

**Impact on economic evaluation**

Several pit optimizations (using Whittle software) are then implemented using as an input to the process either:

>* *The traditional estimates based on the independent ordinary kriging of the accumulated variables panel by panel (leading to what is referred to in the following tables as the 'Traditional Pit'), or

>* *The recoverable resource estimation based on the full multivariate conditional simulation platform (leading to what is referred to in the following tables as the 'Pit based on the CS platform').

To help assess the risk and sensitivity to the uncertainty on the recoverable resources, the reserves are then calculated for both pits using the above resource models as well as the five simulations selected by scenario reduction.

Tables II and III give the ore tonnage for the seven input models (actual = traditional estimates, recov = recoverable resource based on CS platform, simu29, simu49, simu2, simu31, simu32 for the simulations selected by scenario reduction) in the two pits selected. The total tons are split according to the ore circuit they follow (TRADI is ore sent straight for shipment to the mill, UTM is marginal ore that is pre-processed at the plant).

The results highlight the following:

>* *The risk of underachieving actual ore tons mining the current pit geometry

*The higher value, longer LOM within the pit optimized using the conditional co-simulation (CS) platform*

>* *The risk factor based on the selected simulations in the CS pit is acceptable (3% variation on tons; from -1 to +6% on cash flow)

>* *The improved product mix within the CS pit (35% of TRADI (i.e. rich) ore instead of 25% for the actual pit), meaning a reduction in cost and time lost remobilizing material.

**Conclusions**

The objective of the procedure developed for Tiébaghi was to estimate the recoverable proportions of 5×5×3 m^{3} SMUs within each 20×20×3 m^{3} panel when a cut-off is applied on Ni, and to attach to that proportion the grade of all five relevant elements being studied. An original procedure, based on direct block simulations, which is capable of handling the impacts of information and support effects as well as the non-additivity of the input variables has been developed.

The advantage of resorting to a conditional simulation platform is that it allows a full and proper characterization of the variability of key elements at all scales within the orebody. The post-processing proposed allows for a robust single recoverable estimate to be produced that is amenable to treatment by pit optimization. Although this summarizes the information in a way that allows an adapted economic valuation of the project (which properly captures the multivariate nature of the valuation criteria) to be undertaken, valuable information regarding the characterization of variability that is carried by the platform of simulations is lost in the process.

The scenario reduction procedure that is implemented offers a way to circumvent that issue by allowing an efficient selection of a subset of realizations that properly samples the space of uncertainty and paves the way for an industrial characterization of the risk attached to the project due to the uncertainty on the resource.

**References**

Anand, R.R. and Paine, M. 2002. Regolith geology of the Yilgarn Craton, Western Australia: implications for exploration. *Australian Journal of Earth Sciences,* vol. 49, no. 1. pp. 3-162. [ Links ]

Armstrong, M., Ndiaye, A., Razanatsimba, R., and Galli, L. 2013. Scenario reduction applied to geostatistical simulations in mining. *Mathematical Geosciences,* vol. 45, no. 2. pp. 165-182. [ Links ]

Bleines, C., Deraisme, J., Geffroy, F., Perseval, S., Rambert, F., Renard, D., Touffait, Y., and Wagner, L. 2012. ISATIS software manual. Geovariances and Ecole des Mines de Paris. 531 pp. [ Links ]

Cluzel D., Aitchison J.C., and Picard C. 2001. Tectonic accretion and underplating of mafic terranes in the Late Eocene intraoceanic fore-arc of New Caledonia (SW Pacific): geodynamic implications. *Tectonophysics,* vol. 340. pp. 23-59. [ Links ]

Chilès, JP. and Delfiner, P. 1999. Geostatistics - Modeling Spatial Uncertainty. Wiley, New York. pp. 478-508. [ Links ]

]]>Deraisme, J., Rivoirard, J., and Carrasco Castelli, P. 2008. Multivariate uniform conditioning and block simulations with discrete gaussian model: application to Chuquicamata deposit. *GEOSTATS2008, Proceedings of the 8th International Geostatistics Congress,* Santiago, Chile, 1-5 December 2008. pp. 69-78. [ Links ]

Eggleton, R.A. 2001. The Regolith Glossary. Cooperative Research Centre of Landscape Evolution and Mineral Exploration, Perth, Australia. 144 pp. [ Links ]

Emery, X. and Ortiz, J. 2005. Internal consistency and inference of change of support isofactorial models. *Geostatistics Bannf2004,* vol.2. Leuangthong O. and Deutsch, C.V. (eds.). Springer, Dordrecht. pp. 1057-1066. [ Links ]

Journel, A.G. and Huijbregts, Ch. J. 1978. Mining Geostatistics. Academic Press, London. p. 199. [ Links ]

Nahon, D., Colin, F., and Boulange, B. 1992. Metallogeny of weathering: an introduction in weathering and soils and paleosols. *Weathering, Soils and Paleosols.* Martini I. P. and Chesworth W. (eds.). *Developments in Earth Surface Processes,* vol. 2. Elsevier, Amsterdam. pp. 445-471. [ Links ]

Pelletier, B.G. 2003. Les minerais de nickel de Nouvelle-Calédonie. *Revue Géologues,* vol. 138. pp. 30-38. [ Links ]

Rivoirard, J. 1994. Introduction to Disjunctive Kriging and Non-linear Geostatistics. Clarendon Press, Oxford. [ Links ]

Tardy, Y. 1997. Petrology of Laterites and Tropical Soils. Balkema, Amsterdam. [ Links ]

]]>