Stoch Environ Res Risk Assess

DOI 10.1007/s00477-014-0997-x

ORIGINAL PAPER

Simulation of geo-domains accounting for chronology and contact

relationships: application to the Rı́o Blanco copper deposit
Nasser Madani • Xavier Emery

 Springer-Verlag Berlin Heidelberg 2014

Abstract The plurigaussian model is increasingly used for 1 Introduction

simulating geo-domains and quantifying geological uncer-
tainty in the subsurface. However, because they rely on the The spatial modeling of geological domains (hereafter
truncation of only two Gaussian random fields, the current named geo-domains) is important to characterize rock
implementations of this model are often restricted in the properties and geological heterogeneities in the subsurface.
number of geo-domains that can be simulated and in their Examples of applications include groundwater hydrology,
contact relationships. A solution to overcome these restric- mining and petroleum exploration, in which it is of interest
tions is to increase the number of underlying Gaussian random to describe heterogeneous aquifers, deposits and reservoirs
fields. Such an approach yields a very flexible model, able to at multiple scales, in order to predict the distribution of
reproduce the contact relationships between geo-domains in in situ resources and to manage subsurface activities (e.g.,
agreement with their chronology (i.e., such that younger geo- forecast of production performance and remediation pro-
domains cross-cut the older ones), as well as the geo-domain cedures) (De Marsily et al. 2005; Journel and Huijbregts
proportions and spatial correlation structure. The proposed 1978; Deutsch 2002).
approach is applied to a dataset from the Rı́o Blanco copper Since the boundaries of the geo-domains are generally
deposit in the Chilean Central Andes, in which it is of interest uncertain, stochastic spatial models are increasingly used in
to simulate the layout of seven rock units (andesite, granitoid, place of, or in complement to, deterministic models. They
tourmaline breccia, monolithic breccia, magmatic breccia, rely on the use of simulation techniques that no longer
porphyry and pipe). The results are used to map the proba- provide a unique layout of the geo-domains, but a set of
bilities of occurrence of the rock units and to identify the alternative scenarios that mimic the actual distribution of
sectors where the interpreted rock model is uncertain. the geo-domains in space (Chilès and Delfiner 2012).
To date, a variety of simulation approaches have been
Keywords Geological uncertainty  Geological developed for modeling spatial geo-domains, including
heterogeneity  Domaining  Plurigaussian simulation  object-based models (Lantuéjoul 2002), sequential indicator
Hierarchical modeling simulation (Journel and Alabert 1990), simulation based on
multiple-point statistics (Strebelle 2002; Rezaee et al. 2014),
truncated Gaussian simulation (Matheron et al. 1987, 1988)
N. Madani  X. Emery (&)
and plurigaussian simulation (Armstrong et al. 2011). The
Department of Mining Engineering and Advanced Mining
Technology Center, University of Chile, Avenida Tupper 2069, latter offers a flexible framework, as it allows incorporating
8370451 Santiago, Chile topological constraints (contact relationships) on the simu-
e-mail: xemery@ing.uchile.cl lated geo-domains, in addition to reproducing the propor-
N. Madani tions and the spatial correlation structure of these domains. In
e-mail: nasser.madani@ing.uchile.cl the past decade, it has found wide acceptance in the field of
natural resources evaluation, in particular, for modeling
N. Madani
CSIRO-Chile International Center of Excellence in Mining and aquifers, petroleum reservoirs and ore bodies of varied nat-
Mineral Processing, Santiago, Chile ures such as porphyry, epithermal, granite-hosted, carbonate

123
 Stoch Environ Res Risk Assess

hosted, roll front and lateritic deposits (Skvortsova et al. where C(h) is the auto-covariance function of the standard
2001, 2002; Deraisme and Field 2006; Fontaine and Beucher Gaussian random field, G is the standard normal cumulative
2006; Carrasco et al. 2007; Emery and González 2007a, b; distribution function, and ct(h) is the variogram of the indi-
Emery et al. 2008; Mariethoz et al. 2009; Rondon 2009; cator obtained by truncating the Gaussian random field at
Yunsel and Ersoy 2011, 2013; Jeannée et al. 2013; Talebi threshold t. The calculation of the integral in Eq. (1) can be
et al. 2013, 2014; Deutsch and Deutsch 2014). In most of performed by numerical integration or by using expansions
these studies, if not all, it calls the attention that the simu- into series (Kyriakidis et al. 1999; Emery 2007). In practice,
lation is restricted to a few geo-domains (no more than five), there is information about the geo-domain indicators (0 or 1)
which can be explained because of the difficulty in imposing at sampling locations, which allows calculating the indicator
specific contact relationships with a higher number of geo- variograms and fitting a suitable model for the covariance
domains. function of the Gaussian random field.
The aim of this paper is twofold. First, from a methodo- Once the model is specified, the simulation conditional
logical point of view, it is of interest to present an extension to the sampling information can be performed through the
of the currently used plurigaussian model, in order to following steps (Lantuéjoul 2002):
increase its versatility and to allow modeling numerous geo-
(1) Simulation of the Gaussian random field at the
domains with possible complex contact relationships. Sec-
locations where the geo-domains are known. To this
ond, our aim is to apply this model to the simulation of seven
end, an iterative procedure (Gibbs sampler) is used,
rock type domains in the Rı́o Blanco deposit, a giant por-
which allows converting the indicator data (0 or 1)
phyry copper deposit located in the Central Andes of Chile.
into Gaussian data.
(2) Simulation of the Gaussian random field at the target
locations. At this stage, many algorithms can be
2 Methodology
used: turning bands, sequential, spectral, etc. (Lant-
uéjoul 2002; Chilès and Delfiner 2012).
2.1 Truncated Gaussian simulation
(3) Truncation in order to obtain a simulation of the two
geo-domains (geo-domain 1 when the simulated
The truncated Gaussian model was coined in the late 1980s
Gaussian random field is smaller than the truncation
(Matheron et al. 1987, 1988) for simulating geo-domains
threshold t, geo-domain 2 otherwise).
(lithofacies) in oil reservoirs. It relies on a Gaussian ran-
dom field (usually assumed to be second-order stationary,
with zero mean and unit variance) and a set of truncation 2.2 Plurigaussian simulation
thresholds, defining each geo-domain as the set of points in
space where the value of the Gaussian random field The plurigaussian model (Galli et al. 1994; Le Loc’h and Galli
belongs to an interval defined with these thresholds. For 1997) is a generalization of the truncated Gaussian model,
example, to simulate two geo-domains D1 and D2, one aimed at producing a wider range of patterns and allowing for
would define a unique threshold t and associate D1 with the more complicated types of contacts between geo-domains.
points of space where the Gaussian random field is below t, The basic idea is to use several Gaussian random fields and to
and D2 with the points where this field is above t (Fig. 1). define the geo-domains by a multivariate truncation rule. For
Likewise, to simulate N geo-domains, one has to define instance, a three-domain truncation rule can be defined on the
N - 1 thresholds. Once the ordering of the geo-domains basis of two Gaussian random fields and two thresholds: if the
has been chosen, there is a one-to-one relationship between first Gaussian random field is less than the first threshold, then
the thresholds and the proportion of space covered by each one finds the first geo-domain; otherwise, if the second
geo-domain (Armstrong et al. 2011). Gaussian random field is less than the second threshold, one
The auto-covariance function C(h) of the underlying finds the second geo-domain; otherwise, one finds the third
Gaussian random field is related to the spatial correlation geo-domain. Such a truncation rule can be represented
structure of the geo-domains, which can be represented by through a two-dimensional flag, in which each axis represents
the variogram of the geo-domain indicators (Armstrong one Gaussian random field (Fig. 2c). In practice, to restrict the
et al. 2011; Chilès and Delfiner 2012). For example, in the number of parameters and to ease their inference, one usually
case of two geo-domains defined by a single truncation works with two Gaussian random fields, which are further-
threshold t, one has: more assumed to be independent (such an independence
ct ðhÞ ¼ GðtÞ ½1  GðtÞ assumption is sometimes not made, e.g., Dowd et al. 2003;
Z arcsin ½CðhÞ   Armstrong et al. 2011). With this simplification, the plurig-
1 t2
 exp  dh; ð1Þ aussian model requires defining:
2p 0 1 þ sin h

123
Stoch Environ Res Risk Assess

Fig. 1 Simulated Gaussian random field (left) and associated geo-domains (right)

• the truncation rule: it has an impact on the contact (2) One can also control the chronology relationships
relationships between geo-domains (which domains can between geo-domains. For instance, in Fig. 2c, the
be in contact and which domains cannot); first geo-domain (1) appears in the foreground (first
• the truncation thresholds: they have an impact on the layer) and can then be interpreted as a younger domain
geo-domain proportions; that cross-cuts the other two older geo-domains (2 and
• the auto-covariance functions of the Gaussian random 3), which appear in the background (last layer).
fields: they have an impact on the variograms of the
However, by restricting to a two-dimensional flag (as in
geo-domain indicators.
most applications), one is often limited in the number of
One of the current limitations of the plurigaussian model geo-domains that can be modeled, especially when all
is the design of a two-dimensional flag representing the these geo-domains are mutually in contact, a situation that
truncation rule. The choice of this rule has implications on often arises in practice. For instance, Fig. 3 presents a few
the spatial relationships between geo-domains, in particu- examples of flags with four geo-domains (modified from
lar, with respect to permissible and forbidden contacts. For Armstrong et al. 2011): none of them allows each geo-
instance, Fig. 2c shows that the three geo-domains can domain to be in contact with each other, as there is always
touch each other, whereas Fig. 2d shows two geo-domains some prohibited contact between geo-domains.
that cannot have contact together. Most often, the flag To overcome this limitation, following Xu et al. (2006) and
consists of a partition of the two-dimensional space into Emery (2007), the key idea of this work is to increase the
rectangles (Armstrong et al. 2011) or unions of rectangles number of underlying Gaussian random fields, so as to increase
(Emery 2007), although more general designs can be the dimensionality of the flag used to represent the truncation
considered (Allard et al. 2012; Deutsch and Deutsch 2014). rule and to allow the geo-domains to be in contact altogether.
As an example, consider the following four-domain truncation
rule, based on three underlying Gaussian random fields {Y1, Y2,
2.3 Proposal: hierarchical approach Y3} and three truncation thresholds {t1, t2, t3}, which can be
represented by a three-dimensional flag (Fig. 4):
2.3.1 Principle Geo-domain at location x
8
> 1 if Y1 ðxÞ  t1
Defining a truncation rule through a two-dimensional flag >
>
< 2 if Y ðxÞ [ t and Y ðxÞ  t
1 1 2 2
offers the following advantages: ¼
>
> 3 if Y1 ðxÞ [ t1 ; Y2 ðxÞ [ t2 and Y3 ðxÞ  t3
(1) One can visually control the contact relationships >
:
4 if Y1 ðxÞ [ t1 ; Y2 ðxÞ [ t2 and Y3 ðxÞ [ t3 ð2Þ
between geo-domains: two geo-domains are in
contact in space if they are in contact on the flag A realization of this model is shown in Fig. 5. The first
representing the truncation rule. geo-domain (painted in blue) appears in the foreground or

123
 Stoch Environ Res Risk Assess

Fig. 2 Simulated Gaussian

random fields (top) and
associated geo-domains
(bottom)

first ‘‘layer’’ (it may represent a younger domain that cross- geo-domains that are numbered chronologically is pre-
cuts the other ones). The second geo-domain (red) is in a sented in Fig. 6 (geo-domain 1 is the youngest and cross-
second layer, while geo-domains 3 and 4 (purple and cuts the other ones, while geo-domains 6 and 7 are the
green) are located in the background (third layer). oldest). It is based on a truncation rule that involves six
Such a hierarchical approach can easily be generalized Gaussian random fields {Y1, Y2, Y3, Y4, Y5, Y6} and six
to more than four geo-domains. An example with seven thresholds {t1, t2, t3, t4, t5, t6}:

8
>
> 1 if Y1 ðxÞ  t1
>
> 2 if Y1 ðxÞ [ t1 and Y2 ðxÞ  t2
>
>
>
>
<3 if Y1 ðxÞ [ t1 ; Y2 ðxÞ [ t2 and Y3 ðxÞ  t3
Geo-domain at location x ¼ 4 if Y1 ðxÞ [ t1 ; Y2 ðxÞ [ t2 ; Y3 ðxÞ [ t3 and Y4 ðxÞ  t4 ð3Þ
>
>
>5
> if Y1 ðxÞ [ t1 ; Y2 ðxÞ [ t2 ; Y3 ðxÞ [ t3 ; Y4 ðxÞ [ t4 and Y5 ðxÞ  t5
>
>
>
> 6 if Y1 ðxÞ [ t1 ; Y2 ðxÞ [ t2 ; Y3 ðxÞ [ t3 ; Y4 ðxÞ [ t4 ; Y5 ðxÞ [ t5 and Y6 ðxÞ  t6
:
7 if Y1 ðxÞ [ t1 ; Y2 ðxÞ [ t2 ; Y3 ðxÞ [ t3 ; Y4 ðxÞ [ t4 ; Y5 ðxÞ [ t5 and Y6 ðxÞ [ t6

123
Stoch Environ Res Risk Assess

Fig. 3 Examples of truncation

rules with four geo-domains
(modified from Armstrong et al.
2011). The geo-domains
represented in yellow and green
are never in contact

To ease inference, in the following, the Gaussian ran- under consideration, while the design of the truncation rule
dom fields are assumed to be mutually independent and is the same at every location. An example of such a spa-
second-order stationary, with zero mean and unit tially-varying proportion model will be presented in
variance. Sect. 3.

2.3.2 Inference of truncation thresholds 2.3.3 Variogram analysis

The previous truncation rule depends on six scalar As mentioned in Sect. 2.1, the auto-covariance functions
thresholds {t1, t2, t3, t4, t5, t6}, which can be defined in (equivalently, the variograms) of the underlying Gaussian
order to reproduce the proportions of the different geo- random fields are related to the indicator variograms
domains (Armstrong et al. 2011). For instance, the pro- associated with the geo-domains (Eq. (1)). In Eq. (1), it is
portion of the first geo-domain is equal to G-1(t1), where seen that, when the value of C(h) increases from -1 to 1,
G stands for the standard normal cumulative distribution the value of the integral increases, insofar as the integrand
function: threshold t1 can then be defined in order to yield a is positive and the upper bound of the integral increases, so
given proportion for geo-domain 1 (for instance, the pro- that the value of the indicator variogram ct(h) decreases.
portion inferred from an interpreted geological model). Accordingly, there exists a one-to-one relationship between
Likewise, since Y1 and Y2 are independent, the proportion C(h) and ct(h), and one can express the former as a func-
of the second geo-domain is [1 - G-1(t1)] 9 G-1(t2): as t1 tion of the latter (Kyriakidis et al. 1999):
is already known, threshold t2 can now be defined in order
CðhÞ ¼ ut ðct ðhÞÞ; ð4Þ
to yield a given proportion for geo-domain 2. The same
procedure can be applied for the remaining thresholds with ut a function that depends on the chosen truncation
(t3–t6), which relate to the proportions of geo-domains 3–7. threshold t.
Spatial variations in the geo-domain proportions can be In the proposed hierarchical model, the geo-domain data
imposed by using spatially varying thresholds, i.e., by can be converted into indicator data (0, 1, or unknown value)
considering that t1, t2, t3, t4, t5, t6 depend on the location on whether or not the Gaussian random fields are below their

123
 Stoch Environ Res Risk Assess

Fig. 4 Three-dimensional flag

representing a truncation rule
with four geo-domains that are
mutually in contact. Each axis
of this flag corresponds to an
underlying Gaussian random
field

Fig. 6 Plurigaussian realization with seven geo-domains distributed

Fig. 5 Plurigaussian realization with four geo-domains distributed in in six successive layers
three successive layers
Sects. 2.3.1–2.3.3 improves these proposals in two aspects.
The first one is the definition of spatially varying truncation
respective thresholds (see Table 1 for an illustration to the thresholds in order to account for spatially varying pro-
truncation rule presented in Eq. (3)). The indicator vario- portions of the geo-domains. The second one relates to
grams can therefore be experimentally calculated and, using variogram analysis. In former proposals (except for the
Eq. (4), converted into covariance values (equivalently, into simple case of truncated Gaussian simulation, based on a
variograms) for the underlying Gaussian random fields. single Gaussian random field), the variograms of the
These experimental Gaussian covariances/variograms can underlying Gaussian random fields are fitted indirectly, on
finally be fitted by theoretical covariance/variogram models. the basis of the experimental indicator variograms:
accordingly, the user is not aware of which basic nested
2.3.4 Comments on the modeling process structures are the most appropriate and a trial-and-error
process is often used (Armstrong et al. 2011). In our pro-
Although the use of more than two Gaussian random fields posal, the experimental indicator variograms are trans-
for plurigaussian modeling has already been proposed by formed into experimental variograms for the Gaussian
several authors (Xu et al. 2006; Emery 2007; Deutsch and random fields, allowing a direct and more accurate fitting
Deutsch 2014), the hierarchical approach presented in of these variograms.

123
Stoch Environ Res Risk Assess

Table 1 Codification of geo- Geo-domain Y1 \ t1? Y2 \ t2? Y3 \ t3? Y4 \ t4? Y5 \ t5? Y6 \ t6?
domain data into indicator data
1 1 Unknown Unknown Unknown Unknown Unknown
2 0 1 Unknown Unknown Unknown Unknown
3 0 0 1 Unknown Unknown Unknown
4 0 0 0 1 Unknown Unknown
5 0 0 0 0 1 Unknown
6 0 0 0 0 0 1
7 0 0 0 0 0 0

3 Case study: Rı́o Blanco ore deposit molybdenite and minor pyrite. These breccias are dated at
between 5.9 and 5.1 Ma and are found to the west of the
3.1 Geographical and geological settings potassic zone, around the tourmaline-sulfide Donoso
breccia, and to the south-east, around the Sur-Sur tour-
The Rı́o Blanco–Los Bronces ore deposit is a breccia- maline-sulfide-iron oxide breccia.
hosted cogenetic porphyry Cu–Mo deposit located between A detailed geological description of the Rı́o Blanco deposit
latitudes 33070 4500 and 33100 2000 in the Central Chilean can be found in the specialized literature (Skewes and Stern
Andes, about 80 km northeast of Santiago. Los Bronces is 1995, 1996; Serrano et al. 1996; Kay et al. 1999; Vargas et al.
emplaced in the western part and is being mined by Anglo 1999; Kay and Mpodozis 2001, 2002; Skewes et al. 2003;
American, while the Rı́o Blanco deposit is emplaced in the Frikken 2003; Frikken et al. 2005; Hollings et al. 2005).
eastern part, with mineralized bodies (Rı́o Blanco, La
Unión, Central, Don Luis, Sur-Sur and La Americana)
3.2 Lithology
distributed over an area of 6 km from north to south and
2 km from east to west. La Unión, Don Luis and Sur-Sur
Seven main rock units can be recognized within the Rı́o
are mined by Codelco (Andina Division) by open pit and
Blanco deposit. These are, from the oldest to the youngest,
underground operations.
the following (Table 2).
The Rı́o Blanco deposit is partly hosted in the western
margin of the San Francisco batholith, a multiphase intrusive (1) Andesite (AND) The andesite unit consists of volca-
containing quartz–diorite, granodiorite, quartz–monzonite, nic rocks (mainly basaltic andesite, traquiandesite
quartz–monzodiorite and granite, and partly in Miocene and dacite) of the Abanico and Farellones forma-
volcanics of the surrounding Abanico and Farellones for- tions, from late Cretaceous to Oligocene for the
mations, which consist of sub-horizontally dipping andesitic Abanico Formation and Miocene for the Farellones
lavas (Warnaars et al. 1985; Stambuk et al. 1988). Formation. These rocks are affected by quartz-
The earlier alteration phase within the San Francisco sericitic alteration, biotitization and partial chloriti-
batholith is characterized by actinolite ± magnetite ± zation (Stambuk et al. 1988).
titanite ± plagioclase disseminations, narrow actino- (2) Granitoids (GD) The San Francisco Batholith
lite ± magnetite veins and magnetite ± clinopyroxene ± embraces different types of granitoids, the main ones
actinolite ± plagioclase breccias. being the Rı́o Blanco granodiorite in the northern
This phase is postdated by an intense potassic alteration sector, characterized by a coarser grain size and the
of both andesitic lavas and granitoid rocks, centered on the presence of orthoclase feldspar and plagioclase, and
Rı́o Blanco sector and associated with the emplacement of the cascade granodiorite in the southern sector,
quartz monzonite porphyry intrusions, stockwork veins characterized by a finer grain size and enriched in
with biotite ± K feldspar ± quartz ± magnetite ± anhy- oligoclase-andesine. Both granitoids mainly corre-
drite ± sulfides, and multiple mineralized breccias whose spond to granodiorite, quartz-monzonite, tonalite and
matrices are dominated by biotite, tourmaline, quartz, diorite (Stambuk et al. 1988; Frikken 2003).
anhydrite and/or specularite with sulfides (chalcopyrite, (3) Tourmaline breccia (TOB) This rock unit is mainly
bornite, molybdenite and minor pyrite). hosted in granitoids. It outcrops in the Sur-Sur sector
A second mineralization phase comprises a set of over 5 km long striking N10W to N30W and
younger, tourmaline rich breccia pipes with varying pro- presents sharp and steeply dipping contacts with the
portions of quartz, chlorite, sericite, anhydrite, specularite, host rock (Serrano et al. 1996). The matrix is a kind
biotite and sulfides, including chalcopyrite, bornite, of milled rock flour replaced by secondary

123
 Stoch Environ Res Risk Assess

Table 2 Main rock units in the Rı́o Blanco deposit

Code Rock unit Rock types Abbreviation Age

7 Andesite Andesite AND Late Cretaceous to early Pliocene

6 Granitoids Rı́o Blanco granodiorite GD Early Pliocene
Cascade granodiorite
Breccia cascade granodiorite
Diorite
Granodiorite
5 Tourmaline breccia Andesite tourmaline TOB Early Pliocene
Rı́o Blanco granodiorite tourmaline
Cascade granodiorite tourmaline
Tourmaline breccia
Tourmaline breccia with andesite
Tourmaline breccia with cascade granodiorite
Tourmaline breccia with Rı́o Blanco tourmaline
Tourmaline breccia with porphyry
Tourmaline breccia with abundant tourmaline
Quartz sericitic breccia
Tourmaline breccia with granodiorite
4 Monolithic Breccia Breccia porphyry MOB Early Pliocene
Magmatic breccia of granodiorite
Molybdenite breccia
Monolithic breccia
Tourmaline monolithic breccia
Castellana breccia
Don Luis porphyry breccia
Quartz breccia with Don Luis porphyry
Tourmaline breccia of Don Luis porphyry
Biotite metasomatic rocks
Igneous breccia
3 Magmatic breccia Rı́o Blanco granodiorite breccia MAB Late Miocene
Anhydrite breccia
Hematite breccia
Quartz anhydrite breccia
2 Porphyry Quartz Monzonite Porphyry POR Middle Miocene
Feldespathic Porphyry
Monzonite Porphyry
Don Luis Porphyry
1 Pipe Dacite contact breccia PIP Early Miocene
Rhyolitic Pipe
Dacitic Pipe
Rhyolite contact breccia

tourmaline cement containing specular hematite, (4) Monolithic breccia (MOB) This rock unit is present
sulfide, quartz, magnetite, sulfate (anhydrite, gyp- in the western sector of Rı́o Blanco and has low
sum) and biotite, with open spaces filled by tourma- copper contents. It was formed after the tourmaline
line-quartz-sulfide. For lower tourmaline contents, breccia (Frikken 2003) and is characterized by a
the clasts are angular to sub-angular, while they chloritic alteration, which appears in both the clasts
become more rounded for higher tourmaline and the matrix, with a varying intensity (Vargas et al.
contents. 1999).

123
Stoch Environ Res Risk Assess

Fig. 7 Plan view showing a the

drill hole data and b the
interpreted lithological model
for elevation 3348 m a.m.s.l

(5) Magmatic breccia (MAB) This rock unit is charac- The drill holes are in an irregular design along the northwest-
terized by crystalline igneous cement. A heteroge- southeast direction (Fig. 7a); the data consist of information
neous brecciation resulted in variable proportions of on the prevailing rock types logged at drill hole cores and
rock flour matrix and clasts. The clasts are different codified into the previously defined rock units (Table 2). As
in size, from 2 mm to more than 5 m wide, and are for the lithological model, it consists of a discretization of the
altered with primary igneous minerals replaced by region of interest into blocks with a support of
secondary biotite, magnetite, sulfides and K-feld- 15 m 9 15 m 9 16 m, for which the (assumed) prevailing
spar. The matrix is partly substituted for hydrother- rock unit is assigned to each block (Fig. 7b). Such a block
mal minerals such as biotite, sulfates, quartz, model has been constructed on the basis of the available drill
sulfides, magnetite, specularite and tourmaline. hole data, as well as on the interpretation of the geological
(6) Porphyry (POR) The porphyry unit intruded the processes, emplacements and alterations.
magmatic and tourmaline breccias and consists of three
main rock types (Serrano et al. 1996): the Don Luis 3.4 Plurigaussian modeling
Porphyry, located in the center of the deposit and bounded
by the magmatic breccia (this porphyry predominates in 3.4.1 Truncation rule
the region that will be considered in the next sections); the
Quartz Monzonite Porphyry, which cross-cuts the mag- The rock units will be simulated by using the plurigaussian
matic breccia; and the Feldspar Porphyry, less abundant model presented in Sect. 2.2, with the hierarchical
and presenting a clear-cut contact with the magmatic approach defined in Sect. 2.3. To allow every unit to be in
breccia. These porphyries contain quartz, plagioclase, contact with any other unit, six independent Gaussian
K-feldspar and biotite phenocrystals. random fields {Y1, Y2, Y3, Y4, Y5, Y6} will be considered,
(7) Pipe This unit is composed of two main rock types with the same truncation rule as that defined in Fig. 6 and
(Serrano et al. 1996): the Dacitic Pipe, which forms Eq. (3). Such a truncation rule agrees with the rock unit
an inverted cone limited by the magmatic breccia, chronology: rock unit 1 (PIP) appears as the youngest and
containing phenocrystals of quartz and plagioclase rock unit 7 (AND) as the oldest. A younger rock unit is
with minor biotite in a quartz-feldspar matrix; and likely to cross-cut older ones.
the Rhyolitic Pipe, which corresponds to the final
intrusion in the mineralized sectors, characterized by
numerous felsic intrusions. 3.4.2 Truncation thresholds

As explained in Sect. 2.3.2, the truncation thresholds {t1,

t2, t3, t4, t5, t6} are defined in order to reproduce the rock
3.3 Data presentation unit proportions. For a greater realism of the realizations,
local proportions (consequently, local thresholds) can be
In the following, we will focus on the northern sector of the used at this stage. For instance, from the interpreted lith-
Rı́o Blanco deposit, for which two types of data are available: ological model (Fig. 7b), it is seen that pipe (PIP) is more
exploration drill holes and an interpreted lithological model. frequent in the northern part of the region under study,

123
 Stoch Environ Res Risk Assess

while tourmaline breccia (TOB) and porphyry (POR) are very small window, the local proportions will tend to 0 or
more frequent in the southern part. To determine the local 1, depending on whether or not the rock unit is present in
proportions and local thresholds, the following procedure is the interpreted model at the target block or sample; in such
proposed: a case, one is extremely confident in the interpreted model
and the simulated rock units will not depart from this
(1) Based on the interpreted lithological model
model too much. In the present study, a rectangular moving
(Fig. 7b), for each block and each drill hole sample,
window of 135 m 9 135 m 9 144 m (i.e., nine blocks
calculate the local proportion of each rock unit
along each direction) has been chosen for the calculation of
(proportion of the rock unit observed in a moving
the local proportions of rock units; this size is a trade-off
window centered at the target block or target drill
between our uncertainty and our reliance on the interpreted
hole sample). The use of moving windows to infer
model. Such a choice can be validated by checking the
spatially varying proportions is not novel (Arm-
ability of the plurigaussian model to accurately measure the
strong et al. 2011), although this approach is often
geological uncertainty, see Sect. 3.6 on the split-sample
applied to the sole drill hole data and the resulting
method.
proportions are further interpolated to regions with-
out data. Here, there is no need for such an
interpolation, insofar as the local proportions of
3.4.3 Variogram analysis
rock units are derived from the interpreted litholog-
ical model, which covers the entire region of interest.
The drill hole data on the prevailing rock units are codified
(2) Convert the local proportions into thresholds, leav-
into indicators (Table 1) and the experimental indicator
ing unchanged the design of the truncation rule.
variograms are calculated along the three identified main
The size of the moving window used in step (1) depends anisotropy directions: N20W, N70E and vertical. For the
on how much the interpreted model can be trusted. With a first two directions, the lag separation distance is set to
very large window, the calculated proportions will tend to 15 m with a tolerance of 7.5 m, while this distance is set to
the global proportions, meaning that the spatial variations 10 m with a tolerance of 5 m along the vertical. From the
of rock unit proportions that are perceptible in the inter- indicator variograms so obtained, one can derive vario-
preted model tend to be ignored. On the contrary, with a grams for the underlying Gaussian random fields (Eq. (4)).

Fig. 8 Experimental (crosses) and fitted (solid lines) variograms for the six underlying Gaussian random fields, along the main anisotropy
directions: N20W (green), N70E (blue) and vertical (black)

123
Stoch Environ Res Risk Assess

Table 3 Specification of fitted variogram models 3.5 Plurigaussian simulation

Gaussian Basic Sill Range Range Range
random nested along along along 3.5.1 Construction of realizations of rock units
field structure N20W N70E vertical
(m) (m) (m) Provided with the complete model (truncation rule, local
1 Cubic 1.000 1,000 450 900 truncation thresholds and variograms of underlying
2 Cubic 0.713 70 40 80 Gaussian random fields), one can simulate the rock units
Cubic 0.287 500 60 Infinite over the domain of interest, conditionally to the drill hole
3 Cubic 1.000 1,200 280 600
data (i.e., restricting the realizations to reproduce the
4 Cubic 0.084 20 30 20
observed rock units at the drill hole locations). A total of
100 realizations have been constructed, using an adaptation
Cubic 0.403 300 250 400
of a program published by Emery (2007). Table 4 indicates
Cubic 0.513 800 250 Infinite
the main implementation parameters used for simulation.
5 Cubic 0.013 10 60 60
Figure 9 shows a plan view of four realizations, which
Cubic 0.117 60 60 60
globally agree with the interpreted model depicted in
Cubic 0.314 300 200 220
Fig. 7b. Each realization gives an alternative scenario of
Cubic 0.556 Infinite 500 2,500
what could be the distribution of rock units, as it repro-
6 Cubic 1.000 700 280 500
duces the geostatistical properties of the actual rock units
(chronology and contact relationships controlled by the
truncation rule; local proportions controlled by the trun-
Table 4 Parameters for plurigaussian simulation
cation thresholds; spatial continuity controlled by the fitted
Parameter Value variograms) and, furthermore, is consistent with the drill
Number of iterations (sweeps over the dataset) for Gibbs 30 hole conditioning data. The set of realizations additionally
sampler provides an insight into geological uncertainty, as one can
Number of lines for turning bands simulation 500 appreciate whether or not the rock units change from one
Number of neighboring data for conditioning 40 scenario to another.
Radius (m) of search ellipsoid along direction N20W 200
Radius (m) of search ellipsoid along direction N70E 120
3.5.2 Post-processing the realizations
Radius (m) of search ellipsoid along vertical direction 120
The realizations can be used to assess the uncertainty in the
volume of each rock unit within the simulation domain or,
The next step is to fit variogram models (independently)
equivalently, in the proportion of the simulation domain
to the experimental variograms. Owing to the fact that the
covered by each rock unit. This can be done by calculating,
plurigaussian model is based on standard Gaussian random
for each realization, the proportion of blocks belonging to
fields, it requires the variogram sills to be equal to 1. In
each rock unit, then the statistics on the proportions so
order to ease modeling and to ensure a unit sill value, a
obtained over the realizations. As the simulation domain is
semi-automated fitting algorithm has been used (Emery
bounded, the theory of ergodic fluctuations or increasing-
2010). The fitting minimizes the mean squared error
domain asymptotics (Cressie 1993; Chilès and Delfiner
between experimental and modeled variograms (Fig. 8). At
2012) states that these proportions are likely to deviate
this stage, two comments are worth being made. First, the
from the expected proportions (those corresponding to the
cross-variograms between Gaussian random fields are set
lithological model) and vary from one realization to
to zero, since these random fields are assumed to be
another. In particular, one observes a greater uncertainty in
independent. Second, the choice of the basic nested struc-
the proportion of granitoids, which varies between 46.1 and
tures used for modeling the direct variograms has an
53.7 % (Table 5). This can be explained because this rock
impact on the nature of the boundaries between rock units:
unit is under-sampled with respect to the other ones
structures with a linear behavior near the origin (like the
(Fig. 7), so that the layout of its boundaries is more
spherical and exponential) lead to irregular boundaries,
uncertain. The range of uncertainty (7.6 %) lies in-between
while structures with a parabolic behavior near the origin
what would be obtained in the following two extreme
(like the cubic and Gaussian) lead to smooth boundaries
cases:
(Lantuéjoul 2002). Because this last property is deemed
more realistic for the current case study, the sample vari- • No spatial correlation (pure nugget variograms): in this
ograms are fitted with nested anisotropic cubic models case, the proportion of granitoids would be the average
(Table 3). of 23,750 independent Bernoulli variables (one variable

123
 Stoch Environ Res Risk Assess

Fig. 9 Plan views of four

realizations of the rock units for
elevation 3348 m a.m.s.l,
conditioned to the drill hole data
of Fig. 7a

Table 5 Statistics on the proportions of blocks belonging to each • Perfect spatial correlation (variograms increasing infi-
rock unit nitely slowly along all the directions): in this case, the
Rock Minimum Maximum Average over Interpreted proportion of granitoids would be the result of a single
unit over over realizations model (%) Bernoulli variable with mean 0.509, identical for all the
realizations realizations (%) blocks. Over 100 realizations, the expected variation in
(%) (%)
this proportion would be 100 %.
PIP 5.1 9.2 7.3 7.5
One can also assess the uncertainty in the rock units at a
POR 20.1 24.5 22.6 22.2
local (block-by-block) scale, by means of probability maps.
MAB 1.6 4.5 2.8 1.1 The maps are constructed by calculating, for each block,
MOB 3.7 8.0 5.7 6.9 the frequency of occurrence of each rock unit over the 100
TOB 8.1 14.0 10.3 10.2 conditional realizations (Fig. 10). They constitute a com-
GD 46.1 53.7 50.3 50.9 plement to the interpreted lithological model, insofar as
AND 0.2 2.2 1.1 1.2 they show the risk of finding a rock unit different from the
one that has been scheduled. The sectors with little
per block) with mean equal to 0.509 (expected propor- uncertainty are those associated with a high probability for
tion of granitoids, as per the interpreted lithological a given rock unit (painted in red in Fig. 10), indicating that
model). As the standard deviation of such an average is there is little risk of not finding this rock unit, or those
very small (0.0032), the variation in the proportion of associated with a very low probability (painted in dark blue
granitoids over the realizations would be close to zero. in Fig. 10), indicating that one is pretty sure of not finding

123
Stoch Environ Res Risk Assess

Fig. 10 Plan views showing the probabilities of occurrence of each rock unit at elevation 3348 m a.m.s.l, calculated using 100 realizations

this unit, while the other sectors (painted in light blue, unit matches the most probable one, and the blocks for
green or yellow in Fig. 10) are more uncertain. which there is a mismatch and the interpreted model may
From these probability maps, one can construct a lith- be revised (Fig. 11c). The latter blocks are mostly located
ological model by selecting, for each block, the most near the boundaries of pipe, porphyry and granitoids, and
probable rock unit (Fig. 11a). This model can then be correspond to the blocks with a moderate probability
compared to the interpreted lithological model (Fig. 7b) in (0.3–0.6) for the most probable rock unit (in other words,
order to identify the blocks for which the interpreted rock none of the units has a high probability of occurrence)

123
 Stoch Environ Res Risk Assess

Fig. 11 Plan views showing a the most probable rock unit at elevation 3348 m, b the probability of occurrence of the most probable unit, and
c the matches and mismatches with respect to the interpreted model (Fig. 7b)

Table 6 Numbers of blocks per Interpreted model Most probable unit

rock unit, according to
interpreted model and to model PIP POR MAB MOB TOB GD AND Total
based on the most probable rock
unit PIP 1,699 26 0 0 0 57 0 1,782
POR 33 4,182 29 59 206 752 6 5,267
MAB 17 27 142 3 11 52 0 252
MOB 31 443 39 511 181 424 17 1,646
TOB 43 179 28 42 1,679 431 13 2,415
GD 191 285 170 205 307 10,938 0 12,096
AND 0 135 1 24 12 72 48 292
Total 2,014 5,277 409 844 2,396 12,726 84 23,750

Table 7 Transition When leaving a block of rock unit… One finds a block belonging to rock unit…
probabilities calculated by
considering the simulated rock PIP POR MAB MOB TOB GD AND
units on adjacent blocks
(average statistics over 100 PIP 0.9315 0.0180 0.0057 0.0056 0.0045 0.0345 0.0000
realizations) POR 0.0057 0.7930 0.0114 0.0385 0.0439 0.0997 0.0078
MAB 0.0142 0.0924 0.7125 0.0207 0.0403 0.1152 0.0046
MOB 0.0072 0.1532 0.0102 0.6290 0.0513 0.1385 0.0106
TOB 0.0032 0.0960 0.0109 0.0280 0.7703 0.0874 0.0043
GD 0.0050 0.0452 0.0064 0.0157 0.0181 0.9068 0.0028
AND 0.0012 0.1673 0.0130 0.0582 0.0449 0.1362 0.5791

(Fig. 11b). The numbers of matches and mismatches for pipe (PIP) and andesite (AND). Also, around 8–14 % of
the blocks at elevation 3348 m a.m.s.l are indicated in the porphyry (POR), breccia (TOB, MOB, MAB) and
Table 6: one obtains 19,199 matches over a total of 23,750 andesite (AND) blocks are in contact with granitoid (GD),
blocks, i.e., the interpreted rock unit corresponds to the and 9–15 % of breccia (TOB, MOB, MAB) blocks are in
most probable unit for 81 % of the blocks. contact with porphyry (POR).
Other useful information that can be obtained from the
realizations relates to the so-called transition probabilities, 3.6 Model validation
which give the frequencies of simulated rock units at two
adjacent blocks, on average over the simulation domain In order to make sure that plurigaussian simulation pro-
and over the realizations (Table 7). In particular, one vides accurate results, a split-sample validation exercise is
observes that every rock unit is likely to be in contact with carried out. It consists in partitioning the drill hole dataset
every other, except the youngest and oldest ones, namely into two subsets and simulating the rock units at the

123
Stoch Environ Res Risk Assess

Fig. 12 Split-sample validation: calibration plot comparing the probability of occurrence of a rock unit (abscissa) and the proportion of rock unit
actually observed at the testing subset (ordinate). The red line corresponds to the identity

locations of one of these subsets (referred to as the testing of the testing subset. Then, for each rock unit, the data of
subset), using as conditioning information the data of the the testing subset are sorted by increasing probability of
other subset (training subset). In the present case, the drill occurrence of this unit: for a given probability p (up to
hole data belonging to the testing and training subsets have some calculation tolerances), one expects that a proportion
been chosen at random, with half of the data in each subset. p of the testing data actually match this rock unit. The
Using the realizations so obtained, the probability of match or mismatch can be visualized in a calibration plot,
occurrence of each rock unit can be calculated at each data as shown in Fig. 12.

123
 Stoch Environ Res Risk Assess

Fig. 13 Calibration plots associated with pipe (PIP) for different window sizes used for calculating the local rock unit proportions

Arguably, this validation exercise makes use of the always match the interpreted rock units and the testing data
interpreted lithological model (through the calculation of are perfectly reproduced, insofar as these data agree with
local rock unit proportions and local thresholds, see the interpreted lithological model, giving the impression
Sect. 3.4.2), which is constructed on the basis of the entire that the model is validated. In such a case, the calibration
drill hole dataset, and not only the training subset. The plot consists of only two points, corresponding to proba-
testing data are thus indirectly accounted for, which may bilities 0 and 1 (see Figs. 13a and 14a for an illustration to
lead to an overly optimistic assessment of the ability of the two rock units).
plurigaussian model to quantify uncertainty. For instance, Such a drawback is mitigated when the local rock unit
if the local proportions are calculated within a small proportions are calculated using a large moving window, in
moving window, then these proportions will be either 0 or which case the value of any particular testing data is not
1, depending on whether the rock unit at the block under influential any more in the local proportion model. Based
consideration is absent or present in the interpreted litho- on this premise, to avoid the validation results being mis-
logical model. In such a case, the simulated rock units leading, the size of the moving window should be the

123
Stoch Environ Res Risk Assess

Fig. 14 Calibration plots associated with monolithic breccia (MOB) for different window sizes used for calculating the local rock unit
proportions

largest one that yields calibration plots close to the diag- quantify uncertainty and the request of an as large as
onal line. In particular, if a stationary model were adequate, possible moving window. The window size could be
the moving window should be infinitely large and the reduced if our confidence in the interpreted lithological
global rock unit proportions should be used. In the present model increases, but this would be a somewhat subjective
case study however, a stationary model is questionable and decision.
the use of global proportions does not yield a fully satis-
factory validation (Figs. 13d, 14d). When decreasing the
window size, the validation results improve, becoming 4 Conclusions
acceptable with a size of 135 m 9 135 m 9 144 m
(Figs. 13c, 14c) or smaller (Figs. 13b, 14b). This is the The plurigaussian model offers a simple framework for
reason why a moving window of size 135 m 9 simulating geo-domains, but the current implementations
135 m 9 144 m has been used in Sect. 3.4.2, as this option are often restricted in the number of geo-domains and in
fulfills the trade-off between the ability to adequately their contact relationships. The proposed hierarchical

123
 Stoch Environ Res Risk Assess

approach overcomes these limitations and increases the (ed) 6th International mining geology conference, Darwin.
model versatility, providing a set of realistic models that Australasian Institute of Mining and Metallurgy, Melbourne,
pp 193–203
constitute alternative representations of the true (unknown) Deutsch CV (2002) Geostatistical reservoir modeling. Oxford Uni-
distribution of geo-domains in space: versity Press, New York
Deutsch JL, Deutsch CV (2014) A multidimensional scaling approach
• Each realization reproduces the contact relationships to enforce reproduction of transition probabilities in truncated
between geo-domains. These relationships can be plurigaussian simulation. Stoch Environ Res Risk Assess
designed in order to agree with chronology, so that 28(3):707–716
Dowd PA, Pardo-Igúzquiza E, Xu C (2003) Plurigau: a computer
younger geo-domains appear in the foreground and are
program for simulating spatial facies using the truncated
likely to cross-cut older geo-domains, which appear in plurigaussian method. Comput Geosci 29(2):123–141
the background. Emery X (2007) Simulation of geological domains using the
• Each realization also reproduces the local geo-domain plurigaussian model: new developments and computer programs.
Comput Geosci 32(9):1189–1201
proportions, via the choice of local truncation thresh-
Emery X (2010) Iterative algorithms for fitting a linear model of
olds, and the spatial correlation of the geo-domain coregionalization. Comput Geosci 36(9):1150–1160
indicators, via the choice of the variograms of the Emery X, González KE (2007a) Incorporating the uncertainty in
underlying Gaussian random fields. geological boundaries into mineral resources evaluation. J Geol
Soc India 69(1):29–38
• Finally, each realization is consistent with the data used
Emery X, González KE (2007b) Probabilistic modelling of litholog-
as conditioning information in the simulation process. ical domains and its application to resources evaluation. J S Afr
Inst Min Metall 107(12):803–809
In addition to these properties, the construction of Emery X, Ortiz JM, Cáceres AM (2008) Geostatistical modelling of
multiple realizations (as opposed to a single interpreted rock type domains with spatially varying proportions: application
model) allows assessing geological uncertainty, through to a porphyry copper deposit. J S Afr Inst Min Metall
the calculation of probability maps and transition proba- 108(5):285–292
Fontaine L, Beucher H (2006) Simulation of the Muyumkum uranium
bility matrices, and identifying sectors for which the actual roll front deposit by using truncated plurigaussian method. In:
geo-domains are more uncertain. These outputs are helpful Dominy S (ed) 6th International mining geology conference,
for decision-making related to the valorization, planning Darwin. Australasian Institute of Mining and Metallurgy, Mel-
and operation of the mining project. bourne, pp 205–215
Frikken PH (2003) Breccia-hosted copper-molybdenum mineraliza-
tion at Rı́o Blanco, Chile. PhD Dissertation. University of
Acknowledgments This research was funded by the Chilean Com- Tasmania, Australia
mission for Scientific and Technological Research, through Project Frikken PH, Cooke DR, Walshe JL, Archibald D, Skarmeta J, Serrano
CONICYT/FONDECYT/REGULAR/No. 1130085. The authors L, Vargas R (2005) Mineralogical and isotopic zonation in the
acknowledge the support from Claudio Martı́nez from Codelco-Chile Sur-Sur tourmaline breccia, Rı́o Blanco-Los Bronces Cu-Mo
(Andina Division), who provided the dataset used in this work, as well deposit, Chile: implications for ore genesis. Econ Geol
as the comments by Grégoire Mariethoz and another anonymous 100:935–961
reviewer, who helped to improve the manuscript. Galli A, Beucher H, Le Loc’h G, Doligez B, Heresim Group (1994) The
pros and cons of the truncated Gaussian method. In: Armstrong M,
Dowd PA (eds) Geostatistical simulations. Kluwer, Dordrecht,
pp 217–233
References Hollings P, Cooke D, Clark A (2005) Regional geochemistry of
tertiary igneous rocks in Central Chile: implications for the
Allard D, D’Or D, Biver P, Froidevaux R (2012) Non-parametric geodynamic environment of giant porphyry copper and epither-
diagrams for plurigaussian simulations of lithologies. In: Ninth mal gold mineralization. Econ Geol 100:887–904
international geostatistics congress, Oslo. http://geostats2012.nr. Jeannée N, Bardou E, Faucheux C, Ornstein P (2013) Geostatistical
no/pdfs/1746887.pdf. Accessed 23 Nov 2014 assessment of ice content distribution within the Glacier
Armstrong M, Galli A, Beucher H, Le Loc’h G, Renard D, Doligez B, Bonnard. Math Geosci 45:591–599
Eschard R, Geffroy F (2011) Plurigaussian simulations in Journel AG, Alabert FG (1990) New method for reservoir mapping.
geosciences. Springer, Berlin J Pet Technol 42(2):212–218
Carrasco P, Ibarra F, Rojas R, Le Loc’h G, Séguret S (2007) Journel AG, Huijbregts CJ (1978) Mining geostatistics. Academic
Application of the truncated Gaussian simulation method to a Press, London
porphyry copper deposit. In: Magri E (ed) 33rd International Kay SM, Mpodozis C (2001) Central Andes ore deposits linked to
symposium on applications of computers and operations research evolving shallow subduction systems and thickening crust. GSA
in the mineral industry, Gecamin, Santiago, pp 31–39 Today 11:4–9
Chilès JP, Delfiner P (2012) Geostatistics: modeling spatial uncer- Kay SM, Mpodozis C (2002) Magmatism as a probe to the Neogene
tainty, 2nd edn. Wiley, New York shallowing of the Nazca plate beneath the modern Chilean flat-
Cressie NAC (1993) Statistics for spatial data, rev edn. Wiley, New York slab. J S Am Earth Sci 15:39–57
De Marsily G, Delay F, Gonçalves J, Renard P, Teles V, Violette S Kay SM, Mpodozis C, Coira B (1999) Neogene magmatism,
(2005) Dealing with spatial heterogeneity. Hydrogeol J tectonism and mineral deposits of the central Andes. In: Skinner
13(1):161–183 BJ (ed) Geology and ore deposits of the Central Andes. Special
Deraisme J, Field M (2006) Geostatistical simulations of kimberlite Publication No. 7. Society of Economic Geologists, Littleton,
orebodies: application to sampling optimization. In: Dominy S pp 7–59

123
Stoch Environ Res Risk Assess

Kyriakidis PC, Deutsch CV, Grant ML (1999) Calculation of the styles and metallogeny. Special Publication No. 5. Society of
normal scores variogram used for truncated Gaussian lithofacies Economic Geologists, Littleton, pp 119–130
simulation: theory and Fortran code. Comput Geosci 25(2):161– Skewes MA, Holmgren C, Stern CR (2003) The Donoso copper-rich,
169 tourmaline-bearing breccia pipe in central Chile: petrologic, fluid
Lantuéjoul C (2002) Geostatistical simulation: models and algo- inclusion and stable isotope evidence for an origin from
rithms. Springer, Berlin magmatic fluids. Miner Deposita 38(1):2–21
Le Loc’h G, Galli A (1997) Truncated plurigaussian method: Skvortsova T, Armstrong M, Beucher H, Forkes J, Thwaites A,
theoretical and practical points of view. In: Baafi EY, Schofield Turner R (2001) Applying plurigaussian simulations to a granite-
NA (eds) Geostatistics Wollongong’ 96. Kluwer, Dordrecht, hosted orebody. In: Kleingeld WJ, Krige DG (eds) Geostats 200
pp 211–222 Cape Town. Geostatistical Association of Southern Africa,
Mariethoz G, Renard P, Cornaton F, Jaquet O (2009) Truncated Johannesburg, pp 904–911
plurigaussian simulations to characterize aquifer heterogeneity. Skvortsova T, Beucher H, Armstrong M, Forkes J, Thwaites A,
Ground Water 47(1):13–24 Turner R (2002) Simulating the geometry of a granite-hosted
Matheron G, Beucher H, de Fouquet C, Galli A, Guerillot D, Ravenne uranium orebody. In: Armstrong M, Bettini C, Champigny N,
C (1987) Conditional simulation of the geometry of fluvio- Galli A, Remacre A (eds) Geostatistics Rio 2000. Kluwer,
deltaic reservoirs. In: 62nd Annual technical conference and Dordrecht, pp 85–99
exhibition of the society of petroleum engineers, pp 591–599. Stambuk V, Aguilar C, Blondel J, Galeb M, Serrano L, Vargas R
SPE Technical Paper 16753 (1988) Geologı́a del yacimiento Rı́o Blanco. Technical Report.
Matheron G, Beucher H, de Fouquet C, Galli A, Ravenne C (1988) Codelco-Chile, División Andina
Simulation conditionelle à trois facies d’une falaise de la Strebelle S (2002) Conditional simulation of complex geological
formation du Brent. Sciences de la Terre, Série Informatique structures using multiple-point statistics. Math Geol 34(1):1–22
Géologique 28:213–249 Talebi H, Asghari O, Emery X (2013) Application of plurigaussian
Rezaee H, Asghari O, Koneshloo M, Ortiz JM (2014) Multiple-point simulation to delineate the layout of alteration domains in
geostatistical simulation of dykes: application at Sungun Sungun copper deposit. Cent Eur J Geosci 5(4):514–522
porphyry copper system, Iran. Stoch Environ Res Risk Assess Talebi H, Asghari O, Emery X (2014) Simulation of the lately
28(7):1913–1927 injected dykes in an Iranian porphyry copper deposit using the
Rondon O (2009) A look at the plurigaussian simulation for a nickel plurigaussian model. Arab J Geosci 7(7):2771–2780
laterite deposit. In: Dominy S (ed) 7th International mining & Vargas R, Gustafson LB, Vukasovic M, Tidy E, Skewes MA (1999)
geology conference. The Australasian Institute of Mining and Ore breccias in the Rı́o Blanco-Los Bronces copper deposit,
Metallurgy, Melbourne Chile. In: Skinner BJ (ed) Geology and ore deposits of the
Serrano L, Vargas R, Stambuk V, Aguilar C, Galeb M, Holmgren C, Central Andes. Special Publication No. 7. Society of Economic
Contreras A, Godoy S, Vela I, Skewes MA, Stern CR (1996) The Geologists, Littleton, pp 281–297
late Miocene to early Pliocene Rı́o Blanco-Los Bronces copper Warnaars F, Holgrem C, Barassi S (1985) Porphyry copper and
deposit, Central Chilean Andes. In: Camus F, Sillitoe RH, tourmaline breccias at Los Bronces—Rı́o Blanco, Chile. Econ
Petersen R (eds) Andean copper deposits: new discoveries, Geol 80:1544–1565
mineralization, styles and metallogeny. Special Publication No. Xu C, Dowd PA, Mardia KV, Fowell RJ (2006) A flexible true
5. Society of Economic Geologists, Littleton, pp 119–130 plurigaussian code for spatial facies simulations. Comput Geosci
Skewes MA, Stern CR (1995) Genesis of the late Miocene to Pliocene 32(10):1629–1645
copper deposits of central Chile in the context of Andean Yunsel T, Ersoy A (2011) Geological modeling of gold deposit based
magmatic and tectonic evolution. Int Geol Rev 37:893–909 on grade domaining using plurigaussian simulation technique.
Skewes MA, Stern CR (1996) Late Miocene mineralized breccias in Nat Resour Res 20(4):1–19
the Andes of central Chile: Sr and Nd isotopic evidence for Yunsel T, Ersoy A (2013) Geological modeling of rock type domains
multiple magmatic sources. In: Camus F, Sillitoe RH, Petersen R in the Balya (Turkey) lead-zinc deposit using plurigaussian
(eds) Andean copper deposits: new discoveries, mineralization, simulation. Cent Eur J Geosci 5(1):77–89

123