Assessing geologic model uncertainty - a case study comparing methods

Abstract

Evaluating mineral resources requires the prior delimitation of geologically homogeneous stationary domains. The knowledge about the ore genesis and geological processes involved are translated into three dimensional models, essential for planning the production and decision-making. The mineral industry usually considers grade uncertainty for resource evaluation; however, uncertainty related to the geological boundaries are often neglected. This uncertainty, related to the location of the boundary between distinct geological domains can be one of the major sources of uncertainty in a mineral project, and should be assessed due to its potential impact on the ore tonnage, and consequently, on enterprise profitability. This study aims at presenting three different methodologies capable of generating multiple geomodel realizations and thus, assessing uncertainty. A real dataset with high geological complexity is used to illustrate the methodology. The results are compared to a deterministic model used as a reference scenario.

keywords:
geological model; multipoint geostatistics; implicit modeling; uncertainty

1. Introduction

The geological model is one of the main tools used to evaluate mineral resources and their viability as mineral reserves. Comprehending geological formation processes of the mineral deposit and translating it into the model is essential for planning and decision making.

Each subsurface geological unit or stationary domain is characterized by a volume and shape. The contacts between these different geological units are determined by complex discontinuities generated by intrusions, erosions, faults or different depositional processes. There is uncertainty in the limit’s determination in the subsurface, which directly impacts the estimated ore quantity and quality. Mining project success relies on the information about ore volume and content; in this sense, the geological model’s uncertainty must be assessed.

Construction of tridimensional geological models begins with the interpretation of the mineral deposit, studies on the geological formation process and geological data integration. Traditionally, geological models are created explicitly, from interpreted geological sections and surface maps. Different domains are digitized in a time consuming and laborious process, which does not permit prompt model updates as new data is acquired. The explicit methodology transfers into the model the geomodeler’s judgment and expertise. The result is a single subjective model that does not allow reproduction (Silva, 2015). Due to the time constraint, it is rarely possible to create different explicit models to assess the uncertainty associated with the geological model.

Implicit methods emerged to replace manual digitizing with automatic procedures to overcome the disadvantages associated with the traditional methodology. Implicit methods are based on the creation of a continuous mathematical representation of an attribute on a given volume (Cowan et al., 2003). These methods are classified as deterministic, in which only one model is generated, and stochastic, in which multiple possible scenarios for the mineral deposit are generated. Implicit deterministic methods are straightforward, flexible and fast, drastically reducing the time required to construct geological models, along with new information incorporation being fast and easy. Implicit established methods are: discrete smoothed interpolation (Mallet, 2002), potential field method (Chilés et al., 2004, Calcagno et al., 2008) and signed distance interpolation by radial functions (RBF) (Cowan et al., 2002; Cowan et al., 2003). All implicit methodologies share the same mechanics where the only difference lies in the volume function.

Signed distances are by far the most used volume function for implicit modeling. Signed distances are calculated as the Euclidean distance between a sample and the nearest sample coded with an opposite indicator (belonging to a distinct lithotype), calculated distances are signed to indicate if the sample belongs to the domain being modeled. Distance values are interpolated to unsampled locations, and positive or negative values for interpolated distances determine whether a location belongs or not to a given domain. This method is intended to generate binary models in which the contact is defined by the zero isosurface. Silva and Deutsch (2012) presented an extension of this technique to multiple geological domains modelled simultaneously, by retaining the category responsible for the minimum estimated distance.

Scarce data precludes the exact location of the contacts. The exact location is unknown in all its extent, and that is why the uncertainty of the contacts must be assessed. Deterministic models cannot assess contact uncertainty. However stochastic techniques based on geostatistical simulation, such as indicator simulation and truncated Gaussian simulation (Journel, 1989; Matheron et al. 1987), can. These techniques require great computational effort and have been developed to build multiple outcomes for the geological domains.

Distance function can be modified to assess contact uncertainty. Cárceres et al. (2011) present a methodology in which distance values are transformed into normal distributions conditionally simulated. Truncating distributions at locations where distance is equal to zero results in the contact variability. Wilde and Deutsch (2011) proposed the calibration of a parameter C to generate a band of uncertainty. Within this region (uncertainty band), uniform values are simulated, which are then compared against the interpolated values. As a result, multiple realizations of the geological model are generated.

Another alternative to generate multiple geological domain realizations, is based on the so-called Multipoint Statistics Simulation (MPS) that uses a training set to create conditioning probabilities from where scenarios of geological events are drawn. Each scenario aims at reproducing the geological characteristics of the training images. The method was initially proposed by Guardiano and Srivastava (1993), but several variations have been developed based on the initial workflow, such as the SNESIM algorithm developed by Strebelle and Journel (2001) and the method proposed by Silva and Deutsch (2012) that uses multiple training images. In order to access the uncertainty of the contacts, Boucher et al. (2014) proposed the contactsim algorithm that allows to disturb the contacts between domains in a zone obtained from a reference model. Although efficient, the method demands conceptual models or realistic training images that require effort for their production or in some cases, they cannot be built due to the lack of information.

This study investigates three methodologies used to assess geologic model uncertainty, applying them to a real iron ore deposit.

2. Methodologies

Three methodologies for uncertainty assessment are presented. The first one consists of a multipoint statistics simulation algorithm aiming at reproducing lithology proportions, spatial distribution, shape, volume and contact patterns. The other two methodologies are based on distance functions to generate implicit models and the related uncertainty.

3. SNESIM multipoint geostatistics

Models of natural phenomena are usually done through images constructed from interpretation, representation and modeling of these processes. In multipoint geostatistic methods, a training image (TI) provides information to serve as a repository of the spatial patterns and their probabilities found in the phenomenon.

One of the first MPS algorithms was presented by Guardiano and Srivastava (1993), based on the sequential simulation paradigm, in which the inference of the conditional probabilities from which the simulated values are drawn, are obtained through conditional proportions obtained scanning the TI for patterns. This training image scanning process is performed from a template or search neighborhood τ_J , composed of a set of J vectors {h_j, j =1, …, J }, that informs which TI points will be considered relevant in relation to a central point u₀ of the grid node to be simulated. Figure 1 presents an example of a two-dimensional search neighborhood with four vectors. The simulation draws a value from a conditional probability mass function (cpmf) at an unsampled location of the simulation grid. As this method requires scanning the training image for each node of the grid to be sampled during the simulation, it generates a large computational demand.

Figure 1
Example of template or search neighborhood with 4 vectors (Rasera, 2014).

As an alternative, Strebelle (2002) proposed to scan the training image only once, and storing the TI patterns in a file structure known as the search tree. This file allows access to the information needed for simulating each unsampled node. The SNESIM algorithm developed aggregates two parts: construction of the search tree in which all the proportions taken from the training image are stored, and the part of the simulation in which the proportions are read and used to retrieve simulated values.

Considering a random variable S, which can belong to K classes {s_k, k =1, …, K } in which S represents different categories. The training image is the representation of how the values of a given property s_k are distributed and connected in space (Strebelle, 2002). Thus, the training image is considered a non-conditional representation of the random function S(u), a prior conceptual representation of the spatial distribution values of the variable s_k, which does not necessarily need to honor the location of the sample set n(u) of s_k.

This algorithm can be summarized in a series of steps: first scan the training image to generate a search tree at a chosen size (template) τ_J ={h _α, α=1, …, n }. Then, define a random path to visit each location u only once and generate a conditional cumulative distribution function (ccdf) at location u_j using proportions extracted from TI. Randomly draw one value from the ccdf as a possible category at location u_j and use it as a conditioning value for other locations. Repeat the process until all grid nodes are visited. Multiple realizations are obtained by repeating the process using different random paths.

The image on the left in Figure 2 shows a training image with two categories: white and black, and the search template for the central region in gray, referring to the unknown location to be simulated. Adjacent locations (local neighborhood) listed as 1 through 4 refer to known or unknown categorical data, defined as a data event. These numbers are the conventional order to set up the search tree. The modeler informs the maximum number of points and the geometric parameters of an ellipsoid, indirectly defining the template: dimensions in terms of its axes and rotation angles. The training image with two categories is scanned for the construction of the probability distribution, considering the search template.

Figure 2
Training image for two categories and search template shown on the left. On the right the search tree obtained from the training image and the template.

The data structure called search tree is defined within a neighborhood to reduce the size of the object searched (number of pixels), optimizing computational time for assembling cpmf, which is used to draw variables at a given location. The size of the template controls the complexity of the search tree and the memory and time spent during the simulation. The number of categories and the size of the TI also affect the size of the search tree.

Figure 2 on the right shows a search tree constructed from the training image and the template. Note that each level corresponds to a region of the search neighborhood. The search tree begins with the determination of the root, as an empty data event, in which the associated frequency corresponds to the overall proportion of black or white categories found in an eroded TI. This image is considered eroded by being equivalent to the largest subset of the training image, in which the template is fully contained. The variables B and P indicate the number of occurrences of the categories at the unknown location and the frequency of each pattern given from the eroded TI. In this way, the search tree stores the number of occurrences of each category in relation to the gray central node of the template, which results in the proportion of occurrence of the two categories.

4. Implicit metods

Implicit modeling techniques aim at assisting the geological interpretation and the construction of models, by demarcating the boundaries between domains in an objective way, allowing fast and accurate reproduction.

Some methodologies based on distance functions and simulations generate multiple geological models and make possible the uncertainty assessment.

4.1 Distance function

Distance functions used to model geological domains are based in the Euclidian distance between one sample and another sample coded as a different indicator, which makes it necessary to previously code the database as belonging or not to the geological domain being modeled. The value of the distance function is given by the shortest distance between two samples, which are coded as opposite indicators. If the location belongs to the domain being modeled, the value of the distance function is negative, otherwise positive, according to the equation:

$\begin{array}{c}\mathit{df}\left({\mathit{u}}_{\alpha }\right)=+\left({\mathit{u}}_{\alpha }-{\mathit{u}}_{{\alpha }^{\prime }}\right),\mathrm{if}\mathit{i}\left({\mathit{u}}_{\alpha }\right)=0\hfill \\ \mathit{df}\left({\mathit{u}}_{\alpha }\right)=-\left({\mathit{u}}_{\alpha }-{\mathit{u}}_{{\alpha }^{\prime }}\right),\mathrm{if}\mathit{i}\left({\mathit{u}}_{\alpha }\right)=1\hfill \end{array}$

The C parameter controls the uncertainty bandwidth by increasing the difference between positive and negative distance function values. It works by adding C to positive values outside the domain and subtracting from negatives inside the domain:

$\begin{array}{c}\widehat{\mathit{df}}\left({\mathit{u}}_{\alpha }\right)=\mathit{df}\left({\mathit{u}}_{\alpha }\right)+C,\mathrm{if}\mathit{i}\left({\mathit{u}}_{\alpha }\right)=0\hfill \\ \widehat{\mathit{df}}\left({\mathit{u}}_{\alpha }\right)=\mathit{df}\left({\mathit{u}}_{\alpha }\right)-C,\mathrm{if}\mathit{i}\left({\mathit{u}}_{\alpha }\right)=1\hfill \end{array}$

There are three ways to determine the value of the C parameter: empirical determination, partial calibration and full calibration. The empirical calibration is based on sample spacing and on the deposit`s geological knowledge. Interpreting geological sections and comparing interpolated models enables partial calibration. The full calibration, which is time-consuming and laborious, compares the multiple distance function interpolated models (reference models) to determine the C parameter (Karpekov & Deutsch, 2016).

4.2 Distance function interpolation

Once the signed distance is calculated, interpolation permits obtaining the distance function at all grid nodes. Thus, interpolating a modified distance function by the parameter C generates the uncertainty zone between values of -C and +C on the interpolated grid. Any interpolation method could be used such as: inverse of distance, ordinary global kriging, locally varying anisotropy (LVA) kriging or radial basis functions produce satisfactory results. Global interpolation algorithms generate smooth and artifact free models, and in this sense are recommended. However, due to the high computational demand required, global algorithms cannot be applied to all cases (Wilde & Deutsch, 2011). Ordinary kriging could be used as an alternative, which makes it necessary to consider a large search neighborhood and great amount of data within this neighborhood.

The distance function variogram exhibits non-stationary behavior, i.e., it does not reach a sill. The Gaussian model is indicated due to its parabolic behavior near the origin, which contributes to defining smooth boundaries between domains when these types of variograms were used combined with ordinary kriging.

4.3 Contact simulations

There are two main ways to generate multiple boundary realizations: direct signed distance conditional simulations and non-conditional simulations compared against interpolations in the uncertainty zone.

Multiple contacts generated by direct signed distance conditional simulations in fine grids use a large search neighborhood and a great amount of data. The methodology is very sensitive to simulation parameters and can lead into unrealistic models, very softened or with the excessive presence of salt and pepper texture, which is not geologically realistic.

Non-conditional simulations generate domain limits by comparing each realization with the interpolated value within the uncertainty zone. For such purpose, it is necessary that the simulation be performed uniformly between - C and + C, and the standard normal value, y' (u), be simulated and transformed by:

$\mathit{d}{\mathit{f}}^{\prime }\left(\mathit{u}\right)=2*\mathit{C}*{\mathit{G}}^{-1}\left({\mathit{y}}^{\prime }\left(\mathit{u}\right)\right)-\mathit{C}\forall \mathit{u},\forall \mathit{l}$

where df' (u) is the value of the simulated distance function, y' (u) is the standard normal value of the non-conditional simulation, and G^-1 is the determination of the value of the normal standard cumulative distribution corresponding to y' (u). To ensure that the values belong to the determined uncertainty zone, simulated values are multiplied by 2C and subtracted from C.

To classify a region as belonging or not to the domain of interest, a comparison is performed between the interpolated and simulated values. Interpolated values smaller than simulated values represent locations inside the domain, otherwise, the location is outside the domain. The limit is defined were the simulated and interpolates values are the same, as demonstrated in the Figure 3 below:

Figure 3
Sites classification comparing estimated and simulated values (Modified from Wilde e Deutsch, 2011a).

5. Case study

The case under study was conducted in one of the biggest iron ore deposits in the world. The dataset comprehends 6924 samples obtained from 718 core samples and 44236 data points obtained from interpreted geological sections. Different lithotypes present in the deposit were grouped into four categories that represent the main aspects of the mineralization: hematite in category 1, structural laterite in category 2, jaspelite in category 3 and waste in category 0 (Figure 4). The deterministic model was made available by the iron ore mine staff and was used as reference and as the training image (Figure 5). It was constructed from mixed methods combining indicator kriging and explicit modeling.

Figure 4
The figure shows the location of the samples; the image of the upper region shows the plan view of the drilling hole collars and geological sections. The bottom image shows a longitudinal vertical view of the data.

Figure 5
The upper figure shows the reference model of the iron ore mineralization. The zone of uncertainty with 200m is shown in the lower image.

The uncertainty on boundaries between geological domains can be assessed by delimiting a region in which the occurrence of the contact is likely to occur, generating multiple possible scenarios at this marked location. There are different ways of determining the uncertainty zone. Boucher (2004) determines the zone of uncertainty considering the limits between domains in the training image. Souza (2017) determines this region considering a coefficient U of the smallest signed distances, which reflects contact uncertainty. Another alternative suggested by Munroe and Deutsch (2008) is to modify the distance function by a parameter C, interpolate this modified function for all nodes of the grid and select the region between -C and +C signed distance values as the uncertain zone. In this case study, the zone was delimited based on the contacts from the reference model, i.e. a zone within 200m was established to access the model uncertainty and compare methodologies fairly (Figure 5).

The first method applied consists in quantifying the uncertainty by using the SNESIM algorithm within the uncertainty zone. It is necessary to erode the geological reference model in the uncertainty region and use the training image to acquire domain patterns within this region. Multiple contacts are generated simulating the categorical property inside the region using the patterns obtained previously. For each realization of the geological model, it is necessary to merge the SNESIM realization with the frozen blocks, outside the uncertainty zone.

In order to facilitate the comparison between the methodologies, similar parameters were used for all the methodologies. In this sense, the search strategy parameters and grid characteristics are shown in Table 1. As a result of the simulations, 20 visually satisfactory realizations were generated, realistically reproducing the domains of each category.

Table 1
Main parameters used in simulation and kriging: search radius, rotation, number of samples and grid dimensions.

Number of samples			Kriging						Grid Dimension
Max.	Min.	Oct.	Az	Dip	Rake	Radius Max.	Radius Int.	Radius Min.	Dim. X	Dim. Y	Dim. Z
45	3	3	90	0	0	2000	1500	1000	25	25	15
Number of samples			Simulation						Grid Dimension
Max.	Min.	Oct.	Az	Dip	Rake	Radius Max.	Radius Int.	Radius Min.	Dim. X	Dim. Y	Dim. Z
45	3	3	90	0	0	850	500	350	12.5	12.5	12.5

Figure 6 shows the visual validation for realization #15, in which sample points are confronted with three sections from the model evidencing that the data were honored during the simulation. In addition, the fluctuation of the proportions of each geological domain remained similar to the reference model.

Figure 6
Easting slices from realization #15 demonstrating that the data points matched with SNESIM model. Note that the transitional zones also respect the data points. In the images, the dark blue represents waste, the light blue hematite, the yellow superficial laterite and red jaspillite.

Implicit models aim at providing objective, flexible models that allow reproducibility. Both of the next techniques are based on distance functions to generate multiple realizations for the contacts between domains.

The first implicit technique directly simulates the signed distances within the region of interest. First, it is necessary to calculate the signed distances for all sample locations and for each category. A spatial continuity analysis is needed for the calculated distances, and the results obtained are shown in Table 2. The signed distance variogram shows, as expected, nonstationarity behavior. Thus, the variogram tends to increase exponentially with the pair’s separation distance. This feature may make the definition of the variogram distance function model somewhat arbitrary, primarily in the definition of its range.

Table 2
Summary of spatial continuity of data sets. In the upper table, the models adopted for the variography of the signed distances for the direct signed distance simulation methodology are presented. The signed distances calculated for each category are represented by the data sets SD-0, SD-1, SD-2 and SD-3. The table below presents the models adopted for the boundary simulation methodology. The three subsets created by the hierarchical approach are represented by G-I, GII and G-III.

Direct signed distance simulation
Data	Structures			Range
	N. E.	Gauss 1	Gauss 2	R1 90|00	R1 00|00	R1 00|90	R2 90|00	R2 00|00	R2 00|90
SD-0	0.00	0.67	0.33	6686.58	1461.14	1005.81	6710.41	1461.14	1029.30
SD-1	0.00	0.01	0.99	6397.99	1759.88	1181.48	6397.99	1759.88	1181.48
SD-2	0.00	0.22	0.78	10625.10	1809.18	568.64	10710.09	5480.20	568.64
SD-3	0.00	0.34	0.66	6508.65	1366.44	1119.57	6521.02	1366.44	1826.35
Boundary Simulation
Data	Structures			Range
	N. E.	Gauss 1	Gauss 2	R1 90|00	R1 00|00	R1 00|90	R2 90|00	R2 00|00	R2 00|90
G-I	0.00	0.33	0.67	11782.11	1245.42	513.78	11990.13	15168.39	513.95
G-II	0.00	0.99	0.01	5683.50	1200.42	784.80	5723.59	1201.06	779.10
G-III	0.00	0.42	0.58	964.95	1021.68	614.65	22448.82	10181.35	533.14

Conditional simulation is performed within the zone of uncertainty using the calculated distances as conditioning data. The Sequential Gaussian Simulation method (Goovaerts, 1997) was selected for the process, which resulted in twenty distance function realizations for each category in the total of eighty realizations performed using the GSLIB software (Deutsch & Journel, 1992). The parameters used in the simulation are shown in Table 1.

Each realization within the uncertainty zone is constructed by selecting the lowest (most negative) value for the simulated distance functions on each block. Considering, for example, the first realization of each simulated distance function, for each grid node distance function, the simulated value is compared and the category that represents the most negative value assigned to the node results in one realization of the model in the uncertainty zone.

As before, to obtain each realization of the geological model, it is necessary to merge the direct signed distance simulation realization with the pre-defined (frozen) blocks outside the uncertainty zone. Twenty realizations obtained realistically reproduced the domains of each category. Visual validation of realizations demonstrated that the sampling points were honored during the simulations. Realization #10 validation is shown in Figure 7. The fluctuation of the proportions of each geological domain remained similar to the reference model, presenting small difference among models.

Figure 7
Easting slices from realization #10 demons trating that the data points matched with direct signed distance simulation realization. All data points are honored in the model realization. Note that even the isolated data were represented in small structures in the model. In the images, the dark blue represents waste, the light blue hematite, the yellow superficial laterite and red jaspillite.

In the second methodology of boundary simulation, the signed distance interpolation is compared against non-conditional simulations within the uncertainty region, resulting in multiple realizations. In the previous methodology, the calculation of the signed distances is straightforward; however, the unconditional methodology using parameter C requires a hierarchical approach when dealing with multiple lithologies (Amarante et al., 2018). Normally, this hierarchical approach is done geochronologically, that is, the domains are constructed from the oldest to the most recent. However in this case, hierarchy was established by the amount of information and structural patterns. In this way, the categories were divided into subsets according to the hierarchy shown in Figure 8.

Figure 8
The hierarchical division of the database and construction of boundary realizations.

The same 200 m value, selected for the size of the uncertainty zone, was used as parameter C. Consequently, the distance functions previously obtained were modified, adding C on samples outside the domain and subtracting C from samples inside the domain.

Spatial continuity analysis from the modified distance functions were performed individually for each data set. Each data group has different variograms in order to realistically represent the spatial continuity of the domains represented by each group. Table 2 summarizes the directional variograms obtained, adopting the direction of greater continuity N90º / dip 00º.

As the contact is simulated, it is recommended to use a Gaussian model variogram for the non-conditional simulations, as it allows the short-range continuity to be reproduced (Wilde & Deutsch, 2011). Also, it is recommended to use a small nugget effect to stabilize the mathematical calculation. Contact roughness can be controlled by the variogram range.

The modified distances were interpolated by ordinary kriging using the same grid of the reference model. Non-conditional simulations were performed within the zone of uncertainty, using the Sequential Gaussian Simulation (Goovaerts, 1997) implemented at GSLIB (Deutsch & Journel, 1992). Although the simulation is not distance conditioned, the simulation used the same spatial search strategy used in the previous methods. The parameters used in kriging of each set of data and used in the non-conditional simulation are shown in Table 1.

In this last method multiple realizations were obtained by classifying each grid node within the uncertainty zone by comparing interpolated and simulated values. When the interpolated value is less than the simulated value, the location is assigned as internal to the domain, and if the interpolated value is greater than the simulated value, the location is assigned as external to the domain. This classification is performed for all grid nodes in the zone, contemplating all contacts between the lithologies separated in pairs.

The construction of the boundary realizations follows the same criteria established before, grouping lithologies in one model following the hierarchy and generates the realization within the uncertainty zone. Each realization of the geological model is constructed by merging the boundary realization with the frozen blocks outside the uncertainty zone. Realization #15 validation is shown in Figure 9 and demonstrates that the sampling points were honored during the simulation.

Figure 9
Easting slices from realization #15 showing that data points match with boundary realization. Note that even in transition zones, the model generated honors the data. In the images, the dark blue represents waste, the light blue hematite, the yellow superficial laterite and red jaspillite.

The spatial continuity from the geological domains directly affects the evaluation of mineral resources, since it determines the volume of ore available. Thus, the impact is significant in estimating the ore tonnage above cutoff, leading to biased forecasts, errors in mine planning, unexpected costs in the operation and deviations in expected revenues (Srivastava, 2005). The bar graphs (Figure 10) describe the volumetric variability obtained with the several realizations constructed by the three methods: SNESIM, direct signed distance simulation and boundary simulation.

Figure 10
The image graph shows the variability obtained by comparing realizations of each methodology. The smaller and larger volumes of blocks of each category obtained using the different methodologies are compared with the reference model Lito 4.

The volumetric variability obtained in the simulation can be evaluated by comparing the maximum and minimum number of blocks for each category in the zone of uncertainty, as shown in Figure 10. The deterministic geological model is used as reference for model proportions within the uncertainty zone. Table 3 presents the maximum and minimum number of blocks for each category as well as the values of the difference obtained in the simulations.

Table 3
Summary of categories volume variability obtained in each methodology.

	Boundary Simulation			Direct signed distance simulation			SNESIM
	Max.	Min.	Difference	Max.	Min.	Difference	Max.	Min.	Difference
Waste	243889	203575	40314 (17.37%)	226512	226409	103 (0.04%)	223159	219459	3700 (1.59%)
Hematite	173394	145371	28023 (18.5%)	158867	158688	179 (0.12%)	158845	154995	3850 (2.54%)
Laterite	70933	53678	17255 (31.18%)	61188	61096	92 (0.17%)	67123	65593	1530 (2.76%)
Jaspelite	230412	181207	49205 (22.96%)	206867	206705	162 (0.08%)	210759	205556	5203 (2.43%)

The volumetric variability obtained when applying SNESIM is represented by the red bars. Waste and jaspilite show a small relative variability of 1.59% and 2.43%, respectively while hematite and laterite presented higher variability of 2.54% and 2,76%, respectively. A distinct behavior was observed in the waste and laterite, where the variability obtained in the realizations is far from the volume presented in the reference model. The process systematically overestimated waste over laterite, which was underestimated in all realizations.

The green bars represent the values obtained in the direct signed distance simulation methodology, which shows the smallest volumetric difference obtained when comparing the realizations and the reference model. The four categories presented a small variability among realizations, the highest variability obtained belongs to the laterite of only 0.17%. When comparing the difference between realizations with the reference model, it is noted that the difference among realizations is smaller than the difference between these and the reference model.

The results of the boundary simulation methodology are shown by the blue bars, in which the significant difference between max / min values is noted, reaching a variation of 31.18% for laterite and 22.96% for jaspilite, and a smaller variability for waste of 17.37% and hematite of 18.5%. This variability demonstrates the need to evaluate the uncertainty and the possible impact of this variability on the mineral project. It is noteworthy that the volume of each reference model domain belongs to the range of variability of the realizations, which makes the variability consistent for the evaluation of the uncertainty.

The geological sections presented in Figure 11 and Figure 12 allow to compare models generated by the three methods and the Lito 4 reference model. The main objective in the construction of geological models is to obtain meaningful representations from the geobodies and subsurface structures. The realizations generated by all methodologies presented realistic structures, as demonstrated in the easting section slice #50 (Figure 11) and northing sections slice #180 (Figure 12).

Figure 11
The image presents the comparison between easting model sections slice #50 and the reference model, the first image on top. The second image represents the boundary simulation #15 realization, the third image the direct signed distance simulation #15 realization, and the image below the SNESIM #15 realization.

Figure 12
The image presents the comparison between northing model sections slice #180 and the reference model, the first image on top right. The second image on top left represents the boundarysimulation #10 realization, the third image on bottom right the direct signed distance simulation #10 realization, and the SNESIM #10 realization on the bottom left image.

Note in Figure 11 that in region 1, demarcated by the circle, there is the lack of intrusive rocks present in the reference model and absent in the other models. The models generated by signed distances were insensitive to some small structures present in the reference model, such as small intrusions and jaspelite lenses demonstrated in circle 2. In order to reproduce small-scale structures, SNESIM generates innumerable intrusive bodies shown in circle 3, however, many of them are without geological meaning as isolated geobodies and scattered blocks.

Figure 12 shows northing section slice #180 of the reference model and realization #10 of each method. Note in 1 that the intrusive rock structure present in the reference model is not reproduced by any of the methodologies. Most of the geological structures of the reference model were well reproduced, however a few structures generated during the simulation do not present realistic geological characteristics for this deposit, such as the laterite discontinuous on the surface at 2, waste isolated and dispersed blocks at 3.

It is possible to construct maps or probability models from multiple geological domain realizations. In these models, we can identify the probability of a given location, belonging or not to a particular geological domain. Figure 13 shows the probability models generated using realizations from the three methods, one image for each category varying from 0% in blue to 100% in red as the probability of that category.

Figure 13
The figure shows sections of probability maps for geological contacts in which each column represents a given methodology. Each line represents a distinct category, starting from waste in the upper images, followed by the hematite and the laterite, ending with the three images from jaspilites. The color scale represents the associated probability of occurrence of a certain category, varying from 0% in blue to 100% in red.

Note in the upper row of waste on the left image the presence of the intrusive rock in the SNESIM model, a feature that is absent in the realizations from other methodologies. In the second row, the variability of the hematite is depicted. It is worth mentioning the irregular variability within the domain present in the SNESIM methodology in the right image. The images present in the third row show the variability of the laterite, highlighting the high variability shown by the boundary simulation methodology in the right image. The bottom row shows the variability of the jaspilites that makes clear the presence of the intrusive rocks that crosscut the reference model but in the realizations, it only appears evidently in the base of the models. The images of the central column referring to the direct signed distance simulation realizations demonstrate the small difference obtained by the direct simulation of the signed distances.

6. Conclusions

The three methods were able to assess the uncertainty of the geological model and satisfactory results were obtained. The methods tested were valid alternatives for mineral deposit geologic modeling and uncertainty assessment. Each method presents distinct advantages when applied to different stages of the mining project.

The multipoint geostatistical method allows to incorporate the expertise into interpretation when building geological models used as training images. However, the construction of the training image requires a great amount of information and a great computational effort, which is not normally what the geological model specialist is looking for. In this way, the use of this method is indicated for situations where there is a great amount of information about the deposit to adequately condition the simulations.

The implicit methods are indicated for all phases of the mineral deposit modelling, especially for the initial phases in which there is not much information about the deposit. The construction of models by implicit methods allows one to drastically reduce the work time and the subjectivity into the geological model, since it is possible to reproduce the model from the same initial parameters. Moreover, as demonstrated by the case study, it allows to evaluate the uncertainty of the model through construction of multiple scenarios for the geological contacts.

The use of automatic or semi-automatic methods does not completely replace the experienced geomodeler. However, its use can help reduce the time needed to construct the geological model and allow the uncertainty associated with the model to be quantified.

Geological modelling is crucial for the mineral deposit evaluation, since estimates and simulations of mineral content are performed in the pre-established domains. The uncertainty in determining the boundaries between geological domains in subsurface impacts on the quantity and quality of the mineral of interest, as well as its location. In this way, this uncertainty must be evaluated. This article presented three different methodologies that allow uncertainty assessment satisfactorily.

References

• AMARANTE, F. A. N., ROLO, R. M., AND COSTA, J. F. C. L. Incerteza associada ao modelamento implícito e hierárquico. In: CONGRESSO BRASILEIRO DE MINAS A CÉU ABERTO E MINAS SUBTERRÂNEAS (CBMINA), 9 Belo Horizonte, 2018.
• BOUCHER, A. Considering complex training images with search tree partitioning. Computers and Geosciences, n. 35, p. 1151-1158, 2009.
• BOUCHER, A., COSTA, J. F., RASERA, L. G., AND MOTTA, E. Simulation of geological contacts from interpreted geological model using multiple-point statistics. Mathematical Geosciences, v. 46, n. 5, p. 561-572, 2004.
• CÁCERES, A., EMERY, X., AEDO, L., GÁLVEZ, O. Stochastic geological modeling using implicit boundary simulation. In: INTERNATIONAL SEMINAR ON GEOLOGY FOR THE MINING INDUSTRY GEOMIN, 2. Procedings... 2011. p. 21.
• CALCAGNO, P., CHILÈS, J. P., COURRIOUX, G., GUILLEN, A. Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potential-field interpolation and geological rules. Physics of the Earth and Planetary Interiors. v. 171, n. 1-4, p. 147-157., 2008.
• CHILÈS, J. P., AUG, C., GUILLEN, A., LEES, T. Modelling the geometry of geological units and its uncertainty in 3D from structural data: the potential-field method. In: INTERNATIONAL SYMPOSIUM ON OREBODY MODELLING AND STRATEGIC MINE PLANNING. Proceedings... Perth, Australia, p. 313-320, 2004.
• COWAN, E. J., BEATSON, R. K., FRIGHT, W. R., MCLENNAN, T. J., MITCHELL, T. J. Rapid geological modelling. Applied Structural Geology for Mineral Exploration and Mining. In: INTERNATIONAL SYMPOSIUM, p. 23-25. 2002.
• COWAN, E. J., BEATSON, R. K., ROSS, H. J., FRIGHT, W. R., MCLENNAN, T. J., EVANS, T. R., CARR, J. C., LANE, R. G., BRIGHT, D. V., GILLMAN, A. J., OSHUST, P. A., TITLEY, M. Practical implicit geological modeling. In: INTERNATIONAL MINING GEOLOGY CONFERENCE. Proceedings, n. 8, p. 89-99, 2003. (AusIMM Publication Series).
• DEUTSCH, C. V., JOURNEL, A. G. GSLIB: geostatistical library and user’s guide. Oxford University Press, 1992.
• GOOVAERTS, P. Geostatistics for natural resources evaluation. Oxford University Press, 1997.
• GUARDIANO, F. B., SRIVASTAVA, R. M. Multivariate geostatistics: beyond bivariate moments. Geostatistics Troia’92, p. 133-144, 1993.
• JOURNEL, A. G. Fundamentals of Geostatistics in five lessons. In: Short Course in Geology, American Geophysical Union, Wahington DC, v. 8, 1989. 40 p.
• KARPEKOV, T. B., DEUTSCH, C. V. A review of boundary modeling and one way of selectng the C parameter to Model Boundaries with uncertainty. CCG Annual Report 18, p. 11, 2016.
• MATHERON, G., BEUCHER, H., DE FOUQUET, C., GALLI, A., GUERILLOT, D., RAVENNE, C. Conditional simulation of the geometry of fluvio-deltaic reservoirs. In: SPE ANNUAL TECHNICAL CONFERENCE AND EXHIBITION. Society of Petroleum Engineers, 1987.
• MUNROE, M. J., DEUTSCH, C. V. Full calibration of C and b in the framework of vein-type deposit tonnage uncertainty. CCG Annual Report 10, 2008. 16p.
• RASERA, L. G. Geoestatística de múltiplos pontos aplicada à simulação de modelos geológicos em grids estratigráficos. Universidade federal do Rio Grande do Sul, 2014. (Masters Dissertation).
• SILVA, D. A. Enhanced geologic modeling with data-driven training images for improved resources and recoverable reserves. University of Alberta, 2015. p. 222. (PhD Thesis).
• SILVA, D. A., DEUTSCH, C. V. Modeling multiple rock types with distance functions: metodology and software. CCG Annual Report, n. 14. p. 8, 2012.
• SOUZA, R. R. D. Modelagem geológica implícita através de simulações de funções distância assinaladas. Universidade federal do Rio Grande do Sul, 2017. (Masters Dissertation).
• STREBELLE, S. Conditional simulation of complex geological structrures using multiplepoint statistic. Mathematical Geology, v. 34, p. 1-21, 2002.
• STREBELLE, S. B., JOURNEL, A. G. Reservoir modeling using multiple-point statistics. In: SPE ANNUAL TECHNICAL CONFERENCE AND EXHIBITION. Society of Petroleum Engineers, 2001.
• WILDE, B. J., DEUTSCH, C. V. Simulating Boundary Realizations. CCG Annual Report 13, p. 4, 2011.

Publication Dates

• Publication in this collection
  16 Sept 2019
• Date of issue
  Oct-Dec 2019

History

• Received
  21 Mar 2019
• Accepted
  10 May 2019