Effects of High-Order Simulations on the Simultaneous Stochastic Optimization of Mining Complexes

A mining complex is composed of mines, mineral processing streams, stockpiles, and waste facilities, which culminate with generated products that are delivered to customers. The supply uncertainty and variability of materials extracted from the mines, which flow through a mining complex to generate products, can be quantified through geostatistical simulations and can be used as inputs to the simultaneous optimization of mining complexes. A critical aspect to consider is that mineral deposits are characterized by spatially complex, non-Gaussian geological properties and multiple-point connectivity of high-grades, features that are not captured by conventional second-order simulation methods. This paper investigates the benefits of simultaneously optimizing a mining complex where the simulations of the mineral deposit are generated by a high-order, direct-block simulation approach. The optimized life-of-mine (LOM) production schedule is compared to a case in which the same setting is optimized by having the related simulations generated using a second-order simulation method. The comparison shows that the incorporation of simulations that reproduce the spatial connectivity of high-grades results in a more informed LOM production schedule. The sequence of extraction is driven by the spatial connectivity of high-grades, resulting in a mill throughput with better material quality and reduced waste extraction. Furthermore, the discounted cash-flow increases by more than 5% as compared to the case in which the second-order simulations are used.


Introduction
A mining complex can be perceived as a transfer function that entails complex interactions, starting with the extraction of materials from the mines through to their transformation into saleable products, passing through different processes in the mineral value chain [1,2].It can include, for example, multiple mines, various elements and material types, stockpiles, tailings facilities, waste dumps, processing plants and transportation systems [3].The simultaneous stochastic optimization of a mining complex maximizes the global value of the mining operation by integrating all of the components in a single framework, while also including the uncertainty and variability of materials sourced from the mines as a set of stochastically generated realizations [1,2,[4][5][6][7].
The modeling of large mining operations as an integrated mathematical model has advanced over the last three decades.Newmont Mining Corporation pioneered the area, proposing a linear programming approach to the massive gold mining complex in Nevada [8].BHP also contributed significant improvements with the development of the Blasor optimizer, BHP's in-house software that was applied to the Yandi mine complex in Australia [9], optimizing the life-of-mine extraction sequence from multiple pits.Additional applications have extended the software to optimize and assess more complex requirements [10,11].Whittle [12] and Whittle [13] describe an industry standard tool to sequence the extraction of material from multiple deposits in operations with complex blending and processing requirements.Topal and Ramazan [14] present an approach based on linear network flow to optimize a mining complex in Western Australia with more than 100 pits and 13 processing facilities.Khan and Asad [15] propose a mixed integer linear programming model extending the standard cut-off grade strategy used in the industry for the case of multiple processing streams.Although the methods mentioned above represent significant efforts to integrate the multiple components of a mining complex in order to capitalize on synergies, simplifications exist.Conventional approaches optimize the mining complex by solving each step or component sequentially; thus, they do not benefit from the interaction and cooperation of different components, which leads to suboptimal solutions for a value chain as a whole [1,3].Another limitation in the studies to date is the aggregation of mining blocks into larger decision variable units, such as panels, which is done to make the optimization simpler but can misrepresent the mining selectivity.Moreover, the uncertainty in material supply coming from the mines has been long recognized as the primary cause of technical risk in mining operations [16] and, if left unmanaged, leads to unexpected deviations in production targets [17][18][19][20].
The spatial uncertainty and variability of attributes in geosciences can be quantified via geostatistical simulations [21][22][23][24][25][26], which are founded on the concept of random fields.The sequential simulation approach is an alternative to assess these attributes at each unsampled location of a three-dimensional orebody model through Monte Carlo sampling from a conditional distribution function [21][22][23].Traditional simulation methods are based on the second-order statistics, namely, mean and covariance (variogram), where sequential Gaussian simulation (SGS) [23,27,28], sequential indicator simulation (SIS) [23] and sequential direct block simulation [29,30] are some examples.However, geological attributes of spatially distributed phenomena are represented by complex non-Gaussian and non-linear spatial connectivity of low-and high-grades.Using only up to second-order statistical statistics in spatial simulations is not sufficient to fully describe complex geological attributes [31].In addition, Gaussian simulation methods maximize spatial disorder (maximum entropy) beyond the imposed variogram in the realizations generated [32], thus preventing a more realistic quantification of the connectivity of high-grades.
Multiple-point statistics (MPS) simulation methods have been introduced in an attempt to address the above limitations [33][34][35][36][37][38][39][40][41][42].MPS methods replace the random field model with a geological analogue or training image (TI) and facilitate the reproduction of complex curvilinear and other geologic features in the realizations generated, while avoiding distributional assumptions.A consequence of extracting patterns from a TI is that generated realizations are TI-driven and their spatial statistics may be different from those of the available data, particularly when relatively dense drillhole data sets are available [43,44].As a natural extension of second-order methods, the high-order simulation framework [45][46][47][48][49][50] can reproduce very complex non-linear geometries and spatial statistics from data by explicitly calculating high-order spatial cumulants.The generated realizations present a more realistic and structured connectivity of high-grades (lower entropy) compared to traditional methods, as shown in the example in Figure 1.
The appropriate characterization of spatial connectivity and its impact in flow modelling are well-studied subjects in the context of reservoirs and aquifers [32,[51][52][53].In the mining context, some studies have shown that the use of different simulation frameworks impacts the output of transfer functions [54,55].Geostatistical simulations have been effectively incorporated into state-of-the-art simultaneous stochastic optimization of mining complex frameworks [1,2,[4][5][6]56,57].The next step is to investigate the effects of high-order simulation models in this optimization framework.The appropriate characterization of spatial connectivity and its impact in flow modelling are well-studied subjects in the context of reservoirs and aquifers [32,[51][52][53].In the mining context, some studies have shown that the use of different simulation frameworks impacts the output of transfer functions [54,55].Geostatistical simulations have been effectively incorporated into state-of-the-art simultaneous stochastic optimization of mining complex frameworks [1,2,[4][5][6]56,57].The next step is to investigate the effects of high-order simulation models in this optimization framework.
This paper presents an application of the high-order direct block-support simulation method [48] used as an input to the simultaneous stochastic optimization of a simplified mine complex, using the framework proposed by Goodfellow and Dimitrakopoulos [6].As a means of comparison, the same mine complex is optimized using simulations obtained with a second-order method.The differences in life-of-mine production schedules and related production forecasts, optimized for each case, are analyzed and discussed.The following sections provide a brief description of the high-order direct block simulation framework and the optimization model.This is followed by a comprehensive analysis of the different outcomes from the two approaches, and finally, a conclusion section is presented.

Modeling a Mineral Deposit using Geostatistical Simulations
Considering a random function (RF)   (  ),   ∈   , where   represents the location of the point support grid to be simulated in the domain  ∈    Instead of simulating the entire deposit at the point support discretization, the work by de Carvalho et al. [48] presents an alternative method that allows the generation of high-order simulations directly at the block scale.Assume    is a realization of the RF   (  ) at the block support scale, defined in the same domain  ∈   , where   represents the centroid of the block in This paper presents an application of the high-order direct block-support simulation method [48] used as an input to the simultaneous stochastic optimization of a simplified mine complex, using the framework proposed by Goodfellow and Dimitrakopoulos [6].As a means of comparison, the same mine complex is optimized using simulations obtained with a second-order method.The differences in life-of-mine production schedules and related production forecasts, optimized for each case, are analyzed and discussed.The following sections provide a brief description of the high-order direct block simulation framework and the optimization model.This is followed by a comprehensive analysis of the different outcomes from the two approaches, and finally, a conclusion section is presented.

Modeling a Mineral Deposit Using Geostatistical Simulations
Considering a random function (RF) Z P (x i ), x i ∈ R d , where x i represents the location of the point support grid to be simulated in the domain D ∈ R d and z P i a realization of Z P (x i ).The set of initial data is given by d n = {z p (x k ), k = 1, . . ., n} and represents the values obtained from the exploration data.Now consider the set Λ i with the data and previously simulated nodes, i.e., Λ 0 = {d n } and Λ i = {Λ i−1 ∪ Z(x i )}.Thus, under the sequential simulation framework [23,27,52,58,59] the global conditional distribution can be decomposed in the multiplication of the univariate distributions 2.1.1.High-Order Direct Block-Support Simulation Instead of simulating the entire deposit at the point support discretization, the work by de Carvalho et al. [48] presents an alternative method that allows the generation of high-order simulations directly at the block scale.Assume z v i is a realization of the RF Z v (v i ) at the block support scale, defined in the same domain D ∈ R d , where v i represents the centroid of the block in consideration.Additionally, the conditional distribution, Minerals 2019, 9, 210 4 of 15 where λ 0 and λ v i−1 are the set of locations of Λ 0 and Λ v i−1 , respectively.For simplicity, Thus, the above is approximated as where ϕ m (•) belongs to the set of Legendre-like orthogonal splines [49,60], and the coefficient L i...jk...l is approximated experimentally by The method utilizes a TI represented in both support sizes, point and block.Thus, having a spatial template obtained with the block to be simulated and neighboring values, both at the block and point support, this TI is scanned searching for replicates of the template in consideration.The algorithm for block support high-order simulation can be summarized as: 1.
Upscale the TI inputted at point support to block support.

2.
According to the sequential simulation framework, define a random path to visit all the unsampled block locations.

3.
At each block location: a.
Find the closest point and block support values for conditioning.b.
Obtain a spatial template configuration formed by the block to be simulated and related conditioning values.c.
Scan the TI searching for replicates of the above template.d.
e. Derive the conditional cross-support joint probability density function by first calculating the joint distribution in Equation ( 4), then normalizing the distribution in Equation (2).f.Draw a uniform value from [0, 1] to sample z v i from the conditional cumulative distribution derived from the above.g.
Add z v i to the simulation grid at block support at location v i to be used as conditioning value to the simulation of a subsequent block.

4.
Repeat steps 2 and 3 to generate additional realizations.
The method assumes that the conditional distribution f x i ; z p i |Λ i−1 is Gaussian, which facilitates the simulation process, since its approximation only requires the definition of two parameters, namely, mean and variance.For this reason, the original data is transformed into the Gaussian space typically by a graphical transformation; and at every node, the method solves a kriging system to obtain the posterior mean and variance.The simulation is performed at the point support scale with a subsequent re-blocking procedure to generate block-support models.

Mathematical Formulation of the Simultaneous Optimization of Mining Complexes
The current study uses the simultaneous optimization of mining complexes model proposed in [1,6], which is summarized as follows.In this setting, the mine is discretized into mining blocks indexed in b ∈ B. The framework considers that each block has simulated attributes, such as grades and material types, which will denote a stochastic scenario s ∈ S. The extraction of each block b in period t ∈ T incurs a mining cost MC b,t , but its extraction can only happen if the set of predecessors O(b) has already been extracted.Once extracted, the material can flow from the mine to a stockpile or a destination, such as waste pile, leach pad or mill (processors i ∈ P).The cost associated with material transportation is given by TC i,a,t .The amount of property a in location i, period t and scenario s is quantified in v a,i,t,s .Material that flows from the mine to a location, for example, ore tonnage, is indexed in p ∈ P, whereas transformations, potentially non-linear ones such as recovery, are indexed in h ∈ H. Thus, p h,t,s denotes the unitary value of selling the material property h, in period t and scenario s.PC i,p,t is the processing cost of treating material property p, in location i at the period t.The set of production targets is represented by P c , where deviations are quantified in d ± i,a,t,s and penalized by the cost c ± i,a,t .The model also incorporates smoothing and sink constraints [62,63], where deviations of these targets are quantified in and penalized by d smooth b,t , d sink b,t,v and c smooth b,t , c sink b,t , respectively.Additional sets of constraints such as slope, reserve, capacities, destination policy and processing stream constraints are detailed in [6].
There are three types of decision variables in the described model.Extraction sequence (x b,t ∈ {0, 1}) returns 1 if the block b is extracted at the period t, 0 otherwise.Destination decisions z g,j,t ∈ {0, 1} define where to send a group of material g to the destination j at the period t.These groups are defined in a similar manner to what is described in [11], where pre-defined grade bins are inputted, and the optimizer decides the single element optimal cut-off grade boundaries.Processing stream decisions y i,j,t,s ∈ [0, 1] define the proportion of material that flows from location i to destination j, in period t and scenario s.
The objective function in Equation ( 6) maximizes the value of selling the products in the mineral value chain, while minimizing deviation from production targets.
Part I of Equation ( 6) considers the discounted cash flow obtained by selling the products in the value chain.Part II minimizes the processing cost at each processing facility and the transportation costs involved.The third Part is related to the cost of deviating from processing and mining capacities, respectively.Part IV minimizes the deviation from schedule smoothness and mining Minerals 2019, 9, 210 6 of 15 costs, while Part V aims to reduce deviations from sink rate constraints, respectively.Note that geological risk discounting [62] is applied to all penalty costs associated with production target, i.e., where r is the geological discount rate.
The simultaneous optimization of a mining complex framework presented above is very general and allows for the integration of different sources of uncertainty; in addition, it addresses the non-linearity related to materials blending and stockpiling, as well as transformations that are due to refining materials to output products.The above optimization approach leads to instances with multi-million binary decision variables that cannot be solved with commercial solvers, such as CPLEX.Metaheuristic algorithms offer a practical alternative and have proven to be an efficient solving approach for the stochastic optimization of mines and mining complexes [6,56,[64][65][66].The solution approach used in the current paper is from Goodfellow and Dimitrakopoulos [1,6].

Results and Discussion
The first part of the case study relates to the simulation of the grades of a gold mine.The deposit covers an area of approximately 4.5 km 2 and extends to a depth of 400 m.The three-dimensional orebody model is composed of 510,800, 10 × 10 × 10 m 3 mining blocks.The available data come from 2344 drillholes that are spaced at about 35 m apart and include 40,762, 10 m long gold composites.A set of 15 high-order simulations are generated directly at the scale of mining blocks using a training image generated from blasthole data obtained at a 5 m spacing.Each high-order simulation required approximately 6 h on an Intel ® Core TM i7-7700 CPU with 3.60 GHz and 16GB of RAM, running on Windows 7.For comparison, a set of simulations based on the second-order statistics is also generated using traditional sequential Gaussian simulation (SGS) [23,27,61].Realizations generated at the point-support scale are rescaled to block support using a discretization of 25 nodes per block.

Results, Comparisons and Effects of High-Order and Second-Order Simulations
First, to provide a common ground for comparison, Figure 2 shows the grade-tonnage curve for both simulation frameworks with equivalent overall reporting and quantification of related uncertainty.The graph shows very similar proportions regarding tonnages and grades over the deposit.Although the metal quantity is very comparable in both cases, how each method connects these elements in space can be very different, especially at the high-grade values, as noted earlier in Figure 1. Figure 3 displays cross sections of both second-order and high-order realizations of the deposit, where the high-grade zones are highlighted with red circles.It is possible to visualize the effect of the maximum entropy property over the second-order simulations, which is enhanced by the fact that the simulation process was performed in the Gaussian space.The grades displayed by the simulations generated with SGS are visually more dispersed than those generated using high-order simulations.This connectivity can be quantified [52] and it is presented in Figure 4. Please note that, in all figures, P10 and P90 represent the 10th and 90th percentiles of the reported values, respectively.For the connectivity plot, the cut-off applied is 5 g/ton, corresponding to the 99th percentile of the grade distribution, as the focus of the comparison is on the high-grades.In the NE direction, the second-order realizations are consistently less connected than the high-order realizations for all lags.The difference becomes more pronounced in the NE direction and at the 45º dip, with a considerable gap evident between both simulation methods.As the high-grade mineralization drives the mine production schedule, this plays a vital role in the optimization of the mining complex.

LOM Production Schedule Optimization and Forecasting
The mining complex considered in the test case consists of a single gold mine, where ore and waste material can flow to the following destinations: Leach pad, ore stockpile, waste dump, and a mill circuit (Figure 5).Both the mill and the leach pad process ore to generate sellable products, while encountering different costs and recoveries.Since the high-grade materials are typically processed at the mill, having a substantially higher impact on cash flows, this processing stream receives more attention in this paper.Doing so simultaneously allows for investigating the impact of the high-grade spatial connectivity generated by the two simulation methods examined.The critical parameters for the optimization, displayed in Table 1, are kept the same in both cases.The uncertainty in the material properties sourced from the mine is quantified by the set of simulations generated and is used as input for the optimization framework.The mining complex is first optimized using simulations generated by the high-order direct block support method, which is referred to throughout the remainder of the paper as Case 1.The same mining complex setting is then optimized, using the set of simulations generated by the second-order Gaussian simulation method.The result of this optimization is referred to as Case 2.

LOM Production Schedule Optimization and Forecasting
The mining complex considered in the test case consists of a single gold mine, where ore and waste material can flow to the following destinations: Leach pad, ore stockpile, waste dump, and a mill circuit (Figure 5).Both the mill and the leach pad process ore to generate sellable products, while encountering different costs and recoveries.Since the high-grade materials are typically processed at the mill, having a substantially higher impact on cash flows, this processing stream receives more attention in this paper.Doing so simultaneously allows for investigating the impact of the high-grade spatial connectivity generated by the two simulation methods examined.The critical parameters for the optimization, displayed in Table 1, are kept the same in both cases.The uncertainty in the material properties sourced from the mine is quantified by the set of simulations generated and is used as input for the optimization framework.The mining complex is first optimized using simulations generated by the high-order direct block support method, which is referred to throughout the remainder of the paper as Case 1.The same mining complex setting is then optimized, using the set of simulations generated by the second-order Gaussian simulation method.The result of this optimization is referred to as Case 2.

LOM Production Schedule Optimization and Forecasting
The mining complex considered in the test case consists of a single gold mine, where ore and waste material can flow to the following destinations: Leach pad, ore stockpile, waste dump, and a mill circuit (Figure 5).Both the mill and the leach pad process ore to generate sellable products, while encountering different costs and recoveries.Since the high-grade materials are typically processed at the mill, having a substantially higher impact on cash flows, this processing stream receives more attention in this paper.Doing so simultaneously allows for investigating the impact of the high-grade spatial connectivity generated by the two simulation methods examined.The critical parameters for the optimization, displayed in Table 1, are kept the same in both cases.The uncertainty in the material properties sourced from the mine is quantified by the set of simulations generated and is used as input for the optimization framework.The mining complex is first optimized using simulations generated by the high-order direct block support method, which is referred to throughout the remainder of the paper as Case 1.The same mining complex setting is then optimized, using the set of simulations generated by the second-order Gaussian simulation method.The result of this optimization is referred to as Case 2.   The life of mine production schedules for Case 1 and 2 are optimized and results are discussed below.Cross sections of the LOM production schedules optimized for each case are displayed in Figures 6 and 7 along West-East and North-South directions, respectively.The areas with the same color represent the same period of extraction, and they highlight that the sequences of extraction obtained differ considerably, which is not surprising given the differences in the two simulation methods used.Cross sections of Case 1 show that the sequence of extraction follows a clear direction, highlighted by the red arrow.Note that this trend in the extraction sequence is amplified in the direction where the difference in connectivity is more evident, recall from Figure 4.The higher continuity of high-grades drives the schedule towards areas with more connected ore materials so that they can be processed together.The life of mine production schedules for Case 1 and 2 are optimized and results are discussed below.Cross sections of the LOM production schedules optimized for each case are displayed in Figures 6 and 7 along West-East and North-South directions, respectively.The areas with the same color represent the same period of extraction, and they highlight that the sequences of extraction obtained differ considerably, which is not surprising given the differences in the two simulation methods used.Cross sections of Case 1 show that the sequence of extraction follows a clear direction, highlighted by the red arrow.Note that this trend in the extraction sequence is amplified in the direction where the difference in connectivity is more evident, recall from Figure 4.The higher continuity of high-grades drives the schedule towards areas with more connected materials so that they can be processed together.
Figure 8 shows horizontal sections of the LOM production schedules generated and the differences in the sequences of extraction are again evident.Additionally, these sections show variations in the extension of the ultimate pit limits (UPL).The red circles highlight how much larger the UPL is in Case 2. As second-order simulations methods represent high-grade material as being more scattered, it is logical that the pits have to be larger to encompass all of the ore to be processed, resulting in a higher waste extraction, as shown in Figure 9a.These differences are of particular interest, since after the optimization is complete, infrastructure, such as ramps, access points, equipment placement and facility locations, reduces the flexibility to change the schedule.This optimization opportunity can only be achieved if the degree of connectivity of high-grades is correctly modeled.Figure 8 shows horizontal sections of the LOM production schedules generated and the differences in the sequences of extraction are again evident.Additionally, these sections show variations in the extension of the ultimate pit limits (UPL).The red circles highlight how much larger the UPL is in Case 2. As second-order simulations methods represent high-grade material as being more scattered, it is logical that the pits have to be larger to encompass all of the ore to be processed, resulting in a higher waste extraction, as shown in Figure 9a.These differences are of particular interest, since after the optimization is complete, infrastructure, such as ramps, access points, equipment placement and facility locations, reduces the flexibility to change the schedule.This optimization opportunity can only be achieved if the degree of connectivity of high-grades is correctly modeled.Once the production schedules are generated, the risk assessment of achieving production forecasts is calculated by pushing another set of simulations through each LOM schedule and by evaluating related forecasts.The results are presented in terms of P10 and P90, representing the 10th and 90th percentiles of related performance indicators, respectively.Regarding production targets and forecasts, Figure 9a shows the total tonnage mined over the LOM and the cumulative strip ratio for both cases.The production schedule obtained in Case 2 represents, in total, 5% more material mined than in Case 1.This difference reaches 8% at the end of the 10th year to ensure a similar throughput at the mill, Figure 9b.Mining more, in this case, translates to higher waste production, which is quantified by the higher strip ratios presented by Case 2. This can be explained by the spatial disorder (maximum entropy) that Gaussian-based approaches generate with respect to high-grades.Having ore blocks less connected in space forces the optimizer to mine more dispersed high-grade values (mining blocks), so as to provide a consistent feed rate to the mill; this also leads to different total tonnages of materials mined and the differences in UPL shapes, shown in Figure 8.On the other hand, if the optimizer receives more realistic information regarding spatial grade connectivity, the LOM production schedule, such as in Case 1, pursues high-grades more efficiently.Although the mill's throughput is kept reasonably constant throughout the LOM in both instances, the dissimilarities regarding metal content are stressed in Figure 10.Case 1 can feed the mill with a higher head grade for the majority of the LOM, as shown in Figure 10a.As the optimizer encounters better-connected zones of high-grade, it can bring their extraction to the same period so they can be processed together, increasing the average feed grade at the mill and recovering more ounces earlier, as shown in Figure 10b.Case 1 shows an ounces profile that is consistently higher for the first 17 years; this difference reaches 7% after the 10th year.Case 2 produces, after the 20th year, 2% more gold, but this is not significant due to the effect of discounting and the time value of money.Recovering more ounces sooner brings more cash flow earlier to the operation, which positively impacts the net present value (NPV).Summing up the joint effects of correctly meeting production targets, mining less waste and producing more gold earlier results in a considerable increase in NPV, as shown in Figure 10c.By producing more metal and less waste, the LOM production schedule obtained in Case 1 generates in a total of 5% higher NPV than Case 2, and 16% higher in the initial ten years.The difference is substantial and improves financial returns at the early stages of the Once the production schedules are generated, the risk assessment of achieving production forecasts is calculated by pushing another set of simulations through each LOM schedule and by evaluating related forecasts.The results are presented in terms of P10 and P90, representing the 10th and 90th percentiles of related performance indicators, respectively.Regarding production targets and forecasts, Figure 9a shows the total tonnage mined over the LOM and the cumulative strip ratio for both cases.The production schedule obtained in Case 2 represents, in total, 5% more material mined than in Case 1.This difference reaches 8% at the end of the 10th year to ensure a similar throughput at the mill, Figure 9b.Mining more, in this case, translates to higher waste production, which is quantified by the higher strip ratios presented by Case 2. This can be explained by the spatial disorder (maximum entropy) that Gaussian-based approaches generate with respect to high-grades.Having ore blocks less connected in space forces the optimizer to mine more dispersed high-grade values (mining blocks), so as to provide a consistent feed rate to the mill; this also leads to different total tonnages of materials mined and the differences in UPL shapes, shown in Figure 8.On the other hand, if the optimizer receives more realistic information regarding spatial grade connectivity, the LOM production schedule, such as in Case 1, pursues high-grades more efficiently.
Although the mill's throughput is kept reasonably constant throughout the LOM in both instances, the dissimilarities regarding metal content are stressed in Figure 10.Case 1 can feed the mill with a higher head grade for the majority of the LOM, as shown in Figure 10a.As the optimizer encounters better-connected zones of high-grade, it can bring their extraction to the same period so they can be processed together, increasing the average feed grade at the mill and recovering more ounces earlier, as shown in Figure 10b.Case 1 shows an ounces profile that is consistently higher for the first 17 years; this difference reaches 7% after the 10th year.Case 2 produces, after the 20th year, 2% more gold, but this is not significant due to the effect of discounting and the time value of money.Recovering more ounces sooner brings more cash flow earlier to the operation, which positively impacts the net present value (NPV).Summing up the joint effects of correctly meeting production targets, mining less waste and producing more gold earlier results in a considerable increase in NPV, as shown in Figure 10c.By producing more metal and less waste, the LOM production schedule obtained in Case 1 generates in a total of 5% higher NPV than Case 2, and 16% higher in the initial ten years.The difference is substantial and improves financial returns at the early stages of the development of the mine.

Conclusions
This paper investigates the effects of using high-order simulations of an ore deposit in the simultaneous stochastic optimization of a gold mining complex.The high-order simulations are generated directly at the block-support scale and are used as inputs to the simultaneous optimization (Case 1).The optimized LOM production schedule generated is benchmarked against a case where

Figure 1 .
Figure 1.Connectivity of high-grades along X (a) and Y (b) direction, calculated for the 99th percentile from the data distribution.Figure adapted from Minniakhmetov et al. [49].P10 and P90 represent the 10th and 90th percentiles of the reported values, respectively.

Figure 1 .
Figure 1.Connectivity of high-grades along X (a) and Y (b) direction, calculated for the 99th percentile from the data distribution.Figure adapted from Minniakhmetov et al. [49].P10 and P90 represent the 10th and 90th percentiles of the reported values, respectively.

Figure 2 .
Figure 2. The grade-tonnage curve of the gold deposit for SGS (second-order) and high-order simulations.

Figure 3 .
Figure 3. Cross-sections of the simulations, obtained from high-order (a) and second-order (b) based methods, highlighting differences in connectivities of high-grades.The red circles highlight the highgrade zones.

Figure 2 .
Figure 2. The grade-tonnage curve of the gold deposit for SGS (second-order) and high-order simulations.

Figure 2 .
Figure 2. The grade-tonnage curve of the gold deposit for SGS (second-order) and high-order simulations.

Figure 3 .
Figure 3. Cross-sections of the simulations, obtained from high-order (a) and second-order (b) based methods, highlighting differences in connectivities of high-grades.The red circles highlight the highgrade zones.

Figure 3 .
Figure 3. Cross-sections of the simulations, obtained from high-order (a) and second-order (b) based methods, highlighting differences in connectivities of high-grades.The red circles highlight the high-grade zones.

Figure 4 .
Figure 4. Connectivity of simulated realizations at (a) NE direction; (b) NE/45° direction.The arrow highlights the gap between the outcomes from the two methods.

Figure 5 .
Figure 5. Flow diagram of the mine complex configuration.

Figure 4 .
Figure 4. Connectivity of simulated realizations at (a) NE direction; (b) NE/45 • direction.The arrow highlights the gap between the outcomes from the two methods.

Figure 4 .
Figure 4. Connectivity of simulated realizations at (a) NE direction; (b) NE/45° direction.The arrow highlights the gap between the outcomes from the two methods.

Figure 5 .
Figure 5. Flow diagram of the mine complex configuration.

Figure 5 .
Figure 5. Flow diagram of the mine complex configuration.

Figure 6 .
Figure 6.Cross-section of the life-of-mine (LOM) production schedule (plane East-West) obtained in (a) Case 1 and (b) Case 2. The red arrow shows the preferential mining direction along the direction of higher continuity of high-grades.

Figure 6 .
Figure 6.Cross-section of the life-of-mine (LOM) production schedule (plane East-West) obtained in (a) Case 1 and (b) Case 2. The red arrow shows the preferential mining direction along the direction of higher continuity of high-grades.Minerals 2019, 9, x FOR PEER REVIEW 10 of 16

Figure 7 .
Figure 7. Cross-section of the LOM production schedule (plane North-South) obtained in (a) Case 1 and (b) Case 2. The red arrows show the preferential mining direction in the direction of higher continuity of high-grades.

Figure 7 .
Figure 7. Cross-section of the LOM production schedule (plane North-South) obtained in (a) Case 1 and (b) Case 2. The red arrows show the preferential mining direction in the direction of higher continuity of high-grades.

Figure 7 .
Figure 7. Cross-section of the LOM production schedule (plane North-South) obtained in (a) Case 1 and (b) Case 2. The red arrows show the preferential mining direction in the direction of higher continuity of high-grades.

Figure 8 .
Figure 8. Horizontal sections of the mine production schedule at different elevations for (a) Case 1 and (b) Case 2. The red circles highlight the differences in pit sizes.

Figure 8 .
Figure 8. Horizontal sections of the mine production schedule at different elevations for (a) Case 1 and (b) Case 2. The red circles highlight the differences in pit sizes.

Figure 9 .
Figure 9. Production target assessments: (a) Cumulative mine tonnage and strip ratio; (b) mill throughput with upper and lower bounds.P10 and P90 representing the 10th and 90th percentiles, respectively.

Figure 9 .
Figure 9. Production target assessments: (a) Cumulative mine tonnage and strip ratio; (b) mill throughput with upper and lower bounds.P10 and P90 representing the 10th and 90th percentiles, respectively.

Figure 10 .
Figure 10.Risk analysis of production forecasts for (a) head grade at the mill; (b) cumulative gold recovered at the mill and (c) NPV assessment.P10 and P90 representing the 10th and 90th percentiles, respectively.

Figure 10 .
Figure 10.Risk analysis of production forecasts for (a) head grade at the mill; (b) cumulative gold recovered at the mill and (c) NPV assessment.P10 and P90 representing the 10th and 90th percentiles, respectively.
and    a realization of   (  ).

Table 1 .
Main parameter used in the optimization.

Table 1 .
Main parameter used in the optimization.

Table 1 .
Main parameter used in the optimization.