Isatis Case Studies Mining

ISATIS 2011
Mining Case Studies

Published, sold and distributed by GEOVARIANCES
49 bis Av. Franklin Roosevelt, BP 91, 77212 Avon Cedex, France
Web: http://www.geovariances.com
Isatis Release 2011, March 2011
Contributing authors:
Catherine Bleins
Matthieu Bourges
Jacques Deraisme
Franois Geffroy
Nicolas Jeanne
Ophlie Lemarchand
Sbastien Perseval
Jrme Poisson
Frdric Rambert
Didier Renard
Yves Touffait
Laurent Wagner
All Rights Reserved
1993-2011 GEOVARIANCES
No part of the material protected by this copyright notice may be reproduced or utilized in any form
or by any means including photocopying, recording or by any information storage and retrieval
system, without written permission from the copyright owner.
"... There is no probability in itself. There are only probabilistic models. The
only question that really matters, in each particular case, is whether this or
that probabilistic model, in relation to this or that real phenomenon, has or
has not an objective meaning..."
G. Matheron
Estimating and Choosing - An Essay on Probability in Practice
(Springer Berlin, 1989)
1
Table of Contents
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .5
1 About This Manual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .7
2 In Situ 3D Resource Estimation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .11
2.1 Workflow Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .12
2.2 Presentation of the Dataset & Pre-processing. . . . . . . . . . . . . . . . . .16
2.3 Variographic Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .36
2.4 Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .68
2.5 Global Estimation With Change of Support . . . . . . . . . . . . . . . . . . .78
2.6 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .88
2.7 Displaying the Results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .129
3 Non Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .145
3.1 Introduction and overview of the case study. . . . . . . . . . . . . . . . . . .146
3.2 Preparation of the case study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .148
3.3 Global estimation of the recoverable resources . . . . . . . . . . . . . . . .165
3.4 Local estimation of the recoverable resources . . . . . . . . . . . . . . . . .176
3.5 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .203
3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .222
2
5
Introduction
6
1 About This Manual
Note - The present document only contains case studies related to a specific field of application. The full Case
Studies Manual can be downloaded on Geovariances web site.
A set of case studies is developed in this manual. It is mainly designed:
for new users to get familiar with the software and gives some leading lines to carry a study through,
for all users to improve their geostatistical knowledge by following detailed geostatistical workflows.
Basically, each case study describes how to carry out some specific calculations in Isatis as precisely as possi-
ble. The data sets are located on your disk in a sub-directory, called Datasets, of the Isatis installation directory.
You may follow the work flow proposed in the manual (all the main parameters are described) and then com-
pare the results and figures given in the manual with the ones you get from your test.
Most case studies are dedicated to a given field (Mining, Oil & Gas, Environment, Methodology) and therefore
grouped together in appropriate sections. However, new users are advised to run a maximum of case studies,
whatever their field of application. Indeed, each case study describes different functions of the package which
are not necessarily exclusive to one application field but could be useful for other ones.
Several case studies, namely In Situ 3D Resources Estimation (Mining), Property Mapping (Oil & Gas) and
Pollution (Environment) almost cover entire classic geostatistical workflows: exploratory data analysis, data
selections and variography, monovariate or multivariate estimation, simulations.
The other Case Studies are more specific and mainly deal with particular Isatis facilities, as described below:
Non Linear: anamorphosis (with and without information effect), indicator kriging, disjunctive kriging,
uniform conditioning, service variables and simulations.
Non Stationary & Volumetrics: non stationary modeling, external drift kriging and simulations, volume-
tric calculations, spill point calculation, variable editor.
Plurigaussian: an innovative facies simulation technique.
Oil Shale: fault editor.
Isatoil: multi-layer depth conversion with the Isatoil advanced module.
8 Case Studies
Young Fish Survey, Acoustic Fish Survey: polygons editor, global estimation.
Image Filtering: image filtering, grid or line smoothing, grid operator.
Boolean: boolean conditional simulations.
Note - All case studies are not necessarily updated for each Isatis release. Therefore, the last
update and the corresponding Isatis version are systematically given in the introduction.
In Situ 3D Resource Estimation 11
2 In Situ 3D Resource Esti-
mation
This case study is based on a real 3D data set kindly provided by Vale
(Carajs mine, Brazil).

It demonstrates particular features related to the Mining industry:
domaining, processing of three dimensional data, variogram modeling
and kriging. A brief description of global estimation with change of
support and block simulations is also provided. A simple application of
use of local parameters in kriging and simulations is presented.

Reminder: while using Isatis, the on-line help is accessible anytime by
pressing F1 and provides full description of the active application.

Last update: Isatis version 2012
12
2.1 Workflow Overview
This case study aims to give a detailed description of the kriging workflow and a brief introduction
to the grade simulation workflow of iron grades in an iron productive mine. This workflow over-
view lists the sequence of Isatis applications as they are ordered in the case study in order to run
through it. The list is nearly complete but not exhaustive.
Next to each application, two links are provided:
m the first link opens the application description of the Users guide: this allows the user to
have a complete description of the application as it is implemented in the software;
m the second link sends the user to the corresponding practical application example in the case
study.
Applications in bold are the most important for achieving kriging and simulation:
l File/Import Users Guide Case Study
Import the raw drillhole data.
l File/Selection/Macro Users Guide Case Study
Creates a macro-selection variable for each assay of the raw data based on the lithological code.
It is used to define two domains rich ore and poor ore.
l File/Selection/Geographic Users Guide Case Study
Creates a geographic selection to mask 4 drillholes outside of the orebody.
l Tools/Copy Variable/Header to Line Users Guide Case Study
Copy the selection masking the drillholes header to all assays of the drillholes.
l Tools/Regularization Users Guide Case Study
Assays compositing tool. A comparison of regularization by length or by domains is made. This
step is compulsory to make data additive for kriging. The composites regularized by domains
are kept for the rest of the study.
l Statistics / Quick Statistics Users Guide Case Study
Different modes for making statistics are illustrated: numerical statistics by domain, graphic dis-
plays with boxplots or swathplots.
l Statistics/Exploratory Data Analysis Users Guide Case Study
Isatis fundamental tool for QA/QC, 2D data displays, statistical and variographic analysis.
l Statistics/Variogram Fitting Users guide Case Study
Isatis tool for variogram modeling. Different modes are illustrated:
m manual: the user chooses by himself the basic structures (with their types, anisotropy, ranges
and sills) entering the parameters at the keyboard or for ranges/sills interactively in the Fit-
ting Window. This is used for modeling the variogramof the indicator of rich ore,
m automatic: the model is entireley defined (ranges, anisotropy and sills) from the definition of
the types and number of nested structures the user wants to fit. This is used for modeling the
Fe grade of rich ore.
l Statistics/Domaining/Border Effect Users Guide Case Study
Calculates statistical quantities based on domains indicator and grades to visualize the behav-
iour of grades when getting closer to the transition between domains.
l Statistics/Domaining/Contact Analysis Users Guide Case Study
Represents graphically the behaviour of the mean grade as a function of the distance of samples
to the contact between two domains.
l Interpolate/Estimation/(Co-)Kriging Users Guide Case Study
Isatis kriging application. It is applied here to krige (1) the indicator of rich ore and (2) the Fe
grade of rich ore on blocks 75mx75mx15m. In order to take into account the geo-morphology of
the deposit, kriging with Local Parameters is achieved: the main axis of anisotropy and neigh-
borhood ellipsod are changed between the northern and southern part of the deposit.
l Statistics/Gaussian Anamorphosis Modeling Users Guide Case Study
Isatis tool for normal score transform and modeling of histogram on composites support. This
step is compulsory for any non linear application including simulations. It is applied here on Fe
in the rich ore domain.
l Statistics/Support Correction Users Guide Case Study
Isatis tool for modeling grade histograms on block support. Useful for global estimation and for
non linear techniques (see Non Linear case study).
l Tools/Grade Tonnage Curves Users Guide Case Study
Calculates and represent graphically the grade tonnage curves. From the different possible
modes we compare the kriged panels and the distribution of grades on blocks obtained after sup-
port correction.
l File/Create Grid File Users Guide Case Study
Creates a grid of blocks 25mx25mx15m, on which we will simulate the ore type (1 for rich ore,
2 for poor ore) and the grades of Fe-P-SiO
2
.
l Tools/Migrate Grid to Point Users Guide Case Study
Transfers the selection variable defining the orebody from the panels 75mx75mx15m to the
blocks 25mx25mx15m.
14
l Interpolate/Conditional Simulations/Sequential Indicator/Standard Neighborhood Users
Guide Case Study
Simulations of the indicator of rich ore by SIS method.
l Statistics/Gaussian Anamorphosis Modeling Users Guide Case Study
That application is run again, for the purpose of a multivariate grade simulation, to transform
Fe-P-SiO
2
grades of composites. The P grade distribution is modelled differently from Fe and
SiO
2
, because of the presence of many values at the detection limit. The zero-effect distribution
type is then applied. It results that the gaussian value assigned to P has a truncated gaussian
distribution.
The Exploratory Data Analysis is used for calculating the experimental variogram on the gauss-
ian transform of P.
The variogram fitting is used with the Truncation Special Option for modeling the gaussian
experimental variogram of the gaussian transform of P.
l Statistics/Statistics/Gibbs Sampler Users guide Case Study
The Gibbs Sampler algorithm is used to generate the final gaussian transforms of P with a true
Gaussian distribution instead of a truncated one.
The Exploratory Data Analysis is used now for calculating the experimental variogram on the
gaussian transform of Fe-P-SiO
2
.
The variogram fitting is used for modeling the threevariate gaussian experimental variograms of
the gaussian transform of Fe-P-SiO
2
. The Automatic Sill Fitting mode is used: the sills of all
basic structures are automatically calculated using a least square minimization procedure.
l Statistics/Modeling/Variogram Regularization Users guide Case Study
The threevariate variogram model of the gaussian grades is regularized on the block support. A
new experimental variogram is then obtained.
The variogram fitting is used for modeling the threevariate gaussian experimental variograms of
the gaussian transform of Fe-P-SiO
2
on the block support (25mx25mx15m). The Automatic Sill
Fitting mode is used.
l Statistics/Modeling/Gaussian Support Correction Users guide Case Study
Transforms the point anamorphosis and the variogram model referring to the gaussian variables
regularized on the block support. The result is a gaussian anamorphosis on a block support and a
variogram model referring to the block gaussian variables (0-mean, variance 1). These steps are
compulsory for carrying out Direct Block Simulations.
l Interpolate/Conditional Simulations/Direct Block Simulations Users Guide Case Study
Simulations using the Turning Bands technique in the discrete gaussian model framework
(DGM).
l Statistics/Variogram on Grid Users Guide Case Study
Calculates, for QC purpose, the experimental variograms on the simulated gaussian block val-
ues.
l Statistics/Data Transformation/Raw<->Gaussian Transformation Users guide Case Study
Transforms the block gaussian simulations into raw block values.
l Tools/Copy Statistics/ Grid-> Grid Users Guide Case Study
Calculates rich ore tonnage and metal quantities in the panels 75mx75mx15m from the simu-
lated blocks 25mx25mx15m.
l File/Calculator Users Guide Case Study
Transforms the previous results into real ore tonnages and metals.
l Tools/Simulation Post-Processing Users Guide Case Study
Presents examples of Post-Processing of simulations.
l 3D viewer Users Guide Case Study
Some brief description of the 3D viewer module.
16
2.2 Presentation of the Dataset & Pre-processing
The data set is located in the Isatis installation directory (sub-directory Datasets/Mining) and con-
stituted of two different ASCII files:
l borehole measurements are stored in the ASCII file called boreholes.asc;
l a simple 3D geological model resulting from previous geological work (block size: 75 m hori-
zontally and 15 m vertically) is provided in a 3D grid file called block model_75x75x15m.asc.).
Firstly, a new study has to be created using the File / Data File Manager facility; then, it is advised
to verify the consistency of the units defined in the Preferences / Study Environment / Units win-
dow. In particular, it is suggested to use:
l Input Output Length Options:
Default Unit... = Length (m) Default Format...= Decimal (10,2)
l Graphical Axis Units:
X Coordinate = Length (km)
Y Coordinate = Length (km)
Z Coordinate = Length (m)
2.2.1 Borehole data
2.2.1.1 Data import
The boreholes.asc file begins with a header (commented by #) which describes its contents:
#
# structure=line , x_unit=m , y_unit=m , z_unit=m
#
# header_field=1 , type=alpha , name="drillhole ID"
# header_field=2 , type=xb , f_type=Decimal , f_length=8 , f_digits=2 , unit="m"
# header_field=3 , type=yb , f_type=Decimal , f_length=8 , f_digits=2 , unit="m"
# header_field=4 , type=zb , f_type=Decimal , f_length=8 , f_digits=2 , unit="m"
# header_field=5 , type=numeric , name="depth" , ffff=" " , bitlength=32 ;
# f_type=Decimal , f_length=8 , f_digits=2 , unit="m"
# header_field=6 , type=numeric , name="inclination" , ffff=" " ,
bitlength=32 ;
# f_type=Decimal , f_length=8 , f_digits=2 , unit="deg"
# header_field=7 , type=numeric , name="azimuth" , ffff=" " , bitlength=32
;
# f_type=Decimal , f_length=8 , f_digits=2 , unit="deg"
#
# field=1 , type=xe , f_type=Decimal , f_length=8 , f_digits=2 , unit="m"
# field=2 , type=ye , f_type=Decimal , f_length=8 , f_digits=2 , unit="m"
# field=3 , type=ze , f_type=Decimal , f_length=8 , f_digits=2 , unit="m"
#
# field=4 , type=numeric , name="Sample length" , ffff=" " , bitlength=32
;
# f_type=Decimal , f_length=6 , f_digits=2 , unit="m"
# field=5 , type=numeric , name="Fe" , ffff=" " , bitlength=32 ;
# f_type=Decimal , f_length=6 , f_digits=2 , unit="%"
# field=6 , type=numeric , name="P" , ffff=" " , bitlength=32 ;
# field=7 , type=numeric , name="SiO2" , ffff=" " , bitlength=32 ;
# field=8 , type=numeric , name="Al2O3" , ffff=" " , bitlength=32 ;
# field=9 , type=numeric , name="Mn" , ffff=" " , bitlength=32 ;
# field=10 , type=alpha , name="Lithological code ALPHA" , ffff=" "
# field=11 , type=numeric , name="Lithological code INTEGER" , ffff=" "
, bitlength= 8 ;
# f_type=Integer , f_length= 4 , unit=" "
#
# ++++ --------- +++++++++ --------- +++++++++ --------- +++++++++
# ++++++++++ --------- +++++++++ --------- +++++++++ --------- +++++++++ ---
------ +++++++++ --------- ---------
*---- 1 026 1400.00 -195.00 804.21 144.46 90.00 0.00
1 1400.00 -195.00 799.71 4.50 65.90 0.13 0.20
0.90 0.07 6 6
2 1400.00 -195.00 795.32 4.39 66.70 0.12 0.10
0.90 0.08 6 6
3 1400.00 -195.00 791.22 4.10 67.70 0.11 0.20
0.50 0.08 3 3
The samples are organized along lines and the file contains two types of records:
l The header record (for collars), which starts with an asterisk in the first column and introduces a
new line (i.e borehole).
l The regular record which describes one core of a borehole.
The file contains two delimiter lines which define the offsets for both records.
The dataset is read using the File / Import / ASCII procedure and stored in two new files of a new
directory called Mining Case Study:
18
l The file Drillholes Header, which contains the header of each borehole, stored as isolated
points.
l The file Drillholes, which contains the cores measured along the boreholes.
(snap. 2.2-1)
You can check in File / Data File Manager (by pressing s for statistics on the Drillholes file) that
the data set contains 188 boreholes, representing a total of 5954 samples. There are five numeric
variables (heterotopic dataset), whose statistics are given in the next table (using Statistics/Quick
Statistics...):
We will focus mainly on Fe variable. Also note the presence of an alphanumeric variable called
Lithological code Alpha.
Number Minimum Maximum Mean St. Dev.
Al
2
O
3
3591 0.07 44.70 1.77 4.14
Fe 5069 4.80 69.40 60.51 14.19
Mn 5008 0. 30.70 0.58 1.75
P 5069 0. 1. 0.06 0.08
Si O
2
3594 0.05 75.50 1.54 4.32
2.2.1.2 Borehole data visualization without the 3D viewer
Note - To visualize boreholes with the Isatis 3D viewer module, see the dedicated paragraph at the
end of this case study.
All the 2D Display facilities are explained in detail in the Displaying & Editing Graphics chapter
of the Beginner's Guide.
To visualize the lines without the 3D viewer, perform the following steps:
l Click on Display / New Page,
l In the Contents, for the Representation Type, choose Perspective,
l Double-click on Lines. An Item Contents for: Lines window appears:
m In the Data area, select the file Mining Case Study/Drillholes, without selecting any vari-
able as we are looking for a display of the boreholes geometry.
m Click on Display, and OK. The Lines appear in the graphic window.
l To change the View Point, click on the Camera tab and choose for instance:
m Longitude = -46
m Latitude = 20.
l Using the Display Box tab, deselect the toggle Automatic Scales and stretch the vertical dimen-
sion Z by a factor of 3.
l Click on Display.
l You should obtain the following display. You can save this template to automatically reproduce
it later: just click on Application / Store Page as in the graphic window.
(fig. 2.2-1)
20
The data set is contained in the following portion of the space:
Most of the boreholes are vertical and horizontally spaced approximately every 150m. The vertical
dimension is oriented upwards.
2.2.1.3 Creation of domains
In order to demonstrate Isatis capabilities linked to domaining, a simplified approach is presented
here. It consists in splitting the assays into two categories:
m the first one called rich ore corresponds to the lithological codes 1, 3 and 6,
m the second one called poor ore corresponds to the lithological codes 10 and above
A macro-selection final lithology[xxxxx] is created using File / Selection/Macro ...
After asking to create a New Macro Selection Variable and defining its name final lithology in the
Data File, you have to click on New.
(snap. 2.2-2)
Minimum Maximum
X 0.009 km 3.97 km
Y -0.35 km 3.77 km
Z -54.9 m +811.8 m
For creating Rich ore, Poor ore and Undefinedindices, you should give the name you want
(this has to be repeated three times). Then in the bottom part of the window you will define the
rules to apply. For each rule, you will have then to choose which variable it depends to, here Litho-
logical Code Integer, and the criterion to apply among the list you get by clicking on the button
proposing Equals as default:
m in the case of Rich ore you choose Is Lower or Equals to 9
m in the case of Poor ore you choose to match 2 rules (see snap shot on the previous page).
m in the case of Undefined you choose to match any of two rules (see next snap shot).
(snap. 2.2-3)
2.2.1.4 Drillholes selection
From the display of the drillholes, we can see that 4 are outside of the area covered by the other
drillholes. We will mask these drillholes for the rest of the study by using the File / Selection / Geo-
graphic menu.
The procedure "File / Selection / Geographic" is used to visualize and to perform a masking opera-
tion based on complete boreholes or more selectively on composites within a borehole.
We create the selection mask drillholes outside in the Drillholes header file.
22
(snap. 2.2-4)
When pressing the "Display as Points" button, the following graphic window opens representing by
a + symbol in green (according to the menu Preferences / Miscellaneous). the headers of all the
boreholes in a 2D XOY projection.
(snap. 2.2-5)
By picking with the mouse left button the 4 boreholes, their symbols are blinking, they can then be
masked by using the menu button of the mouse and clicking on Mask, the 4 masked boreholes are
then represented with the red square (according to the menu Preferences / Miscellaneous).
In the Geographic Selection window the number of selected samples (i.e.boreholes) is appearing
(184 from 188). To store the selection you must click on Run.
0 1000 2000 3000 4000
X (m)
0
1000
2000
3000
4000
Y

(
m
)
24
(snap. 2.2-6)
This selection is defined on the drillhole collars. In order to apply this selection to all samples of the
drillholes, a possible solution is to use the menu Tools / Copy Variable / Header Point -> Line.
(snap. 2.2-7)
2.2.1.5 Borehole data compositing
The compositing (or regularization) is an essential phase of a study using 3D data, especially in the
mining industry, although the principle is much more general. The idea is that geostatistics will
consider each datum with the same importance (prior to assigning a weight in the kriging process
0 1000 2000 3000 4000
X (m)
0
1000
2000
3000
4000
Y

(
m
)
for example) as it does not make sense to combine data that does not represent the same amount of
material.
Therefore, if data is measured on different support sizes, a first, essential task is to convert the
information into composites of the same dimension. This dimension is usually a multiple of the size
of the smallest sample, and is related to the height of the benches, which is in this case 15m.
l This operation can be achieved in different ways:
m the boreholes are cut into intervals of same length from the borehole collar, or in intervals
intersecting the boreholes and a regular system of horizontal benches. It is performed with
the Tools / Regularization by Benches or by Length facility, consists in creating a replica of
the initial data set where all the variables of interest in the input file are converted into com-
posites.
m the boreholes are cut into intervals of same length, determined on the basis of domain defini-
tion. Each time the domain assigned to the assay is changed a new composite is created. The
advantage of that method is to get more homogeneous composites. It is performed with the
Tools / Regularization by Domains facility.
m We will work on the 5 numerical variables Al
2
0
3
, Fe, Mn, P and SiO
2
.
m The regularization by length is performed on 5 numerical variables Al
2
0
3
, Fe, Mn, P and
SiO
2
and on the lithological code, in order to keep for each composite the information on the
most abundant lithology and the corresponding proportion. The new files are called:
- Composites 15m by length header for the header information (collars).
- Composites 15m by length for the composite information.
m Regularization mode: By Length measured along the borehole: this is the selected option as
some boreholes are inclined, with a constant length of 15m.
m Minimum Length: 7.5 m. It may happen that the first composite, or the last composite (or
both) do not have the requested dimension. Keeping too many of those incomplete samples
will lead us back to the initial problem of having samples of different dimensions being con-
sidered with the same importance: this is why the minimum length is set to 7.5 m (i.e. half of
the composite size).
26
(snap. 2.2-8)
m Three boreholes are not reproduced in the composite file as their total length is too small
(less than 7.5m): boreholes 93, 163 and 171. There are 1282 composites in the new output
file.
l The regularization by domain will calculate composites for two domains rich ore and poor
ore. The macro selection defining the domains in the input file is created with the same indices
in the output composites file. The selection mask drillholes outside is activated to regularize
only the boreholes within the orebody envelope. Only Fe, P, SiO2 are regularized. The new files
are called:
m Composites 15m header for the header information (collars).
m Composites 15m for the composite information.
m The Undefined Domain is assigned to the Undefined index. It means that when a sample is
in the Undefined Domain the composition procedure keeps on going (see on-line Help for
more information).
m The Analysed Length is kept for each grade element.
m The option Merge Residual is chosen, which means that the last composite is merged with
the previous one if its length is less than 50% of the composite length.
(snap. 2.2-9)
There are 1556 composites on the 184 boreholes in the new output file. From now on all geostatis-
tical processes will be applied on that regularized by domains composites file.
Using Statistics / Quick Statistics we can obtain different types of statistics, as for example:
The statistics on the Fe grades by domains. You note that after compositing there are no more
Undefined composites.
28
(snap. 2.2-10)
(snap. 2.2-11)
l Graphic representations with Boxplots by slicing according the main axes of the space.
(snap. 2.2-12)
30
32
l Swathplots by slicing according the main axes of the space.
(snap. 2.2-13)
(snap. 2.2-14)
The swathplots along OY shows for Fe rich ore a trend to decrease from South to North.
2.2.2 Block model
2.2.2.6 Grid import
The block model_75x75x15m.asc file begins with a header (Isatis format, commented by #) which
describes its contents:
#
# structure=grid, x_unit="m", y_unit="m", z_unit="m";
# sorting=+Z +Y +X ;
# x0= 150.00 , y0= -450.00 , z0= 310.00 ;
# dx= 75.00 , dy= 75.00 , dz= 15.00 ;
# nx= 28 , ny= 47 , nz= 31 ;
# theta= 0 , phi= 0 , psi= 0
# field=1, type=numeric, name="geographic domain", bitlength=32;
# ffff="N/A", unit="";
# f_type=Integer, f_length=9, f_digits=0;
# description="Creation Date: Mar 21 2006 15:13:15"
#
#+++++++++
0
0
0

The file contains only one numeric variable named geographic domain which equals 0, 1 or 2:
l 0 means the grid node lies outside the orebody,
l 1 means the grid node lies in the southern part of the orebody,
l 2 means the grid node lies in the northern part of the orebody.
Launch File/Import/ASCII... to import the grid in the Mining Case Study directory and call it 3D
Grid 75x75x15 m.
You have now to create a selection variable, called orebody, for all blocks where the geographic
code is either 1 or 2, by using the menu File / Selection / Intervals.
34
(snap. 2.2-15)
2.2.2.7 Visualization without the 3D viewer
Note - To visualize with the Isatis 3D viewer module, see the dedicated paragraph at the end of this
case study.
Click on Display / New Page in the Isatis main window. In the Contents window:
l In the Contents list, double click on the Raster item. A new Item contents for: Raster window
appears, in order to let you specify which variable you want to display and with which color
scale:
m Grid File...: select orebody variable from the 3D Grid 75x75x15 m file,
m In the Grid Contents area, enter 16 for the rank of the section XOY to display.
m In the Graphic Parameters area below, the default color scale is Rainbow.
m In the Item contents for: Raster window, click on Display.
m Click on OK.
l Your final graphic window should be similar to the one displayed hereafter.
(fig. 2.2-2)
The orebody lies approximately north-South, with a curve towards the southwestern part. The
northern part thins out along the northern direction and has a dipping plane striking North with a
western dip of 15 approximately. This particular geometry will be taken into account during vario-
graphic analysis.
500 1000 1500 2000
X (m)
0
1000
2000
3000
Y

(
m
)
36
2.3 Variographic Analysis
This step describes the structural analysis performed on 3D data set. In a first stage we consider the
Fe grade only of the rich ore (univariate analysis) on the 15 m composites. The estimation requires
to estimate for each block the proportion of rich ore and its grade. The analysis has then to be made:
l on the indicator of rich ore variable, which is defined on all composites
l and on the rich ore Fe grade, which is defined on rich ore composites.
The Exploratory Data Analysis (EDA) will be used in order to perform Quality Control, check sta-
tistical characteristics and establish the experimental variograms. Then variogram models will be
fitted.
2.3.1 Variographic analysis of rich ore indicator
The workflow that has been applied illustrates some important capabilities of Exploratory Data
Analysis, the decisions that are taken would probably require more detailed analysis in a real study.
The main steps of the workflow, that will be detailed in the next pages are:
l Calculation of the rich ore indicator.
l Variogram map in horizontal slices to confirm the existence of anisotropy.
l Calculations of directional variograms in horizontal plane. For simplification we keep 2 orthog-
onal directions East-West (N90) and North-South (N0).
l Check that the main directions of anisotropy are swapped when looking to northern or southern
boreholes.
l Save the Indicator variogram in the northern part (where are most of the data), with the idea
that the variogram in the Southern part is the same as in the North by inverting N0 and N90
directions of the anisotropy. In practice this will be realized at the kriging/simulation stage by
the use of Local Parameters for the variogram structures.
l Variogram Fitting using a combination of Automatic and Manual mode.
2.3.1.1 Calculation of the indicator
Use File / Calculator to assign the macro-selection index corresponding to rich ore to a float vari-
able Indicator rich ore.
(snap. 2.3-1)

2.3.1.2 Experimental Variogram of the Indicator
Launch Statistics/Exploratory Data Analysis... to start the analysis on the variable Indicator rich
ore:
38
(snap. 2.3-2)
Highlight the Indicator rich ore variable in the main EDA window and open the Base Map and His-
togram:
(fig. 2.3-1)
The mean value gives the proportion of rich ore samples.
The variogram map allows to check potential anisotropy. After clicking on the variogram map, the
Define Parameters Before Initial Calculations being on, you should choose the parameters as
shown in the next figure. You define parameters for horizontal slices, i.e. Ref.Plane UV with No
rotation.
Switch off the button Define the Calculations in the UW Plane and in the VW Plane, using the cor-
responding tabs.
With 18 directions each direction makes an angle of 10 with the previoius one. By asking a Toler-
ance on Directions of 2 sectors, the variograms are calculated from pairs in a given direction +/-
25.
0 1000 2000
X (m)
0
1000
2000
3000
4000
Y

(
m
)
0.0 0.5 1.0
Indicator rich ore
0.0
0.1
0.2
0.3
0.4
0.5
0.6
F
r
e
q
u
e
n
c
i
e
s
Nb Samples: 1556
Minimum: 0.000
Maximum: 1.000
Mean: 0.627
Std. Dev.: 0.484
40
(snap. 2.3-3)
(snap. 2.3-4)
After pressing OK you get the representation of the Variogram Map. In the Application Menu ask
Invert View Order to have variogram map and extracted experimental variograms in a landscape
view.
In the Application Menu ask Graphic Specific Parameters and change the Color Scale to Rain-
bow Reversed.
In the variogram map representation drag with the mouse a zone containing all directions. With the
menu button ask Activate Direction. You will then visualize the experimental variograms in the 18
directions of the horizontal plane. It exhibits clearly anisotropic behaviour.
42
(snap. 2.3-5)
We will now calculate the experimental variograms directly from the main EDA window by click-
ing on the Variogram bitmap at the bottom of the window. In the next figure we can see the param-
eters used for the calculation of 4 directional variograms in the horizontal plane and the vertical
variogram.
(snap. 2.3-6)
(snap. 2.3-7)
(snap. 2.3-8)
For sake of simplicity we decide to keep only 2 directions N0, showing more continuity and the
perpendicular direction N90.
The procedure to follow is:
44
l In the List of Options, change from Omnidirectional to Directional.
l In Regular Direction choose Number of Regular Directions 2 and switch on Activate Direction
Normal to the Reference Plane. Click Ok and go back to the Variogram Calculation Parameters
window.
(snap. 2.3-9)
You have then to define the parameters for each direction. Click the parameter table to edit:
l You have then to define the parameters for each direction. Click the parameter table to edit. For
applying the same parameters on the 2 horizontal directions, you must highlight these directions
in the Directions list of the Directions Definition window.
l The two regular directions choose the following parameters:
m Label for direction 1: N90 (default name)
m Label for direction 2: N0
m Tolerance on direction: 45 (in order to consider all samples without overlapping)
m Lag value: 90 m (i.e. approximately the distance between boreholes)
m Number of lags: 15(so that the variogram will be calculated over 1350 m distance)
m Tolerance on Distance (proportion of the lag): 0.5
m Slicing Height: 7.55 m (adapted to the height of composites)
m Number of Lags Refined: 1
m Lag Subdivision: 45m (so that we can have the variogram at short distance from the drill-
holes closely spaced).
l The normal direction with the following parameters:
m Label for direction 1: Vertical
m Tolerance on angle: 22.5
m Lag value: 15 m
m Number of lags: 10
m Tolerance on lags (proportion of the lag): 0.5
l In the Application Menu ask for Graphic Specific Parameters and click on the toggle button
for the display of the Histogram of Pairs.
(snap. 2.3-10)
Because the general shape of the orebody is anisotropic, we will calculate the variogram restricted
to the northern part and to the southern part of the orebody.
To do so you will use capabilities of the linked windows of EDA, by masking samples in the Base
Map. Automatically the variograms will be recalculated with only the selected samples.
For instance in the Base Map you drag a box around data in the Southern part (as shown on the fig-
ure) and with the menu button of the mouse you ask Mask. You will then get the variogram calcu-
lated from the northern data.
46
(snap. 2.3-11)
In the next figure we compare the variograms calculated from the northern and the southern data.
The main directions of anisotropy are swapped between North and South.
(snap. 2.3-12)
48
(snap. 2.3-13)
We decide now to fit a variogram model on the northern variogram, which is calculated with the
most abundant data. Then we will apply the same variogram to the southern data by making the
main axes of anisotropy swapped. This will be realized by means of local parameters attached to the
variogram model and to the neighborhood.
In the graphic window containing the experimental variogram in the northern zone, click on Appli-
cation / Save in Parameter File and save the variogram under the name Indicator rich ore North.
2.3.1.3 Variogram Modeling of the Indicator rich ore
You must now define a Model which fits the experimental variogram calculated previously. In the
Statistics / Variogram Fitting application, define:
l the Parameter File containing the set of experimental variograms: Indicator rich ore North.
l the Parameter File in which you wish to save the resulting model: Indicator rich ore
Click on Show Advanced Parameters.

(snap. 2.3-14)
50
l Set the toggles Fitting Window and Global Window ON; the program displays automatically
one default spherical model. The Fitting window displays one direction at a time (you may
choose the direction to display through Application/Variable & Direction Selection...), and the
Global window displays every variable (if several) and direction in one graphic.
l To display each direction in separate views, click in the Global Window on Application /
Graphic Specific Parameters and choose the Manual mode. Choose for Nb of Columns 3,
then Add, in turn for each Current Column, in the Selection by picking in the View Contents
area the First Variable, the Second Variable and the Direction.
(snap. 2.3-15)
l when pressing the Edit button next to the variogram model, the Model Definition sub-window
opens and the user can choose the basic structures. The model must reflect:
m the variability at short distances, with a consistent nugget effect,
m the main directions of anisotropy,
m the general increase of the variogram.
The model is automatically defined with the same rotation definition as the experimental vario-
gram. Three different structures have been defined (in the Model Definition window, use the Add
button to add a structure, and define its characteristics below, for each structure):
l Nugget effect,
l Anisotropic Exponential model with the following respective ranges along U, V and W: 700 m,
550 m and 70 m,
l Anisotropic Exponential model with the following respective ranges along U, V and W: 500
m, 5000 m and nothing (which means that it is a zonal component with no contribution in the
vertical direction).
Do not specify the sill for each structure at this stage, instead:
l click Nugget effect in the main Variogram Fitting window, set the toggle button Lock the Nug-
get Effect Components During Automatic Sill Fitting ON and enter the value .065.
l set the toggle Automatic Sill Fitting ON. The program automatically computes the sills and dis-
plays the results in the graphic windows.
l A final adjustement is necessary, particularly to get a total sill of 0.25, which is the maximum
admissible for a stationary indicator variogram. Set the toggle Automatic Sill Fitting OFF from
the main Variogram Fitting window, then in the Model Definition window set the sill for the
first exponential to 0.14 and the sill for the second exponential to 0.045.
The final model is saved in the parameter file by clicking Run in the Variogram Fitting window.
(snap. 2.3-16)
2.3.2 Variographic Analysis of Fe rich ore
2.3.2.4 Experimental Variogram of Fe rich ore
Launch Statistics/Exploratory Data Analysis... to start the analysis on the variable Fe using the
selection for the rich ore composites.
52
(snap. 2.3-17)
You will calculate the variograms in 2 directions of dipping plane striking North with a western dip
of 15. In the Calculation Parameters you will ask in List of Options a Directional. Click then Reg-
ular Directions a new window Directions pops up where you will define the Reference Direction
and switch on Activate Direction Normal to the Reference Plane.
(snap. 2.3-18)
Click Reference Direction, in 3D Direction Definition window set the convention to User Defined
and define the rotation parameters as shown in the next figure.
(snap. 2.3-19)
The reference direction U (in red) correspond to the N121 main direction of anisotropy.
The calculation parameters are then chosen as shown in the next figure.
54
(snap. 2.3-20)
The next figure shows the experimental variograms.
Two points may be noted:
l the anisotropy is not really marked, we will recalculate isotropic variogram in the horizontal
plane,
l the second point of the variogram for the direction N121, calculated with 42 pairs, shows a peak
that we can explain by using the Exploratory Data Analysis linked windows.
(snap. 2.3-21)
For using the linked windows the following actions have to be made:
56
l ask to display the histogram (accept the default parameters),
l in the Graphic Specific Parameters of the graphic page containing the experimental variogram,
set the toggle button Variogram Cloud (if calculated) OFF, and click on the radio button Pick
from Experimental Variogram.
l in the Calculation Parameters of the graphic page containing the experimental variogram, set
the toggle button Calculate the Variogram Cloud ON.
l In the graphic page click on the experimental point with 43 pairs and ask in the menu of the
mouse Highlight. The variogram is then represented as a blue square, and all data making the
pairs represented the part painted in blue in the histogram.
(snap. 2.3-22)
The high variability due to pairs made of the samples with low values is responsible of the peak in
the variogram. It can be proved by clicking in the histogram on the bar of the minimum values and
clicking with the menu of the mouse on Mask, the variograms are automatically calculated and
dont show anymore the anomalous point as shown on the next figure.
(snap. 2.3-23)
l We now re-calculate the variograms with 2 directions, omni-directional in the horizontal plane
and vertical, with the parameters shown hereafter you enter by clicking Regular Directions....
(snap. 2.3-24)
58
(snap. 2.3-25)
In the graphic containing this last variogram ask for the Application->Save in Parameter File to
save the variogram with the name Fe rich ore.
2.3.2.5 Variogram Modeling of Fe rich ore
In the Statistics / Variogram Fitting application, define:
l the Parameter File containing the set of experimental variograms: Fe rich ore
l the Parameter File in which you wish to save the resulting model: Fe rich ore
Open the Model Intialization window in order to make an automatic model with a nugget effect and
2 spherical (short and long range )
N0
41
157
472
688
1120
1373
1195
1196
900
1108
1222
1155
D-90
6
78 392
325
266
223
183
148
117
0 500 1000 1500
Distance (m)
0
5
10
15
V
a
r
i
o
g
r
a
m

:

F
e
(snap. 2.3-26)
In the Global window, you represent the variograms in two columns, the automatic variogram
looks satisfactory, so you click Run in the Variogram Fitting window to save it.
(fig. 2.3-2)
2.3.3 Analysis of border effects
This chapter may be skipped in a first reading as it does not change anything in the Isatis study. It
helps to decide whether kriging/ simulation will be made using hard or soft boundary.
In order to understand the behaviour of Fe grades when the samples are close to the border between
rich and poor ore, we can use two applications:
60
l Statistics / Domaining / Border effect calculates bi-point statistics from pairs of samples belong-
ing to different domains. The pairs are chosen in the same way as for experimental variogram
calculations.
l Statistics / Domaining / Contact Analysis calculates the mean values of samples of 2 domains
as a function of the distance to the contact between these domains along the drillholes.
2.3.3.6 Statistics on Border effect
Launch Statistics / Domaining / Border effect and choose in the file Composites 15m, the Macro
Selection Variable final lithology[xxxxx], that contains the definition of all domains, and the vari-
able of interest Fe.
In the list of Domains you may pick only some of these, in this case Rich ore and Poor ore, while
you ask to Mask Samples from Domain choosing Undefined.
In the Calculation Parameters sub-window we define the parameters for 3 directions by pressing
the corresponding tabs in turn and switching on the toggle Activate Direction. For the 3 directions
the parameters are:
Switch on the three toggle buttons for the Graphic Parameters and click on Run.
(snap. 2.3-27)
Three graphic pages corresponding to the three statistics are then displayed:
62
l Transition Probability, that, in the case of only 2 domains, is not very informative.
(snap. 2.3-28)
l Mean [Z(x+h)|Z(x)], that shows that when going from Rich ore to Poor ore there is a border
effect (the grade of the new domain, i.e. Poor ore, is higher than the mean Poor ore grade which
means it is influenced at short distance by the proximity to Rich ore samples. Conversely when
going from Poor ore to Rich ore there is no border effect.
(snap. 2.3-29)
Dir
Dir
Dir
0 500 1000 1500
Distance (m)
0
10
20
30
40
50
60
70
F
e

e
n
t
e
r
i
n
g

i
n

R
i
c
h

o
r
e
Dir
Dir
Dir
0 500 1000 1500
Distance (m)
0
10
20
30
40
50
60
70
F
e

x
+
h

i
n

R
i
c
h

o
r
e

|

x

i
n

P
o
o
r

o
r
e
Dir
Dir
Dir
0 500 1000 1500
Distance (m)
0
10
20
30
40
50
60
70
F
e

e
n
t
e
r
i
n
g

i
n

P
o
o
r

o
r
e
Dir
Dir
Dir
0 500 1000 1500
Distance (m)
0
10
20
30
40
50
60
70
F
e

x
+
h

i
n

P
o
o
r

o
r
e

|

x

i
n

R
i
c
h

o
r
e
64
l Mean Diff[Z(x+h)-Z(x)], that shows that when going from Rich ore to Poor ore as well as
going from Poor ore to Rich ore the grade difference is influenced by the proximity of both
domains.
(snap. 2.3-30)
2.3.3.7 Contact Analysis
Launch Statistics / Domaining / Contact Analysis and choose in the file Composites 15m, the
Macro Selection Variable final lithology[xxxxx], that contains the definition of all domains, and
the variable of interest Fe. You set the variables Direct Distance Variable and Indirect Distance
Variable to None, which means that the contact point is determined when the domain changes
down the boreholes.
In the list of Domains you pick Rich ore for Domain 1 and Poor ore for Domain 2, while you let
Use Undefined Domain Variable to Off.
The statistics are calculated as a function of the distance to the contact along the drillhole, you have
the possibility to select only some of the drillholes according to a specific direction with an angular
tolerance. In this case, as most of the drillholes are vertical, we select all drillholes by choosing a
tolerance of 90 on the vertical direction defined by thre rotation angles Az=0, Ay=90, Ax=0 (Math-
ematician Convention). The samples are regrouped by Distance Classes of 15m.
Dir
Dir
Dir
0 500 1000 1500
Distance (m)
-40
-30
-20
-10
0
10
20
30
40
D
i
f
f

F
e
,

x
+
h

i
n

R
i
c
h

o
r
e
,

x

N
O
T
Dir
Dir
Dir
0 500 1000 1500
Distance (m)
-40
-30
-20
-10
0
10
20
30
40
D
i
f
f

F
e
,

x
+
h

i
n

R
i
c
h

o
r
e

|

x

i
n

P
o
o
r

o
r
Dir
Dir
Dir
0 500 1000 1500
Distance (m)
-40
-30
-20
-10
0
10
20
30
40
D
i
f
f

F
e
,

x
+
h

i
n

P
o
o
r

o
r
e
,

x

N
O
T
Dir
Dir
Dir
0 500 1000 1500
Distance (m)
-40
-30
-20
-10
0
10
20
30
40
D
i
f
f

F
e

x
+
h

i
n

P
o
o
r

o
r
e

|

x

i
n

R
i
c
h

o
r
e
(snap. 2.3-31)
Two graphic pages are then displayed:
l Contact Analysis (Oriented) contains two views:
m Direct for statistics calculated in the Reference Direction
m Indirect for statistics calculated in the opposite of the Reference Direction
In the Application menu of the graphic pages we ask the Graphical Parameters, as shown
below, to display the Number of Points and the Mean per Domain.
(snap. 2.3-32)
66
(snap. 2.3-33)
l Contact Analysis (Non-Oriented) displays the average of the two previous ones.
(snap. 2.3-34)
From these graphs it appears that the poor grades are influenced by the proximity to rich grades.
In conclusion we decide for the kriging and simulations steps to apply hard boundary when dealing
with rich ore.
68
2.4 Kriging
We are now going to estimate on blocks 75mx75mx15m the tonnage and Fe grades of Rich ore.
Therefore, we will perform two steps:
l Kriging of the Indicator of Rich ore to get the estimated proportion of rich ore, from which the
tonnage can be deduced.
l Kriging of the Fe grade of rich ore using only the rich ore samples. Each block is then estimated
as if it would be entirely in rich ore, by applying the estimated tonnage, we can then obtain an
estimate of the Fe metal content.
2.4.1 Kriging of indicator of rich ore with local parameters
After the variographic analysis it was found that the variogram model has an horizontal anisotropy
that has a different orientation in the northern and southern part of the orebody. We will then use
that orientation as local parameter recovered from the grid file in a variable called RotZ. As a first
attempt, that should be sufficient in this case because of the orebody shape, we will use two values
90 for blocks in the southern area and 0 for the northern area, both areas being defined by means
of the geographic code variable (respectively 1 and 2). These values are stored in the grid file by
using File / Calculator.
(snap. 2.4-1)
Then you launch Interpolate / Estimation / (Co)Kriging.
(snap. 2.4-2)
You need to specify the type of calculation to Block and the number of variables to 1, then:
l Input File: Indicator rich ore (Composites on 15m with the selection None).
l The names of the variables in the output file (3D Grid 75 x 75 x 15 m), with the orebody selec-
tion active:
m Kriging indicator rich ore for the estimation of Indicator rich ore
m Kriging indicator rich ore std dev for the kriging standard deviation
70
l The variogram model contained in the Parameter File called Indicator rich ore.
l The neighborhood: open the Neighborhood... definition window and specify the name (Indica-
tor rich ore for instance) of the new parameter file which will contain the following parameters,
to be defined from the Edit... button nearby. The neighborhood type is set by default to moving:
(snap. 2.4-3)
m The moving neighborhood is an ellipsoid with No rotation, which means that U,V,W axes
are the original X,Y,Z axes;
m Set the dimensions of the ellipsoid to 800 m, 600 m and 60 m along the vertical direction;
m Switch ON the Use Anisotropic Distances button.
m Minimum number of samples: 4;
m Number of angular sectors: 12
m Optimum Number of Samples per Sector: 5
m Block discretization: as we chose to perform Block kriging, the block discretization has to be
defined. The default settings for discretization are 5 x 5 x 1, meaning each block is sub-
divided by 5 in each X and Y direction, but is not divided in Z direction. The Block Discret-
ization sub-window may be used to change these settings, and check how different discreti-
zations influence the block covariance C
vv
. In this case study, the default parameters 5x5x1
will be kept.
m Press OK for the Neighborhood Definition.
l The Local Parameters: open the Local Parameters Loading... window and specify the name of
the Local Parameters File (3D Grid 75x75x15m). Fore the Model All Structures and Neighbor-
hood tabs switch ON Use Local Rotation (Mathematician convention) then 2D and define as
Rotation/Z the variable Rot Z.
(snap. 2.4-4)
72
It is possible to check both the model and the neighborhood performances when processing on a
grid node, and to display the results graphically: this is the purpose of the Test option at the bottom
of the (Co-)Kriging main window. When pressing it, a graphic page opens where:
l The Indicator rich ore variable is represented with proportional symbols,
l The neighborhood ellipsoid is drawn on a 2D section.
By pressing once on the left button of the mouse, the target grid is shown (in fact a XOY section of
it, you may select different sections through Application/Selection For Display...). The user can
then move the cursor to a target grid node: click once more to initiate kriging. The samples selected
in the neighborhood are highlighted and the weights are displayed. We can see here that the nearest
samples get the higher weights. It is also important to check that the negative weights due to screen
effect are not too important. The neighborhood can be changed sometimes to avoid this kind of
problem (more sectors and less points by sector...).
You can also select the target grid node by giving the indices along X, Y and Z with the Application
menu Target Selection (for instance 6, 11, 16). You can figure out how the local parameters used
for the neighborhood are applied.
(snap. 2.4-5)
(snap. 2.4-6)
Note - From Application/Link to 3D viewer, you may ask for a 3D representation of the search
ellipsoid if the 3D viewer application is already running (see the end of this case study).
Close the Test Window and press RUN.
7814 grid nodes have been estimated. Basic statistics of the variables are displayed below.
(fig. 2.4-1)
The kriging standard deviation is an indicator of the estimation error, and depends only on the geo-
metrical configuration of the data around the target grid node and on the variogram model. Basi-
cally, the standard deviation decrease as an estimated grid node is closer to data.
Some blocks have the kriged indicator above 1. These values will be changed into 1 by means of
File / Calculator.
74
(snap. 2.4-7)
Note - In the main Kriging window, the optional toggle Full set of Output Variables allows to
store in the Output File other kriging parameters: slope of regression, weight of the mean,
estimated dispersion variance of estimates etc...
2.4.2 Kriging of Fe rich ore
In the Standard (Co)Kriging menu specify the type of calculation to Block and the number of
variables to 1, then enter the following parameters:
l Input File: Fe (Composites on 15m with the selection final lithology{rich ore}).
l The names of the variables in the output file (3D Grid 75 x 75 x 15 m), with the orebody selec-
tion active:
m Kriging Fe rich ore for the estimation of Fe;
m Kriging Fe rich ore std dev for the kriging standard deviation.
l The variogram model contained in the Parameter File called Fe rich ore.
l The neighborhood: open the Neighborhood... definition window and specify the name (Fe rich
ore for instance) of the new parameter file which will contain the following parameters, to be
defined from the Edit... button nearby. The neighborhood type is set by default to moving:
m The moving neighborhood is an ellipsoid with No rotation, which means that U,V,W axes
are the original X,Y,Z axes;
m Set the dimensions of the ellipsoid to 800 m, 300 m and 50 m along the vertical direction;
m Switch ON the Use Anisotropic Distances button.
m Minimum number of samples: 4;
m Optimum Number of Samples per Sector: 3
m Block discretization: as we chose to perform Block kriging, the block discretization is kept to
the default 5 x 5 x 1.
l Apply Local Parameters but only for the Neighborhood, where you use Rot Z variable for 2D
Rotation /Z.
(snap. 2.4-8)
76
After Run you can calculate the statistics of the kriged estimate by asking in Statistics / Quick Sta-
tistics to apply as Weight the weight variable Kriging indicator rich ore. 7561 blocks from 7814
have been kriged. By using a weight variable you will obtain the statistics weighted by the propor-
tion of the block in rich ore.
(snap. 2.4-9)
(fig. 2.4-2)
The mean grade is close to the average of the composites grade (65.84). Therefore in the next steps,
carrying out non linear methods which require the modeling of the distribution, we will not apply
any declustering weights.
78
2.5 Global Estimation With Change of Support
The support is the geometrical volume on which the grade is defined.
Assuming the data sampling is representative of the deposit, it is possible to fit a histogram model
on the experimental histogram of the composites. But at the mining stage, the cut-off will be
applied on blocks, not on composites. Therefore, it is necessary to apply a support correction to the
composite histogram model in order to estimate an histogram model on the block support.
Note - When kriging too small blocks with a high error level, applying a cut-off to the kriged
grades will induce biased tonnage estimates due to the high smoothing effect. It is then
recommended to use non-linear estimation techniques, or simulations (see the Non Linear case
study). For global estimation, an other alternative is to use the Gaussian anamorphosis modeling,
as described here below.
2.5.1 Gaussian anamorphosis modeling
Gaussian anamorphosis is a mathematical technique which allows to model histograms, taking the
change of support from composites to blocks into account.
Note - From a support size point of view, composites will be considered as points compared to
blocks.
The technique will not be mathematically detailed here: the reader is referred to the Isatis on-line
help and technical references. Basically, the anamorphosis transforms an experimental dataset to a
gaussian dataset (i.e. having a gaussian histogram). The anamorphosis is bijective, so it is possible
to back transform gaussian values to raw values. A gaussian histogram is often a pre-requisite for
using non linear and simulation techniques. The anamorphosis function may be modelled in two
ways:
l by a discretization with n points between a negative gaussian value of -5 and a positive gaussian
value of +5.
l by using a decomposition into Hermite polynomials up to a degree N. This was the only possi-
bility until the Isatis release V10.0. It is still compulsory for some applications, as will be
explained later on.
Open the Statistics/Gaussian Anamorphosis Modeling window.
(snap. 2.5-1)
l In Input... choose the Composites 15 m file with the selection final lithology{Rich ore};
choose Fe for the raw variable.
l Do NOT ask for a Gaussian Transform.
l Name the anamorphosis function Fe rich ore.
l In Interactive Fitting... choose the Type Standard and switch ON the toggle button Dispersion
with the Dispersion Law set to Log-Normal Distribution. In this mode the histogram will be
modelled by assigning to each datum a dispersion, that accounts for some uncertainty that is
80
globally reflected by an error on the mean value. The variability of the dispersion is controlled
by the Variance Increase parameter, related to the estimation variance of the mean. By default
that variance is set to the statistical variance of the data divided by the number of data.
(snap. 2.5-2)
l Click on the Anamorphosis and Histogram bitmaps. You will visualize the anamorphosis func-
tion and how the experimental histogram is modelled (black bars are for the experimental histo-
gram and the blue bars for the modelled histogram).
(snap. 2.5-3)
Close the Fitting Parameters window.
l Press RUN in the Gaussian Anamorphosis window: because you have not asked for Hermite
Polynomials, the following error message window is displayed to advise you on the applications
requiring these polynomials.
(snap. 2.5-4)
82
2.5.2 Block anamorphosis on SMU support
Using the composite histogram and variogram models, we are now going to take the change of sup-
port into account using Statistics/Support Correction...:
(snap. 2.5-5)
The Selective Mining Unit (SMU) size has been fixed to 25 x 25 x 15 m. Therefore, the correction
will be calculated for a block support of 25 x 25 x 15 m. Each block is discretized by default in 3x3
for the X and Y direction (NX = 3 and NY = 3); no discretization is needed for the vertical direction
(NZ = 1) as the composites are regularized accordingly to the bench height (15 m). Changing the
discretization along X and Y may allow to study the sensitivity on change of support coefficients.
Switch ON the toggle button Normalize Variogram Sill. As the variogram sill is higher than the
variance, the consequence is to reduce a little bit the support correction (r coefficient a bit higher
than without normalization).
Press Calculate at the bottom of the window. The block support correction calculations are dis-
played in the message window:
(snap. 2.5-6)
The block variogram value Gamma (v,v) is calculated and is the base for calculating the real
block variance and the real block support correction coefficient r. We can see that the support
correction is not very important (r not very far from 1), it is because of the variogram model whose
ranges are rather large compared to the smu size. The calculation is made at random, so different
calculations will give similar results, but different. If the differences in the real block variance are
too large, the block discretization should be refined by increasing NX and NY. By pressing Calcu-
late... several times, we statistically check if the discretization is fine enough to represent the vari-
ability inside the blocks. Press OK.
Save the Block Anamorphosis under the name Fe rich ore block 25x25x15 and press RUN.
2.5.3 Grade Tonnage Curves
Launch Tools / Grade Tonnage Curves. You will ask to display two types of curves, calculated
from:
84
l Kriged Fe rich ore on the panels 75mx75mx15m, the Histogram modelled after support correc-
tion on blocks 25mx25mx15m
. (snap. 2.5-7)
For each curve you have to click Edit and Fill the parameters.
For the first curve on kriged panels:
(snap. 2.5-8)
For the second curve, on blocks histogram:
86
(snap. 2.5-9)
After clicking the bitmaps at the bottom of the Grade Tonnage Curves window (M vs. z, T vs z, Q
vs. z, Q vs.T, B vs z) you get the graphics like for instance T(z), M(z):
(snap. 2.5-10)
(snap. 2.5-11)
These curves show as expected that the selectivity is better from true blocks 25x25x15 than from
kriged panels 75x75x15, that have a lower dispersion variance.
The legend is displayed in a Separate Window as was asked in the Grade Tonange Curves win-
dow. By clicking Define Axes you switch OFF Automatic Bounds to change the Axis Minimum and
Axis Maximum for Mean Grade to 60 and 70 respectively.
50 55 60 65
Cutoff
0
10
20
30
40
50
60
70
80
90
100
T
o
t
a
l

T
o
n
n
a
g
e
50 55 60 65
Cutoff
60
61
62
63
64
65
66
67
68
69
70
M
e
a
n

G
r
a
d
e
88
(snap. 2.5-12)
(snap. 2.5-13)
2.6 Simulations
This chapter aims at giving a quick example of conditional block simulations in a multivariate case.
Simulations allow to reproduce the real variability of the variable.
We will focus on the Fe-P-SiO
2
grades of rich ore of blocks 25mx25mx15m. Two steps will then be
achieved:
l simulation of the rich ore indicator. Sequential Indicator method will be applied to generate sim-
ulated model where each block has a simulated code 1 for rich ore blocks and 2 for poor ore
blocks. A finer grid would be required to be more realistic, for sake of simplicity we will make
the indicator simulation on the same blocks 25mx25mx15m.
l simulation of rich ore Fe grade, as if each block would be entirely in rich ore. By intersecting
with the indicator simulation, we will get the final picture.
2.6.1 Simulation of the indicator rich ore
You must first create the grid of blocks 25x25x15 with File / Create Grid File.
(snap. 2.6-1)
90
To create in the grid file the orebody selection we use the migration capability (Tools/Migrate/Grid
to Point...) from the 3D Grid 75x75x15 m file to 3D Grid 25x25x15 with maximum migration dis-
tance of 55 m.
(snap. 2.6-2)
Open the menu Interpolate / Conditional Simulations / Sequential Indicator / Standard Neighbor-
hood.
(snap. 2.6-3)
For defining the two facies 1 for rich ore and 2 for the complementary you have to click on
Facies Definition and enter the parameters as shown below.
92
(snap. 2.6-4)
You may use the same variogram model, the same neighborhood and the same local parameters
as used for the kriging. The only additional parameter is the Optimum Number of Already Simulated
Nodes, you can fix to 30 (the total number being 5 for 12 sectors, i.e. 60). Save the simulation in
SIS indicator rich ore.
You ask 100 simulations, then press on Run.
2.6.2 Block simulations of Fe-P-SiO
2
rich ore
The direct block simulation method, based on the discrete gaussian model (DGM), will be used.
The workflow is the following:
m transform the raw data to gaussian values by anamorphosis. For the case of P grade the
anamorphosis will take into account the fact that many samples are at the detection limit,
that produces an histogram with a significant zero effect.
m do a multivariate variographic analysis on the gaussian data in order to have a gaussian vari-
ogram
m model these gaussian variograms with a linear model of coregionalisation;
m regularize these variograms on the block support;
m perform a support correction on the gaussian transforms;
m perform the simulations using the discrete gaussian model framework, that allows to condi-
tion block simulated values to gaussian point data.
2.6.2.1 Gaussian Anamorphosis
We will perform the gaussian anamorphosis on the three grades of the rich ore domain in one go.
and independently. Note that the three anamorphosis functions must be stored together in the same
Parameter file called Fe-SiO
2
-P rich ore. Note in this case that we also ask to store the Gaussian
transforms in the composites file with the names Gaussian Fe/P/SiO
2
rich ore, ...
94
(snap. 2.6-5)
By clicking on Interactive Fitting, the Fitting Parameters window pops up, you will have to choose
parameters for the three variables in turn, by clicking on the arrow on the side of the area displaying
Parameters for Fe/P/SiO
2
. For Fe and SiO
2
you choose the Standard Type with a Dispersion
using a Log Normal Distribution and the default Variance Increase (as was made before for Fe
alone).
For P many samples have values equal to the detection limit of 0.01. The histogram shows a spike
at the origin, that will be modelled by a zero-effect. You must choose the type Zero-effect and click
on Advanced Parameters to enter the parameters defining the zero effect. In particular we will put
in the atom all values equal to 0.01 with a precision of 0.01, i.e. all samples between 0 and 0.02.
(snap. 2.6-6)
After Run the transformed values of Fe and SiO2 have a gaussian distribution, while for P the gaus-
sian transform has a truncated gaussian distribution. The gaussian values assigned to the samples
concerned by the zero effect are all equal to the same value (gaussian value corresponding to the
frequency of the zero effect).
2.6.2.2 Gaussian transform of P rich ore
The next steps consist of making the gaussian transform of P a true gaussian distribution. This is
achieved by using a Gibbs Sampler algorithm that will generate for all samples of the zero effect a
gaussian value consistent with the structure of spatial correlation with all gaussian values. Practi-
cally 3 steps must be carried out:
l calculation of the experimental variogram of the truncated gaussian values;
l variogram modelling of the gaussian transform using the truncation option;
l Gibbs Sampler to generate the gaussian transform with a true distribution and honouring the
spatial correlation.
Using EDA we calculate the histogram and the experimental variogram on the variable Gaussian
P rich ore (activating the selection final lithology{Rich ore}). In the Application menu of the his-
togram you ask the Calculation Parameters and switch off the Automatic mode to the values shown
below:
(snap. 2.6-7)
96
For the variogram you choose the same parameters as used for Fe (omnidirectional in the horizontal
plane and vertical), by asking in the Application Menu / Calculation Parameters, in the Variogram
Calculation Parameters window click Load Parameters from Standard Parameter File and select
the experimental variogram Fe rich ore.
On the graphic display you see the truncated distribution with about 35% of samples concerned by
the zero effect, the gaussian truncated value is -0.393. The variance displayed as the dotted line on
the variograms is about 0.5. In the Application / Save in Parameter File menu of the graphic con-
taining the variogram you save it under the name Gaussian P rich ore zero effect.
(snap. 2.6-8)
(snap. 2.6-9)
In the Variogram Fitting window you choose the Experimental Variograms Gaussian P rich ore
zero effect and you create a New Variogram Model, called Gaussian P rich ore. Note that the var-
iogram model refers to the gaussian transform (with the true gaussian distribution), it is transformed
by means of the truncation to match the experimental variogram of the truncated gaussian variable.
(snap. 2.6-10)
Click Edit, in the Model Definition window you must first click Truncation.
98
(snap. 2.6-11)
In the Truncation window, switch ON Truncation and click Anamorphosis V1 to select the anamor-
phosis Fe-SiO
2
-P[P].
(snap. 2.6-12)
(snap. 2.6-13)
Coming back to the Model Definition window you enter the parameters of the variogram model as
shown below. It is important to choose sill coefficients summing up to 1 (dispersion variance of the
true gaussian) and not 0.5 the dispersion variance of the truncated gaussian.
(snap. 2.6-14)
N0
1
157
472
688
1120
1373
1195
1196 900
1108
1222
1155
0 500 1000 1500
Distance (m)
0.0
0.1
0.2
0.3
0.4
0.5
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

P

r
i
c
h

o
r
e
D-9
606
478
392
325
266
223
183
148
117
0 50 100 150
Distance (m)
0.0
0.1
0.2
0.3
0.4
0.5
0.6
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

P

r
i
c
h

o
r
e
100
You will now generate gaussian values for the zero effect on P rich ore by using Statistics / Statis-
tics / Gibbs Sampler. Note that the gaussian values not concerned by the zero effect are kept
unchanged.
l The Input Data are the variogram model you just fitted Gaussian P rich ore and the Gaussian
P rich ore variable stored after the GaussainAnamorphosis Modelling.
l The Output Data are a new variogram model Gaussian P rich ore no truncation (which is in
fact the same as the input one without the truncation option) and a new variable in the Compos-
ites 15m file Gaussian P rich ore (Gibbs).
l You ask to perform 1000 iterations.
(snap. 2.6-15)
You can check how the Gibbs Sampler has reproduced the gaussian distribution and the input vari-
ogram. You just have to recalculate the histogram and the variograms on the variable Gaussian P
rich ore (Gibbs). After saving in the Parameter File that experimental variogram, you can superim-
pose to it the variogram model with no truncation using Variogram Fitting menu. For the first dis-
tance the fit is acceptable.
(snap. 2.6-16)
(snap. 2.6-17)
N0
1
157
472688
1120
1373
1195
1196
900
1108
1222
1155
D-9
6
78
92
325
266
223
183
148
117
0 500 1000 1500
Distance (m)
0.0
0.5
1.0
1.5
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

P

r
i
c
h

o
r
e

(
G
i
b
b
s
)
102
2.6.2.3 Multivariate Gaussian variogram modeling
In Statistics / Exploratory Data Analysis you calculate the variograms with the same parameters as
before (one monidirectional horizontal direction and one vertical direction) on the 3 gaussian trans-
forms.
In the graphic window you use Application / Save in Parameter File to save these variograms under
the name Gaussian Fe-SiO2-P rich ore.
(snap. 2.6-18)
In Statistics/Variogram Fitting..., choose the experimental variogram you just saved. Create the
new variogram model with the same name Gaussian Fe-SiO2-P rich ore. Set the toggles Global
Window and ask to display the number of pairs in the graphic window (Application/Graphic
Parameters...).
(snap. 2.6-19)
The model is made using the following method:
104
l enter the name of the new variogram model Gaussian Fe-SiO2-P rich ore and Edit it.
l in the Model Definition window click on Load Model and choose the model made for Gaussian
P rich ore no truncation. The following window pops up:*
(snap. 2.6-20)
Clck on Clear button, then move the mouse to the second line Gaussian P rich ore, click on Link
and on OK in the Selector window to put the variogram made on Gaussian P alone for the same
variable in the three variate variogram. Then you click on OK in the Model Loading window.
l in the Variogram Fitting window click on Automatic Sill Fitting. The Global Window shows
the model that has been fitted. Press Run to save it in the parameter file.
(snap. 2.6-21)
2.6.2.4 Variogram regularization
In order to perform the direct block simulation you have to model the three variate variogram on the
support of the blocks 25x25x15.
N0
1
157
472
688
11201373
1195
1196
900
1108
1222
1155
D-9 6
78
92
325 266 223
183
148
117
0 500 1000 1500
Distance (m)
0.00
0.25
0.50
0.75
1.00
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

F
e

r
i
c
h

o
r
e
N0
1
157
472 688
11201373
11951196900
1108
1222
1155
D-9
6
78
92
325
266
223
183
148
117
0 500 1000 1500
Distance (m)
-0.5
0.0
0.5
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

P

r
i
c
h

o
r
e

(
G
i
b
b
s
)
N0
1
157
472 688
1120
1373
1195
1196
900
11081222
1155
D-9
6
78
92
325
266
223
183
148
117
0 500 1000 1500
Distance (m)
0.0
0.5
1.0
1.5
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

P

r
i
c
h

o
r
e

(
G
i
b
b
s
)
N0
0
56
311
417763
896
663 701
561
683
722
632
D-9
6
75
11
262
216 182 151 122
97
0 500 1000 1500
Distance (m)
-0.5
0.0
0.5
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

S
i
O
2

r
i
c
h

o
r
e

&

G
a
N0
0
56311417
763
896
663 701
561
683
722
632
D-9
6
75
11
262 216
182
151
122
97
0 500 1000 1500
Distance (m)
-0.5
0.0
0.5
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

S
i
O
2

r
i
c
h

o
r
e

&

G
a
N0
0
56
311
417
763 896
663
701
561 683
722 632
D-9
6
75
11
262
216
182 151
122
97
0 500 1000 1500
Distance (m)
0.00
0.25
0.50
0.75
1.00
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

S
i
O
2

r
i
c
h

o
r
e
106
l You first have to launch Statistics / Modeling / Variogram Regularization. You will store in a
new experimental variogram Gaussian Fe-SiO
2
-P rich ore block 25x25x15 3 directional vari-
ograms using a discretization of 5x5x1. You will also ask to Normalize the Input Point Vario-
gram.
(snap. 2.6-22)
l Then you model the regularized variogram using Variogram Fitting and the Automatic Sill Fit-
ting mode, after having loaded the model made on the point samples Gaussian Fe-SiO2-P rich
ore. You note that the Nugget effect is put to zero. When you save the variogram model the
Nugget effect is not stored in the Parameter file
(snap. 2.6-23)
108
(snap. 2.6-24)
2.6.2.5 Gaussian Support Correction
The point gaussian anamorphosis and the regularized variogram model have to be transformed in
gaussian anamorphosis and variogram model related to the gaussian block variable Yv (gaussian
zero-mean, variance-1 variable).
This is achieved by running Statistics / Modeling / Gaussian Support Correction.
N0
N27
D-9
0 500 1000 1500
Distance (m)
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

F
e

r
i
c
h

o
r
e

(
B
l
o
c
k
N0
N27
D-9
0 500 1000 1500
Distance (m)
-0.5
0.0
0.5
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

P

r
i
c
h

o
r
e

(
G
i
b
b
s
)
N0
N27
D-9
0 500 1000 1500
Distance (m)
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

P

r
i
c
h

o
r
e

(
G
i
b
b
s
)
N0
N27
D-9
0 500 1000 1500
Distance (m)
-0.5
0.0
0.5
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

S
i
O
2

r
i
c
h

o
r
e

(
B
l
o
N0
N27
D-9
0 500 1000 1500
Distance (m)
-0.5
0.0
0.5
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

S
i
O
2

r
i
c
h

o
r
e

(
B
l
o
N0
N27
D-9
0 500 1000 1500
Distance (m)
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

S
i
O
2

r
i
c
h

o
r
e

(
B
l
o
(snap. 2.6-25)
2.6.2.6 Direct Block Simulation
It is achieved by running the menu Interpolate / Conditional Simulations / Direct Block Simulation.
It takes some time to get 100 simulations. Depending on the computer it may be more than an hour.
110
l The simulated variables are created with the following names Simu block Gaussian Fe rich
ore ...in the 3D Grid 25x25x15. We store the gaussian values before transform to allow a check
of the experimental variograms on gaussian simulated values with the input variogram model,
that is defined on the gaussian variables.
l The Block Anamorphosis and the Block Gaussian Model are those obtained from the Gaussian
Support Correction.
l The Neighborhood used for kriging Fe rich ore is modified into a new one called Fe rich ore
simulation changing the radius along V to 800m. The reason is just because the Local Parame-
ters for the neighborhood are not implemented in the application Direct Block Simulation.
l Number of simulations: 100 for instance .
l We ask to not Perform a Gaussian Back Transformation, for the reason explained above. The
back transform will be achieved afterwards.
l The turning bands algorithm is used with 1000 Turning Bands.
(snap. 2.6-26)
You can compare the experimental variograms calculated from the 100 simulations in up to 3 direc-
tions with the input variogram model. The directions are entered by giving the increments (number
of grid mesh) of the unit directional lag along X, Y, Z. For instance for the direction 1, the incre-
ments are respectively 1, 0, 0, which makes the unit lag 25m East-West.
112
(snap. 2.6-27)
Three graphic pages (one per direction) are then displayed. The average experimental variograms
are displayed with a single line, the variogram model with a double line. On the next figure the var-
iograms in the direction 3 show a good match up to 100m. For the cross-variogram P-SiO2 where
the correlation is very low, some simulations look anomalous, further analysis could be made to
exclude these simulations for the next post processing steps.
(snap. 2.6-28)
It is then necessary to transform the simulated gaussian values into raw values, using Statistics /
Data Tranformation / Raw Gaussian Transformation. For transforming the three grade you will
have to run that menu three times. You should choose as Transformation Gaussian to Raw Trans-
formation. The New Raw Variable will be created with the same number of indices with names
like Simu block Fe rich ore...
The transform is achieved by means of the block anamorphosis Fe-SiO2-P rich ore block
25x25x15, do not forget to choose on the right side of the Anamorphosis window the right variable.
0 25 50 75 100 125
Distance (m)
0.00
0.25
0.50
0.75
1.00
1.25
V
a
r
i
o
g
r
a
m

:

S
i
m
u

b
l
o
c
k

G
a
u
s
s
i
a
n

F
e

r
i
c
h
0 25 50 75 100 125
Distance (m)
-0.5
-0.4
-0.3
-0.2
-0.1
0.0
V
a
r
i
o
g
r
a
m

:

S
i
m
u

b
l
o
c
k

G
a
u
s
s
i
a
n

P

r
i
c
h

0 25 50 75 100 125
Distance (m)
0.00
0.25
0.50
0.75
1.00
V
a
r
i
o
g
r
a
m

:

S
i
m
u

b
l
o
c
k

G
a
u
s
s
i
a
n

P

r
i
c
h

0 25 50 75 100 125
Distance (m)
-0.4
-0.3
-0.2
-0.1
0.0
V
a
r
i
o
g
r
a
m

:

S
i
m
u

b
l
o
c
k

G
a
u
s
s
i
a
n

S
i
O
2

r
i
0 25 50 75 100 125
Distance (m)
-0.03
-0.02
-0.01
0.00
0.01
0.02
0.03
0.04
0.05
V
a
r
i
o
g
r
a
m

:

S
i
m
u

b
l
o
c
k

G
a
u
s
s
i
a
n

S
i
O
2

r
i
0 25 50 75 100 125
Distance (m)
0.00
0.25
0.50
0.75
1.00
1.25
V
a
r
i
o
g
r
a
m

:

S
i
m
u

b
l
o
c
k

G
a
u
s
s
i
a
n

S
i
O
2

r
i
114
(snap. 2.6-29)
We can now combine the simulations of the rich ore indicator and the grades simulations, by chang-
ing to undefined (N/A) the grades when the block is simulated as poor ore (simulated code 2).
These transformations have to be applied on the 100 simulations using File / Calculator. It is com-
pulsory to create beforehand new macro variables, with 100 indices, called Simu block Fe ... with
Tools / Create Special Variable.
(snap. 2.6-30)
116
(snap. 2.6-31)
If you complete this Case Study by simulating also the grades of poor ore, you will get valuated
grades for all blocks in the orebody. The displays will be presented in the last chapter.
2.6.3 Simulations post-processing
One main advantage of simulations is the possibility to apply non linear calculations (for example
applying different cut-off grades simultaneously, or calculation of the probability for a grade to be
above a threshold etc.) for local reserves estimation. The post-processing may be applied on the
simulated blocks, but in the present case it is more interesting to first regroup the simulated blocks
in the blocks 75x75x15 (called panels) and illustrate some basic post-processing on the tonnage and
metals of rich ore within those panels.
2.6.3.7 Regrouping blocks into panels
We will calculate for each panel the mean grade, tonnage and metal quantitiy of rich ore and the
quantity of rich ore Fe-P-SiO
2
by using Tools / Processing / Grade Simulation Post-Processing, that
applies directly on the macro-variables. The Grade Simulation Post-processing is designed to calcu-
late local grade tonnage curves on panel grid (Q,T,M variables) from simulated grade variables on
block grid. The grade variables can be simulated using the panel Turning bands, Sequential Gauss-
ian Simulation or any kind of simula-tion that generates continuous variables.
The Block Grid usually corresponds to the S.M.U. (Selective Mining Unit). It has to be consistent
with the Panels, in other words the Block Grid must make a partition of this Panel Grid.This
appli-cation handles multivariable cases with a cuttof on the main variable.
(snap. 2.6-32)
118
(snap. 2.6-33)
2.6.3.8 Examples of Post Processing
The menu Tools / Simulation Post-processing offers different options, illustrated hereafter on the
Tonnage and Metal variables stored on the 3D Grid 75x75x15m file:
l Statistical Maps to calculate the average of 100 simulated tonnages
(snap. 2.6-34)
(snap. 2.6-35)
The mean tonnage may be compared to the kriged indicator (after multiplication by the panel ton-
nage).
120
l Iso-Frequency Maps to calculate the quantile at the frequencies of 25%-50%-75% of the Ton-
nage of rich ore. In the previous Simulation Post-Processing window, click the Toggle button
Iso-Frequency Maps, the following window pops up and you define a New Macro Variable
Quantile Tonnage rich ore[xxxxx].
(snap. 2.6-36)
then click Quantiles and choose for Step Between Frequencies 25%. You get a macro-variable with
3 indices, one per frequency: for each panel the tonnage such that 25%, 50%, 75% of the simula-
tions is lower than the corresponding quantile value.
(snap. 2.6-37)
l Iso-Cutoff Maps to calculate the probability for the Metal P rich ore to be above 0, 50, 100,
150, 200.
(snap. 2.6-38)
In the previous Simulation Post-Processing window, click the Toggle button Iso-Cutoff Maps, the
following window pops up and you define a New Macro Variable for Probability to be Above Cut-
off (T), i.e. Proba P rich ore above[xxxxx].
122
(snap. 2.6-39)
then click Cutoff and click Regular Cutoff Definition and choose the parameters as shown below.
You get a macro-variable with 4 indices, one per cutoff: for each panel the probability to be above
0.02,0.03 ...
(snap. 2.6-40)
l Risk Curves to calculate the distribution of 100 simulations of Fe metal quantities of rich ore
over the orebody.
(snap. 2.6-41)
Click Risk Curves then Edit and fill the parameters in the Risk Curves & Printing Format window,
as shown. Only the Accumulations are interesting. For a given simulation the accumulation is
obtained by multiplying the simulated block value (here the Fe metal in tons) by the volume of the
block. It means that the average grade of the block is multiplied twice by the block volume. That is
why in order to get the metal in MTons we have to apply a scaling factor of 75x75x15 (84375) and
multiply it by 10
6
. That scaling is entered in the box just on the left of m3*V_unit of the Accumula-
tions sub-window. By asking Print Statistics the 100 accumulations will be output in the Isatis mes-
sage window. The order of the printout depends of the option Sorts Results by, here we ask
Accumulations.
124
(snap. 2.6-42)
Coming back to the Simulation Post-processing window and press Run. The following graphic is
then displayed.
(snap. 2.6-43)
With the Application / Graphic Parameters you may Highlight Quantiles with the Simulation Value
on Graphic.
(snap. 2.6-44)
The graphic page is refreshed as shown.
126
(snap. 2.6-45)
In the message window we get the 100 simulated metal quantities by increasing order. The column
Macro gives the index of the simulation for each outcome: for instance the minimum metal is
obtained for the simulation #72, the next one for the simulation 97 ...
Rank Macro Frequency Accumulation Volume
1 72 1.00 1140.90MT 3442162500.00m3
2 97 2.00 1156.65MT 3442162500.00m3
3 38 3.00 1171.82MT 3442162500.00m3
4 15 4.00 1179.91MT 3442162500.00m3
5 91 5.00 1181.25MT 3442162500.00m3
6 41 6.00 1185.01MT 3442162500.00m3
7 30 7.00 1191.53MT 3442162500.00m3
8 45 8.00 1191.71MT 3442162500.00m3
9 57 9.00 1194.86MT 3442162500.00m3
10 59 10.00 1195.80MT 3442162500.00m3
11 35 11.00 1196.15MT 3442162500.00m3
12 6 12.00 1196.37MT 3442162500.00m3
13 48 13.00 1197.58MT 3442162500.00m3
14 62 14.00 1199.70MT 3442162500.00m3
15 40 15.00 1201.25MT 3442162500.00m3
16 1 16.00 1201.90MT 3442162500.00m3
17 86 17.00 1204.47MT 3442162500.00m3
18 33 18.00 1206.65MT 3442162500.00m3
19 93 19.00 1206.83MT 3442162500.00m3
20 11 20.00 1210.44MT 3442162500.00m3
...
(snap. 2.6-46)
128
2.7 Displaying the Results
The last chapter consists in visualizing the different result in the 3D grids, through the 2D Display
facility then through the 3D viewer.
2.7.1 Using the 2D Display
2.7.1.1 Display of the Kriged block model
We are going to create a new Display template (Display/New Page...), that consists in an overlay of
a grid raster and isolines. All the Display facilities are explained in detail in the "Displaying & Edit-
ing Graphics" chapter of the Beginner's Guide.
Click on Display / New Page in the Isatis main window. A blank graphic page pops up, together
with a Contents window. You have to specify in this window the contents of your graphic. To
achieve that:
l First, give a name to the template you are creating: Kriging Fe rich ore. This will allow you to
easily display again this template later.
l In the Contents list, double click the Raster item. A new window appears, in order to let you
specify which variable you want to display and the color scale:
m Select the Grid file, 3D Grid 75x75x15m with selection orebody active, select the variable
Kriging Fe rich ore
m Specify the title for the Raster part of the legend, for instance Kriging Fe rich ore
m In the Grid Contents area, enter 16 for the rank of the section XOY to display
m In the Graphic Parameters area, specify the Color Scale you want to use for the raster dis-
play. You may use an automatic default color scale, or create a new one specifically dedi-
cated to the Fe variable. To create a new color scale: click the Color Scale button, double-
click on New Color Scale and enter a name: Fe, and press OK. Click the Edit button. In the
Color Scale Definition window:
- In the Bounds Definition, choose User Defined Classes.
- Choose Number of Classes 22,
- Click on the Bounds... button, enter 60 and 71 as the Minimum and Maximum values.
Press OK.
- Switch on the Invert Color Order toggle in order to affect the red colors to the large Fe
values.
- Click Undefined Values button and select Transparent.
- In the Legend area, switch off the Display all tick marks button, enter 60 as the reference
tickmark and 2 as the step between the tickmarks. Then, specify that you do not want
your final color scale to exceed 7 cm. Switch off the Automatic Format button, and spec-
ify that you want to use integer values of Length 7. Ask to display the Extreme Classes.
Click OK.
(snap. 2.7-1)
m In the Item contents for: Raster window, click Display current item to display the result.
m Click OK.
l Double-click on the Isolines item. A new Item contents window appears. In the Data area, select
the Kriging Fe rich ore variable from the 3D Grid file with the same selection. In the Grid
Contents area, select the rank 16 for the XOY section. In the Data Related Parameters area,
switch on the C1 line, enter 60 and 71 as lower and upper bounds and choose a step equal to 2.
130
Switch off the Visibility button. Click on Display Current Item to check your parameters, then
on Display to see all the previously defined components of your graphic. Click on OK to close
the Item contents window.
l In the Item list, you can select any item and decide whether or not you want to display its leg-
end, by setting the toggle Legend ON. Use the Move Front and Move Back buttons to modify
the order of the items in the final Display.
l Close the Contents window. Your final graphic window should be similar to the one displayed
hereafter.
(fig. 2.7-1)
You can also visualize your 3D grid in perspective. Open again the Contents window of the previ-
ous graphic display (Application/Contents...). Switch the Representation Type from Projection to
Perspective:
500 1000 1500 2000
X (m)
0
1000
2000
3000
Y

(
m
)
Kriging Fe rich ore
Kriging Fe rich ore
N/A
70
68
66
64
62
60
132
l just click on Display: the previous section is represented within the 3D volume. Because of the
extension of the grid, set the vertical axis factor to 3 in the Display Box tab (switch the toggle
Automatic Scales OFF). In the Camera tab, modify the Perspective Parameters: longitude=60,
latitude=40.
(fig. 2.7-2)
335
435
535
635
735
335
435
535
635
735

1
6
3

1
1
6
3

2
7
5

1
2
7
5

2
2
7
5

Kriging Fe rich ore
Kriging Fe rich ore
N/A
70
68
66
64
62
60
l Representing the whole grid as a solid: this is obtained by setting the 3D Grid contents to 3D
Box, both in the Raster and Isolines item contents windows.
l Representing the 3G grid as a solid and penetrating into the solid by digging a portion of the
grid. For each item content window (for raster and isolines), set the 3D Grid contents to Exca-
vated Box. Then define the indices of the excavation corner (for instance: cell=17, 21, 15).
(fig. 2.7-3)
In the Contents window, the Camera tab allows you to animate (animate tab from the main contents
window) the graphic in several ways:
l by animating the entire graphic along the longitude or latitude definition,
l by animating one item property at a time, for instance the grid raster section. To interrupt the
animation, press the STOP button in the main Isatis window.
2.7.1.2 Display of the simulated block model
l Fe grade
m Create a raster image of the Fe simulated macro variable: choose the first simulation (index
1). Display rank 16 of the 25x25x15 m 3D grid file, so you can compare simulations with
the kriging) and choose the grade Fe color scale. Ask to display the legend.
m Create a Base map of the composite data from the Composites 15 m with the selection final
lithology{Rich ore} active and no variable in order to use the same Default Symbol a full
circle of 0.15cm.
335
435
535
635
735
335
435
535
635
735

1
6
3

1
1
6
3

2
7
5

1
2
7
5

2
2
7
5

Kriging Fe rich ore
Kriging Fe rich ore
N/A
70
68
66
64
62
60
134
(snap. 2.7-2)
In the Display Box tab from the contents window, set the mode to Containing a set of items and
click the Raster item: set the toggle Box Defined as Slice around Section ON and set the Slice
Thickness to 45 m.
(snap. 2.7-3)
Press Display:
136
(fig. 2.7-4)
From the Animate tab, select the raster item and choose to animate on the macro index. Set the
Delay to 1s and press Animate. The different simulations appear consecutively: the animation
allows to sense the differences between the simulations. Check that the simulations tend to be simi-
lar around boreholes.
l Display of the probability for the Metal P of rich ore in panels to be above cut-off = 50T:
m Create a new page and display the macro variable Proba P rich ore above from the 3D
Grid 75x75x15m file: choose the macro index n 2 (i.e. cutoff = 50)
m Legend title: probability
m Ask to display rank 16 (horizontal section 16)
m =Make a New Color scale named Proportion as explained before for Fe, but with 20 classes
between 0 and 1.
m press OK
500 1000 1500 2000
X (m)
0
1000
2000
3000
Y

(
m
)
Simu block Fe[00001]
Fe rich ore
N/A
70
68
66
64
62
60
l Ask for the legend and press Display:
(fig. 2.7-5)
2.7.2 Using the 3D Viewer
Launch the 3D Viewer (Display/3D Viewer...).
2.7.2.3 Borehole visualization
l Display the Fe composites:
m Drag the Fe variable from the Composites 15 m file in the Study Contents and drop it in the
display window;
m Magnify by a factor of 2 the scale along Z by clicking the Z Scale button at the top of the
graphic page.
m Click Toggle the Axes in the menu bar on the left of the graphic area.
m From the Page contents, click right on the 3D Lines object to open the 3D Lines properties
window. In the 3D Lines tab
500 1000 1500 2000
X (m)
0
1000
2000
3000
Y

(
m
)
Proba P rich ore above{50.000000}
Probability
N/A
1.00
0.90
0.80
0.70
0.60
0.50
0.40
0.30
0.20
0.10
0.00
138
- select the Tube mode;
- switch on the toggle Selection and choose the final lithology{Rich ore} macro index;
- switch off the toggle Allow Clipping
(snap. 2.7-4)
m In the Color tab, choose the same Fe Isatis color scale;
m In the Radius tab, set the mode to constant with a radius of 20 m
m Press Display and close the 3D Lines Properties window
m In the File menu click Save Page as and give a name (composites rich ore) in order to be
able to recover it later if you wish.
(snap. 2.7-5)
2.7.2.4 Display of the kriged 3D Block model
As an example we will display the kriged indicator of rich ore. In order to make a New Page click
Close Page in the File menu.
l Click Compass in the menu bar on the left of the graphic area.
l Drag the Kriging indicator rich ore variable from the 3D Grid 75 x 75 x 15 m file in the Study
Contents and drop it in the display window;
l Click right on the 3D Grid 75x75x15m file in the Page Contents to open the 3D Grid Proper-
ties:
140
m In the 3D Grid tab, tick the selection toggle, choose the orebody selection;
m in the color tab:
- set the color mode to variable and change the variable to Kriging Indicator rich ore;
- apply the Rainbow reversed Isatis color scale;
- Press Display and close the 3D Grid properties window
(fig. 2.7-6)
l Investigate inside the kriged block model:
m open the clipping plane facility from Toggle the Clipping Plane in the menu bar on the left
of the graphic area: the clipping plane appears across the block model;
m Go in select mode by pressing the arrow button in the function bar;
m Click the clipping plane rectangle and drag it next by the block model for better visibility;
m Click one of the clipping planes axis to change its orientation (be careful to target precisely
the axis itself in dark grey, not its squared extremity nor the center tube in white)
m Add the drill holes (Fe rich ore) as you did for the previous graphic page
m Open the Line Properties window of the Composites 15 m file: set the Allow Clipping tog-
gle ON;
m Click on the clipping planes center white tube and drag it in order to translate the clipping
plane along the axis: choose a convenient cross section, approximately in the middle of the
block model. You may also benefit from the clipping controls parameters available on the
right of the graphic window in order to clip a slice with a fixed width and along the main
grid axes.
m Click on one block of particular interest: its information is displayed in the top right corner:
(snap. 2.7-6)
You may click also on boreholes to display composite data.
142
l Slicing (before hand, click on Toggle the Clipping Plane)
m Edit the 3D Grid 75x75x15m attributes, go in the Slicing tab and set the properties as fol-
low:
(snap. 2.7-7)
Set the toggle Automatic Apply ON, and move the slices to visualize interactively the slicing.
l Save the graphic as a New Page with the name Composites and kriged indicator rich ore.
2.7.2.5 Display of the search ellipsod
From the kriging application (the definition parameters of the 3D kriging of Fe should be kept),
launch the Test window. From Application/Target Selection, select the grid node (20,19,14) for
instance and press Apply. Then, make sure that the 3D viewer is running and, from the same Appli-
cation menu of the Test window, ask to Link to 3D Viewer: a 3D representation of the search ellip-
sod neighborhood is represented, and the samples used for the estimation of this particular node are
highlighted. A new graphic object neighborhood appears in the Page Contents from which you
may change the graphic properties (color, size of the samples for coding the weights or the Fe val-
ues etc...)
(fig. 2.7-7)
144
Non Linear 145
3 Non Linear
This case study, dedicated to advanced users, is based on the Walker
Lake data set, which has been first introduced and analyzed by Mohan
SRIVASTAVA and Edwards H. ISAAKS in their book Applied Geosta-
tistics (1989, Oxford University Press).

Geostatistical methods applicable to perform global and local estima-
tion of recoverable resources in a mining industry context are
described through this case study:
Non linear methods, including four methods used to estimate local
recoverable resources: indicator kriging, disjunctive kriging, uniform
conditioning and service variables.
Conditional simulations of grades, using the two main methods appli-
cable: turning bands and sequential gaussian.

The efficiency of these methods will be evaluated by comparison to
the reality, which can be considered as known in this case because
of the origin of the data set.
Reminder: while using Isatis, the on-line help is accessible anytime by
pressing F1 and provides full description of the active application.

Last update: Isatis version 2012
146
3.1 Introduction and overview of the case study
This case study is dedicated to advanced users who feel comfortable with linear geostatistics and
Isatis.
3.1.1 Why non linear geostatistics?
Non linear geostatistics are used for estimating the recoverable resources. At the difference of the
estimation of in situ resources by conventional kriging (linear geostatistics), the estimation of the
recoverable resources considers the mining aspects of the question. Three points can effectively be
taken into account by non linear geostatistics:
l the support effect, that makes the recovered ore depending on the volume on which the ore/
waste decision is made. In this case the size of the selective mining unit (SMU or blocks) has
been fixed to 5m x 5m. When performing the local estimations we will calculate the ore tonnage
and grade after cut-off in panels of 20m x 20m. It is important to keep these terms of block for
the selective unit and panel for the estimated unit (e.g.: tonnage within the panel of the ore con-
sisting of blocks with a grade above the cut-off). These terms are systematically used in the Isa-
tis interface.
l the information effect, that makes the mis-classification between selected ore and waste depend-
ing on the amount of information used in estimating the blocks. At this stage two notions are
important. Firstly the recovered ore is made of true grades contained in blocks whose estimated
grade is above the cut-off. Secondly the decision between ore and waste will be made with an
additional information (blast-holes...) in the future of the production. The question is then what
can we expect to recover tomorrow, if we assume a future pattern blast-holes for instance.
l the constraint effect, that leads for any technical/economical reason to ore dilution or ore left in
place. The two previously mentioned effects are assuming a free selection of blocks within the
panels, only the distribution of block grades is of importance. When their spatial distribution has
to be considered (the recovered ore will be different if rich blocks are contiguous or spread
throughout the panel), only geostatistical simulations provide an answer.
3.1.2 Organization of the case study
This case study is divided in several parts: the first part 3.2 Preparation of the case study
rehearses geostatistical concepts and Isatis manipulation already described in the In Situ 3D
Resource Estimation case study, consisting of declustering, grid manipulations, variography, ordi-
nary kriging with neighborhood creation. These topics will not be detailed here and the user is
invited to have a look at the previous case study for an extensive description. The following of the
case study describes several different methods for the estimation of recoverable resources; it is also
recommended that the user reads 3.3 Global estimation of recoverable resources before start-
ing any method described in 3.4 Local estimation of the recoverable resources or in 3.5 Sim-
ulations. The dataset allows to compare estimations with real measurements: this will be done
exhaustively in 3.6 Conclusions.
Non Linear 147
3.1.2.1 Global Estimation of Recoverable Resources (developed in 3.3)
The global estimation makes use of the raw data histogram (possibly weighed by declustering coef-
ficients): each grade is attached to a frequency, i.e the global proportion relative to the global ton-
nage of the deposit assuming a perfect sampling. This is a direct statistical approach. Geostatistics
appears as soon as the variogram is used to correct this histogram, i.e the proportion, to reflect the
support effect and/or the information effect. Thus, a histogram model is needed in order to perform
these corrections: the modeling and the corrections are done through the Gaussian Anamorphosis
Modeling and Support Effect panels in Isatis, widely used through the whole case study. Compari-
son to reality and kriging will be done through global grade-tonnage curves.
3.1.2.2 Local estimation of recoverable resources
The local estimation of recoverable resources makes use of non linear estimation or simulation
techniques, involving gaussian anamorphosis. The aim is to estimate the proportion of ore blocks
within larger panels (assuming free selection of blocks within each panel), and the corresponding
metal tonnage and mean grade above cut-off:
l by non linear kriging techniques (developed in 3.4): the main advantage of these methods is
their swiftness, but no information on the location of the ore blocks within the panels is given.
Four methods will be described: Indicator kriging, Disjunctive kriging, Service variables and
Uniform Conditioning.
l by simulation techniques (developed in 3.5): the main advantages of simulations is the possi-
bility to derive simulated histograms and estimate the constraint effect, but the method is quite
heavy and time consuming for big block models. Two methods will be described: Turning
Bands (TB) and Sequential Gaussian Simulations (SGS).
Comparison to reality through a specific analysis of the 600 ppm cut-off will be done through
graphic displays and cross plots of the ore tonnage and mean grade above cut-off.
Note - If you wish to compare the local estimates with reality you will need first to calculate the
real tonnage variables from the real grades for the specific cut-off 600 (this is done in 3.4.1
Calculation of the true QTM variables based on the panels).
148
3.2 Preparation of the case study
The dataset is derived from an elevation model from the western United States, the Walker Lake
area in Nevada. It has been transformed in order to represent measures of concentration in some
elements (economic grades in the deposit we are going to evaluate). From the original data set we
will use only the variable V, considered as the grade of an ore mineral measured in ppm: the multi-
variate aspect of this data set will not be considered, as the non linear estimation methods available
in Isatis are currently univariate (unlike simulations). The data set is two fold, the exhaustive data
set, containing 78 000 measurements points on a 1m x 1m grid, and the sample set resulting from
successive sampling campaigns and containing 470 data locations. Several methods for the estima-
tion of recoverable resources are proposed in Isatis: this case study aims to describe them all and
compare them to the reality issued from the exhaustive set.
3.2.1 Data import and declustering
The data is stored in the Isatis installation directory (sub-directory Datasets/Walker_Lake). Load
the data from ASCII file by using File / Import / ASCII. The ASCII files are Sample_set.hd for the
sample set and Exhaustive_set.hd for the exhaustive data set. The files are imported into two sepa-
rate directories Sample set and Exhaustive set respectively, and files are called Data.
(snap. 3.2-1)
By visualizing the Sample set data (using Display / Basemap/ Proportional), we immediately see
the preferential sampling pattern of high grade zones:
Non Linear 149
(fig. 3.2-1)
In order to correct the bias of preferential sampling of high grade zones, it is necessary to declus-
ter the data. To do so you can use Tools / Declustering: it performs a cell declustering with a mov-
ing window centered on each sample. We store the resulting weights in a variable Weight of the
sample data set: this variable will be used later to weight statistics for the variographic analysis in
the EDA and the gaussian anamorphosis modeling. The moving window size for declustering has
been fixed here to 20m x 20m, accordingly to the approximative sampling loose mesh size outside
the clusters.
Note - A possible guide for choosing the moving window dimensions is to compare the value of the
resulting declustered mean to the mean of kriged estimates (kriging has natural declustering
capabilities).
The statistics before and after declustering are the following:
0
0
100
100
200
200
X (m)
X (m)
0
100
200
300
Y

(
m
)
V
150
(snap. 3.2-2)
Mean: 436.35 -> 279.68
Std dev: 299.92 -> 251.44
The next graphics correspond to the histograms of the Sample set, Exhaustive set and Declustered
sample set; they have been calculated using Statistics / Exploratory Data Analysis (EDA). The his-
togram of the Declustered sample set has been calculated with the Compute Using the Weight Vari-
able option toggle ON, using the Weight variable.
Non Linear 151
(snap. 3.2-3)
152
(fig. 3.2-2)
From these three histograms we clearly see that the declustering process will allow to better repre-
sent the statistical behavior of the phenomenon.
3.2.2 Variographic analysis of the sample grades
We first focus on possible anisotropies of the sample set data. From the Statistics / Exploratory
Data Analysis panel, activate the option Compute using the Weight Variable: we will calculate a
weighted 2D variogram map on the V variable from the sample dataset. By default, the Reference
Direction is set to an azimuth equal to the North (Azimuth = N0.00). The parameters related to
the directions, lags and tolerance may be tuned for a detailed variographic analysis but here we will
base ourselves directly on common parameters: ask for 18 directions (10 degrees each), and we
will define 11 lags of 15 m. Generally, the variogram is calculated with a tolerance on distance set
to 50% of the lag which corresponds to a Tolerance on Lags equal to 0 lag; besides, calculations
are often made with an angular tolerance of 45 (in order to consider all samples once with two
directions) which corresponds to a Tolerance on Directions equal to 4 sectors (4 sectors of 10 +
half sector 5 = 45 ).
If the focus is on short scale, one may decide to calculate a bi-directional variogram along N70 and
N160, considering that N160 is a direction of maximum continuity.
Note - This short scale anisotropy is not clearly visible on the variogram map below: to better
visualize it, you may re-calculate the variogram map on 5 lags only and create a customized color
scale through Application / Graphic Specific Parameters...
Non Linear 153
In the variogram map area you can activate a direction using the mouse buttons, the left one to
select a direction, and the right one for selecting Activate Direction in the menu. Activating both
principal axes (perpendicular directions N160 and N70) displays the corresponding experimental
variograms below. When selecting the variogram, click right and ask for Modify Label... to change
N250 to N70:
(snap. 3.2-4)
The short scale anisotropy is visible on the experimental variogram; it is then saved in a parameter
file Raw V from the graphic window (Application / Save in Parameter File...).
We now have to fit a model based on these experimental variograms using the Statistics / Variogram
Fitting facility. We fit the model from the Manual Fitting tab.
154
(snap. 3.2-5)
Non Linear 155
(snap. 3.2-6)
156
(snap. 3.2-7)
Press Print to check the output variogram and then save the variogram model in the parameter file
under the name of Raw V. It should be noted that the total sill of the variogram is slightly above the
dispersion variance and the low nugget value has been chosen.
3.2.3 Calculation of the true block and panel values
In this case study, during the mine exploitation period, a 5m x 5m block will be the selective mining
unit (SMU). The recoverable resource estimation will be based on this 5m x 5m block support; but
first, the in-situ resource estimation will be done on 20m x 20m panels for more robust estimation.
As we have access to an exhaustive data set of the whole area to be mined, we can assume that we
know the true values for any size of support, just by averaging the real values of the exhaustive
set on the wanted block or panel support.
3.2.3.1 Calculation of the true grade values for 5 m x 5 m SMU blocks
To store this average value on a 5m x 5m block support, we need to create a new grid (new file
called Grid 5*5 in a new directory Grids, using the File / Create Grid File facility) and choose the
coordinates of the origin (center of the block at the lowest left corner) in order to match exactly the
data. The Graphic Check, in Block mode, will help to achieve this task. Enter the following grid
parameters:
m X and Y origin: 3m,
m X and Y mesh: 5m,
m 52 nodes along X, 60 nodes along Y.
Non Linear 157
(snap. 3.2-8)
Using this configuration we have exactly 25 samples from the exhaustive data set for each block of
the new grid. Edit the graphic parameters to display the auxiliary file.
158
(snap. 3.2-9)
(fig. 3.2-3)
Now we need to average the real values on this Grid 5*5 file, using Tools / Copy Statistics / Points
-> Grid. We will call this new variable True V.
Note - Using a moving window equal to zero for all the axes, we constrain the new Mean variable
to a calculation area of 5m x 5m (1 block).
Non Linear 159
(snap. 3.2-10)
(fig. 3.2-4)
Display of the true block grade values (5m x 5m blocks)
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
True V
N/A
1000
900
800
700
600
500
400
300
200
100
0
160
The above figure is a result of two basic actions of the Display Menu: a display grid raster of the
true block grade is performed, then isolines are overlaid. Isolines range from 0 to 1500 by steps of
250 ppm, 1000 ppm isoline has been represented with a bold line type. The color scale has been
customized to cover grades between 0 and 1000 ppm, even if there are values greater than this
upper bound. Each class has a width of 62.5 ppm, the extreme values are represented using the
extreme colors.
Note - Keep in mind that V variable has primarily been deduced from elevation data: we clearly
see on the above map a NW-SE valley responsible of the anisotropy detected during variography.
The Walker Lake itself (consequently with zero values...) is in this valley. One could raise
stationarity issues, as the statistical behavior of elevation data differs from valleys (with a lake) to
nearby ranges. This is not the subject of this case study.
3.2.3.2 Calculation of the true grade values for 20 m x 20 m panels
Create a new grid file Grid 20*20 in the Grids directory with the following parameters:
m X and Y origin: 10.5 m,
m X and Y mesh: 20 m,
m 13 nodes along X, 15 nodes along Y
Non Linear 161
(snap. 3.2-11)
The graphic check with the Grid 5*5 shows that the 5m x 5m blocks describe a perfect partition of
the 20m x 20m panels. This allows to use the specific Tools / Copy Statistics / Grid to grid... for cal-
culating the true panel values True V for the Mean Name:
162
(snap. 3.2-12)
Non Linear 163
3.2.4 Ordinary Kriging - In situ resource estimation
The in-situ resource estimation will be done on the 20 m x 20 m panels through Interpolate / Esti-
mation / (Co)-Kriging...:
l Type of calculation: Block
l Input file: Sample Set/Data/V
l Output file: Grids/Grid 20*20 /Kriging V
l Model: Raw V
l Neighborhood: create a moving neighborhood named octants without any rotation and a con-
stant radius of 70 m, made of 8 sectors with a minimum of 5 samples and the optimum num-
ber of samples by sector set to 2. This neighborhood will be used extensively throughout the
case study.
(snap. 3.2-13)
164
(snap. 3.2-14)
For comparison purposes, it is interesting to do also the same kriging on the small blocks (Grid 5*5)
to quantify the smoothing effect of linear kriging.
3.2.5 Preliminary conclusions
Basic statistics may be done through different runs of Statistics / Quick Statistics...; the results are
summarized below. Interpolation by Inverse Distance ID2 with a power equal to 2 and the same
neighborhood has been done for comparison (through Interpolate / Interpolation / Quick Interpola-
tion...):
VARIABLE Count Minimum Maximum Mean Variance
True V punctual 78000 0.0 1631.2 278.0 62422
Sampled V punctual (declus.) 470 0.0 1528.1 279.7 63221
True V blocks 5x5 3120 0.0 1378.1 278.0 52287
ID2 V blocks 5x5 3120 1.6 1279.3 299.1 39031
Kriging V blocks 5x5 3120 -50.9 1361.1 275.4 44013
True V panels 20x20 195 2.2 997.8 278.0 37617
ID2 V panels 20x20 195 0.7 945.0 279.7 53539
Kriging V panels 20x20 195 -4.4 1011.3 275.8 35973
Non Linear 165
Comparing the true V values for the three different supports (punctual, block 5x5 and panel 20x20):
l as expected, the mean remains exactly identical
l the variance decreases with the support size: this is the support effect
Comparing estimated values vs. true values for one same support:
l punctual: the estimation by declustering is satisfactory because the mean and the variance are
comparable. The bias (279.7 compared to 278.0) is negligible
l block 5x5: ID2 shows an overestimation. For kriging, the bias is negligible and, as expected, the
variance of the kriged blocks (44013) is smaller than the real block variance (52287); this is the
smoothing effect caused by linear interpolation. Beside, there are some negative estimates; the
5m x 5m blocks are too small for a robust in situ estimation.
l panel 20x20: The bias of ID2 is less pronounced, but the variance is not realistic; this is because
strong local overestimation of high grade zones. The variance of the kriged panels is smaller
than the real panel variance, but the difference is less pronounced. Moreover, there is only one
negative panel estimate.
Note - 72 SMU blocks have negative estimates indicating that the 5 m x 5 m block size is too small
in this case.
166
3.3 Global estimation of the recoverable resources
3.3.1 Punctual histogram modeling
Using Statistics / Gaussian Anamorphosis Modeling we model the anamorphosis function linking
the raw values of V (called Z in Isatis) and their normal score transform (called Y in Isatis), i.e the
associated gaussian values. In order to reproduce correctly the underlying distribution we have to
apply the Weight variable previously calculated by the Declustering tool. The Gaussian variable
will be stored under Gaussian V:
(snap. 3.3-1)
Non Linear 167
(snap. 3.3-2)
The Interactive Fitting... gives access to specific parameters for the anamorphosis (intervals on the
raw values to be transformed, intervals on the gaussian values, number of polynomials etc...): the
default parameters will be kept. The distribution function is modeled by specific polynomials called
Hermite polynomials; the more polynomials, the more precise is the fit. There are also QC graphic
windows allowing to check the fit between experimental (raw) and model histograms:
168
(fig. 3.3-1)
Punctual anamorphosis function.
Experimental data is in black, the anamorphosis is in blue.
Save the anamorphosis in a new parameter file called Point and to perform the gaussian transform
with the default Frequency inversion method. This will write the Gaussian V variable on disk and
will be used for the Disjunctive Kriging, Service Variable estimations and for the simulations.
The Point Anamorphosis is equivalent to a histogram model of the declustered raw values V; it may
be used to derive global estimation as an overall view of the potential of an orebody (Grade-Ton-
nage curves are available in the Interactive Fitting... parameters), but it does not take the support
effect nor the information effect into account. This is done hereafter.
3.3.2 Support effect correction
We are going now to quantify the support effect for 5 m x 5 m blocks; that is, how much does the 5
m x 5 m block distribution differ from the punctual grades calculated above. The following is
required:
l a model of the distribution, defined by means of a gaussian anamorphosis function
l the block variance, which can be calculated using the Krige's relationship giving the dispersion
variance as a function of the variogram.
The gaussian discrete model provides then a consistent change of support model.
Use the Statistics/Support Correction... panel with the Point anamorphosis and the Raw V vario-
gram model as input. The 5mx5m block will be discretized in 4x4. At this stage no information
effect is considered, so the corresponding toggle is not activated.
-3 -2 -1 0 1 2 3 4 5
Gaussian values
0
500
1000
1500
V
Non Linear 169
(snap. 3.3-3)
Press Calculate to calculate the Gamma(v,v), and the corresponding Real Block Variance and Cor-
rection are displayed in the message window:
_________________________________________________
_________________________________________________
| | |
| | V |
|--------------------------------------|----------|
| Punctual Variance (Anamorphosis) | 63167.25 |
| Variogram Sill | 66500.00 |
| Gamma(v,v) | 8452.79 |
| Real Block Variance | 54714.47 |
| Real Block Support Correction (r) | 0.9370 |
| Kriged Block Support Correction (s) | 0.9370 |
| Kriged-Real Block Support Correction | 1.0000 |
| Zmin Block | -0.02 |
170
| Zmax Block | 1528.12 |
|______________________________________|__________|

Note - Gamma (v,v) is calculated using random procedures; hence, different results are generated
when pressing the Calculate button. Gamma (v,v) and the resulting Real Block Variance should not
vary too much between different calculations.
By clicking on the anamorphosis and on the histogram bitmaps we can check that, after the support
effect correction, the histogram of blocks is smoother (smaller variance) than the punctual histo-
gram model:
(fig. 3.3-2)
Histograms (punctual in blue and block in red): the block histogram model is smoother
The anamorphosis function will be saved under the name Block 5m * 5m and press RUN to save it.
3.3.3 Support & information effects correction
The grade tonnage curves obtained at this stage consider that the mining selection is based upon
true SMU grade. In reality, the SMU grades will be estimated using the ultimate information from
the blast-holes. The consequence is that the grade tonnage curve is deteriorated as it ignores the
uncertainty of the estimation: this is called the information effect. Knowing the future sampling
pattern, it is possible to consider this information effect.
We suppose that, at the mining stage, there will be one blast-hole at the centre of each block. The
blocks will then be estimated from blast-holes spread on a regular grid of 5m x 5m: we will use the
grid nodes of the Grid 5*5 file to simulate this future blast-hole sampling pattern. In order to calcu-
late the grade tonnage curves taking into account the information effect from this blast-hole pattern
(i.e. the selection between ore and waste is made on the future estimated grades, and not on the real
grades), we should calculate 2 coefficients:
0 500 1000 150
V
0.0
2.5
5.0
7.5
10.0
12.5
F
r
e
q
u
e
n
c
i
e
s

(
%
)
Non Linear 171
l a coefficient that transforms the point anamorphosis in the kriged block one.
l a coefficient that allows to calculate the covariance between true and kriged blocks.
Therefore, the variance of the kriged block and the covariance between real and kriged blocks are
needed: they can be automatically calculated in the same Support Correction panel through the
Information Effect optional calculation sub-panel (... selector next to the toggle):
(snap. 3.3-4)
The final sampling mesh corresponds to the final sampling pattern to be considered: 5x5 m. Press
OK and create a new anamorphosis function Block 5m*5m with information effect. Click on Run
button; two extra support correction coefficients are calculated and are displayed when pressing
RUN from the main panel:
Block Support Correction Calculation:
_________________________________________________
| | |
| | V |
|--------------------------------------|----------|
| Punctual Variance (Anamorphosis) | 63167.25 |
| Variogram Sill | 66500.00 |
| Gamma(v,v) | 9431.85 |
| Real Block Variance | 53735.40 |
| Real Block Support Correction (r) | 0.9293 |
| Kriged Block Support Correction (s) | 0.9117 |
| Kriged-Real Block Support Correction | 0.9859 |
|______________________________________|__________|
3.3.4 Analysis of the results for the global estimation
Open Tools / grade Tonnage Curves... and activate 5 data toggles. This tool allows to compare his-
tograms from different kind of data (histogram models, grade variables, tonnage variables) and
derive grade-tonnage curves for the following QTM key variables:
172
Press Edit... for the first one and then ask for a histogram model kind of data. Choose the Point
anamorphosis function and specify 21 cut-offs from 0 to 1000:
(snap. 3.3-5)
Non Linear 173
(snap. 3.3-6)
174
(snap. 3.3-7)
Press OK then repeat the procedure for the other 4 data with the same cut-off definition and specify-
ing different curve parameters for distinguishing them:
m curve 2: choose histogram model and the Block 5m * 5m anamorphosis function
m curve 3: choose histogram model and the Block 5m * 5m with information effect anamor-
phosis
m curve 4: choose grade variable and select the True V variable from the Grid 5*5 file
m curve 5: choose grade variable and select the Kriging V variable from the Grid 5*5 file
Once the 5 curves have been edited, click on the graphic bitmaps to display the Total tonnage vs.
cut-off and the Mean grade vs. cut-off curves:
Non Linear 175
(fig. 3.3-3)
Total tonnage vs. cut-off - the block histograms are close to the true tonnages.
The ordinary kriging curve under-estimates the total tonnage for high cut-offs, showing the danger
to apply cut-offs on linear estimates for recoverable resources.
176
(fig. 3.3-4)
Mean grade vs. Cut-off
Pressing Print from the main Grade Tonnage Curves window prints the numeric values for each
cut-off. The QTM variables for the particular cut-off 600 are obtained by pressing Print (the total
tonnage T is expressed in %):
| Q | T | M
True block 5x5 | 77.954| 10.385 | 750.67
Point model | 87.738| 11.351| 772.934
Block 5*5 (no info) | 76.103| 10.084| 754.699
Kriged blocks 5x5 | 61.082| 8.077| 756.258
In 3.2.5 we have seen that linear kriging is well adapted to in situ resource estimation on panels.
But when mining constraints are involved (i.e applying the 600 cut-off on small blocks), the kriging
predicts a tonnage of 8.08% instead of 10.38%: the mine will have to deal with a 29% over-produc-
tion compared to the prediction.
On the other hand, the global estimation using the point model over-estimates the reality. The glo-
bal estimation with change of support (block 5*5 no info) gives a prediction of good quality.
Because we know the reality from the exhaustive dataset, it is possible to calculate the true block
grades taking the true information effect into account and compare it to the Block 5x5 with infor-
Non Linear 177
mation effect anamorphosis. The detailed workflow to calculate the true information effect will not
be detailed here, only the general idea is presented below:
l Sample one true value at the center of each block from the exhaustive set (representing the
blasthole sampling pattern with real sampled grades V)
l krige the blocks with these samples: this is the ultimate estimated block grades on which the
ultimate selection will be based
l select blocks where ultimate estimates > 600 and derive the tonnage
l calculate the associated QTM variables based on the true grades
We can now compare the Block 5x5 with info to the real QTM variables calculated with the true
information effect (info):
| Q | T | M
True block 5x5 | 77.95 | 10.38 | 750.67
True block 5x5 (info) | 67.92 | 9.01 | 754.11
Block 5*5 with info | 71.83 | 9.66 | 743.40

As expected, the information effect on the true grades deteriorates the real recovered tonnage and
metal quantity because the ore/waste mis-classification is taken into account: the real tonnage
decreases from 10.38% to 9.01%. The estimation from the Block 5x5 with info anamorphosis
(9.69%) is closer to this reality.
178
3.4 Local Estimation of the Recoverable Resources
We want now to perform the local estimation of the recoverable resources, i.e. the ore and metal
tonnage contained in selective 5m x 5m SMU blocks within 20 m x 20 m panels.
Four main estimation methods will be reviewed: Indicator kriging, Disjunctive kriging, Uniform
conditioning and Services variables. For a set of given cut-offs, these methods will issue the follow-
ing QTM variables:
l the total Tonnage T: the total tonnage is expressed as the percentage or the proportion of SMU
blocks that have a grade above the given cut-off in the panel. Each panel is a partition of 16
SMU blocks, i.e when T is expressed as a proportion, T = 1 means that all the 16 SMU blocks of
the panel have an estimated grade above the cut-off.
l the metal Quantity Q (also referred sometimes as the metal tonnage) is the quantity of metal
relative to the tonnage proportion T for a given cut-off (according to the grade unit);
l the Mean grade M is the mean grade above the given cut-off.
In Isatis, QTM variables for local estimations are calculated and stored in macro-variables (1 index
for each cut-off) with a fixed terminology:
l base name_Q[xxxxx] for the metal Quantity variable
l base name_T[xxxxx] for the Tonnage variable
l base name_M[xxxxx] for the Mean grade above cut-off variable
All three variables are linked by the following relation:
Q = T x M
In order to be able to compare the different methods with the reality, we need first to calculate the
real QTM variables on the panel 20 x 20 support; the cut-off is defined at 600 ppm and each
method is locally compared to reality through this particular cut-off. The global grade tonnage
curves of all methods will be displayed and commented later in the final conclusion ( 3.6).
Non Linear 179
3.4.1 Calculation of the true QTM variables based on the panels
l In Grid 5*5, create a constant 600 ppm variable named Cut-off 600 ppm: this is done through
File / Calculator window:
(snap. 3.4-1)
l Tools / Copy Statistics / Grid -> Grid: in the input area we will select the true block grades True
V from the Grid 5*5 file and the Cut-off 600 ppm as the Minimum Bound Name, i.e only cells
for which the grade is above 600 will be considered. In the output area we will store the true
tonnage above 600 under Number Name and the true grade above 600 under Mean Name in
the Grid 20*20 file. If inside a specific panel no SMU block has a grade greater than 600, then
the true tonnage of this panel will be 0 and its true grade will be undefined:
180
(snap. 3.4-2)
In order to get the true total tonnage T relevant for future comparisons (i.e the ore proportion above
the cut-off 600), we have to normalize the number of blocks contained in each panel by the total
number of blocks in one panel (16):
(snap. 3.4-3)
Non Linear 181
l The metal quantity Q is calculated as Q = T x M. When the true grade above 600 is defined, the
metal quantity is equal to M x T otherwise it is null. A specific ifelse syntax is needed to reflect
this:
(snap. 3.4-4)
if this specific ifelse syntax was not used, the metal quantity in the waste would be undefined
instead of being null.
Now, we have the true tonnage, the true mean and the true metal quantity above 600 ppm to base
our comparisons in the Grid 20*20 file.
Note - Beware that the true grade above 600 is not additive as it refers to different tonnages.
Therefore, it is necessary to use the true tonnage above 600 as weights for computing the global
182
mean of the grade over the whole deposit. Another way to compute the global mean of the grade
above 600 is to divide the global metal quantity by the global tonnage after averaging on the whole
deposit.
3.4.2 Indicator kriging
Indicator kriging is a distribution free method. It is based on the kriging of indicators defined on a
series of cut-off grades. The different kriged indicators are assumed to provide the possible distribu-
tion of block grades (after a block support correction) within each panel, given the neighboring
samples. Indicator kriging can be applied in two ways:
l Multiple indicator (co-)kriging: performs the kriging of the indicator variables with their own
variograms, independently or not, for the different cut-offs.
l Median indicator kriging: supposes that all the indicator variables have the same variogram; that
is, the variogram of the indicator based on the median value of the grade.
Multiple indicator kriging is preferable because of the de-structuring of the spatial correlation with
increasing cut-offs (the assumption of an unique variogram for all cut-offs does not hold for the
whole grade spectrum), but problems of consistency must be corrected afterwards. Besides it has
the disadvantage to be quite tedious because it requires a specific variographic analysis for each
cut-off. Incidentally it is the reason why median indicator kriging has been proposed as an alterna-
tive. In this case study we will use the median indicator kriging of the panels 20m x 20m; using Sta-
tistics / Quick Statistics..., with the declustering weights, the median of the declustered histogram is
found to be 223.9.
3.4.2.1 Calculation of the median indicator variogram
We have first to generate a Macro Indicator variable Indicator V[xxxxx] in the Sample set data file
and in the output grid, by using the Statistics / Indicator Pre Processing panel, with 20 cut-offs from
50 by step of 50.
(snap. 3.4-5)
Non Linear 183
We then calculate the experimental variogram of this macro indicator variable Indicator V [xxxxx]
with the EDA (make sure that the Weight variable is activated). When selecting the Indica-
torV[xxxxx] macro variable from the EDA, you will be asked to specify the index corresponding to
the median indicator: we have chosen the index 5 corresponding to the cut-off 250 which is close
enough to 223.9. If the same calculations parameters of the Raw V variogram are used, the anisot-
ropy is no more visible; hence, the experimental variogram will be omnidirectional and calculated
with 33 lags of 5 m. It is stored in a parameter file Model Indicator, and used through Statistics /
Variogram fitting... to fit a variogram model with the following parameters detailed below the
graphic:
(fig. 3.4-1)
13
244
863
928
1520
1208
1774
1411
2053
1875
2140
1941
2659
2222
2742
2346
2882
2546
2829
2237
3243
2405
2912
2596
3093
2549
3016
2496
2999
2661
2717
2530
2914
0 50 100 150
Distance (m)
0.0
0.1
0.2
0.3
V
a
r
i
o
g
r
a
m

:

I
n
d
i
c
a
t
o
r

V
{
2
5
0
.
0
0
0
0
0
0
}
Isatis
Sample set/Data
- Variable #1 : Indicator V{250.000000}
Experimental Variogram : in 1 direction(s)
D1 :
Angular tolerance = 90.00
Lag = 5.00m, Count = 33 lags, Tolerance = 50.00%
Model : 2 basic structure(s)
Global rotation = (Az=-70.00, Ay= 0.00, Ax= 0.00)
S1 - Nugget effect, Sill = 0.035
S2 - Exponential - Scale = 45.00m, Sill = 0.21
184
It should be noted that the total sill is close to 0.25, which is the maximum authorized value for an
indicator variogram. The model is fit using the tab Manual Fitting. The variogram is saved in the
parameter file under the name Model Indicator.
(snap. 3.4-6)
Non Linear 185
(snap. 3.4-7)
186
3.4.2.2 Kriging of the indicators
We now perform the kriging of the indicators keeping the same variogram whatever the cut-off, by
using Interpolate / Estimation / Bundled Indicator Kriging:
(snap. 3.4-8)
l We ask to calculate a Block estimate: we are estimating the proportion of points above the cut-
offs within the panel.
l As Indicator Definition we define the same cut-offs as previously. In the Cut-off definition win-
dow, by clicking on Calculate proportions we get the experimental probabilities of the grade
being above the different cut-offs. These values correspond to the mean of the indicators and are
used if we perform a simple kriging. In this case because a strict stationarity is not likely, we
prefer to run an ordinary kriging, which is the default option.
l Output panels: Grid 20*20 / Indicator V[xxxxx]
l Model: Model Indicator
l The same moving neighborhood octants will be used.
3.4.2.3 Calculation of the final grade tonnage curves
At the moment we only have 20m * 20m panel estimates of probabilities for a restricted set of spec-
ified cut-offs. Now we need to perform two actions:
Non Linear 187
l rebuild the cumulative density function (cdf) of tonnage, metal and grades above cut-off for
each panel,
l Apply a volume correction (support effect) to take into account the fact that the recoverable
resources will be based on 5m * 5m blocks.
These two actions are done through Statistics / Indicator Post-processing... with the Indicator
V[xxxxx] variable from the panels as input:
(snap. 3.4-9)
l Basename for Q.T.M variables: IK. As the cut-offs used for kriging the indicators and the cut-
offs used here for representing the final grade tonnage relationships may differ (an interpolation
is needed), three different macro-variables will be created:
m IK_T{cut-off} for the ore total Tonnage T above cut-off.
m IK_Q{cut-off} for the metal Quantity Q above cut-off
m IK_M{cut-off} for the Mean grade M above cut-off.
188
l Cut-off Definition... for the QTM variables: 50 cut-offs from 0 by a step of 25.
l Volume correction: a preliminary calculation of the dispersion variance of the blocks within the
deposit is required. A simple way to achieve this consists in using the real block variance calcu-
lated by Statistics/support Correction... choosing the block size as 5 m x 5 m (cf. 3.3.2). The
Volume Variance Reduction Factor of the affine correction is calculated by dividing the Real
Block Variance (53842) by the Punctual Variance (63167). But the real block variance is calcu-
lated from the variogram sill (66500), which is superior to the punctual variance, the difference
being 3333; the real block variance needs to be corrected according to this value:
Corrected Real Block Variance = Real Block Variance - 3333 = 53842 - 3333 = 50509
Thus, the Volume Variance Reduction Factor is:
Volume Variance Reduction Factor = 50509 / 63167 = 0.802
Therefore, enter 0.802 for the Volume Variance Reduction Factor.
l two volume corrections may be applied: affine or indirect lognormal correction. As the original
distribution is clearly not lognormal we prefer to apply the affine correction, which is just
requiring the variance ratio between the 5m * 5m blocks and the points
l Parameters for Local Histogram Interpolation: we keep the default parameters for interpolating
the different parts of the histogram (linear interpolation) including for the upper tail, which is
generally the most problematic. A few tests made with other parameters (hyperbolic model with
exponent varying from 1 to 3) showed great impact on the resources. We need now to define the
maximum and minimum block values of the local block histograms: the Minimum Value
Allowed is 0; the Maximum Value Allowed may be simply approximated by applying the affine
correction by hand on the maximum value from the weighted point histogram and transpose it to
the block histogram with the Volume Variance Reduction Factor (0.8) calculated above: the
obtained value is 1391.
3.4.2.4 Analysis of the results
The Grade-Tonnage curves of the IK estimates will be displayed in 3.6 Conclusions, as for the
other following methods. Here, we will focus on the cut-off V = 600 ppm only, and compare the
results with the true values for this specific cut-off.
Below are displayed the calculated tonnage IK_T{600} compared to the true tonnage:
Non Linear 189
(fig. 3.4-2)
Tonnage T calculated by IK (SMU proportion) compared to the true tonnage.
The color scale is a regular 16-class grey palette between 0 and 1: panels for which
there is strictly less than 1 block (i.e 0 <= proportion < 0.0625) are white.
Below are displayed the calculated mean grade compared to the true grade of panels:
(fig. 3.4-3)
Mean grade calculated by IK compared to the true grades.
The color scale is a regular 16-class grey palette between 600 and 1000 and
undefined values are black: panels for which the tonnage is strictly 0 are black.
Hereafter we show the scatter diagrams of the real panel values and IK estimates for the 600 ppm
cut-off:
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
true tonnage above 600
N/A
1.000
0.875
0.750
0.625
0.500
0.375
0.250
0.125
0.000
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
IK_T{600.000000}
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
true grade above 600
ppm
N/A
1000
950
900
850
800
750
700
650
600
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
IK_M{600.000000}
190
(fig. 3.4-4)
Scatter diagram of the IK estimates vs. the true panel values above 600 ppm
(the black line is the first bisector)
At this stage of the case study we can consider that globally the indicator kriging gives satisfactory
results. At the local scale noticeable differences exist with a tendency to overestimate the grade,
especially in the upper tail of the histogram.
3.4.2.5 Disjunctive kriging
An argument against Indicator Kriging is that it ignores the relationship existing between different
cut-offs. This argument would not hold anymore, if an indicator Co-kriging was performed instead
of an univariate kriging; in practice, it is difficult to establish a model of corregionalization accept-
able for a large number of cut-offs. Disjunctive Kriging solves this problem by transforming the
cokriging problem into N krigings performed independently. One model offering this possibility is
the gaussian anamorphosis model using the Hermite polynomials where the change of support is
just explained by a coefficient (r coefficient of change of support).
In order to achieve the Disjunctive Kriging we have to provide:
l the gaussian data values Gaussian V
l the anamorphosis function on the block support Block 5m * 5m.
l the variogram model of the block gaussian variable. To determine this model we need first to
calculate an experimental block gaussian variogram using the Raw V variogram model and the
block anamorphosis. For mathematical reasons, the sill of Raw V should not exceed the punc-
tual variance of the anamorphosis, which is unfortunately the case here. Therefore, we need first
to compute another block anamorphosis including a sill normalization (cf. 3.3.2 With support
0.0 0.5 1.0
IK T{600.000000}
0.0
0.5
1.0
t
r
u
e

t
o
n
n
a
g
e

a
b
o
v
e

6
0
0
rho=0.906
600 700 800 900 1000
IK M{600.000000}
600
700
800
900
1000
t
r
u
e

g
r
a
d
e

a
b
o
v
e

6
0
0
rho=0.683
Non Linear 191
effect correction) using Statistics / Support Correction... and ask for Normalize Variogram
Sill. Store the anamorphosis in a new parameter file Block 5m * 5m (normalized) to avoid
overwriting the existing block anamorphosis Block 5m * 5m.
Open Statistics / Block Gaussian Variogram... to calculate the experimental block gaussian vario-
gram:
(snap. 3.4-10)
m Variogram model: Raw V
m Block anamorphosis: Block 5m * 5m (normalized)
m Number of directions: 2. It is convenient to make these directions coincide with the main
directions of anisotropy of the raw variogram (N160E and N70E) by setting a rotation of -
70 around positive z axis
m 20 lags of 5 m for each direction
m New experimental variogram: Block Gaussian V
We fit this variogram in Statistics / Variogram Fitting...; as expected the nugget effect has disap-
peared. Two anisotropic structures (cubic + spherical, details below the graphic) combine to a total
sill of 1, and we store the resulting model in a parameter file Block Gaussian V:
192
(fig. 3.4-5)
We are now ready to perform the Disjunctive Kriging with Interpolate / Estimation / Disjunctive
Kriging...:
N16
N70
0 25 50 75 100 125
Distance (m)
0.00
0.25
0.50
0.75
1.00
V
a
r
i
o
g
r
a
m

:

V

(
B
l
o
c
k

G
a
u
s
s
i
a
n
)
Isatis
S1 - Cubic - Range = 42.00m, Sill = 0.4
Directional Scales = ( 42.00m, 60.00m)
S2 - Spherical - Range = 40.00m, Sill = 0.6
Non Linear 193
(snap. 3.4-11)
194
l Input: Gaussian V
l Block anamorphosis...: Block 5m * 5m (normalized)
l Number of Kriged Polynomials: we use the same number as during the modeling of the anamor-
phosis function, i.e. 30.
l Cut-off definition...: we choose 21 cut-offs from 0 by steps of 50. It is compulsory to include the
zero cut-off, which should give the in situ grade estimate.
l we ask to perform Tonnage Corrections with a minimum tonnage of 0.5%.
l the Auxiliary Polynomial File will contain experimental values of the different Hermite polyno-
mials for the data points, that will be also put at the center of the closest block 5m x 5m. They
are calculated before the RUN, as soon as the output grid is defined (it may take a little time).
l Output Grid File...: in the panels Grid 20*20, store the error DK variable
l in the Panel Grid file we will also store the Q.T.M. values for each cut-off from the Basename
DK.
l we use the neighborhood octants as before.
l we choose for the Block Gaussian Variogram Model the variogram model previously fitted
Block Gaussian V.
Graphic displays of the panels for comparison with reality (proportion of SMU above 600 ppm):
(fig. 3.4-6)
Tonnage T calculated by DK (SMU proportion) compared to the true tonnage.
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
N/A
1.000
0.875
0.750
0.625
0.500
0.375
0.250
0.125
0.000
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
DK_T{600.000000}
Non Linear 195
Graphic displays of the panels for comparison with reality (grade above 600 ppm):
(fig. 3.4-7)
Mean grade calculated by DK compared to the true grades.
(fig. 3.4-8)
Scatter diagram of the DK estimates vs. the true panel values above 600 ppm
The results on tonnage look very comparable to those obtained with indicator kriging; but the
grades show a better correlation between Disjunctive kriging estimates and true values.
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
ppm
N/A
1000
950
900
850
800
750
700
650
600
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
DK_M{600.000000}
500 600 700 800 900 1000
DK M{600.000000}
500
600
700
800
900
1000
t
r
u
e

g
r
a
d
e

a
b
o
v
e

6
0
0
rho=0.753
0.0 0.5 1.0
DK T{600.000000}
0.0
0.5
1.0
t
r
u
e

t
o
n
n
a
g
e

a
b
o
v
e

6
0
0
rho=0.925
196
3.4.3 Uniform conditioning
This method aims to calculate directly the distribution of the blocks 5m x 5m within each panel, by
using the panel estimate and the anamorphosis functions to take the change of support into account.
To achieve the Uniform Conditioning we have to provide:
l the kriged 20m x 20m panel grades,
l two anamorphosis functions, one for the panel and one for the block support (Block 5m * 5m).
The calculation of the panel anamorphosis requires the value of the kriged panel dispersion vari-
ance. The two anamorphosis models must be consistent, that is, created from the same samples.
3.4.3.1 Kriging of panels
The panels have already been kriged during the in situ resource estimation (cf 3.2.4) but we need
to calculate the local dispersion variance of these estimates. In Interpolate / Estimation / (Co-)krig-
ing.:
(snap. 3.4-12)
Non Linear 197
m Set to Block mode and activate the Full set of Output Variables option
m Input: Sample set / Data / V
m Output: in Grids / Grid 20*20. Because we have asked for the Full set of Output Variables,
we are able to store the local estimated dispersion variance Variance of Z* for V under a
new variable Local dispersion Var Z*
m variogram model: Raw V
m Neighborhood: octants
Below are displayed the panel estimates:
(fig. 3.4-9)
Map of the kriged panels 20m x 20m
The Uniform Conditioning recreates a local gaussian histogram of the SMU in each panel, the mean
of this histogram being the gaussian equivalent of the kriged estimate. The panel dispersion vari-
ance (Local dispersion var Z*, estimated at the kriging step above) is also needed to re-construct
these histograms.
3.4.3.2 Uniform Conditioning
We then run Interpolate/Estimation/ Uniform Conditioning as shown below. The Block 5m * 5m
anamorphosis will be chosen for the block anamorphosis and a Tonnage correction of 0.5% will be
performed. The Basename for Output Variables is UC_no info, as the block anamorphosis has no
information effect. The same set of cut-offs as for the disjunctive kriging (21 cut-offs ranging from
0 to 1000) will be defined:
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
Kriging V
ppm
N/A
1000
900
800
700
600
500
400
300
200
100
0
198
(snap. 3.4-13)
Graphic displays of the panels for comparison with reality:
(fig. 3.4-10)
Tonnage T calculated by UC (SMU proportion) compared to the true tonnage.
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
N/A
1.000
0.875
0.750
0.625
0.500
0.375
0.250
0.125
0.000
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
UC_no info_T{600.000000}
Non Linear 199
(fig. 3.4-11)
Mean grade calculated by UC compared to the true grades.
(fig. 3.4-12)
Scatter diagram of the UC estimates vs. the true panel values above 600 ppm
The quality of local estimation is satisfying.
Moreover, UC allows to take the information effect into account by changing the block anamorpho-
sis to the Block 5*5 with information effect instead of block 5*5.
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
ppm
N/A
1000
950
900
850
800
750
700
650
600
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
UC_no info_M{600.000000}
0.0 0.5 1.0
UC_no info_T{600.000000}
0.0
0.5
1.0
t
r
u
e

t
o
n
n
a
g
e

a
b
o
v
e

6
0
0
rho=0.928
600 700 800 900 1000
UC_no info_M{600.000000}
600
700
800
900
1000
T
r
u
e

G
r
a
d
e

a
b
o
v
e

6
0
0
rho=0.785
200
Note - Some grade inconsistencies may appear when taking the information effect into account,
because the cut-off have to be applied on a histogram of kriged values. These grade inconsistencies
affect low grades for small tonnages, therefore it may be corrected by suppressing the lowest
tonnage values (as done here with a minimum tonnage fixed at 0.5%).
Do not forget to change the Basename for Output Variables to UC_with info and press RUN:
(snap. 3.4-14)
The statistical results are presented in 3.6.
In conclusion, Disjunctive kriging and Uniform Conditioning give both good results; in practice, on
real datasets, Uniform Conditioning is often preferred because it is less sensitive to stationarity
hypothesis.
3.4.3.3 Localized Uniform Conditioning
A criticism addressed to non linear techniques, including Uniform Conditioning, is that the outputs
are probability of smus grades within bigger units. We dont have a representation of the spatial dis-
tribution of smu grades, like for instance with simulations.
One way to get such a representation is to apply the Localized Uniform Conditioning methodology
(see Abzalov, M.Z. Localized Uniform Conditioning (LUC): A New Approach to Direct Modelling
of Small Blocks, Mathematical Geology 38(4) pp 393-411).
The principle is the following: the tonnage and metal at different cutoffs contained in each panel are
distributed over the smus according to a preference based on the ranking of smus kriged grade. The
metal for higher cutoff is first assigned to the smus whose kriged grades is the highest, and so on.
Non Linear 201
As there are enough data to get a realistic estimate of the kriged smus, we can apply that method to
the results of Uniform Conditioning (without information effect for instance).
As the kriging of smus has already been achieved (see 3.2.4) you just have to run Statistics / Pro-
cessing / Localized Uniform Conditioning.
Note: the same method can be used in the multivariate case, the metal of other elements are
assigned according to the ranking of the main variable kriged smus.
After Run we get the following Error message:
It is due to the fact that it is compulsory that for the highest cutoff the tonnage represents less than
the tonnage of one smu.
The solution consists in Re-running Uniform Conditioning with 41 cutoffs from 0 with a step of 50.
Then running Localized Uniform Conditioning does not produce anymore error message.
The statistics and the displays show that after Localized Uniform Conditioning the variability of
actual block grades is much better reproduced compared to the true smu grades.
202
With Tools / Grade Tonnage Curve we can also check that the QTM values obtained from Uniform
Conditioning (with Tonnage Variables option) are the same as those obtained from grades estimated
using Localized Uniform Conditioning method.
Variable Count Minimum Maximum Mean Std. Dev. Variance
True V 5x5 3120 0 1378.1 278.0 228.7 52287
KrigingV 5x5 3120 -51.06 1361.57 275.29 210.13 44153.00
LUC V 5x5 3120 0 1438.57 276.02 228.7 52374.87
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
Kriging V
V
N/A
2000
1900
1800
1700
1600
1500
1400
1300
1200
1100
1000
900
800
700
600
500
400
300
200
100
0
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
LUC V
Non Linear 203
3.4.4 Service variables
The Service Variables method is based on the transformation of grades into two variables represent-
ing the ore and metal tonnage above a given cut-off for a block centered around the data point. This
transformation requires a change of support model. Each variable is then kriged by ordinary krig-
ing. We can apply this technique for the cut-off 600 ppm (Tools / Service Variables...):
(snap. 3.4-15)
The scatter diagram between the Ore and the Metal above 600 ppm shows a very strong (non linear)
correlation.
(fig. 3.4-13)
0.0 0.5 1.0
Ore Tonnage T above 600 ppm
0
500
1000
M
e
t
a
l

Q
u
a
n
t
i
t
y

Q

a
b
o
v
e

6
0
0

p
p
m
rho=0.987
204
Consequently, we will perform independently the kriging of both variables. The experimental vari-
ograms are omnidirectional and calculated with 16 lags of 10 m (with the declustering weights
active). They have been fitted as shown below:
(fig. 3.4-14)
The declustering weights have great impact on the short scale structure; the variograms at short
scale are not satisfactory.
Then, the kriging of Ore and Metal is performed, with the usual octants neighborhood; the variables
Service Var Ore Tonnage T > 600 and Service var Metal Q > 600 are created.
91
1524
2572
3124
3696
3972
4885
5035
5319
5224
5537
5390
5578
5579
5627
5254
0 50 100 150
Distance (m)
0
10000
20000
30000
40000
50000
60000
V
a
r
i
o
g
r
a
m

:

M
e
t
a
l

Q
u
a
n
t
i
t
y

Q

a
b
o
v
e

6
0
0

Isatis
S1 - Nugget effect, Sill = 8100
S2 - Spherical - Range = 53.00m, Sill = 2.876e+004
91
1524
2572
3124
3696
3972
4885
5035
5319
5224
5537
5390
5578
5579
5627
5254
0 50 100 150
Distance (m)
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.10
V
a
r
i
o
g
r
a
m

:

O
r
e

T
o
n
n
a
g
e

T

a
b
o
v
e

6
0
0

p
p
m
Isatis
Non Linear 205
(snap. 3.4-16)
206
(snap. 3.4-17)
Because a linear kriging is performed, some panels have negative or unacceptable low Tonnage T
values: for all panels having a tonnage T < 0.02 (i.e 2%), T and Q are set to 0 (this is done using
File / Calculator...). Using the Calculator once more, we derive from the kriged variables
Service var Metal Q > 600 and Service Var Ore Tonnage T > 600 the variable Service var grade
M > 600 using the same relation M = Q / T.
Non Linear 207
(snap. 3.4-18)
208
(fig. 3.4-15)
The scatter diagrams show that some grades overestimation, and a slight under-estimation of high
tonnage values.
0.0 0.5 1.0
Service Var Ore Tonnage T>600
0.0
0.5
1.0
T
r
u
e

T
o
n
n
a
g
e

a
b
o
v
e

6
0
0
rho=0.924
600 700 800 900 1000
Service Var grade M>600
600
700
800
900
1000
T
r
u
e

G
r
a
d
e

a
b
o
v
e

6
0
0
rho=0.644
Non Linear 209
3.5 Simulations
After having reviewed the non linear estimation techniques, we can also perform simulations to
answer the same questions on the recoverable resources. Because we are in a 2D framework, we
can perform 100 simulations within a reasonable computation time.
Two techniques, both working under multigaussian hypothesis, will be described: Turning Bands
(TB) and Sequential Gaussian (SGS). This multigaussian hypothesis requires that the input variable
is gaussian: the Gaussian V variable, calculated previously ( 3.3.1 Punctual Histogram Model-
ing), will be used.
Simulations will be performed on the SMU blocks of 5 m x 5 m (Grid 5*5): this will allow to com-
pare results with the non linear estimation techniques. Therefore, block simulations require a gaus-
sian back transformation and a change of support from point to block: this implies specific remarks
discussed hereafter.
3.5.1 Before starting... important comments on block simulations
3.5.1.1 Block discretization optimization
In the standard version of Isatis, only points may be simulated and the change of support from point
to block is done by averaging simulated points. In practice, each block is discretized in n sub-cells
and each sub-cell is approximated as a point: the number n has to be large enough to validate this
approximation. But if n increases, the CPU time calculation increases, as each block will require n
simulation process (basically the CPU time is proportional to n). Thus, the choice of the block dis-
cretization is the result of a compromise between performance and precision.
The block discretization is defined through the neighborhood definition panels, and Isatis gives
some guidance to the best compromise by calculating the mean block covariance C
vv
. The block
covariance depends only on the variogram model and the block geometry. Theoretically, if n was
infinite the mean block covariance would tend to its true value.
Note - In Isatis the default block discretization is 5 x 5 and may be optimized, as explained later (
3.5.4.1).
3.5.1.2 Gaussian back transformation
When simulating in Block mode, Isatis performs automatically the following workflow:
l from the input gaussian data, simulate gaussian point grades according to the block discretiza-
tion parameters as discussed above;
l gaussian back transformation (gaussian -> raw) of the point grades using a point anamorphosis
l block grade = averaging of the raw point grades
the averaging is done automatically at the end of the simulation run. Hence the required anamor-
phosis function to perform the gaussian back transformation is the Point anamorphosis based on the
sample (point) support, which has already been calculated during the 3.3.1 Punctual Histogram
210
Modeling. The block anamorphosis Block 5m*5m (which includes a change of support correction)
should not be used here.
3.5.2 Simulations workflow summary
The aim is to simulate 5 m x 5 m block grades and to calculate the ore Tonnage T, the metal Quan-
tity Q and the mean grade M above 600 ppm for 20 m x 20 m panels. The workflow will consist in:
l Variographic analysis of the gaussian sample grades (the name of the variogram model will be
Point Gaussian V)
l Simulate the SMU grades (5 m x 5 m blocks) with Turning Bands (TB) or Sequential Gaussian
(SGS) method with the following parameters:
m Block mode
m input data: Sample Set / Data / Gaussian V
m output macro-variables to be created: Grids / Grid 5*5 / Simu V TB or Simu V SGS
m Number of simulations: 100
m Starting index: 1
m Gaussian back transformation enabled using the Point anamorphosis
m Model...: Point Gaussian V defined at the previous step
m Seed for Random Number Generator: leave the default number 423141. This seed is sup-
posed to be a large prime number; the same seed allows reproducibility of realizations.
The neighborhood and other parameters specific to each method will be detailed in the relevant
paragraph.
l Calculation of the QTM variables for both techniques (described for TB): ore Tonnage T (i.e
SMU proportion within each panel), metal Quantity Q, and mean grade M of blocks above 600
ppm among each 20 m x 20 m panel (M = Q / T). The panel mean grades can not be averaged
directly on the 100 simulations: the mean grade is not additive because it refers to different ton-
nages (the tonnage may differ between different simulations). Therefore it has to be weighed by
the ore proportion T. One way to do this is to use an accumulation variable for each panel:
m calculate the ore proportion T and the metal quantity Q (the metal quantity is the accumula-
tion variable: Q = T x M) for each simulation
m calculate the average (T) and average (Q) of the 100 simulations
m calculate the average mean grade: average (M) = average (Q) / average (T)
3.5.3 Variographic analysis of gaussian sample grades
The experimental variogram of gaussian variables often show more visible structures and make
their interpretation easier; the analysis of anisotropy using the variogram map gives similar infor-
mation about the main directions of continuity. In Statistics / Exploratory Data Analysis..., the
experimental variogram Point Gaussian V is calculated with the same rotation parameters than
Non Linear 211
Raw V. A variogram model using 3 structures has been fitted and saved under the name Point
Gaussian V:
(fig. 3.5-1)
3.5.4 Simulation with the Turning Bands method
3.5.4.1 Simulations
We run Interpolate / Conditional Simulations / Turning Bands... with the parameters defined in the
workflow summary ( 3.5.2):
N160
N250
0 50 100 150
Distance (m)
0.00
0.25
0.50
0.75
1.00
1.25
V
a
r
i
o
g
r
a
m

:

G
a
u
s
s
i
a
n

V
Isatis
212
(snap. 3.5-1)
l Gaussian back transformation... enabled: select the Point anamorphosis.
l Neighborhood...: create a new neighborhood parameter file named octants for TB. Press Edit...
and from the Load... button reload the parameters from the octants neighborhood. We are now
going to optimize the block discretization: press the ... button next to Block Discretization: the
Discretization Parameters window pops up where the number of discretization points along the
x,y,z directions may be defined. These numbers are set to their default value (5 x 5 x 1). Press
Calculate C
vv
, the following appears in the message window (values differ at each run due to the
randomization process):
Non Linear 213

Regular discretization: 5 x 5 x 1
In order to account for the randomization, 11 trials are performed
(the first value will be kept for the Kriging step)
Variables Gaussian V
Cvv = 0.811792
Cvv = 0.809978
Cvv = 0.812136
Cvv = 0.811752
Cvv = 0.810842
Cvv = 0.812900
Cvv = 0.808768
Cvv = 0.811977
Cvv = 0.810781
Cvv = 0.810921
Cvv = 0.812400
11 mean block covariances have been calculated with 11 different randomizations. The mini-
mum value is 0.808768 and the maximum is 0.812900; the maximum relative variability is
approximately 0.5% which is more than acceptable: the 5 x 5 discretization is a very good
approximation of the punctual support and may be optimized.
Note - Note that, for reproducibility purposes, the first value of C
vv
will be kept for the simulations
calculation
For optimization, we decrease the number of discretization points to 3x3:
214
(snap. 3.5-2)
Press Calculate C
vv
:

Regular discretization: 3 x 3 x 1
In order to account for the randomization, 11 trials are performed
(the first value will be kept for the Kriging step)
Variables Gaussian V
Cvv = 0.809870
Cvv = 0.814197
Cvv = 0.808329
Cvv = 0.812451
Cvv = 0.819093
Cvv = 0.809922
Cvv = 0.814171
Cvv = 0.811332
Cvv = 0.805993
Cvv = 0.806053
Cvv = 0.807459
Non Linear 215
The minimum value is 0.805993 and the maximum value is 0.819093: the maximum relative
variability is approximately 1.6%. As expected, it has increased but remains acceptable: there-
fore, the 3 x 3 discretization is a good compromise and will be kept for the simulations (i.e each
simulated block value will be the average of 3 x 3 = 9 simulated points). Press Close then OK
for the neighborhood definition window.
l Number of Turning Bands: 300. The more turning bands, the more precise are the realizations
but CPU time increases. Too few turning bands would create visible 1D-line artefacts.
Press RUN: calculations may take a few minutes.
We represent in the next figure five simulations, compared to the true map:
216
(fig. 3.5-2)
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
True V
ppm
N/A
1000
900
800
700
600
500
400
300
200
100
0
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
Simu V TB[00020]
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
Simu V TB[00030]
50 100 150 200 250
X ( )
50
100
150
200
250
300
Y

(
m
)
Simu V TB[00040]
50 100 150 200 250
50
100
150
200
250
300
Y

(
m
)
Simu V TB[00050]
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
Simu V TB[00002]
Non Linear 217
3.5.4.2 Calculation of the QTM variables
From Statistics/Processing/Grade simulation Post-processing compute the metal quantity, mean
grade and tonnage on the 20*20 grid from the 5*5 grid simulation.
(snap. 3.5-3)
We can then display the ore Tonnage T and mean grade M above 600 ppm calculated by Turning
Bands and compere them to the true values:
218
(eq. 3.5-1)
Tonnage T calculated by TB (SMU proportion) compared to the true tonnage. The color scale is
a regular 16-class grey palette between 0 and 1: panels for which there is strictly less than 1
block (i.e 0 <= proportion < 0.0625) are white.
(fig. 3.5-3)
Mean grade calculated by TB compared to the true grades.
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
N/A
1.000
0.875
0.750
0.625
0.500
0.375
0.250
0.125
0.000
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
TB_ mean ore tonnage above 600
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
ppm
N/A
1000
950
900
850
800
750
700
650
600
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
TB_mean (mean grade above 600)
Non Linear 219
(fig. 3.5-4)
Scatter diagrams of ore tonnage and mean grade above 600 ppm between
the mean of 100 TB simulations and the true values of panels.
3.5.5 Simulation with the Sequential Gaussian method
Two different algorithms are available for SGS in Isatis, using two different kinds of neighborhood:
l Interpolate / Conditional Simulation / Sequential Gaussian / Standard Neighborhood...: a stan-
dard elliptical neighborhood is used taking the point data & the previously simulated grid nodes
into account.
l Interpolate / Conditional Simulation / Sequential Gaussian / Sequential Neighborhood...: the
sequential neighborhood performs first a migration of point data on the nearest grid node; the
neighborhood is then defined by a moving window made of x blocks around the target block.
We will use the standard neighborhood option because it is more accurate from a theoretical point
of view, and moreover the Block simulation is possible (automatic averaging of point values).
3.5.5.4 Simulations
Open Interpolate / Conditional Simulations / Sequential Gaussian / Standard neighborhood.... and
enter the same parameters described in the workflow summary ( 3.5.2):
600 700 800 900 1000
TB_mean (mean grade above 600)
600
700
800
900
1000
T
r
u
e

G
r
a
d
e

a
b
o
v
e

6
0
0
rho=0.869
0.0 0.5 1.0
TB_mean ore tonnage above 600
0.0
0.5
1.0
T
r
u
e

T
o
n
n
a
g
e

a
b
o
v
e

6
0
0
rho=0.936
220
(snap. 3.5-4)
Non Linear 221
l The Gaussian Back Transformation is enabled with the Point anamorphosis function
l Special Model Options...: by default, a Simple Kriging (SK) is performed using a constant
mean equal to zero
l Neighborhood...: create a new neighborhood named octants for SGS with the following param-
eters (you may load the parameters from the octants for TB parameter file):
(snap. 3.5-5)
m The search ellipsoid is maintained to 70 m.
m minimum number of samples: 5
m Optimum Number of Samples per Sector: 4, which adds to a maximum of 32 samples. The-
oretically, the SGS technique would require a unique neighborhood and use all the previ-
ously simulated grid nodes to reproduce exactly the variogram; in practice, it is impossible,
so it is recommended to increase the Optimum Number in respect to the Optimum Number of
Already Simulated Node (to be defined below in the main SGS window) and the capacity of
the computer.
222
m in the Advanced tab, set the Minimum distance between two samples to 2 m; as two different
sets of data are used to condition the simulations (i.e the actual data points combined with
the previously simulated grid nodes), this minimum distance criterion avoids fictitious
duplicates between original data points and simulated grid nodes. It allows also to spread
conditioning data for a better reproducibility of the variogram.
m The same Block Discretization of 3 x 3 will be used.
l Optimum Number of Already Simulated Node: 16. This means that the software will load all the
real samples and the 16 closest already simulated nodes in memory for the search neighborhood
algorithm. The maximum number of samples being 32, there will be 16 real samples used for
each node simulation, as for the Turning Bands method. The TEST window allows to evaluate
the impact of these different parameters on the neighborhood.
l Leave the other parameters to their default values and press RUN
Note - Isatis offers the possibility to perform the different simulations with independent paths
(optional toggle in the main SGS window). By default, this toggle is set OFF, meaning that the same
random path is used for all simulations: the independency is no more certain, but the algorithm is
much quicker. If the toggle is set ON, the CPU time will approximately be multiplied by the number
of simulations. Here, it has been checked that both options show negligible differences in the final
results.
The resulting outcomes are very similar to the TB method.
Non Linear 223
3.5.5.5 Calculation of the QTM variables
From Statistics/Processing/Grade simulation Post-processing compute the metal quantity, mean
grade and tonnage on the 20*20 grid from the 5*5 grid simulation.
(snap. 3.5-6)
224
(fig. 3.5-5)
Tonnage T calculated by SGS (SMU proportion) compared to the true tonnage.
(fig. 3.5-6)
Mean grade calculated by SGS compared to the true grades.
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
N/A
1.000
0.875
0.750
0.625
0.500
0.375
0.250
0.125
0.000
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
SGS_ mean ore tonnage above 600
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
ppm
N/A
1000
950
900
850
800
750
700
650
600
50 100 150 200 250
X (m)
50
100
150
200
250
300
Y

(
m
)
SGS_ mean (mean grade above 600)
Non Linear 225
(fig. 3.5-7)
Scatter diagrams of ore tonnage and mean grade above 600 ppm between
the mean of 100 SGS and the true values of panels
We observe that SGS simulations give very similar results to TB and are also well correlated to the
reality.
0.0 0.5 1.0
SGS_ mean ore tonnage above 600
0.0
0.5
1.0
t
r
u
e

t
o
n
n
a
g
e

a
b
o
v
e

6
0
0
rho=0.938
600 700 800 900 1000
SGS_ mean (mean grade above 600)
600
700
800
900
1000
t
r
u
e

g
r
a
d
e

a
b
o
v
e

6
0
0
rho=0.870
226
3.6 Conclusions
The objective of the case study was to illustrate several non linear methods (global and local) to
estimate recoverable resources, and compare them to linear kriging. All methods take the same sup-
port effect for 5 m x 5 m blocks into account, but only a few take the information effect into
account. Therefore, we will first focus on results without information effect.
3.6.1 Global estimation
3.6.1.1 Without information effect
Grade Tonnage curves
The following methods will be compared to the true values (True): Ordinary Kriging (OK), block
anamorphosis (block 5x5), Indicator Kriging (IK), Disjunctive Kriging (DK) and Uniform Condi-
tioning (UC). The grade-tonnage curves for all these methods will be presented; Service Variables
(SV) and simulations (TB and SGS) have been calculated only for one particular cut-off V = 600
ppm so we can not display G-T curves for these methods.
Open Tools / Grade Tonnage Curves...: Activate 6 curves. For IK, DK and UC outcomes, we need
to ask for Tonnage Variables. For instance, for the Indicator Kriging (IK): press Edit..., choose the
Tonnage Variables option then IK_Q[xxxxx] for the Metal Quantity and IK_T[xxxxx] for the Total
Tonnage:
Non Linear 227
(snap. 3.6-1)
Repeat the same for DK and UC, and change the curve parameters and labels for optimal visibility.
By clicking on the graphic windows below, ask for the following Grade Tonnage curves: Mean
grade vs. cut-off, Total tonnage vs. cut-off, Metal tonnage vs. cut-off and Metal tonnage vs. Total
tonnage. The graphics are presented here below:
228
(fig. 3.6-1)
Mean Grade vs. Cutoff
(fig. 3.6-2)
Total Tonnage vs. Cutoff
0 250 500 750 1000
Cutoff
0
250
500
750
1000
1250
M
e
a
n

G
r
a
d
e
True
OK
Block 5*5
IK
DK
UC
0 250 500 750 1000
Cutoff
0
10
20
30
40
50
60
70
80
90
100
T
o
t
a
l

T
o
n
n
a
g
e
True
OK
Block 5*5
IK
DK
UC
Non Linear 229
(fig. 3.6-3)
Metal Tonnage vs. Cutoff
(fig. 3.6-4)
Metal Tonnage vs. Total Tonnage
0 250 500 750 1000
Cutoff
0
50
100
150
200
250
M
e
t
a
l

T
o
n
n
a
g
e
True
OK
Block 5*5
IK
DK
UC
0 10 20 30 40 50 60 70 80 90 10
Total Tonnage
0
50
100
150
200
250
M
e
t
a
l

T
o
n
n
a
g
e
True
OK
Block 5*5
IK
DK
UC
230
The True curve is black and represented with a bold line type. We clearly see that the OK tonnage
curves are shifted compared to others: linear kriging induces a significant smoothing effect despite
a refined sampling and a good coverage of the domain.
All non linear methods provide similar and suitable results; a zoom centered on V = 600 allows
a more precise comparison around this particular cut-off:
(fig. 3.6-5)
Grade-Tonnage curves with a zoom on the 600 ppm cutoff of interest (same legend)
Little differences are noticeable: IK overestimates the grades whereas DK overestimates the ton-
nages.
570 580 590 600 610 620 630 64
Cutoff
720
730
740
750
760
770
780
790
800
M
e
a
n

G
r
a
d
e
570 580 590 600 610 620 630
Cutoff
8
9
10
11
12
13
T
o
t
a
l

T
o
n
n
a
g
e
525 550 575 600 625 650
Cutoff
60
70
80
90
M
e
t
a
l

T
o
n
n
a
g
e
8.5 9.0 9.5 10.0 10.5 11.0 11.
Total Tonnage
73
74
75
76
77
78
79
80
81
M
e
t
a
l

T
o
n
n
a
g
e
True
OK
Block 5*5
IK
DK
UC
Non Linear 231
As we had to choose a particular cut-off for comparing these methods with SV and simulations, we
have chosen V = 600 and the global results according to this cut-off are presented hereafter.
Global statistics on cut-off V = 600 ppm
The following tables give the statistics on ore tonnage, metal quantity and grade above 600 for the
different methods on the 195 panels. The true values are compared to the following methods (using
Statistics / Quick Statistics...): Turning Bands (TB), Sequential Gaussian Simulations (SGS), Indi-
cator Kriging (IK), Disjunctive Kriging (DK), Uniform conditioning (UC), Service Variables (SV),
global estimation with support effect (Block 5x5 without information effect, results already shown
in 3.3.4 Analysis of the results for the global estimation p.94) and ordinary kriging (OK):
Statistics on Ore Tonnage above 600 (proportion)
Statistics on Metal Quantity above 600
As the Mean grade M defined on the panels refers to different tonnages, it is not additive so the cal-
culation of the mean and the standard deviation needs to be weighed by the tonnages. Therefore,
VARIABLE Count Minimum Maximum Mean Std. Dev. Variance
True 195 0 1 0,104 0,21 0,044
TB 195 0 0,994 0,104 0,185 0,034
SGS 195 0 0,996 0,1 0,185 0,034
IK 195 0 0,99 0,101 0,2 0,04
DK 195 0 1 0,115 0,192 0,037
UC 195 0 0,987 0,096 0,189 0,036
SV 195 0 0,884 0,098 0,163 0,027
OK 195 0 1 0,081 0,205 0,042
Block 5x5 0,101
True 195 0 997,8 78,0 169,7 28796,8
TB 195 0 996,5 79,4 155,4 24148,8
SGS 195 0 1002,4 75,9 156,4 24453,6
IK 195 0 982,3 78,2 165,9 27533,9
DK 195 0 1002,5 86,0 159,7 25488,7
UC 195 0 1004,6 71,7 154,6 23893,9
SV 195 0 804,2 74,3 131,9 17231,0
OK 195 0 1005,7 61,1 165,6 27418,7
Block 5x5 76,0
232
use Statistics / Quick statistics 8 times on each grade variable of each method with the relevant ton-
nage as the Weight variable:
Statistics on Mean Grade above 600
These statistics are attached to the specific cut-off 600: no global conclusion on the performances of
the methods may be assessed here. Besides, the dataset may not be compared to a realistic explora-
tion campaign.
3.6.1.2 With information effect
Comparisons will be made for the anamorphosis Block 5*5 with information effect and the Uni-
form Conditioning (UC_with info[xxxxx]). Results for the block anamorphosis have already been
discussed (cf. 3.3.4 Analysis of the results for the global estimation p.94). Only global statistics
for the cut-off V = 600 ppm have been made:
| Q | T | M
True block 5x5 | 77.95 | 10.38 | 750.67
True block 5x5 (info) | 67.92 | 9.01 | 754.11
Block 5*5 with info | 72.03 | 9.69 | 743.05
UC_with info | 69.20 | 9.17 | 754.60
For the cut-off V = 600 ppm, UC has correctly quantified the information effect.
3.6.2 Local estimation
For each local estimation method, a scatter diagram of the panel estimates with true values (ton-
nages and grades) with the correlation coefficients has already been done (cf. relevant paragraphs).
Here, the error for each panel has been calculated and reported:
error = estimate - true value
Therefore, positive error values reveal overestimation.
True 66 603,0 997,8 700,0 79,0 9225,2
TB 173 600,3 1009,7 689,0 52,0 2706,7
SGS 166 606,2 1015,7 684,0 53,0 2815,5
IK 91 604,5 992,2 722,0 87,0 7560,3
DK 116 132,6 1002,5 659,5 118,6 14066,0
UC 115 653,4 1017,2 691,5 50,0 2496,2
SV 103 607,3 951,3 735,0 63,3 4000,5
OK 44 601,0 1005,7 756,7 102,0 10411,6
Block5x5 754,5
Non Linear 233
The table below summarizes the main results for the error on tonnages:
Local statistics of error on tonnages estimates and correlation
with true tonnage values (for cut-off = 600 ppm)
The true global tonnage is 0.104; the bias for all non linear methods remains acceptable.
The table below summarizes the main results for the error on mean grades above 600:
Local statistics of error on mean grades above 600 and correlation
with true values (for cut-off = 600 ppm)
IK and SV methods show a global overestimation of the grades and a lower correlation with reality.
VARIABLE Count Minimum Maximum Mean Std. Dev. Correlation
TB 195 -0,397 0,248 -0,001 0,075 0,94
SGS 195 -0,404 0,237 -0,004 0,074 0,94
IK 195 -0,395 0,39 -0,003 0,089 0,91
DK 195 -0,399 0,363 0,01 0,08 0,93
UC 195 -0,314 0,219 -0,009 0,079 0,93
SV 195 -0,411 0,277 -0,006 0,086 0,93
OK 195 -0,5 0,375 -0,023 0,085 0,92
ID2 195 -0,563 0,188 -0,024 0,08 0,93
TB 66 -88,8 98,6 18,1 39,0 0,87
SGS 66 -81,9 98,3 12,7 38,9 0,87
IK 57 -93,4 275,3 33,8 69,2 0,68
DK 65 -485,6 161,2 0,8 82,9 0,75
UC 65 -126,4 94,0 7,9 49,3 0,79
SV 65 -113,1 188,2 41,6 61,2 0,67
OK 40 -100,5 67,9 -22,6 40,6 0,89
ID2 44 -130,9 134,6 -21,0 50,1 0,82
234
The table below summarizes the main results for metal quantity:
Local statistics of error on metal quantity and correlation
with true values (for cut-off = 600 ppm)
All non linear methods give consistent results for the metal quantity.
3.6.3 Final conclusions
The conclusions based on these numerical results only concern this particular dataset and should
not be interpreted as a straightforward classification of the methods.
Despite a refined sampling, linear interpolation methods (linear kriging, inverse distance...) induce
a smoothing effect that has a significant impact on recoverable resources. Non linear geostatistics
provide practical solutions and this case study shows that all methods are globally consistent;
though some little differences appear at the local scale.
Global estimation techniques, based on anamorphosis functions, showed satisfying results and are
quick to proceed.
Simulations techniques (TB and SGS) showed good results but these techniques are time consum-
ing and quite heavy to proceed. Indicator Kriging showed some little differences at the local scale
(as service variables), and requires some specific pre/post-processing. Disjunctive Kriging and Uni-
form Conditioning both make use of anamorphosis functions, but Uniform Conditioning has the
advantage to base itself on ordinary kriging estimates instead of the global mean for Disjunctive
Kriging, which requires a stronger stationarity hypothesis. Besides, Uniform Conditioning is
straightforward to the global estimation techniques and allows to take the information effect into
account.
TB 195 -276,5 174,9 -1,4 51,9 0,95
SGS 195 -281,1 166,8 -2,1 51,3 0,95
IK 195 -266,8 260,5 0,2 59,1 0,94
DK 195 -277,0 253,1 8,0 55,8 0,94
UC 195 -213,7 153,2 -6,2 54,8 0,95
SV 195 -279,0 192,2 -3,6 62,9 0,94
OK 195 -350,3 242,3 -16,8 57,5 0,94
ID2 195 -389,0 120,8 -18,5 54,7 0,95

Isatis Case Studies Mining

Uploaded by

Copyright:

Available Formats

Isatis Case Studies Mining

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Isatis Case Studies Mining

Uploaded by

Copyright:

Available Formats

ISATIS 2011

Mining Case Studies

You might also like