A critical task of any mining projects is to construct a three-dimensional block model mainly representing the tonnage and grade distribution of the mineralised deposit. It is constituted by arrays of blocks holding necessary geological attributes, i.e. mineral and metal grades, lithology codes, density and tonnage. The quality of a resource block model is highly dependent on how attributes of blocks are estimated from exploration structures. Geostatistical estimation methods, such as ordinary kriging and simple kriging, have widely been used since the 60 s of the last century given their superiority in considering spatial distribution of samples in the calculations (David 1977; Isaaks and Srivastava 1989; Ravenscroft and Armstrong 1990). Nevertheless, there is still no available algorithm which can provide 100% accurate estimates. Thus, estimation errors are inherent in all resource block models and this phenomenon is generally referred to as geological uncertainty. Considering this drawback, the true grade and tonnage of a deposit can be vastly different from the estimated figures and the mine plan constructed based on this resource model tends to fail to achieve production targets and revenue expectation (Baker and Giacomo 1998; Vallee 2000; Dimitrakopoulos et al. 2002; Groeneveld and Topal 2011; Groeneveld et al. 2012).

To fulfil the inability of estimation techniques in characterising variability, conditional simulation is developed to generate a series of plausible possibilities of an orebody, which are commonly termed as realisations. Each realisation has an equal chance to be real and this provides an appropriate platform to analyse the risk-associated with the geological conditions (Goovaerts 1997; Vann et al. 2002; Dimitrakopoulos 2011; Topal and Ramazan 2012). To deal with multivariate deposits, where at least two attributes of interest are correlated, multivariate conditional simulation techniques have been developed (Matheron 1979; Carr and Myers 1985; Gómez-Hernández and Journel 1993; Goovaerts 1993; Verly 1993; Desbarats and Dimitrakopoulos 2000; Leuangthong 2003), as discussed in the next section.

Once the resource block model is available, strategic mine planning is the next critical step to the success of a mining project as it decides on the economic output. Although one of those realisations can be true, detecting the correct one is not possible for the time being. Therefore, deterministic estimation techniques, such as ordinary kriging, are still preferable in the industry and will not be replaced soon in the role of providing resource models for mine planning. Conditional simulation, however, can be applied in two promising areas: (1) to characterise the geological uncertainty of the kriged resource model and (2) to analyse risk of the kriging-based mine plan against the geological conditions. This information is crucial for decision makers and shareholders as a new dimension for analysing the potential of mining projects. In this paper, a new systematic framework to quantify the geological risk in mining projects is proposed. Then, it is followed by the implementation of the prosed framework onto an iron deposit to demonstrate its practical aspects.

The remainder of the paper is outlined as follows, in Sect. 2 methodology to characterise geological uncertainty of a multivariate deposit is developed using a new multivariate conditional simulation framework and subsequently implemented to an iron ore deposit in Western Australia. In Sect. 3, strategic mine planning is implemented using an in-house mine planning tool and the risk profiles associated with that mine plan are constructed using simulations. Conclusions and recommendations follow at the end.

It is common in mining industry that the deposits contain multiple correlated attributes of interest. A typical example is iron ore mining, where up to six variables need to be considered in resource estimation and mine planning, namely iron (Fe), silica (SiO₂), alumina (Al₂O₃), phosphor (P), loss on ignition (LOI) and Ochreous Goethite (GOL), and they are cross-correlated. As demonstrated in Table 1, iron and silica contents have a strong negative correlation of −0.72 which means that if an ore block has a high iron grade, the silica grade must be low and vice versa. Because the nature of simulation is random, performing simulation on these two attributes without considering their correlation with other attributes will lead to an unrealistic result (Mai et al. 2016).

2.1 A new framework for multivariate conditional simulation of iron ore deposits

Several methods for multivariate conditional simulation have been developed in the past decades, including co-simulation (Matheron 1979; Carr and Myers 1985; Verly 1993), Stepwise Conditional Transformation (SCT) (Leuangthong and Deutsch 2003), Principal Component Analysis (PCA) (Switzer and Green 1984; Goovaerts 1993) and Minimum/maximum Autocorrelation Factors (MAF) (Desbarats and Dimitrakopoulos 2000). Each of these above mentioned approaches has their own strengths and weaknesses. Co-simulation can be applied directly on the correlated variables but it becomes computationally inefficient to work with more than three variables. SCT and PCA are capable of de-correlating variables only at zero or small lag distance, in addition, SCT requires an intensive amount of samples to ensure the accuracy of its transformation procedure. Finally, MAF is the most suitable approach for doing simulation of iron ore deposits, where MAF can be implemented in any lag distances and has no specific requirements for the number of samples. Moreover, MAF is available in ISATIS™, a commercial geostatistical software package (Bleines and de Paris 2000).

As MAF works in normal score space, a Gaussian-based simulation technique is preferable to perform simulation on MAF factors. In this study, Sequential Gaussian simulation (SGS) was selected because its application in mining industry has been widely recognised (Johnson 1987; Dowd 1994; Journel 1994).

A six-step framework is proposed to perform the joint simulation using MAF and SGS:

Step 1: Primary Gaussian transformation: Running Gaussian anamorphosis to transform the sample data \(Z\left( x \right)\) into normal (Gaussian) scores \(Y\left( x \right)\);

Step 2: MAF transformation: Transform \(Y\left( x \right)\) into uncorrelated MAF factors \(M\left( x \right)\);

Step 3: Secondary Gaussian transformation: Run secondary Gaussian anamorphosis if the Gaussianity of MAF factors are not adequate. After this transformation, the output data is called normal score MAF factors \(N\left( x \right)\);

Step 4: Continuity analysis and variography of normal score MAF factors \(N\left( x \right)\) individually;

Step 5: Perform SGS on each normal score MAF factors \(N\left( x \right)\) individually;

Step 6: Back transformation of step 3, 2, and 1 sequentially.

where \(Z\left( x \right)\) is the original data; \(Y\left( x \right)\) is the normal scores; \(M\left( x \right)\) is the MAF factors; \(N\left( x \right)\) is the normal score MAF factors; \(N^{*} \left( x \right)\) is the simulated result of \(N\left( x \right)\); \(\phi\) is the primary Gaussian anamorphosis; \(\phi^{\prime}\) is the secondary Gaussian anamorphosis; A ^T is the MAF transformation matrix; \(y\) is the sequential Gaussian simulation.

The simulation procedure is schematically demonstrated in Fig. 1.

Schematic presentation of the proposed simulation framework. NS normal score, Transf. transformation, anam. anamorphosis

2.2 Implementation of the proposed framework: Kriging versus Simulation

The proposed framework is implemented on an iron ore block model consisting of 985,088 blocks with a block size of 25 m × 12.5 m × 10 m. 20 orebody realisations were generated using the proposed simulation framework. Ordinary kriging was also performed on the same block model. The grade/tonnage curves of the realisations and kriged model are presented in Fig. 2.

Grade/tonnage curves of regular kriging versus 20 SGS realisations

By comparing the kriged resource model with 20 orebody realisations, the geological uncertainty of the deterministic resource model can be identified as:

• The estimated Fe grade would be considerably lower than the reality by approximately 1%;

• The kriged model overestimates total ore tonnage when Fe grade is lower than the mean of the composite Fe (59.23%) and underestimates when Fe grade is higher. This is an evidence of smoothing effect of ordinary kriging;

• The in situ tonnage variability of iron ore is approximately 3 million tonnes;

• At cut-off grade of 57.5% Fe, there is a risk that the actual iron ore tonnage would be approximately 10–13 million tonnes less than expectation;

• The number of 20 realisations generated for the given iron ore deposit is adequate to capture its uncertainty characteristics according to the reasonable similarity between realisations.

Plan views of the kriged and example of simulated orebody models are presented in Figs. 3 and 4.

Kriged resource model for Fe element

SGS realisation 1 for Fe element

3.1 Utilisation of an in-house strategic mine planning tool

In this study, an in-house mine planning tool was deployed to find the optimal production schedule over the life of mine upon the kriged resource model generated in Sect. 2.2. The mine planning process consists of two phases: In phase 1, blocks are aggregated using TopCone Algorithm (TCA) to create TopCones (TCs). The main purpose of this phase is to significantly reduce the amount of data being processed in the production scheduling progress. Then in phase 2, TCs are fed into an integer programming (IP)-based model to optimise the long-term production schedule, where the project’s NPV is maximised subject to various operational constraints.

The basic ideas of our in-house strategic mine planning tool can be described as follows:

3.1.1 Phase 1. Block aggregation

TCA formulates and solves a series of linear programming models to combine blocks into TCs. The main features of TCA are:

• TCs have positive economic value and can be mined in a certain order without violating the slope safety

• The number of TCs generated can be controlled by setting the minimum number of blocks per TC

• The combination of all TCs forms an ultimate pit

The performance of TCA is demonstrated via a hypothetical two-dimensional deposit of 15 blocks, of which blocks 7, 8, 9 and 13 are ore, as presented in Fig. 5. Implementing TCA with a minimum number of block per TC is one, three TCs and an ultimate pit were generated. Joint support ability of TCA can be observed at TC 1 as either ore block 7 or 8 is strong enough to support their overlying waste blocks. In addition, the extractions of TC 1 then 2 and 3 sequentially always secure the slope safety.

Three TopCones were generated with minimum size of one block

To demonstrate the ability of TCA in controlling the number of TCs generated from the aggregation process, the algorithm was also implemented with a minimum size of three blocks per TC, the number of TCs reduced to two and still inside the same ultimate pit, as shown in Fig. 6.

Two TopCones were generated with minimum size of three blocks

Similarly, keep increasing minimum size of TC and the whole ultimate pit only contained a single TC. This ability of TCA to control number of TCs is critical for the application of mathematical programming on open pit mine planning, as the solution time and computational intensiveness of solving mathematical models are exponentially related to the amount of data imported. Generally, a standard computing system is not able to solve an IP-based production scheduling model for 5 years life of mine over 10,000 blocks within a practical timeframe, meanwhile industry-standard resource models normally contain much more than that. In other words, data scale is the greatest obstacle of optimising mine plans using operations research techniques and TCA was developed to directly tackle this challenge. In addition, the precise mining sequence between TCs allows to significantly reduce the number of sequencing constraints in the downstream mathematical model and contribute to cut the solution time further. In addition, the ability to find an ultimate pit limit of TCA eliminates the implementation of external pit optimisation algorithms in the proposed mine planning procedure.

Running TCA on the given iron ore deposit with the minimum size of 200 blocks per TC, 500 TCs were generated in 10 min using an office-standard computer having an Intel(R) Core(TM) i7 with 3.4 GHz CPU processor and 8 Gb of RAM. The visualisation of 500 TCs and the ultimate pit limit are presented in Fig. 7.

Visualisation of 500 TCs and the ultimate pit limit

3.1.2 Phase 2. Production scheduling using integer programming

Once TCs are generated, each of them has a specific economic value, mineral/metal grade, ore tonnage, waste tonnage and mining sequence with other TCs. The formulation of the IP model followed traditional format with an objective function of maximising discounted cash flow, operational constraints and slope safety. Similar IP models can be found at the works of Caccetta and Hill (2003); Gleixner (2009); Ramazan and Dimitrakopoulos (2004). As the number of TCs is relatively small, the IP model was efficiently solved by CPLEX™ (CPLEX 2009) in less than 20 min. A set of hypothetical scheduling parameters, as presented in Table 2, for a mine plan of 6 years was used to prove the concepts. The model, however, can easily be adjusted for real input values.

To prove the practical mining sequence of the result, some typical cross-sections of the mine plan are presented in Fig. 8.

Typical cross-sections of mining sequence generated by TCA-based IP model

3.2 Quantifying risk

To analyse the risk associated with the mine plan which was determined based on the kriged resource model, its mining sequence, as presented in Fig. 8, was sequentially applied on the 20 orebody realisations generated in Sect. 2. The purpose of this task is to find all possible production scheduling outcomes of the mining project under the impact of geological uncertainty. The results are presented in Fig. 9. The red lines are the scheduling results of IP model implemented on kriged resource model while blue dots are the associated risk profiles.

Risk profiles of the mine plan using TCA-based IP and kriged resource model

From the risk profiles of material production, it is likely that there are less ore tonnage and more waste tonnage in all periods over the life of mine. The blending grades of 6 attributes are also highly deviated from predicted values by the kriging-based mine plan. For example, Fe grade is considerably higher than expectation whereas those of SiO₂, Al₂O₃, P and LOI are lower. A 60% chance of violating upper bound of the GOL grade can be observed at the 2nd period, which does not exist in the deterministic mine plan.

In addition, calculating the risk associated with project’s NPV, as shown in Fig. 10, points out the likelihood that the project will not achieve the expected NPV. Indeed, the NPV deficit due to the impact of geological uncertainty is between 6.70% and 7.34%, equal to 135 M$ and 151 M$. This great loss of NPV again emphasises the detrimental impact of geological uncertainty on project’s valuation.

NPV risk profiles of the mine plan using TCA-based IP and kriged resource model

In this study, the authors successfully applied joint simulation using MAF and SGS on a multivariate deposit consisting of six attributes of interest. 20 realisations and a kriged resource model were generated to facilitate the analysis of grade/tonnage variability and smoothing effect of ordinary kriging.

The authors deployed an in-house mine planning tool where block aggregation and integer programming were used to find the optimal production schedule for the kriged resource model. By testing the mining sequence against the 20 orebody realisations, the risk profiles of the deterministic mine plan were constructed. The results proved that under the impact of geological uncertainty, the project’s NPV can be considerably less than expectation at a magnitude of 6.70% to 7.34%, as well as some strong deviations of blending grade and tonnage.

Given that kriging estimation and deterministic mine planning techniques will well remain their popularity in the near future, it is recommended that any mine plan built based on deterministic resource model should be validated against geological uncertainty. Ignoring this factor could turn a sound mine plan in feasibility study into a project failure in the future.

Similar to iron ore deposits, the nature of coal deposits is typically multivariate, where attributes of interest, such as calorific value, sulphur, nitrogen, and ash content, are commonly cross-correlated. To this extent, the application of multivariate conditional simulation and mining planning methodology as outlined in this paper can be easily adapted to coal mining.