Natural Resources Research, Vol. 28, No.

3, July 2019 ( 2018)

https://doi.org/10.1007/s11053-018-9444-x

Original Paper

Nonparametric Geostatistical Simulation of Subsurface

Facies: Tools for Validating the Reproduction of,
and Uncertainty in, Facies Geometry

Nasser Madani,1,4 Mohammad Maleki,2,3 and Xavier Emery2,3

Received 5 July 2018; accepted 7 December 2018

Published online: 15 December 2018

Delineation of facies in the subsurface and quantification of uncertainty in their boundaries

are significant steps in mineral resource evaluation and reservoir modeling, which impact
downstream analyses of a mining or petroleum project. This paper investigates the ability of
nonparametric geostatistical simulation algorithms (sequential indicator, single normal
equation and filter-based simulation) to construct realizations that reproduce some expected
statistical and spatial features, namely facies proportions, boundary regularity, contact
relationships and spatial correlation structure, as well as the expected fluctuations of these
features across the realizations. The investigation is held through a synthetic case study and a
real case study, in which a pluri-Gaussian model is considered as the reference for comparing
the simulation results. Sequential indicator simulation and single normal equation simulation
based on over-restricted neighborhood implementations yield the poorest results, followed
by filter-based simulation, whereas single normal equation simulation with a large neigh-
borhood implementation provides results that are closest to the reference pluri-Gaussian
model. However, some biases and inaccurate fluctuations in the realization statistics (facies
proportions, indicator direct and cross-variograms) still arise, which can be explained by the
use of a single finite-size training image to construct the realizations.
KEY WORDS: Geological uncertainty, Pluri-Gaussian model, Sequential indicator simulation, Single
normal equation simulation, Filter-based simulation, Statistical fluctuations.

INTRODUCTION alization zones or rock-type domains in mineral

deposits, litho-facies in reservoirs or hydro-facies in
Context and Motivation aquifers) is often required to characterize quantita-
tive regionalized variables of interest, such as rock
In subsurface modeling, delimiting the spatial permeability, porosity or metal grades (Armstrong
layout of facies or geological domains (e.g., miner- et al. 2011; Rossi and Deutsch 2014). A limited
number of boreholes across the region under study
1
Department of Mining Engineering, School of Mining and are usually available, but the exact locations of the
Geosciences, Nazarbayev University, Astana, Kazakhstan. boundaries between facies remain uncertain. Prob-
2
Department of Mining Engineering, University of Chile, abilistic modeling, however, allows quantifying the
Santiago, Chile. uncertainty in these boundaries through the con-
3
Advanced Mining Technology Center, University of Chile,
struction of multiple realizations of facies that can
Santiago, Chile.
4
To whom correspondence should be addressed; e-mail: be viewed as equiprobable (equally likely to be
nasser.madani@nu.edu.kz simulated) outcomes (Chilès and Delfiner 2012).

1163
1520-7439/19/0700-1163/0  2018 International Association for Mathematical Geosciences
1164 Madani, Maleki, and Emery

Stochastic simulation approaches to construct tained through surface, satellite or object-based

facies outcomes can be classified into two main modeling (Al-Mudhafar 2018). The choice of a TI
families according to a parametric or nonparametric and its representativeness is a challenging issue that
statistical framework. The former rely on the spec- is subject to ongoing developments (Boisvert et al.
ification of a parent spatial random field, charac- 2007; Emery and Lantuéjoul 2014; Mariethoz and
terized by a set of all its finite-dimensional Caers 2014; Dubrule 2017).
distributions. For example, truncated Gaussian Several qualitative and quantitative tools and
simulation (Matheron et al. 1987; Galli et al. 1994) validation techniques have been proposed to assess
defines facies by truncating a Gaussian random field the underlying algorithms and the quality of the
and can be applied when facies obey an ordered realizations they produce. Leuangthong et al. (2004)
sequence that can be confirmed through the geo- suggested generic criteria to check the approximate
logical setting, like strata in sedimentary deposits. reproduction of the distribution and spatial corre-
Pluri-Gaussian simulation is an extension of trun- lation structure of the realizations. In the specific
cated Gaussian simulation that allows more complex case of facies simulation, De Iaco and Maggio (2011)
contact relationships between facies and can adapt and de Iaco (2013) analyzed the reproduction of
to a broader range of geological settings. It relies on geometrical properties of facies and of high-order
the specification of several Gaussian random fields cumulants that give an insight into the spatial con-
and on a truncation rule to convert these random nectivity of facies. Boisvert et al. (2010) presented a
fields into facies. The reader is referred to the lit- set of metrics based on first- and second-order
erature (Le LocÕh and Galli 1997; Lantuéjoul 2002; statistics, global and local trends, multi-scale his-
Armstrong et al. 2011; Yunsel and Ersoy 2011; tograms and multiple-point histograms, which allow
Chilès and Delfiner 2012; Madani and Emery 2015; inquiring into the properties of simulated facies at
Madani and Emery 2017; Madani et al. 2018) for different spatial scales. Emery and Ortiz (2011)
details on model parameters and algorithms for showed, however, that the same or very similar
conditional simulation. multiple-point distributions can yield large differ-
In contrast, nonparametric approaches rely on ences in up-scaled properties and responses to
statistics and spatial continuity measures that can be transfer functions, such as spatial connectivity and
directly inferred from borehole data or from a fluid flow simulation. For an enhanced analysis of
training image (TI), which are deemed representa- spatial connectivity, several authors (Renard and
tive of the spatial arrangement of facies in the Allard 2013; Oriani and Renard 2014; Rongier et al.
subsurface. For instance, sequential indicator simu- 2016) proposed indicators relevant to the spatial
lation (SISIM) is based on coding of facies into arrangement, shape and topology of connected
indicators and on modeling of indicator variograms components (also called geo-bodies) of simulated
or covariances, which constitute two-point spatial facies. In all cases, the assessment of a simulation
continuity measures (Alabert 1987; Journel and algorithm relies on the comparison between the
Alabert 1990; Journel and Gómez-Hernández statistical parameters of the TI as the reference in-
1993). However, previous researches have shown put and the simulation results as the output. Dis-
that SISIM may not reproduce correctly the desired tance-based approaches also provide metrics for
spatial distribution and the contact relationships selecting the optimum simulation algorithm based
between facies (Emery 2004; Modis and Sideri on pattern and spatial uncertainty reproduction
2013). To circumvent that shortcoming, several (Tan et al. 2014).
algorithms based on multiple-point statistics have However, many of the aforementioned metrics
emerged in the past decades, such as single normal require gridded realizations or training images and
equation simulation (SNESIM, Strebelle 2002) or cannot be calculated when only a set of borehole
filter-based simulation (FILTERSIM, Zhang et al. data scattered in space are available. In such a sit-
2006), which enable the reproduction of more uation, the true values of the metrics remain un-
complex geometries and spatial patterns. These known. However, as will be shown in this paper,
multiple-point simulation algorithms rely on a TI simple one- and two-point metrics that can be in-
that gives a representative layout of the facies to be ferred from scattered borehole data are sufficient to
simulated. The TI can be hand-drawn by an expert capture three essential characteristics of facies.
who is knowledgeable of the geological setting of These are: (1) facies proportions; (2) regularity of
the region under study (De Iaco 2013) or it ob- their boundaries; and (3) contact relationships be-
Nonparametric Geostatistical Simulation of Subsurface Facies 1165

tween facies, which give insight into the geometry TI in the case of a multiple-point simulation algo-
and spatial arrangement of facies. rithm.
As guideline for best practices, this paper fo-
cuses on common nonparametric simulation algo-
Aims, Contribution and Scope of this Work rithms in order to detect whether or not they provide
an accurate response to modeling the layout of fa-
Therefore, this research focuses on the above- cies in the subsurface. The algorithms are selected
mentioned three characteristics and propose simple on account of their popularity among the practi-
qualitative and semi-quantitative tools to check their tioners in academia and industry and their avail-
reproduction with a given geostatistical simulation ability in commercial software. The comparison is
algorithm. Of particular interest is the examination held in a synthetic case study and a real case study,
of direct and cross-variograms of facies indicators, in which a parametric approach (pluri-Gaussian
which, in practice, are often overlooked although simulation, PLURISIM) is used as a reference rep-
they provide valuable information on the charac- resenting real facies distribution, so that the non-
teristics of interest (Emery and Lantuéjoul 2011; parametric simulation outputs could be assessed in
Maleki et al. 2017). The proposed tools are appli- comparison with such reference.
cable for determining the suitability of an algorithm
to be used for simulating facies in the subsurface, or
for comparing the performances of several simula- BACKGROUND
tion algorithms, when the exact values of multiple-
point metrics or of up-scaled properties are Essential Statistical and Geometrical Characteristics
unavailable. of Facies
This paper has two contributions. First, as
mentioned above, the recommended metrics for In this section, it is of interest to show some
validation can be calculated even if data represent- characteristics of facies through six examples, in
ing ground truth are scarce and scattered in space which the spatial layout exhibits regular boundaries
(like borehole data) and provide valuable informa- and there may exist forbidden contacts between fa-
tion on facies spatial distribution, particularly cies. These examples are collected from geological
geometry of boundaries and contact relationships. basins in static modeling of reservoir projects that
Second, the validation considers not only the have been carried out in Iran during the past decade.
reproduction of the target average over of a large However, for confidentiality reasons, the names of
number of realizations, but also the reproduction of the reservoirs are not disclosed.
their spreads or fluctuations across the realizations,
which are key elements in geostatistical modeling of
uncertainty. The analysis of statistical fluctuations is Example 1: Carbonate Platform A
too often ignored for validating geostatistical simu-
lation algorithms, with the notable exception of a This setting consists of limestone (Lst), marl
few works related to the simulation of Gaussian (Ml), argillite limestone (argLst) and shale (Sh). In
random fields (Lantuéjoul 1994; Emery and Lan- this sedimentary sequence, the Lst has regular con-
tuéjoul 2006; Emery 2008). To the best of the au- tacts with the other three facies (Fig. 1, top left).
thorsÕ knowledge, this criterion has never been
suggested for the simulation of categorical random
fields representing the spatial layout of facies in the Example 2: Carbonate Platform B
subsurface. However, it can enlighten on the ability
of an algorithm to correctly reproduce uncertainty In this example, anhydrite (An), Sh, argLst and
and, thus, on the appropriateness of the chosen in- Lst have regular boundaries. The An can touch the
puts and implementation parameters, e.g., relevant Lst, the Lst is in contact with Sh and argLst, but one
to the size of the search template (moving neigh- cannot expect to see a contact between An and Sh or
borhood) or to the size and representativeness of the argLst (Fig. 1, middle left).
1166 Madani, Maleki, and Emery

Figure 1. Maps of facies and their contact relationships for varied geological settings.

Example 3: Siliciclastic A Example 5: Siliciclastic C

This setting is composed of Sand, Sand shale This setting is formed by Sand, Lst and Sh. All
(SandSh), Sh, Lst and argLst. In addition to regular the facies can be in contact altogether while showing
boundaries, the facies exhibit forbidden contacts. In regular boundaries (Fig. 1, middle right).
particular, Sand is in contact with SandSh and
SandSh can be in contact with Sh, while Sh is in
contact with both Lst and argLst. This gradually Example 6: Siliciclastic D
converts the Sand into Lst and argLst (Fig. 1, bot-
tom left). The facies in this setting are Sand, SandSh, Sh
and argLst. Owing to the sedimentation process,
argLst is surrounded by SandSh and Sh but never
Example 4: Siliciclastic B touches Sand (Fig. 1, bottom right).

The Sand, Shale sand (ShSand), SandSh and Sh

are the facies in this example. In this case, Sand has Geostatistical Tools to Describe Facies Boundaries
contact with ShSand and the latter has contact with
both SandSh and Sh. However, Sand cannot have From a mathematical viewpoint, the concept of
direct contact with SandSh and Sh (Fig. 1, top right). regularity of facies boundary can be defined as fol-
Nonparametric Geostatistical Simulation of Subsurface Facies 1167

lows. A facies in d-dimensional Euclidean space has The preceding discussion leave questions on
a regular boundary if the specific (d-1)-volume (i.e., nonparametric simulation approaches, in which the
the average perimeter per unit surface when d = 2, random field model is replaced by a set of statistics
or the average surface area per unit volume when derived from sampling data or from a TI. Even if
d = 3) of the facies boundary is finite (Lantuéjoul they are based on statistics that are fully coherent
2002). When the facies can be modeled by a sta- with a given parametric random field model, these
tionary random set, the boundary regularity is nonparametric simulation approaches may not be
reached when the direct variogram of the facies able to reproduce the correct (parametric) distribu-
indicator along any direction of space has a finite tion of uncertainty (Emery and Lantuéjoul 2014).
slope at the origin (Lantuéjoul 2002; Emery and Instead, they may provide a biased distribution or
Lantuéjoul 2011; Dubrule 2017). In other words, the too small fluctuations for the quantity of interest
average slope at the origin of the facies indicator which is a common complaint among practitioners
variogram is a sign of whether the facies has regular when the realizations bear too much resemblance to
(finite slope) or irregular (infinite slope) boundary. all or part of the TI (Mariethoz and Caers 2014).
Likewise, the average slope at the origin of the
cross-variogram between two facies indicators con-
veys information on how much contact is shared by SYNTHETIC CASE STUDY
these facies: If the two facies have no contact, their
indicator cross-variogram has zero slope at the ori- Geological Setting
gin; otherwise, the slope at the origin is strictly
negative and increases in absolute value when the This synthetic case study is inspired by the sixth
facies are preferentially in contact (Maleki et al. oil reservoir example (siliciclastic D) above, with
2017). Accordingly, the indicator direct and cross- four facies such that facies 1 can have contact with
variograms turn out to be simple and informative facies 2 and 3, but not with facies 4. This setting has
tools to portray facies boundaries and contact rela- some relevance in mining applications, in which fa-
tionships between facies, reason for which they cies stand for alteration, mineralization or rock-type
among the preferred tools proposed in this research. domains. For instance, according to the conceptual
models of porphyry-Cu deposits (Lowell and Guil-
bert 1970), facies 1 can be the inner potassic alter-
Statistical Fluctuations and Uncertainty Modeling ation domain of a deposit, which is surrounded by
propylitic (facies 2) and sericitic (facies 3) alteration
Following Matheron (1989), a regionalized domains; in turn, an argillic alteration (facies 4)
property (here, facies in the subsurface) can be domain can be in contact with the latter two do-
modeled as a realization of a spatial random field. mains, but not with the potassic alteration. Alter-
Any quantity of interest about the regionalized natively, facies 1 to 4 can represent a superficial
property (e.g., the proportion of a region covered by gravel domain, a leached capping, an oxide and a
a given facies or the amount of contact between two mixed (oxide-sulfide) mineralization domains,
facies, reflected through their indicator cross-vari- respectively, where the last domain is found at depth
ogram at the origin) can be calculated from a large and does not touch the gravel. As another example,
set of realizations constrained to conditioning data, facies 1 may stand for an overburden, facies 2 for an
in case that sampling information is available. The intrusive dyke, facies 3 for a friable (with fine
fluctuations observed across these realizations allow granulometry) rock type, and facies 4 for a compact
appraising the uncertainty in the true but unknown (with coarse granulometry) rock type located un-
value of the original quantity (Chilès and Delfiner derneath, with no contact with the overburden.
2012). Accordingly, having specified a parametric
random field model and the region of space where
this random field exists, any simulation algorithm Simulation Settings
should accurately reproduce these fluctuations: Too
small or too large fluctuations across the realizations Pluri-Gaussian simulation has proven to be an
would indicate an incorrect assessment of the effective parametric approach to construct realistic
uncertainty (Lantuéjoul 1994, 2002). numerical models of oil reservoirs and mineral de-
1168 Madani, Maleki, and Emery

Table 1. Variogram parameters of the Gaussian random fields practical implementation considers the use of indi-
used in pluri-Gaussian simulation cator simple cokriging, with moving neighborhood
Gaussian random field Theoretical variogram Range Sill containing up to 100 previously simulated nodes for
the indicator cokriging used to model the successive
Y1 Isotropic cubic model 300 1
conditional cumulative distribution functions. The
Y2 Isotropic cubic model 300 1
target grid nodes are visited according to a random
sequence with three multiple grids.
To run the multiple-point simulation algorithms
(SNESIM and FILTERSIM), one needs to define a
TI. This TI is chosen as an additional pluri-Gaussian
realization with the same size 1000 9 1000, so as to
be perfectly compatible with the previous pluri-
Gaussian model. Three cases are considered for
SNESIM, which use a moving neighborhood with
only 1 previously simulated node (case 1), 2 previ-
ously simulated nodes (case 2) and 10 previously
simulated nodes (case 3, which agrees more with the
current practice of SNESIM), respectively, in order
Figure 2. Pluri-Gaussian truncation rule to sensitize the results to the neighborhood size. For
giving the contact relationships between FILTERSIM, a moving neighborhood with 10 pre-
four facies. Each axis represents a viously simulated nodes is considered to define the
Gaussian random field (Y1, Y2). The
data events. For each algorithm and case, one hun-
horizontal and vertical lines indicate
truncation thresholds (t1, t2, t3) dred non-conditional realizations are drawn on the
partitioning the bivariate Gaussian target grid of size 1000 9 1000.
space into facies.

posits (Armstrong et al. 2011, and references Analysis of Results

therein), and so we adopt it for this synthetic case
study. Thus, two independent stationary Gaussian Facies Proportions and their Fluctuations
random fields, each with an isotropic cubic vari-
ogram model (Table 1), are used to define the four In all the algorithms trialed in this synthetic
facies according to the truncation rule presented in case study, the prior probability of occurrence of
Figure 2. Facies 1 is in contact only with facies 2 and each facies has been set to 0.25, i.e., each facies is
3, and its contact with facies 4 is forbidden. expected (on average over many realizations) to
One hundred non-conditional realizations of the cover one quarter of the region under study. This
four facies are generated by pluri-Gaussian simula- way, the facies are equally representative and one
tion on a two-dimensional region of 1000 9 1000 avoids conditions (low proportions and lack of spa-
nodes. The two underlying Gaussian random fields tial continuity of some facies) that may complicate
are simulated by means of a spectral-turning bands the interpretation of results. Figure 3 displays the
algorithm with 5000 turning lines, which accurately cumulative facies proportions over 100 realizations
reproduces the desired distributions and variograms for all the trials. Apart from PLURISIM, only SI-
(Emery et al. 2016). Hence, the pluri-Gaussian real- SIM and SNESIM with a 1-point or 2-point neigh-
izations are assumed error-free and will constitute a borhood provide unbiased results (i.e., proportions
benchmark for comparing the nonparametric simu- that fluctuate around 0.25 for each facies), whereas
lation algorithms (SISIM, SNESIM and FILTER- the other algorithms produce realizations that
SIM), particularly in reference to the accuracy and underestimate the expected proportion of facies 1
fluctuations of the simulated statistics. and overestimate the expected proportion of facies
For non-conditional sequential indicator simu- 4. Figure 3 also shows that the proportions on all the
lation (SISIM), the mean values and the direct and nonparametric algorithms, especially the ones based
cross-variograms of the facies indicators are derived on a training image (SNESIM and FILTERSIM),
from the previous pluri-Gaussian model (Lantuéjoul fluctuate in a smaller range in comparison with the
2002; Emery 2007; Armstrong et al. 2011). The pluri-Gaussian approach, which means that the
Nonparametric Geostatistical Simulation of Subsurface Facies 1169

Figure 3. Cumulative facies proportions over 100 realizations obtained with each simulation algorithm under consideration
(synthetic case study).
1170 Madani, Maleki, and Emery

uncertainty in the facies proportions is underesti- cies. In Figure 5, the indicator cross-variograms be-
mated by all these algorithms. This can be explained tween facies 1 and 4 are displayed for all the
because the TI provides information on the spatial realizations and simulation algorithms. For SISIM and
distribution of the facies (variability of the pluri- SNESIM with a neighborhood of 1 or 2 points, the
Gaussian random field in space), but lacks infor- cross-variogram has a strictly negative slope at the
mation on the variability across realizations of the origin. This is an indication that the realizations allow
pluri-Gaussian random field. the two facies to be in contact (Beucher and Renard
2016; Maleki et al. 2017), although they should not be.
In contrast, for SNESIM with a neighborhood
Indicator Direct Variograms and their Fluctuations of 10 points and FILTERSIM, the cross-variogram is
equal to zero near the origin, up to a distance of
One may also be interested in seeing how the about 50 that indicates the minimum separation
spatial continuity varies through the realizations. As distance between facies 1 and 4. On average over the
example, the experimental indicator variogram for 100 realizations, the experimental cross-variogram
facies 1 is shown in Figure 4 for all the realizations matches the true theoretical cross-variogram derived
and simulation algorithms. The average and confi- from the pluri-Gaussian model. The fluctuations
dence limits (± 2 standard deviations around the around the average are, however, overestimated
average) of the indicator variograms calculated over when using FILTERSIM, whereas SNESIM yields
the 100 realizations are also depicted. Since the pluri- fluctuations that are comparable to PLURISIM, al-
Gaussian model is taken as the reference for setting though differences are still perceptible.
the parameters of all the other algorithms, the indi-
cator variograms are expected to bear resemblance to
those of the pluri-Gaussian case. It can be seen, Facies Maps
however, that SISIM and SNESIM with a moving
neighborhood of 1 or 2 points fail to reproduce the Two factors are deemed critical for assessing the
overall value, shape and spread (fluctuations) of the quality of the realization maps. The first one is regu-
indicator variograms. In particular, the slope of the larity of facies boundaries, and the second one is
indicator variogram at the origin turns out to be much contact relationships between facies (in particular, the
larger than the one provided by PLURISIM, which forbidden contact between facies 1 and 4). As seen in
implies that the boundaries of the facies obtained Figure 6, the pluri-Gaussian realizations reproduce the
with SISIM and SNESIM are more irregular than contact relationships defined by the truncation rule
those obtained with PLURISIM (Lantuéjoul 2002; (Fig. 2) and yield regular boundaries between the fa-
Emery and Lantuéjoul 2011; Dubrule 2017; Maleki cies. The latter can be explained by the use of cubic
et al. 2017). FILTERSIM and SNESIM with a mov- variogram models for the underlying Gaussian random
ing neighborhood of 10 points yield indicator vari- fields, which have a regular behavior near the origin
ograms that are closer to those of PLURISIM, in (Lantuéjoul 2002). In contrast, SISIM and SNESIM
particular with a comparable slope at the origin (thus, with a moving neighborhood of 1 or 2 points fail to
a comparable regularity of the facies boundaries). reproduce both the regular boundaries and the contact
However, on average over the 100 realizations, the relationships and provide patchy maps that bear little
variogram sill appears to be underestimated, which resemblance with the pluri-Gaussian maps. Only
relates to the underestimation of the facies propor- FILTERSIM and SNESIM with a moving neighbor-
tion observed in the previous subsection (recall that hood of 10 points correctly reproduce the regular
the sill of an indicator variogram is nothing but the boundaries and the contact relationships between fa-
indicator variance, which is equal to the facies pro- cies, although, in both approaches, spotted occurrences
portion multiplied by the complementary proportion) of facies 2 are overrepresented with respect to those in
(Dubrule 2017). the TI (Fig. 6).

Indicator Cross-variograms and their Fluctuations Contact Probabilities

The indicator cross-variograms are useful tools It is also of interest to calculate the probabilities
for assessing the spatial cross correlation between fa- of contact between two different facies. In Table 2,
Nonparametric Geostatistical Simulation of Subsurface Facies 1171

Figure 4. Indicator direct variograms for facies 1 obtained for each simulation algorithm. Green: experimental direct
variograms of 100 realizations. Red: average of 100 experimental direct variograms. Blue: average ± 2 standard deviations.
Black: theoretical direct variogram derived from the pluri-Gaussian model (synthetic case study).

these probabilities have been calculated by consid- forbidden contacts, the probability of contact should
ering the facies simulated at two adjacent grid nodes be zero in the corresponding entry of the table,
over the 100 realizations. Since facies 1 and 4 have which is the actually case for PLURISIM, FIL-
1172 Madani, Maleki, and Emery

Figure 5. Indicator cross-variograms between facies 1 and 4 obtained for each simulation algorithm. Green: experimental
cross-variograms of 100 realizations. Red: average of 100 experimental cross-variograms. Blue: average ± 2 standard
deviations. Black: theoretical cross-variogram derived from the pluri-Gaussian model (synthetic case study).

TERSIM and, up to a small approximation, for contact probabilities that are closest to those of
SNESIM with a moving neighborhood of 10 points. PLURISIM. In contrast, when using only 1 or 2
Overall, the last algorithm is the one that yields points in the moving neighborhood, the SNESIM
Nonparametric Geostatistical Simulation of Subsurface Facies 1173

Figure 6. Training image (right) and first four non-conditional realizations obtained with each simulation algorithm under
consideration (left and middle) (synthetic case study).
1174 Madani, Maleki, and Emery

Table 2. Contact probabilities between facies, calculated as the frequencies of the transitions between the corresponding facies for two
adjacent grid nodes over 100 realizations (synthetic case study). For each contact, the probability closest to that of PLURISIM is
highlighted in bold

Contact PLURISIM SISIM SNESIM (1 point) SNESIM (2 points) SNESIM (10 points) FILTERSIM

1 fi 2 0.329 0.275 0.432 0.407 0.342 0.306

1 fi 3 0.670 0.414 0.495 0.543 0.653 0.692
1 fi 4 0.000 0.310 0.071 0.049 0.003 0.001
2 fi 1 0.273 0.266 0.349 0.289 0.215 0.227
2 fi 3 0.368 0.402 0.321 0.341 0.344 0.344
2 fi 4 0.358 0.331 0.328 0.362 0.440 0.428
3 fi 1 0.341 0.317 0.326 0.292 0.258 0.258
3 fi 2 0.225 0.319 0.261 0.258 0.217 0.173
3 fi 4 0.432 0.362 0.412 0.449 0.524 0.568
4 fi 1 0.000 0.275 0.064 0.035 0.001 0.000
4 fi 2 0.337 0.304 0.368 0.368 0.345 0.274
4 fi 3 0.662 0.419 0.566 0.595 0.652 0.724

Figure 7. Three-dimensional view of borehole data (left) and interpretive geological model (right). Color
indicate lithological facies (real case study).

realizations fail to reproduce the contact relation- Santiago and is mined by open pit and underground
ships between facies (in particular, the absence of methods. Here, we focus on a sector of the deposit
contacts between facies 1 and 4); the same with SI- with a volume of 600 m 9 775 m 9 450 m that has
SIM. been recognized by diamond drilling. A set of 3085
samples composited at a length of 16 m are available
with information on predominant rock type, which
REAL CASE STUDY comprises three main lithological facies (Emery and
Silva 2009; Madani and Emery 2015): granitoids
Background Information (GD), tourmaline breccia (TOB) and porphyry
(POR) (Fig. 7, left). The spatial distribution of the
To further illustrate the applicability of the facies was modeled by the mining geologists in the
proposed tools, we turn to a real case study in sector of interest (Fig. 7, right). This model corre-
mineral resources evaluation. The Rı́o Blanco–Los sponds to an interpretation, with a rather rough
Bronces porphyry-Cu deposit is located in the Chi- resolution (blocks of 15 m 9 15 m 9 16 m), and a
lean Central Andes at about 70 km northeast of single facies is assigned to each block.
Nonparametric Geostatistical Simulation of Subsurface Facies 1175

SNESIM Simulation Based on Interpretive consistent with the geology of the deposit insofar as
Geological Model as TI the POR is a post-mineralization intrusive body,
younger than the GD and the tourmaline breccia.
Nonparametric simulation is used to assess the The truncation thresholds (t1 = -0.80 and t2 = 0.03)
uncertainty in the facies boundaries. For concise- are chosen in order to reproduce the facies propor-
ness, we restrict the analysis to the SNESIM algo- tions over the sector of interest (38.5% for GD,
rithm and use implementation parameters close to 40.3% for TOB, 21.2% for POR), as per the inter-
the standard practice of this algorithm, in particular pretive geological model. The variograms of the
a neighborhood of 50 points and, for conditioning Gaussian random fields are fitted in order to repro-
the realizations to the borehole data, a servo-system duce the direct and cross-variograms of the facies
weight factor of 0.5 (Mariethoz and Caers 2014). indicators. The modeled variograms consist of nes-
The interpretive geological model is held to a reso- ted spherical models with a geometric anisotropy
lution of 2.5 m 9 2.5 m 9 16 m and is used as a TI with main directions N20W, N70E and vertical.
to construct 100 realizations over a horizontal cross One SNESIM realization and one PLURISIM
section within the region of interest. realization, restricted to the area surrounding the
To obtain a reference against which to compare sampling data, are mapped and are shown in Fig-
the SNESIM realizations, a pluri-Gaussian model is ure 9, together with the interpretive geological
considered, based on two independent Gaussian model for the same sector. The facies proportions,
random fields (Y1, Y2) and two truncation thresholds the confidence limits (average ± 2 standard devia-
(t1, t2). The truncation rule (Fig. 8) is defined so that tions) of the indicator direct and cross-variograms,
the POR crosscuts the other two facies, which is and the contact probabilities obtained with 100
SNESIM and 100 PLURISIM realizations are shown
in Figures 10, 11 and Table 3, respectively.
Strong differences can be observed between the
SNESIM and PLURISIM realizations. The latter
exhibits higher short-scale variability, which is
noticeable when examining the maps, the contact
probabilities and the facies indicator direct and
cross-variograms. In particular, the steeper increase
in these variograms near the origin (Fig. 11) reflects
more irregular facies boundaries (Maleki et al.
Figure 8. Truncation rule used for pluri-
2017). With both algorithms, the proportion of GD
Gaussian modeling of facies (real case fluctuates around the ‘‘true’’ proportion calculated
study). from the interpretive geological model in the sector

Figure 9. Maps of a SNESIM realization (left), of the interpretive geological model (center) and of a pluri-Gaussian
realization (right) in the simulation area (cross section with elevation 3650 m). The interpretive geological model was used
as the training image for SNESIM. The surrounding sampling data are superimposed on each map (real case study).
1176 Madani, Maleki, and Emery

Figure 10. Cumulative facies proportions over 100 realizations obtained with PLURISIM and SNESIM (real case study;
interpretive geological model used as TI). Dotted lines indicate the proportions calculated from the interpretive geological
model in the sector of simulation (cross section with elevation 3650 m).

of simulation (33.5%); however, a slight bias can be shown in Figure 12. The facies proportions, indica-
observed in the SNESIM realizations concerning the tor direct and cross-variogram confidence limits and
proportions of TOB and POR when compared with contact probabilities obtained with 100 realizations
the interpretive geological model in the sector of are shown in Figures 13, 14 and Table 3, respec-
simulation (56.6% for TOB and 9.9% for POR). In tively.
contrast, no bias is perceptible with the PLURISIM The map of the simulated facies exhibits the
realizations. Thus, one doubts the quality of the same texture as the pluri-Gaussian TI and is con-
SNESIM realizations, which seem to provide an sistent with the facies chronology. (The younger
inaccurate reproduction of the facies statistics and POR crosscuts the other two facies.) However, a
spatial correlation (proportions, probabilities of bias is still observed in the reproduction of the facies
contacts and indicator direct and cross-variograms). proportions, indicator direct and cross-variograms
The explanation mostly lies in the TI (i.e., inter- and contact probabilities when comparing SNESIM
pretive geological model) used to construct the with the PLURISIM realizations. The bias in the
SNESIM realizations, which lacks resolution and is proportions can be explained by the differences
over-smoothed with respect to reality, as reflected between the proportions calculated within the whole
by the borehole data. volume of the TI (600 m 9 775 m 9 450 m) and
those calculated within the sector of simulation (a
2D region enclosing the sampling data). The use of a
SNESIM Simulation Based on a Pluri-Gaussian servo-system weight factor (0.5) aims to reproduce
Realization Model as TI the former proportions, but the sampling data
coerces the realizations to reproduce the latter
To make a fair comparison between SNESIM proportions. Similar arguments can be made about
and the reference model (PLURISIM), we repeat the contact probabilities and indicator variograms,
the preceding exercise by changing the TI used in for which SNESIM is torn between the TI statistics
SNESIM: Instead of the interpretive geological within the whole sector and the data statistics within
model, we use an extra pluri-Gaussian simulation the area of simulation.
calculated on the entire sector (600 m 9 775 m 9 One more point worthy to mention concerns
450 m) at a fine resolution (2.5 m 9 2.5 m 9 16 the fluctuations of the indicator direct and cross-
m). Such a TI and a resulting SNESIM realization variograms across the realizations, as reflected by
over the area surrounding the sampling data are the range of the interval between the lower and
Nonparametric Geostatistical Simulation of Subsurface Facies 1177

Figure 11. Confidence limits for the indicator direct and cross-variograms over 100 realizations obtained with PLURISIM
and SNESIM (real case study; interpretive geological model used as TI).
1178 Madani, Maleki, and Emery

Table 3. Contact probabilities between facies, calculated as the frequencies of the transitions between the corresponding facies for two
adjacent grid nodes over 100 realizations (real case study). For each contact, the probability closest to that of PLURISIM is highlighted in
bold

Contact PLURISIM SNESIM (TI = interpreted model) SNESIM (TI = pluri-Gaussian realization)

GD fi TOB 0.748 0.806 0.754

GD fi POR 0.251 0.193 0.245
TOB fi GD 0.430 0.431 0.450
TOB fi POR 0.569 0.568 0.549
POR fi GD 0.202 0.154 0.210
POR fi TOB 0.797 0.845 0.789

Figure 12. Maps of a SNESIM realization (left) and of the training image (pluri-Gaussian realization, right) on
the cross section with elevation 3650 m. The sampling data are superimposed on the SNESIM map (real case
study).

upper confidence limits (Fig. 14). At short scales, the  The indicator direct variogram provides
ranges observed with SNESIM turn out to be less information on the regularity of boundaries
than the ranges observed with PLURISIM, a situa- (variogram slope at the origin), proportion
tion that could be explained by the finite sizes of the (variogram sill) and spatial correlation struc-
moving neighborhood and of the TI (as discussed ture (variogram shape and correlation range)
below). of facies. It also allows measurement of the
fluctuations of these features across the real-
izations.
DISCUSSION  The indicator cross-variogram gives infor-
mation on whether two facies are in contact
As shown in the preceding sections, different (negative variogram slope at the origin) or
statistical and graphical tools can be used to gain not (cross-variogram equal to zero until the
insight into the properties of simulated facies: minimum separation distance between fa-
cies), and on the joint correlation structure of
 The distribution of facies proportions across these facies (variogram shape and correlation
the realizations allows assessment of biases range). Through the direct variogram, one
and inconsistent fluctuations across the real- can also measure fluctuations in these fea-
izations. tures across the realizations.
Nonparametric Geostatistical Simulation of Subsurface Facies 1179

tative of the ground truth. An over-smoothed TI

corresponding to an interpretive geological model
can provoke strong biases in the reproduction of the
facies indicator variograms and contact probabili-
ties. In any case, some biases still arise when the
statistics of the TI differ, locally, from the statistics
of the conditioning data. Furthermore, the fluctua-
tions in the facies proportions and/or in the indicator
direct and cross-variograms obtained with SNESIM
may not coincide with those obtained with the ref-
erence pluri-Gaussian model, indicating an inexact
reproduction of the uncertainty.
The discrepancies observable between SNESIM
and the reference model can be explained by the
restrictions that affect the implementation of SNE-
SIM, mainly the finite size of the moving neighbor-
hood and the finite size of the TI. Increasing the
latter so that it contains as many different patterns
Figure 13. Cumulative facies proportions over 100 realizations as possible, with a repeatability as high as possible,
obtained with SNESIM (real case study; pluri-Gaussian may prevent the simulated statistics to be globally
realization used as TI). Dotted lines indicate the proportions biased with respect to the model statistics. The TI
calculated from the interpretive geological model in the sector
should also be large in comparison with the moving
of simulation.
neighborhood in order to avoid understating statis-
tical fluctuations. To better understand this limita-
 Mapping of realizations allows qualitative
tion, imagine that the neighborhood used in
assessment of regularity of boundaries and
SNESIM has no restriction on the number of points,
contact relationships between facies.
which would be the ideal implementation of the
 The contact probabilities give quantitative
algorithm. In such a case, the SNESIM realization
information on the contact relationships (in
would necessarily be a piece of the TI so that only a
particular, about forbidden contacts) between
finite number of realizations could be drawn (Emery
facies.
and Lantuéjoul 2014; Maleki et al. 2017) (i.e., only a
The experiments in the synthetic case study consider finite number of pieces can be extracted from the
the pluri-Gaussian model as the reference case and TI) and fluctuations across the realizations would
attempt to determine whether or not a given consequently be understated. Accordingly, a trade-
nonparametric algorithm (SISIM, SNESIM or FIL- off in the SNESIM implementation emerges:
TERSIM) correctly reproduces the main features Statistics are poorly reproduced when the neigh-
imposed to the pluri-Gaussian realizations. Based on borhood contains too few conditioning points, while
the aforementioned tools, the results are shown in fluctuations across realizations may be incorrect
Table 4. SNESIM implemented with a large moving when the neighborhood contains too many points
neighborhood turns out to provide the closest and/or when the TI is (too) small. An alternative
features to that of the pluri-Gaussian model, fol- solution could be the use of several TI (ideally, as
lowed by FILTERSIM, although biases in some many TIs as desired realizations).
simulation statistics or inexact fluctuations across Other statistical tools may be used to comple-
the realizations are perceptible with both algo- ment the presented analysis, involving multiple-
rithms. In contrast, SISIM and SNESIM imple- point metrics (Boisvert et al. 2010; De Iaco and
mented with an over-restricted neighborhood yield Maggio 2011; De Iaco 2013) or related to the
the poorest results so that their use should be geometry of the connected components (Renard and
discouraged. Allard 2013; Oriani and Renard 2014; Rongier et al.
The real case study corroborates the finding 2016) of the simulated facies. These tools provide a
that SNESIM with a large moving neighborhood richer description of the realizations but, as men-
globally reproduces the features of the pluri-Gaus- tioned in the Introduction, cannot be calculated
sian realizations, provided that the TI is represen- based on scattered borehole data only; their use
1180 Madani, Maleki, and Emery

Figure 14. Confidence limits for the indicator direct and cross-variograms over 100 SNESIM realizations (real case study;
pluri-Gaussian realization used as TI).
Nonparametric Geostatistical Simulation of Subsurface Facies 1181

Table 4. Summary of whether or not the nonparametric simulation algorithms under study reproduce the features of the reference pluri-
Gaussian model (synthetic case study)

SISIM SNESIM SNESIM SNESIM FILTERSIM

(1 point) (2 points) (10 points)

Contact relationships No No No Yes Yes

Boundary regularity No No No Yes Yes
Average facies proportions Yes Yes Yes No No
Fluctuations in facies proportions No No No No No
Average spatial correlation structure No No No Mixed Mixed
(indicator direct and cross-variograms)
Fluctuations in spatial correlation structure No No No Mixed No

require other sources of information, such as geo- particular, multiple-point simulation algorithms like
physical observations, flow and transport responses SNESIM may fail to reproduce the expected statis-
and dynamic data (e.g., analysis of borehole tests, tical fluctuations if a too small TI or a too large
production data, pressure for detecting the presence neighborhood is selected for its practical imple-
of permeable pathways between boreholes) (Renard mentation. A solution for this could be the use of a
and Allard 2013; Chautru et al. 2015). very large TI containing as many different patterns
as possible, or of several TI (ideally, as many TIs as
desired realizations).
CONCLUSIONS

This research argues in favor of two simulation

algorithms based on multiple-point statistics, SNE- ACKNOWLEDGMENTS
SIM and FILTERSIM, with a moving neighborhood
containing a large number of points (at least 10), as The first author acknowledges the Nazarbayev
the ones that are more capable of realistically University for financial supporting. The second and
reproducing the regularity of facies boundaries, the third authors acknowledge the support of the Chi-
contact relationships and the spatial correlation lean Commission for Scientific and Technological
structure of facies indicators; however, they may still Research (CONICYT), through projects CONIC
yield slightly biased reproduction of facies propor- YT/FONDECYT/POSTDOCTORADO/N3180655
tions and/or indicator variograms. In contrast, the and CONICYT PIA Anillo ACT1407, respectively.
two-point simulation algorithm (SISIM) and the The examples presented in ‘‘Essential Statistical and
multiple-point (SNESIM) simulation algorithms Geometrical Characteristics of Facies’’ section were
based on an over-restricted neighborhood are cap- provided by Dr. Bijan Biranvand from the Depart-
able of reproducing the facies proportions but not ment of Petroleum Geology, Research Institute of
geometry of the facies boundaries and the spatial Petroleum Industry, Tehran, Iran. The data for the
correlation structure of facies indicators. case study were provided by CODELCO—Chile.
A warning should, however, be issued for all the Constructive comments by two anonymous review-
algorithms regarding possible inconsistencies about ers are gratefully acknowledged. Also, we deeply
the fluctuations of realization statistics across the thank Dr. Carranza for the valuable comments on
realizations. Even if the statistics are correctly the final version of manuscript.
reproduced on average over the realizations, the
previous global and local experimental analyses over
the simulation results through both synthetic and
real case studies suggest that their statistical fluctu- REFERENCES
ations may not be the correct ones. This is true with
SISIM, SNESIM and FILTERSIM in regard to fa- Alabert, F. (1987). Stochastic imaging of spatial distributions
cies proportions and to spatial correlation structure using hard and soft information. MasterÕs thesis, Department
of Applied Earth Sciences, Stanford University, Stanford, pp.
(direct and cross-variograms) of facies indicators. In 332.
1182 Madani, Maleki, and Emery

Al-Mudhafar, W. J. (2018). Multiple-point geostatistical lithofa- Lantuéjoul, C. (1994). Non conditional simulation of stationary
cies simulation of fluvial sand-rich depositional environment: isotropic multigaussian random functions. In M. Armstrong
a case study from Zubair formation/South Rumaila oil field. & P. A. Dowd (Eds.), Geostatistical Simulations (pp. 147–
SPE Reservoir Evaluation & Engineering, 21(1), 39–53. 177). Dordrecht: Kluwer Academic.
Armstrong, M., Galli, A., Beucher, H., Le LocÕh, G., Renard, D., Lantuéjoul, C. (2002). Geostatistical simulation, models and
Renard, B., et al. (2011). Plurigaussian Simulations in Geo- algorithms (p. 256). Berlin: Springer.
sciences (p. 187). Berlin: Springer. Le LocÕh, G., & Galli, A. (1997). Truncated plurigaussian method:
Beucher, H., & Renard, D. (2016). Truncated Gaussian and de- Theoretical and practical points of view. In E. Y. Baafi & N.
rived methods. Comptes Rendus Géoscience, 348, 510–519. A. Schofield (Eds.), Geostatistics WollongongÕ96 (pp. 211–
Boisvert, J. B., Pyrcz, M. J., & Deutsch, C. V. (2007). Multiple- 222). Dordrecht: Kluwer Academic.
point statistics for training image selection. Natural Resources Leuangthong, O., McLennan, J. A., & Deutsch, C. V. (2004).
Research, 16(4), 313–321. Minimum acceptance criteria for geostatistical realizations.
Boisvert, J. B., Pyrcz, M. J., & Deutsch, C. V. (2010). Multiple Natural Resources Research, 13(3), 131–141.
point metrics to assess categorical variable models. Natural Lowell, J. D., & Guilbert, J. M. (1970). Lateral and vertical
Resources Research, 19, 165–175. alteration-mineralization zoning in porphyry ore deposits.
Chautru, J. M., Meunier, R., Binet, H., & Bourges, M. (2015). Economic Geology, 65, 373–408.
Geobodies stochastic analysis for geological model parame- Madani, N., & Emery, X. (2015). Simulation of geo-domains
ter inference. Petroleum Geostatistics 2015 (pp. 293–297). accounting for chronology and contact relationships: Appli-
Houten: European Association of Geoscientists & Engineers. cation to the Rı́o Blanco copper deposit. Stochastic Envi-
Chilès, J. P., & Delfiner, P. (2012). Geostatistics: Modeling spatial ronmental Research and Risk Assessment, 29(8), 2173–2191.
uncertainty (p. 699). New York: Wiley. Madani, N., & Emery, X. (2017). Plurigaussian modeling of geo-
De Iaco, S. (2013). On the use of different metrics for assessing logical domains based on the truncation of non-stationary
complex pattern reproduction. Journal of Applied Statistics, Gaussian random fields. Stochastic Environmental Research
40(4), 808–822. and Risk Assessment, 31(4), 893–913.
De Iaco, S., & Maggio, S. (2011). Validation techniques for geo- Madani, N., Naderi, A., Biranvand, B., & Keshavarz, N. (2018).
logical patterns simulations based on variogram and multiple- Lithofacies uncertainty modeling in a siliciclastic reservoir
point statistics. Mathematical Geosciences, 43, 483–500. setting by incorporating geological contacts and seismic
Dubrule, O. (2017). Indicator variogram models: do we have information. Journal of Petroleum Exploration and Produc-
much choice? Mathematical Geosciences, 49(4), 441–465. tion Technology. https://doi.org/10.1007/s13202-018-0531-7.
Emery, X. (2004). Properties and limitations of sequential indi- Maleki, M., Emery, X., & Mery, N. (2017). Indicator variograms
cator simulation. Stochastic Environmental Research and Risk as an aid for geological interpretation and modeling of ore
Assessment, 18(6), 414–424. deposits. Minerals, 7(12), 241.
Emery, X. (2007). Simulation of geological domains using the Mariethoz, G., & Caers, J. (2014). Multiple-point geostatistics:
plurigaussian model: New developments and computer pro- Stochastic modeling with training images (p. 376). New York:
grams. Computers & Geosciences, 33(9), 1189–1201. Wiley.
Emery, X. (2008). Statistical tests for validating geostatistical Matheron, G. (1989). Estimating and choosing: An essay on
simulation algorithms. Computers & Geosciences, 34(11), probability in practice (p. 141). Berlin: Springer.
1610–1620. Matheron, G., Beucher, H., Galli, A., Guérillot, D., & Ravenne,
Emery, X., Arroyo, D., & Porcu, E. (2016). An improved spectral C. (1987). Conditional simulation of the geometry of fluvio-
turning-bands algorithm for simulating stationary vector deltaic reservoirs. In: 62nd Annual technical conference and
Gaussian random fields. Stochastic Environmental Research exhibition of the society of petroleum engineers, pp. 591–599.
and Risk Assessment, 30, 1863–1873. SPE Paper 16753, Dallas.
Emery, X., & Lantuéjoul, C. (2006). TBSIM: A computer pro- Modis, K., & Sideri, D. (2013). Geostatistical simulation of
gram for conditional simulation of three-dimensional Gaus- hydrofacies heterogeneity of the West Thessaly aquifer sys-
sian random fields via the turning bands method. Computers tems in Greece. Natural Resources Research, 22(2), 123–138.
& Geosciences, 32(10), 1615–1628. Oriani, F., & Renard, P. (2014). Binary upscaling on complex
Emery, X., & Lantuéjoul, C. (2011). Geometric covariograms, heterogeneities: The role of geometry and connectivity. Ad-
indicator variograms and boundaries of planar closed sets. vances in Water Resources, 64, 47–61.
Mathematical Geosciences, 43(8), 905–927. Renard, P., & Allard, D. (2013). Connectivity metrics for sub-
Emery, X., & Lantuéjoul, C. (2014). Can a training image be a surface flow and transport. Advances in Water Resources, 51,
substitute for a random field model? Mathematical Geo- 168–196.
sciences, 46(2), 133–147. Rongier, G., Collon, P., Renard, P., Straubhaar, J., & Sausse, J.
Emery, X., & Ortiz, J. M. (2011). A comparison of random field (2016). Comparing connected structures in ensemble of ran-
models beyond bivariate distributions. Mathematical Geo- dom fields. Advances in Water Resources, 96, 145–169.
sciences, 43(2), 183–202. Rossi, M. E., & Deutsch, C. V. (2014). Mineral resource estimation
Emery, X., & Silva, D. A. (2009). Conditional co-simulation of (p. 332). New York: Springer.
continuous and categorical variables for geostatistical appli- Strebelle, S. (2002). Conditional simulation of complex geological
cations. Computers & Geosciences, 35(6), 1234–1246. structures using multiple-point statistics. Mathematical
Galli, A., Beucher, H., Le LocÕh, G., & Doligez, B. (1994). The Geology, 34(1), 1–22.
pros and cons of the truncated Gaussian method. In M. Tan, X., Tahmasebi, P., & Caers, J. (2014). Comparing training-
Armstrong & P. A. Dowd (Eds.), Geostatistical simulations image based algorithm using an analysis of distance. Mathe-
(pp. 217–233). Dordrecht: Kluwer. matical Geosciences, 46(2), 149–169.
Journel, A. G., & Alabert, F. G. (1990). New method for reservoir Yunsel, T., & Ersoy, A. (2011). Geological modeling of gold de-
mapping. Journal of Petroleum Technology, 42(2), 212–218. posit based on grade domaining using plurigaussian simula-
Journel, A. G., & Gómez-Hernández, J. J. (1993). Stochastic tion technique. Natural Resources Research, 20(4), 1–19.
imaging of the Wilmington clastic sequence. SPE Formation Zhang, T., Switzer, P., & Journel, A. (2006). Filter-based classi-
Evaluation, 8(1), 33–40. fication of training image patterns for spatial simulation.
Mathematical Geology, 38(1), 63–80.