Chapter 30

Multiple Point Statistics: A Review

Pejman Tahmasebi

Abstract Geostatistical modeling is one of the most important tools for building an
ensemble of probable realizations in earth science. Among them, multiple-point
statistics (MPS) has recently gone under a remarkable progress in handling complex
and more realistic phenomenon that can produce large amount of the expected
uncertainty and variability. Such progresses are mostly due to the recent increase in
more advanced computational techniques/power. In this review chapter, the recent
important developments in MPS are thoroughly reviewed. Furthermore, the
advantages and disadvantages of each method are discussed as well. Finally, this
chapter provides a brief review on the current challenges and paths that might be
considered as future research.

30.1 Introduction

Characterization and modeling of geological structures have been investigated for

several years in geosciences. Geostatistics is one of the such methods that can be
used to analyze the data effectively. Such analysis can be performed both spatially
and temporally. Lack of data is one of the intrinsic issues in the earth science
applications, which causes a significant uncertainty and ambiguity in these prob-
lems. Kriging, as one of the most widespread geostatistical tools, was developed for
dealing with such problems. The basic mathematically equations of Kriging, after
developing by Daniel Krige, was further advanced by Matheron (Journel and
Huijbregts 1978; Matheron 1973). Kriging is a deterministic method, meaning that
it only produces one outcome from the available sparse data, which intrinsically
cannot be used to effectively quantify the uncertainty. This method requires a prior
model of variability and correlation between the variables, known as the variogram

P. Tahmasebi (✉)
Department of Petroleum Engineering, University of Wyoming, Laramie,
WY 82071, USA
e-mail: [email protected]

© The Author(s) 2018 613

B. S. Daya Sagar et al. (eds.), Handbook of Mathematical Geosciences,
https://doi.org/10.1007/978-3-319-78999-6_30
614 P. Tahmasebi

(Chiles and Delfiner 2011; Cressie and Wikle 2011; Deutsch and Journel 1998;
Goovaerts 1997; Kitanidis 1997).
It has been shown that Kriging produces excessively smooth results (Deutsch
and Journel 1998; Journel and Zhang 2006) and it cannot represent the hetero-
geneity and non-smooth phenomena. One consequence is the underestimation and
overestimation for low and high values, respectively. This problem becomes evi-
dent when important parameters such as water breakthrough is intended to be
predicted. Thus, the results of Kriging cannot be used for these situations as they
ignore the connectivity and variability.
Stochastic simulation can be used to overcome the limitations of Kriging
(Goovaerts 1997; Journel and Huijbregts 1978). Several simulation methods have
been proposed that can produce various equi-probable realizations. Methods such as
sequential Gaussian simulation (SGSIM) and sequential indicator simulation
(SISIM) have become popular among different fields of earth sciences. These
methods, give a number of “realizations” or interpolation scenarios, which allow
assessing the uncertainty and quantifying it more accurately. It should be noted that
Kriging is still the main algorithm used in the above stochastic methods. An example
for the application of Kriging and stochastic modeling is provided in Fig. 30.1.
Due to relying on variogram (i.e. covariance), kriging-based geostatistical sim-
ulations are not able to reproduce complex patterns. Clearly, considering only two
points is not sufficient for reproducing complex and heterogeneous models. Thus,
several attempts in the recent years in the context of multiple point geostatistics
(MPS) have been made that can use more than two points simultaneously. Using the
information from multiple points require a big source of data, which is not usually
available in the earth science problems as they come with sparse and incomplete
data. Such data, instead, can be browsed in the form a conceptual image, called
training image (TI).

Fig. 30.1 Comparison between the results of Kriging (b) and stochastic simulation (c) using
conditioning point data in (a)
30 Multiple Point Statistics: A Review 615

Technically, geostatistical methods can be divided into three main groups.

Object-based (or Boolean) simulation methods are in the first group (Kleingeld
et al. 1997). These methods consider the medium as a group of stochastic objects
that are defined based on a specific statistical distribution (Deutsch and Wang 1996;
Haldorsen and Damsleth 1990; Holden et al. 1998; Skorstad et al. 1999).
Pixel-based methods are considered in the second group. These methods are
based a set of points/pixels that represent various properties of a phenomenon.
Mathematically speaking, such methods vary from the LU decomposition of the
covariance matrix (Davis 1987), sequential Gaussian simulation (Dimitrakopoulos
and Luo 2004), frequency- domain simulation (Borgman et al. 1984; Chu and
Journel 1994), simulated annealing (Hamzehpour and Sahimi 2006), and the
genetic algorithm. The last two methods, namely optimization techniques, also
belong to this group as they gradually change an earth model in a pixel-by-pixel
manner.
Each of the above methods has some advantages and limitations. For example,
geological structures can be reproduced accurately using the object- based simu-
lations. However, conditioning in these methods to well and soft data require
intensive computation.
Pixel-based methods simulation on one pixel at a time. Such techniques produce
the conditioning point data exactly. One drawback of these methods is that they are
based on variograms that represent two-point statistics and, thus, they cannot
reproduce the complex and realistic geological structures. Consequently, the gen-
erated models using these techniques cannot represent an accurate representation of
any physics-based simulations (e.g. flow, grade distribution, contaminate fore-
casting and etc.).
In the MPS methods, the spatial statistics are not either extracted using vari-
ogram, but a conceptual tool named training image (TI), which is an example of the
spatial structure to be reproduced, is provided that can represent the necessary data.
During the recent years, several MPS methods have been developed to address
issues related to CPU time and improved graphical representation of the models
produced. This chapter, thus, reviews the existing concepts in MPS and discusses
the available methods. The main two-point based stochastic simulation methods are
first reviewed. Then, the basic terminologies and concepts of MPS are demon-
strated. Next, different MPS methods are explained and the advantages and dis-
advantages associated with each method are demonstrated. Finally, some avenues
for future research are discussed.

30.2 Two-Point Based Stochastic Simulation

The smoothing effect of Kriging can be avoided using the sequential simulation,
which helps to quantify the uncertainty accurately. Consider a set of N random
variables Z ðuα Þ, α = 1, . . . , N defined at locations uα . The aim of sequential sim-
ulation is to produce realizations fzðuα Þ, α = 1, . . . , N g, conditioned to n available
616 P. Tahmasebi

data and reproducing a given multivariate distribution. For this aim, the multivariate
distribution is decomposed into a set of N univariate conditional cumulative dis-
tribution functions (ccdfs):

Fðu1 , . . . , uN ; z1 , . . . , zN jðnÞÞ = Fðu1 ; z1 jðnÞÞ ×

Fðu2 ; z2 jðn + 1ÞÞ × ⋯ ×
ð30:1Þ
FðuN − 1 ; zN − 1jðn + N − 2ÞÞ ×
FðuN ; zN jðn + N − 1ÞÞ

where FðuN ; zN jðn + N − 1ÞÞ = Prob fZ ðuN Þ ≤ zN jðn + N − 1Þg is the conditional
ccdf of Z ðuN Þ conditioned to a set of n original data and ðN − 1Þ previously sim-
ulated values.

30.2.1 Sequential Gaussian Simulation (SGSIM)

In this method, the multivariate distribution and the higher order are constructed
based on the lower order statistical such as histogram and variogram. In other
words, the mean and covariance matrix are used to build a Gaussian function.
Therefore, along a random path, the mean and variance of the Gaussian distribution
is estimated via Kriging and Kriging variance. The overall algorithm of SGSIM can
be summarized as follows. First, a random path is defined over all visiting points on
the simulation grid. Then, the ccdf at each node based on the hard data and pre-
viously simulated data are considered in Kriging. Then, a random value from the
obtained Gaussian ccdf is drawn and added to the simulation grid. Next, based on
the predefined random path, another node is chosen and simulated. Finally, another
realization can be generated using a different random path.
It is worth noting that the conditioning data should be normally distributed. If it
is not the case, it entails transforming them into a Gaussian distribution in order to
be useable for SGSIM. Finally, the results must be back-transferred at the end of
simulation. Such transformations can be accomplished using normal-score trans-
forms or histogram anamorphous through Hermite polynomials.

30.2.2 Sequential Indicator Simulation (SISIM)

Indicator simulation follows the same principle as SGSIM. This method, however,
is suited for categorical data, which do not have an order relationship. Typical
examples in earth science are rock type, lithology codes and some other categorical
properties. The similar sequential procedure based on the estimation of the ccdf
conditioning to neighboring data is applied here as well. This algorithm is based on
two-point indicator variograms, which represent the spatial variability of each
30 Multiple Point Statistics: A Review 617

category. An indicator variable is defined for each variable, equal to 1 if at location

u a particular category is found, and zero otherwise. Also, E fI ðuÞg = p is the
stationary proportion of a given category. The indicator variogram can be expressed
as:

ProbfI ðuÞ = 1, I ðu + hÞ = 1g
= E fI ðuÞI ðu + hÞg ð30:2Þ
= Prob fI ðuÞ = 1jI ðu + hÞ = 1g

Usually, the categorical variables expressed as a set of K discrete categories that

zðuÞf0, . . . , k − 1g. Therefore, the indicator value for each of the defined classes
can be expressed as:


1 Z ðuÞ = k
I ðu, kÞ = ð30:3Þ
0 otherwise

The aim of the indicator formulation is to estimate the probability of Z ðuÞ to be

less than the predefined threshold for a category conditional to the data (n) retained:

I * ðu, zk Þ = E* ðI ðu, zk ÞjðnÞÞ

ð30:4Þ
= Prob* ðZ ðuÞ < zk ðnÞÞ

We can rewrite the above equation for categorical variables by using simple
Kriging as:
n
I * ðu, zk Þ − E ðI ðu, kÞÞ = ∑ λα ðuÞðI ðuα , kÞ − EfI ðuα , kÞgÞ
α=1
K
ð30:5Þ
∑ I * ðu, k Þ = 1
k=1

where E fI ðu, kÞg is the marginal probability for category k.

The above formulation can be applied within the sequential scheme which
known as SISIM. Indicator Kriging (IK) is used to estimate the probability of each
category. This algorithm can be described as follow. Similarly, as SGSIM, a ran-
dom path is defined by which all of the nodes are visited. Then, using Simple
Kriging, the indicator random variable for each category is estimated for each node
on the random path based on the neighboring data. Next, the conditional probability
density function (cpdf) is obtained and a value is randomly drawn from that cpdf
and assigned to the simulated node. This procedure is repeated sequentially for all
the visiting nodes until the simulation grid is completed. By choosing another
random path, one can generate another realization. More information on this
method can be found in Goovaerts (1997).
618 P. Tahmasebi

30.3 Multiple Point Geostatistics (MPS)

One of the bottlenecks in the two-point based geostatistical simulations is their

inability in dealing with complex and heterogeneous spatial structures. Such
methods cannot fully reproduce the existing physics and most of their parameters
usually do not have an equivalent in the reality. In particular, these methods cannot
convey the connectivity and variability when the considered phenomenon contains
definite patterns or structures. For example, models containing regular structures
cannot be reproduced using the SGSIM method. Thus, increasing the number of
points can help reproducing the connectivities and complex features. The MPS
methods, indeed, intend to reproduce the physics in natural phenomena and they all
are based on a set of training images. Below, some preliminary concepts are first
reviewed.

30.3.1 Training Image

Training image (TI) is one of the most important inputs in the MPS techniques.
Thus, providing a representative TI, or a set of TIs, is the biggest challenge in the
MPS applications. In general, TIs can be generated using the physics derived from
process-based methods or statistical methods or by using the extracted and observed
rules for each geological system. The TI can be of any type, ranging from an image
to statistical properties in space and time. In fact, TIs let us to include subjectivity in
the geological modeling, as they are difficult to be taken into account in the tra-
ditional statistical methods. In a broader sense, TI can be constructed based on the
traditional statistical methods. These outcomes, however, do not represent
the deterministic aspects of geological models, as they usually tend to signify the
randomness fragment. Geologically speaking, most of the images in natural sci-
ences represent some degree of complexity and uniqueness. Some examples of the
available TIs are shown in Fig. 30.2.
The available methods for constructing the TIs are divided into three main
groups:
• Outcrop Data: An example of TI is the outcrop images, which are one of the
preliminary sources of information at the first step of geological modeling. They
provide a unique and direct representation of geological features. They also
provide a clear illustration of the geometry and spatial continuity that allow
visual inspection of the existing structures in 2D sections.
• Object-based Methods: An alternative for constructing structured categorical
models is the object-based (or Boolean) method (Deutsch and Wang 1996;
Haldorsen and Damsleth 1990; Holden et al. 1998; Lantuéjoul 2002; Skorstad
et al. 1999). These methods are defined based on some shape parameters (e.g.
size, direction, and sinuosity). The results can be used within an iterative
30 Multiple Point Statistics: A Review 619

Fig. 30.2 a Wagon Rock Caves outcrop (Anderson et al. 1999), b digitized outcrop driven from
(a), c Herten gravel pit (Bayer et al. 2011), d litho and hydrofacies distribution extracted from (c),
e a 3D object-based model (Tahmasebi and Sahimi 2016a), f some 2D section of the 3D model
shown in (e), g a 2D model generated using the process-based techniques (Tahmasebi and Sahimi
2016a), h a 3D model generated by the process-based methods (Tahmasebi and Sahimi 2016a)
620 P. Tahmasebi

algorithm to provide any further alterations. The results of object-based simu-

lation methods are one of the best and most accessible sources for TIs.
• Process-based Methods: Process-based methods (Biswal et al. 2007, 1999;
Bryant and Blunt 1992; Gross and Small 1998; Lancaster and Bras 2002; Pyrcz
et al. 2009; Seminara 2006) try to develop 3D models by mimicking the
physical processes that form the porous medium. Though realistic, such meth-
ods are, however, computationally expensive and require considerable calibra-
tions. Moreover, they are not general enough, because each of them is
developed for a specific type of formation, as each type is the outcome of some
specific physical processes.

30.4 Simulation Path

Geostatistical techniques are conducted on a simulation grid G, which is con-

structed on several cells. These cells are visited in diverse ways on a predefined
path, either in random or in structural manner (i.e. raster path).

30.4.1 Random Path

Random path is one of the most commonly used visiting path in sequential sim-
ulation algorithms. In this particular path, a series of random number equal to the
number of unknown cells, based on a random seed, is generated for each realization
and the unvisited points on G are simulated accordingly. Clearly, the number of
simulated (i.e. known) points increase as the simulation proceeds. Each realization
is generated using a simulation path. These paths commonly come with unbi-
asedness around the conditioning point data.

30.4.2 Raster Path

Algorithms based on raster path are popular in the stochastic modeling. These paths
are constructed based on structural 1D path, meaning that the simulation cells are
visited systematically and one can predict the future visiting points. Daly (2005)
presented a Monte Carlo algorithm that utilized raster path. Then, patch-based
algorithm was used based on this path by El Ouassini et al. (2008). Next, Parra and
Ortiz (2011) used a similar path in their study. Finally, Tahmasebi et al. (2014,
2012a, b) implemented a raster path along a fast similarity computation and
achieved high-quality realizations. Such paths usually produced high quality real-
izations that can barely be produced using the random path algorithms.
30 Multiple Point Statistics: A Review 621

One of the advantages in using such paths is the small number of constraints that
help the algorithms to better identify the matching data (or patterns). For instance,
one only deals with 1–2 overlap regions in 2D simulations, which is much more
efficient when four overlaps are used in the random path algorithms. Thus, one
should expect more discontinuities and artefacts when then number of overlaps are
increased. Indeed, identifying a pattern from TI based on four constraints is very
difficult, if not impossible. Therefore, using small number of overlaps is desirable as
they result in high-quality realizations. Raster path algorithms offer such a prospect
and one can achieve realizations with higher quality.
Dealing with conditioning data (e.g. point and secondary data) is one of the
crucial issues in these paths. They, in fact, cannot account for the conditioning data
that are ahead of them. Therefore, some biases have been observed in these algo-
rithms, particularly around the conditioning point data. Some complementary
methods such as template splitting (Tahmasebi et al. 2012a) and co-template (Parra
and Ortiz 2011; Tahmasebi et al. 2014) have addressed this issue partially.

30.4.3 Some Other Definitions

Simulation Grid (G): a 2D/3D computation grid on which the geostatistical mod-
eling is performed and is composed of several cells, depending on the size of
domain and simulation. It contains no information for unconditional simulation,
while the hard data are distributed in their corresponding cells.
Data-Event: a set of points that are characterized by a distance, namely lag,
which are considered around a visiting point (cell) on G.
Template: a set of points that are organized systematically and used for finding
similar patterns in TI.

30.5 Current Multiple Point Geostatistical Algorithms

Generally, the MPS methods have been developed in both pixel- and pattern-based
states, each of which, as discussed, have similar pros and cons. For example, the
pixel-based MPS methods can perfectly match the well data, whereas, these
methods, in some complex geological models, produce unrealistic structures. On
the other hand, pattern-based techniques bring a more accurate representation of the
subsurface model, while they usually miss the conditioning data. The pattern-based
methods simulate a group of points at a time. Currently, these techniques are under
different progress, due to their ability for simultaneous reproduction of conditioning
data and geologically realistic structures. As mentioned, conditioning to well data is
one of the critical issues in the pattern-based techniques. Thus, taking advantage of
the capabilities of both pixel- and pattern-based techniques in the MPS methods
622 P. Tahmasebi

through the hybrid frameworks will result in an efficient algorithm. Such a com-
bination is reviewed thoroughly in this chapter as well.
Most of the available MPS methods can be used with non-stationary systems, the
ones in which the statistical properties of a region is different from other parts
(Chugunova and Hu 2008; Honarkhah and Caers 2012; Mariethoz et al. 2010;
Strebelle 2012; Tahmasebi and Sahimi 2015a; Wu et al. 2008).

30.5.1 Pixel-Based Algorithms

i. Extended Normal Equation Simulation (ENESIM)

The ENESIM is the first method wherein the idea of MPS was raised (Guardiano
and Srivastava 1993). This method is based on an extended concept of indicator
kriging, which allows reproduction of multiple-event inferred from a TI. It first
finds the data even at each visiting point and then scans the TI for identifying all
occurrences. Then, a conditional distribution for all the identified occurrences is
constructed. Next, a sample from the generated histogram is drawn and placed in
the visiting point on G. One of the main drawbacks of this algorithm is scanning the
TI for each visiting point, which makes it unpractical for large G and TI. This
algorithm was later redesigned in the SNESIM algorithm by aid of search tree so
one does not need to rescan the TI for each visiting point, but it can be done once
before the simulation begins. Some of the results of this algorithm are presented in
Fig. 30.3.

ii. Simulated Annealing

Simulated annealing (SA) is one of the popular methods in optimization that is used
to the global minima. Suppose E represent the energy:
n      2
E = ∑ f O j − f Sj ð30:6Þ
j=1

where Oj and Sj represent the observed (or measured) and the corresponding
simulated (calculated) properties of a porous medium, respectively, with n being the
number of data points. If there are more than one set of data for distinct properties
of the medium, the energy E is generalized to
m
E = ∑ ωi E i ð30:7Þ
i=1
30 Multiple Point Statistics: A Review 623

Fig. 30.3 The results of the ENESIM algorithm. a cross-bedded sand, which is used as TI, b one
realization generated using ENESIM, c TI: fractured model generated using physical rock
propagation, d one conditional realization based on the TI in (c) and 200 conditioning data points.
The results are browed from Guardiano and Srivastava (1993)

where Ei is the total energy for the data set i, and ωi the corresponding weight, as
two distinct set of data for the same porous medium do not usually have the same
weight or significance.
An initial guess is usually considered as the structure of medium by which the
algorithm can start. Then, a small perturbation is made on the initial model and the
new energy E′ and the difference ΔE = E′ − E, are then computed. Based on a
probability this interchange is then accepted. The interchange is then accepted with
a probability pðΔE Þ. Then, according to the Metropolis algorithm,


1, ΔE ≤ 0,
pðΔE Þ = ð30:8Þ
exp ð − ΔE ̸ T Þ ΔE > 0,

where T is a fictitious temperature.

Based on statistical mechanics, it is well known that the equilibrium state of
ground state can be achieved when it is heated up to a high temperature T and then
slowly cooled down to absolute zero. The cooling is usually considered slow to
allow the system to reach its true equilibrium state. It, indeed, allows the systems to
not trap in a local energy minimum. At each annealing step i the system is allowed
to evolve long enough to “thermalize” at T ðiÞ. Then, the temperature T is decreased
624 P. Tahmasebi

and continues until the true ground state of the system is reached. This process stops
when the E is deemed to be small enough. This method has been used widely for
reconstruction of fine-scale porous media by Yeong and Torquato (1998a, b),
Manwart et al. (2000) and Sheehan and Torquato (2001), as well as large-scale oil
reservoirs.
Using this framework, the algorithm starts on a spatially random distribution for
the simulation grid G. It should be noted that the hard data are placed in the G at the
same time. Afterwards, the simulation cells are visited and the energy of realization
(i.e. global energy) is calculated using the terms considered in the objective func-
tion. Then, the probability of acceptance/rejection is calculated and the new value
will be used/ignored accordingly. Consequently, the objective function will be
updated and next T can be defined afterwards. This process continues until the
predefined stopping criteria are meet.
Deutsch (1992) used this algorithm to reproduce some MPS properties. He
considered an objective function that satisfies some constrains such as histogram
and variogram. Furthermore, some researchers applied simulated annealing for
simulation of continuous variables (Fang and Wang 1997). However, simulated
annealing has drawbacks, a major one being CPU time. Therefore, one can only
consider a limited number of statistics as constrains, because increasing the number
of constrains has a strong effect on CPU time. In addition, this algorithm has many
parameters which should be tuned and therefore need a large amount of trial and
error to achieve optimal values. Peredo and Ortiz (2011) used speculative parallel
computing to accelerate the simulated annealing; however the computation times
are still far from what is obtained with sequential simulation methods (Deutsch and
Wen 2000). A overall comparison between the SA algorithm and the traditional
algorithms is presented in Fig. 30.4. In a similar fashion, the multiple point sta-
tistical methods have also used the effect of iterations on removing the artifacts
(Pourfard et al. 2017; Tahmasebi and Sahimi 2016a). It should be noted that the
sequential algorithms can be parallel using different strategies which are not dis-
cussed in this review chapter (Rasera et al. 2015; Tahmasebi et al. 2012b).
iii. Markov Random Field (MRF)

These models incorporate constraints by formulating high-order spatial statistics

and enforcing them on the simulated domain using a Metropolis-Hastings algo-
rithm. In this case, the computational problem of the previous methods remains
because the Metropolis-Hastings algorithm, although always converging in theory,
may not converge in a reasonable time. The model parameters are inferred from the
available data, namely TI.
The Markovian properties are usually expressed as a conditional probability:

pðZjall previous ZÞ = pðz1 Þpðz2 jz1 Þ . . . pðzN jzN − 1 , zN − 2 , . . . , z2 , z1 Þ ð30:9Þ

|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
pðzN jZΦN Þ
30 Multiple Point Statistics: A Review 625

TI SA realization

SGESIM SISIM

SA

TI SA

Fig. 30.4 A comparison between the results of the SA algorithm and traditional two-point based
geostatistical simulations (Deutsch 1992). It should be noted that the results of the SGESIM,
SISIM and SA algorithms are generated based on the TI shown in Fig. 30.3a. The last raw is
browed from Peredo and Ortiz (2011)

where ZΦN indicates the conditional probability of zN and pðzN Þ > 0∀zN .
Fully utilizing the MRF algorithm for large 3D simulation grids in earth science
is not practical. Thus, researchers have focused on less computationally demanding
626 P. Tahmasebi

Fig. 30.5 An illustration of the MMM method. The gray cells represent the unvisited points. The
neighborhood is shown in a red polygon. This figure is taken from Stien and Kolbjørnsen (2011)

algorithms such as Markov Mesh Models (MMM) (Daly 2005; Stien and
Kolbjørnsen 2011; Toftaker and Tjelmeland 2013). In this algorithm, the simulation
is only restricted to a reasonable small window around the visiting point, see
Fig. 30.5. Thus, Eq. (9) can be shorten as:

pðZÞ = pðz1 Þpðz2 jz1 Þ . . . pðzn jzn − 1 , zn − 2 , . . . , z2 , z1 Þ ð30:10Þ

|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
pðzn jZΦn Þ

where n ≪ N.
Tjelmeland and Eidsvik (2005) used a sampling algorithm that incorporates an
auxiliary random variable. These methods suffer from extensive CPU demand and
instability in convergence. Besides, the large structures cannot be reproduced finely,
a series of factors that make them difficult to use for 3D applications. Some of the
results of this method are shown in Fig. 30.6.
iv. Single Normal Equation Simulation (SNESIM)

The single normal equation simulation (SNESIM) is an improved version of the

original algorithm proposed by Guardiano and Srivastava (1993). The SNESIM
algorithm scans the input TI for once and then stores the frequency/probability of
all pattern occurrences in a search tree (Boucher 2009; Strebelle 2002), which
reduces the computational time significantly. Then, the probabilities are retrieved
from the constructed search-tree based on the existing data in the data-event.
The SNESIM algorithm is a pixel-based algorithm, which can perfectly reproduce
the conditioning point data.
The SNESIM algorithm is a sequential algorithm and, thus, each cell S can take
k possible states fsk , k = 1, . . . , K g, which usually represents facies unit. This
algorithm, like any other conditional techniques, calculates the joint probability
over n discrete points using:
30 Multiple Point Statistics: A Review 627

Fig. 30.6 A demonstration of the results of MRF (Daly 2005; Stien and Kolbjørnsen 2011). a,
c TI and b, d realizations

n

Φðh1 , . . . , hn ; k1 , . . . , kn Þ = E ∏ I ð u + hα ; k α Þ ð30:11Þ
α=1

where h, k, u and E represent separation vector (lag), state value, visiting location
and expected value, respectively. I ðu; k Þ also denotes the indicator value at location
u. This equation, thus, gives the probability of having n values ðk1 , . . . , kn Þ at the
locations sðu + h1 Þ, . . . , sðu + hn Þ. The above probability is replaced with the
following equation in SNESIM:

cðdn Þ
Φ ð h1 , . . . , hn ; k 1 , . . . , k n Þ ≅ ð30:12Þ
Nn

where Nn and cðdn Þ denote the total number of patterns in the TD and number of
replicates for the data event dn = fsðu + hn Þ = skα , α = 1, . . . , ng.
628 P. Tahmasebi

Fig. 30.7 Demonstration of

multiple-grid approach in
SNESIM. The figure is taken
from Wu et al. (2008)

This algorithm benefits from multiple-grid by which the large structures are first
captured using a smaller number of nodes and then the details are added. This
concept is illustrated in Fig. 30.7.
One of the limitations in the SNESIM algorithm is lack of producing realistic,
highly connected and large-scale geological features. This algorithm, however, can
be used only on categorical TIs. The SNESIM algorithm is still inefficient for the
real multimillion cells applications (Tahmasebi et al. 2014). Several other methods
were latter proposed to improve the efficiency and quality of the SNESIM algorithm
(Cordua et al. 2015; Straubhaar et al. 2013). A new technique has recently been
presented that can take the realizations and perform a structural adjustment to match
the well data (Tahmasebi 2017) (Fig. 30.8).
V. Direct Sampling

Direct sampling method is very similar to SIMPAT algorithm (see below) in that
sense it only scans a part of TI and pastes one single pixel (Mariethoz et al. 2010).
Since the TI is scanned in each loop of the simulation, thus, there is no need to
make any database and less RAM is required. Like the pattern-based techniques,
this algorithm uses a distance function for finding the closest patterns in TI. This
method can be used for both categorical and continuous variables.
The DS algorithm selects the known data at each visiting point. Then, the
similarity of the data-event with the TI is calculated based on a predefined searching
portion. As soon as the first occurrence of a matching data event in the TI is found
(corresponding to a distance under a given threshold acceptance), the value of the
central node of the data event in the TI is accepted and pasted in the simulation. It
30 Multiple Point Statistics: A Review 629

Fig. 30.8 The results of SNESIM. The realizations shown in (a, b) are generated using the TI in
Fig. 30.6c. c TI and d a realization based on the TI in (c)

should be noted that the searching phase stops if the algorithm finds a pattern that is
similar up to a given threshold. If not, the most similar pattern found in the pre-
defined portion of TI is selected and its central node is pasted on the simulation grid
(Fig. 30.9).
vi. Cumulants

More information beyond two-point statistics can be inferred using the cumulants
approach. This method, indeed, can extract such a higher information directly from
the existing data, rather than the TI. Dimitrakopoulos et al. (2010) first used this
method to simulate geological structures. The geological process, anisotropy and
pattern redundancy are the important factors that should be considered in selecting
the necessary cumulants (Mustapha and Dimitrakopoulos 2010). The conditional
probability is first calculated based on the available data. Then, the TI is only
researched if not sufficient replicates cannot be found in the data. One requires
selecting appropriate spatial cumulants for each geological scenario and there is no
specific strategy on this. Some of the results of this method are shown in Fig. 30.10.
630 P. Tahmasebi

TI Realization

TI Realization

Fig. 30.9 The results of the DS algorithm for modeling of a hydraulic conductivity field (upper
row) and a continuous property (b). The results are taken from Mariethoz et al. (2010) and Rezaee
et al. (2013)

30.5.2 Pattern-Based Algorithms

Pixel-based algorithms can have problems to preserve the continuity of the geo-
logical structures. To palliate this, some pattern-based methods have been devel-
oped which briefly are introduced bellow. Their commonality is that they do not
simulate one pixel at a time, but they paste an entire “patch” in the simulation. One
of the main aims of using pattern based simulation methods is their ability to
preserve the continuity and overall structure observed in TI.
i. Simulation of Pattern (SIMPAT)

The algorithm of simulation of patterns was first introduced to address some of the
limitations in the SNESIM algorithm, namely the CPU time and connectivity of
patterns (Arpat and Caers 2007). This method replaces the probability with a dis-
tance for finding most similar pattern. The algorithm can be summarized as follows.
30 Multiple Point Statistics: A Review 631

TI Realization

TI Realization

Fig. 30.10 The results of cumulants for modeling of two complex channelized systems. The
results are taken from Mustapha and Dimitrakopoulos (2010, 2011)

The TI is first scanned using a predefined template T and all the extracted patterns
are stored in a pattern database. Then simulation points are visited based on the
given random path and the corresponding data-event is extracted accordingly. One
of the patterns in pattern database is selected randomly if the data-event at the
visiting point contains no data. Otherwise, the most similar pattern is selected based
on the similarity between the data-event the patterns in pattern database. The above
steps are repeated for all visiting points. The results of SIMPAT algorithm are
realistic. However, it requires an extensive CPU time and encounters various
serious issues in 3D modeling. Furthermore, the produced results manifest a con-
siderable similarity with TI as this algorithm seeking for the best matching pattern.
Thus, this method seems to underestimate the spatial uncertainty. Some of the
results of SIMPAT are shown in Fig. 30.11.
In a similar fashion, the pattern-based techniques can be used within a Bayesian
framework (Abdollahifard and Faez 2013). This process, however, can be very
632 P. Tahmasebi

TI Realization

TI Realization

Fig. 30.11 The results of SIMPAT. These results are taken from Arpat (2005)

CPU demanding. Other enhancements on SIMPAT was also considered later by

incorporating wavelet decomposition (Gardet et al. 2016).
ii. Filter-based Simulation (FILTERSIM)

As pointed out, SIMPAT suffers from its computational cost, as it requires calcu-
lating the distance between the data-event the entire patterns in pattern database.
One possible solution is to summarize both the TI and data-event. Zhang et al.
(2006) proposed a new method, FLITERSIM, in which various filters (6 and 9
filters in 2D and 3D) have been used in order to reduce the spatial complexity and
dimensions. This allows reducing the complexity and computation time. Thus, the
patterns are first filtered using the pre-defined linear filters. Then, the outputs are
clustered based the similarity of the filtered patterns. Next, a prototype pattern is
computed for each cluster that represents the average of all the patterns in the
cluster. Afterwards, similar to SIMPAT, the most similar prototype is identified
using a distance function and one of its patterns is selected randomly. These steps
are continued until the simulation grid is filled. Due to using a limited number of
filters, this algorithm requires less computational time compared to SIMPAT. The
distance function in FILTERSIM was later replace with wavelet (Gloaguen and
Dimitrakopoulos 2009). The drawbacks of the wavelet-based method are that it has
a lot of parameters (e.g. wavelet decomposition level) that can effect on both quality
30 Multiple Point Statistics: A Review 633

and CPU time. Such parameters require an extensive tuning in order to achieve
good results in a reasonable time.
In a similar fashion, Eskandari and Srinivasan (2010) proposed Growthism to
integrate the dynamic data in the simulation. This method begins with the locations
of data and grows gradually and completes the simulation grid.
The most important shortcoming of FILTERSIM is that, it uses a limited set of
linear filters that cannot always convey all the information and variability in the TI.
Moreover, selecting the appropriate filters and several user-dependent parameters
for each TI is an issue that is that common among many MPS methods. Some of the
generated realizations using the FILTERSIM algorithm are shown in Fig. 30.12.
iii. Cross-Correlation based Simulation (CCSIM)

One of the recent algorithms of MPS is the cross correlation-based simulation

(CCSIM) algorithm that utilizes a cross-correlation function (CCF) along a
1D-raster path (Tahmasebi et al. 2012a). The CCF, which represents a multi-point

Realization #1 Realization #2

TI Realization

Fig. 30.12 The results of FILTERSIM. These results are taken from Zhang et al. (2006)
634 P. Tahmasebi

characteristic function, is used to match the patterns in a realization with those in

the TI. This algorithm has been adopted for different scenarios and computational
grids. For example, multi-scale CCSIM (Tahmasebi et al. 2014) can be used when
the simulation grid is very larger. This algorithm, similar to the SNESIM algorithm,
is also based on calculating the joint probability. The SNESIM algorithm calculates
the probability using a search-tree algorithm. However, calculating the above
conditional probability for every single point in the simulation grid, in the presence
of a large 2D/3D TD, is computationally prohibitive. Unless a very small neigh-
borhood is used, which leads to poorly connected features in the outcome
realizations.
In the CCSIM algorithm, the above limitation is addressed differently. First, the
probability function is replaced with a similarity function, called cross-correlation
function (CCF), which is much more efficient than drawing the probability. Sec-
ondly, based on the Markov Chain theory, the CCSIM algorithm uses a similar
search template (i.e. radius). However, unlike the previous algorithms where all the
data in the search template are used for simulating the visiting point, only a small
data-event located in the boundaries, called overlap region OL, is considered in the
calculations. Furthermore, except the point fallen in the OL region, the rest of the
points are removed from the visiting points and they would not be simulated again.
Thus, instead of simulating each single point, this algorithm ignores some of them
and partitions the grid into several blocks. The CCF can be calculated as follow:

Dx − 1 Dy − 1
CTD, DT ði, jÞ = ∑ ∑ TIðx + i, y + jÞDT ðx, yÞ, ð30:13Þ
x=0 y=0

with i ∈ ½0Tx + Dx − 1Þ and j ∈ ½0 Ty + Dy − 1Þ and i, j ∈ Z. The i and j represent the

shift steps in the x and y directions. TIðx, yÞ represents the location at point

ðx, yÞ of
TD of size Lx × Ly , with x ∈ f0, . . . , Dx − 1g and y ∈ 0, . . . , Dy − 1 . An OL
region of size Dx × Dy and a data event DT are used to match the pattern in the TI.
T represents the size of template used in CCSIM.
The CCSIM algorithm can realistically reproduce the large-scale structures in
diminutive time. These techniques, however, do not fully match the well data and
some artefacts are generated around the point conditioning data. Recently, this
techniques has been used within an iterative framework along with boundary cut-
ting methods by which the efficiency and conditioning data reproduction have been
increased significantly (Gardet et al. 2016; Kalantari and Abdollahifard 2016;
Mahmud et al. 2014; Moura et al. 2017; Scheidt et al. 2015; Tahmasebi and Sahimi
2016a, b; Yang et al. 2016). Some of the results of CCSIM are shown in Fig. 30.13.
Furthermore, this method has been successfully implemented for fine-scale mod-
eling in digital rock physics (Karimpouli et al. 2017; Karimpouli and Tahmasebi
2015; Tahmasebi et al. 2016a, b, c, 2017a, b).
30 Multiple Point Statistics: A Review 635

Realization #1 Realization #2

TI Realization

Fig. 30.13 The results of the CCSIM algorithm

30.5.3 Hybrid Algorithms

Each of the current MPS has some specific limitations. For example, the
pixel-based techniques are good in conditioning the point data, while they barely
can produce long-range connectivities. Similarly, the pattern-based methods can
produce such structures, but they are unable to preserve the hard data. Thus, the
idea of hybrid MPS method can be interesting if one uses the strength of both group
effectively. Following, the available hybrid methods are reviewed and their
advantages and disadvantages are discussed as well.
i. Hybrid Sequential Gaussian/Indicator Simulation and TI

Ortiz and Deutsch (2004), under an assumption of independence of the different

data sources, integrate the indicator method with MPS. Hence, instead of using a TI,
the MPS properties are obtained directly from the available hard data (variogram)
and integrated with the results of indicator kriging. Finally, a value is drawn from
this new distribution. These methods were further investigated by Ortiz and Emery
(2005). However, in most cases, the initial results of indicator kriging highly
influence final realization.
636 P. Tahmasebi

ii. Hybrid Pixel- and Pattern-based Simulation (HYPPS)

The strength of both the pixel- and pattern-based algorithms can be combined and
make a hybrid algorithm. Tahmasebi (2017) has combined these two algorithm and
proposed a new hybrid algorithm, called HYPPS. This algorithm discretizes the
simulation grid into regions with/without the conditioning data. One needs to
consider more attention the location containing the well data as providing them
require considerable cost. Thus, the SNESIM algorithm, as a pixel-based method
can be used around such locations. Regardless of the type of any geostatistical
methods, reproducing of patterns and well data are the most important factors.
Producing realistic models, without taking into account the conditioning data, or
vise versa, is not deasirable. Any successful algorithm should be able to manintin
both of the above crietra at the same time.
In the HYPPS algorithm the simulation grid is divided into two regions. Then,
the geostatistical methods are applied on each of them. Following this step, the
HYPPS algorithm uses the CCSIM, as a pattern-based algorithm, for the location
where no HD exist and, similarly, the SNESIM algorithm is used around the well
data, which can precisely reproduce the conditioning data. Thus, the hybrid state of
the pixel-based and pattern-based techniques can be written as follow:

n
 n
Φ ð h1 , . . . , h n ; k 1 , . . . , k n Þ = E ∏ I ðu + hα ; kα Þ + ∏ Φðhα jhΦn Þ ð30:14Þ
α=1 α=1

which implies that the joint event over a template where both methods are working
simultaneously can be expressed as the summation of the two probability distri-
butions defined earlier (see the SNESIM algorithm). Thus, a normalization terms,
namely nx and np , should be included such that nx + np = 1. Note that nx and np
represent normalized number (or percentage) of the simulated points used in the
pixel- and pattern-based methods. An equivalent form of the above probability can
be expressed as:
0   1
cðdn Þ
nx +
B Nn
!C
Φðh1 , . . . , hn ; k1 , . . . , kn Þ ≅ B
@ Dx − 1 Dy − 1 C ̸ nx + np
A
np ∑ ∑ TIðx + i, y + jÞDT ðx, yÞ
x=0 y=0

ð30:15Þ

The second term in the above equation on the right side is used for the areas
where the CCSIM algorithm is utilized. While, the visiting points around the well
data are evaluated jointly.
It is worth mentioning that co-template (Tahmasebi et al. 2014) can be used with
the CCSIM to give the priority to the patterns that contain the conditioning data
ahead of the raster path. Therefore, long-range connectivity structures are taken into
30 Multiple Point Statistics: A Review 637

Fig. 30.14 The application of the SNESIM algorithm for simulating a grid when the borders are
conditioned (a): the TI, b: boundary data, c generated realizations using the boundary data, d:
boundary data along with well data, e generated realizations using the boundary and well data.
Note, the sizes of the TI and simulation grid in (a) and (b) are 250 × 250 and 100 × 100,
respectively. The shale and sand facies are indicated by blue and red colors, individually

account even in the blocks with no conditioning data. An example is provided in

Fig. 30.13e. Although the density of conditioning data, compared to the previous
scenario, is increased, but the produced realizations represent the real heterogeneity
represented in the TI. The HYPPS algorithm can be used to integrate data at
different scales as well (Fig. 30.14).
638 P. Tahmasebi

30.6 Current Challenges

The MPS techniques have been developed extensively to deal with complex and
more realistic problems in which the geology and other source of data such as well
and secondary data are reproduced. There are still some critical challenges in MPS
that require more research. Some of them are listed below:
• Conditioning the model to dense point/secondary datasets (e.g. well and
seismic data): This issue is one of the main challenges in the pattern-based
techniques. In fact, such methods require calculating distance function between
the patterns of heterogeneities in the TI and the model and between the con-
ditioning data that must be honored exactly in conditional simulation (Hashemi
et al. 2014a, b; Rezaee and Marcotte 2016; Tahmasebi and Sahimi 2015b).
• Lack of similarity between the internal patterns in the TI: Producing a
comprehensive and diverse TI is very challenging. Thus, serious discontinuity
and unrealistic structures are generated using the current MPS techniques when
the TI is not enough large. On the other hand, using very large or multiple TIs
are costly for most of such methods.
• Deficiency of the current similarity functions for quantitative modeling:
Most of the current distance functions are based on some simple and two-point
criteria by which the optimal pattern cannot be identified. Such distance func-
tions are very limited in conveying the information. Thus, more informative
similarity functions are required.
• Better and more realistic validation methods: There is not so many number of
methods that can be used to evaluate the performance of new developed MPS
algorithms (Dimitrakopoulos et al. 2010; Tan et al. 2014). For example, the
realizations can show a considerable different with each other and TI, while such
a variability cannot be quantified using the current methods. Thus, visual
comparison is still one of the popular method for verifying the performance of
the MPS methods.
Many issues must be addressed yet. For example, the current MPS methods are
designed for stationary TIs, whereas the properties of many large-scale porous
media exhibit non-stationary features. Some progress has recently been made in this
direction (Chugunova and Hu 2008; Honarkhah and Caers 2012; Tahmasebi and
Sahimi 2015a). In addition, associated with every TI is large uncertainties. Thus, if
several TIs are available, it is necessary to design methods that can determine which
TI(s) to use in a given context.
30 Multiple Point Statistics: A Review 639

References

Abdollahifard MJ, Faez K (2013) Stochastic simulation of patterns using Bayesian pattern
modeling. Comput Geosci 17:99–116. https://doi.org/10.1007/s10596-012-9319-x
Anderson KS, Hickson TA, Crider JG, Graham SA (1999) Integrating teaching with field research in
the wagon rock project. J Geosci Educ 47:227–235. https://doi.org/10.5408/1089-9995-47.3.227
Arpat B (2005) Sequential simulation with patterns. Stanford University
Arpat B, Caers J (2007) Stochastic simulation with patterns. Math Geol 39:177–203
Bayer P, Huggenberger P, Renard P, Comunian A (2011) Three-dimensional high resolution
fluvio-glacial aquifer analog: Part 1: field study. J Hydrol 405:1–9. https://doi.org/10.1016/j.
jhydrol.2011.03.038
Biswal B, Manwart C, Hilfer R, Bakke S, Øren PE (1999) Quantitative analysis of experimental
and synthetic microstructures for sedimentary rock. Phys Stat Mech Appl 273:452–475. https://
doi.org/10.1016/S0378-4371(99)00248-4
Biswal B, Øren P-E, Held RJ, Bakke S, Hilfer R (2007) Stochastic multiscale model for carbonate
rocks. Phys Rev E 75:61303. https://doi.org/10.1103/PhysRevE.75.061303
Borgman L, Taheri M, Hagan R (1984) Three-dimensional, frequency-domain simulations of
geological variables. In: Geostatistics for natural resources characterization. Springer,
Netherlands, Dordrecht, pp 517–541. https://doi.org/10.1007/978-94-009-3699-7_30
Boucher A (2009) Considering complex training images with search tree partitioning. Comput
Geosci 35:1151–1158
Bryant S, Blunt M (1992) Prediction of relative permeability in simple porous media. Phys Rev A
46:2004–2011. https://doi.org/10.1103/PhysRevA.46.2004
Chiles J-P, Delfiner P (2011) Geostatistics : modeling spatial uncertainty. Wiley-Blackwell
Chu J, Journel AG (1994) Conditional fBm simulation with dual kriging. Springer Netherlands,
pp. 407–421. https://doi.org/10.1007/978-94-011-0824-9_44
Chugunova TL, Hu LY (2008) Multiple-point simulations constrained by continuous auxiliary
data. Math Geosci 40:133–146. https://doi.org/10.1007/s11004-007-9142-4
Cordua KS, Hansen TM, Mosegaard K (2015) Improving the pattern reproducibility of
multiple-point-based prior models using frequency matching. Math Geosci 47:317–343.
https://doi.org/10.1007/s11004-014-9531-4
Cressie NAC, Wikle CK (2011) Statistics for spatio-temporal data. Wiley
Daly C (2005) Higher order models using entropy, markov random fields and sequential
simulation. Springer, Netherlands, pp 215–224. https://doi.org/10.1007/978-1-4020-3610-1_22
Davis MW (1987) Production of conditional simulations via the LU triangular decomposition of
the covariance matrix. Math Geol 19:91–98. https://doi.org/10.1007/bf00898189
Deutsch CV (1992) Annealing techniques applied to reservoir modeling and the integration of
geological and engineering (well test) data/. Stanford University
Deutsch CV, Journel AG (1998) GSLIB: geostatistical software library and user’s guide, 2nd edn.
Oxford University Press, New York
Deutsch CV, Wang L (1996) Hierarchical object-based stochastic modeling of fluvial reservoirs.
Math Geol 28:857–880. https://doi.org/10.1007/BF02066005
Deutsch CV, Wen XH (2000) Integrating large-scale soft data by simulated annealing and
probability constraints. Math Geol 32:49–67. https://doi.org/10.1023/A:1007502817679
Dimitrakopoulos R, Luo X (2004) Generalized sequential gaussian simulation on group size and
screen-effect approximations for large field simulations. Math Geol 36:567–591. https://doi.
org/10.1023/B:MATG.0000037737.11615.df
Dimitrakopoulos R, Mustapha H, Gloaguen E (2010) High-order statistics of spatial random fields:
exploring spatial cumulants for modeling complex non-gaussian and non-linear phenomena.
Math Geosci 42:65–99. https://doi.org/10.1007/s11004-009-9258-9
El Ouassini A, Saucier A, Marcotte D, Favis BD (2008) A patchwork approach to stochastic
simulation: a route towards the analysis of morphology in multiphase systems. Chaos Solitons
Fractals 36:418–436. https://doi.org/10.1016/j.chaos.2006.06.100
640 P. Tahmasebi

Eskandari K, Srinivasan S (2010) Reservoir modelling of complex geological systems–a

multiple-point perspective. J Can Pet Technol 49:59–69. https://doi.org/10.2118/139917-PA
Fang JH, Wang PP (1997) Random field generation using simulated annealing vs. fractal-based
stochastic interpolation. Math Geol 29:849–858. https://doi.org/10.1007/BF02768905
Gardet C, Le Ravalec M, Gloaguen E (2016) Pattern-based conditional simulation with a raster
path: a few techniques to make it more efficient. Stoch Environ Res Risk Assess 30:429–446.
https://doi.org/10.1007/s00477-015-1207-1
Gloaguen E, Dimitrakopoulos R (2009) Two-dimensional conditional simulations based
on the wavelet decomposition of training images. Math Geosci 41:679–701. https://doi.org/
10.1007/s11004-009-9235-3
Goovaerts P (1997) Geostatistics for natural resources evaluation. Oxford University Press
Gross LJ, Small MJ (1998) River and floodplain process simulation for subsurface characteri-
zation. Water Resour Res 34:2365–2376. https://doi.org/10.1029/98WR00777
Guardiano FB, Srivastava RM (1993) Multivariate geostatistics: beyond bivariate moments. In:
Geostatistics Troia’92. Springer, pp 133–144
Haldorsen HH, Damsleth E (1990) Stochastic Modeling (includes associated papers 21255 and
21299). J Pet Technol 42:404–412. https://doi.org/10.2118/20321-PA
Hamzehpour H, Sahimi M (2006) Development of optimal models of porous media by combining
static and dynamic data: the porosity distribution. Phys Rev E 74:26308. https://doi.org/10.
1103/PhysRevE.74.026308
Hashemi S, Javaherian A, Ataee-pour M, Tahmasebi P, Khoshdel H (2014a) Channel
characterization using multiple-point geostatistics, neural network, and modern analogy: a
case study from a carbonate reservoir, southwest Iran. J Appl Geophys 111. https://doi.org/10.
1016/j.jappgeo.2014.09.015
Hashemi S, Javaherian A, Ataee-pour M, Tahmasebi P, Khoshdel H (2014b) Channel
characterization using multiple-point geostatistics, neural network, and modern analogy: a
case study from a carbonate reservoir, southwest Iran. J Appl Geophys 111:47–58. https://doi.
org/10.1016/j.jappgeo.2014.09.015
Holden L, Hauge R, Skare Ø, Skorstad A (1998) Modeling of fluvial reservoirs with object
models. Math Geol 30:473–496. https://doi.org/10.1023/A:1021769526425
Honarkhah M, Caers J (2012) Direct pattern-based simulation of non-stationary geostatistical
models. Math Geosci 44:651–672. https://doi.org/10.1007/s11004-012-9413-6
Journel A, Zhang T (2006) The necessity of a multiple-point prior model. Math Geol 38:591–610
Journel AG, Huijbregts CJ (1978) Mining geostatistics. Academic Press
Kalantari S, Abdollahifard MJ (2016) Optimization-based multiple-point geostatistics: a sparse
way. Comput Geosci 95:85–98. https://doi.org/10.1016/j.cageo.2016.07.006
Karimpouli S, Tahmasebi P (2015) Conditional reconstruction: an alternative strategy in digital
rock physics. Geophysics 81. https://doi.org/10.1190/geo2015-0260.1
Karimpouli S, Tahmasebi P, Ramandi HL, Mostaghimi P, Saadatfar M (2017) Stochastic modeling
of coal fracture network by direct use of micro-computed tomography images. Int J Coal Geol
179. https://doi.org/10.1016/j.coal.2017.06.002
Kitanidis PK, Peter K (1997) Introduction to geostatistics: applications to hydrogeology.
Cambridge University Press
Kleingeld WJ, Thurston ML, Prins CF, Lantuéjoul C (1997) The conditional simulation of a Cox
process with application to deposits with discrete particles. Geostat Wollongong 96:683–694
Lancaster ST, Bras RL (2002) A simple model of river meandering and its comparison to natural
channels. Hydrol Process 16:1–26. https://doi.org/10.1002/hyp.273
Lantuéjoul C (2002) Geostatistical simulation: models and algorithms, in: geostatistical
simulation: models and algorithms. Springer, Berlin, Heidelberg, pp 1–6. https://doi.org/10.
1007/978-3-662-04808-5_1
Mahmud K, Mariethoz G, Caers J, Tahmasebi P, Baker A (2014) Simulation of earth textures by
conditional image quilting. Water Resour Res 50:3088–3107. https://doi.org/10.1002/
2013WR015069
30 Multiple Point Statistics: A Review 641

Manwart C, Torquato S, Hilfer R (2000) Stochastic reconstruction of sandstones. Phys Rev E

62:893–899. https://doi.org/10.1103/PhysRevE.62.893
Mariethoz G, Renard P, Straubhaar J (2010) The direct sampling method to perform multiple‐point
geostatistical simulations. Water Resour Res 46
Matheron G (1973) The intrinsic random functions and their applications. Adv Appl Probab
5:439–468. https://doi.org/10.1017/S0001867800039379
Moura P, Laber E, Lopes H, Mesejo D, Pavanelli L, Jardim J, Thiesen F, Pujol G (2017) LSHSIM:
a locality sensitive hashing based method for multiple-point geostatistics. Comput Geosci
https://doi.org/10.1016/j.cageo.2017.06.013
Mustapha H, Dimitrakopoulos R (2010) High-order stochastic simulation of complex spatially
distributed natural phenomena. Math Geosci 42:457–485
Mustapha H, Dimitrakopoulos R (2011) HOSIM: a high-order stochastic simulation algorithm for
generating three-dimensional complex geological patterns. Comput Geosci 37:1242–1253.
https://doi.org/10.1016/j.cageo.2010.09.007
Ortiz JM, Deutsch CV (2004) Indicator simulation accounting for multiple-point statistics. Math
Geol 36:545–565. https://doi.org/10.1023/B:MATG.0000037736.00489.b5
Ortiz JM, Emery X (2005) Integrating multiple-point statistics into sequential simulation algorithms.
Springer, Netherlands, pp 969–978. https://doi.org/10.1007/978-1-4020-3610-1_101
Parra Á, Ortiz JM (2011) Adapting a texture synthesis algorithm for conditional multiple point
geostatistical simulation. Stoch Environ Res Risk Assess 25:1101–1111. https://doi.org/10.
1007/s00477-011-0489-1
Peredo O, Ortiz JM (2011) Parallel implementation of simulated annealing to reproduce
multiple-point statistics. Comput Geosci 37:1110–1121. https://doi.org/10.1016/j.cageo.2010.
10.015
Pourfard M, Abdollahifard MJ, Faez K, Motamedi SA, Hosseinian T (2017) PCTO-SIM:
Multiple-point geostatistical modeling using parallel conditional texture optimization. Comput
Geosci 102:116–138. https://doi.org/10.1016/j.cageo.2016.12.012
Pyrcz MJ, Boisvert JB, Deutsch CV (2009) ALLUVSIM: a program for event-based stochastic
modeling of fluvial depositional systems. Comput Geosci 35:1671–1685. https://doi.org/10.
1016/j.cageo.2008.09.012
Rasera LG, Machado PL, Costa JFCL (2015) A conflict-free, path-level parallelization approach
for sequential simulation algorithms. Comput Geosci 80:49–61. https://doi.org/10.1016/j.
cageo.2015.03.016
Rezaee H, Marcotte D (2016) Integration of multiple soft data sets in MPS thru multinomial
logistic regression: a case study of gas hydrates. Stoch Environ Res Risk Assess 1–19. https://
doi.org/10.1007/s00477-016-1277-8
Rezaee H, Mariethoz G, Koneshloo M, Asghari O (2013) Multiple-point geostatistical simulation
using the bunch-pasting direct sampling method. Comput Geosci 54:293–308. https://doi.org/
10.1016/j.cageo.2013.01.020
Scheidt C, Tahmasebi P, Pontiggia M, Da Pra A, Caers J (2015) Updating joint uncertainty in
trend and depositional scenario for reservoir exploration and early appraisal. Comput Geosci
19. https://doi.org/10.1007/s10596-015-9491-x
Seminara G (2006) Meanders. J Fluid Mech 554:271. https://doi.org/10.1017/S002211200
6008925
Sheehan N, Torquato S (2001) Generating microstructures with specified correlation functions.
J Appl Phys 89:53–60. https://doi.org/10.1063/1.1327609
Skorstad A, Hauge R, Holden L (1999) Well conditioning in a fluvial reservoir model. Math Geol
31:857–872. https://doi.org/10.1023/A:1007576801266
Stien M, Kolbjørnsen O (2011) Facies modeling using a markov mesh model specification. Math
Geosci 43:611–624. https://doi.org/10.1007/s11004-011-9350-9
Straubhaar J, Walgenwitz A, Renard P (2013) Parallel multiple-point statistics algorithm based on
list and tree structures. Math Geosci 45:131–147. https://doi.org/10.1007/s11004-012-9437-y
Strebelle S (2012) Multiple-point geostatistics: from theory to practice. Ninth international
geostatistics congress. Springer, Oslo, Norway, pp 11–15
642 P. Tahmasebi

Strebelle S (2002) Conditional simulation of complex geological structures using multiple-point

statistics. Math Geol 34:1–21
Tahmasebi P (2017) Structural adjustment for accurate conditioning in large-scale subsurface
systems. Adv Water Resour 101. https://doi.org/10.1016/j.advwatres.2017.01.009
Tahmasebi P, Hezarkhani A, Sahimi M (2012a) Multiple-point geostatistical modeling based on
the cross-correlation functions. Comput Geosci 16:779–797. https://doi.org/10.1007/s10596-
012-9287-1
Tahmasebi P, Sahimi M, Mariethoz G, Hezarkhani A (2012b) Accelerating geostatistical
simulations using graphics processing units (GPU). Comput Geosci 46:51–59. https://doi.org/
10.1016/j.cageo.2012.03.028
Tahmasebi P, Javadpour F, Sahimi M (2016a) Stochastic shale permeability matching:
three-dimensional characterization and modeling. Int J Coal Geol 165:231–242. https://doi.
org/10.1016/j.coal.2016.08.024
Tahmasebi P, Javadpour F, Sahimi M, Piri M (2016b) Multiscale study for stochastic
characterization of shale samples. Adv Water Resour 89:91–103. https://doi.org/10.1016/j.
advwatres.2016.01.008
Tahmasebi P, Sahimi M (2015a) Reconstruction of nonstationary disordered materials and media:
watershed transform and cross-correlation function. Phys Rev E 91:32401. https://doi.org/10.
1103/PhysRevE.91.032401
Tahmasebi P, Sahimi M (2015b) Geostatistical simulation and reconstruction of porous media by a
cross-correlation function and integration of hard and soft data. Trans Porous Media 107:871–
905. https://doi.org/10.1007/s11242-015-0471-3
Tahmasebi P, Sahimi M (2016a) Enhancing multiple-point geostatistical modeling: 2. Iterative
simulation and multiple distance function. Water Resour Res 52:2099–2122. https://doi.org/10.
1002/2015WR017807
Tahmasebi P, Sahimi M (2016b) Enhancing multiple-point geostatistical modeling: 1. Graph
theory and pattern adjustment. Water Resour Res 52:2074–2098. https://doi.org/10.1002/
2015WR017806
Tahmasebi P, Sahimi M, Andrade J (2017a) Direct modeling of granular materials. In:
Poromechanics VI. American society of civil engineers, Reston, VA, pp. 1436–1442. https://
doi.org/10.1061/9780784480779.178
Tahmasebi P, Sahimi M, Andrade JE (2017b) Image-based modeling of granular porous media.
Geophys Res Lett https://doi.org/10.1002/2017gl073938
Tahmasebi P, Sahimi M, Caers J (2014) MS-CCSIM: accelerating pattern-based geostatistical
simulation of categorical variables using a multi-scale search in fourier space. Comput Geosci
67:75–88. https://doi.org/10.1016/j.cageo.2014.03.009
Tahmasebi P, Sahimi M, Kohanpur AH, Valocchi A (2016) Pore-scale simulation of flow of CO2
and brine in reconstructed and actual 3D rock cores. J Pet Sci Eng. https://doi.org/10.1016/j.
petrol.2016.12.031
Tan X, Tahmasebi P, Caers J (2014) Comparing training-image based algorithms using an analysis
of distance. Math Geosci 46. https://doi.org/10.1007/s11004-013-9482-1
Tjelmeland H, Eidsvik J (2005) Directional metropolis : hastings updates for posteriors with
nonlinear likelihoods. Springer, Netherlands, pp. 95–104. https://doi.org/10.1007/978-1-4020-
3610-1_10
Toftaker H, Tjelmeland H (2013) Construction of binary multi-grid markov random field prior
models from training images. Math Geosci 45:383–409. https://doi.org/10.1007/s11004-013-
9456-3
Wu J, Boucher A, Zhang T (2008) A SGeMS code for pattern simulation of continuous and
categorical variables: FILTERSIM. Comput Geosci 34:1863–1876. https://doi.org/10.1016/j.
cageo.2007.08.008
Yang L, Hou W, Cui C, Cui J (2016) GOSIM: a multi-scale iterative multiple-point statistics
algorithm with global optimization. Comput Geosci 89:57–70. https://doi.org/10.1016/j.cageo.
2015.12.020
30 Multiple Point Statistics: A Review 643

Yeong CLY, Torquato S (1998a) Reconstructing random media. II. Three-dimensional media from
two-dimensional cuts. Phys Rev E 58:224–233. https://doi.org/10.1103/PhysRevE.58.224
Yeong CLY, Torquato S (1998b) Reconstructing random media. Phys Rev E 57:495–506. https://
doi.org/10.1103/PhysRevE.57.495
Zhang T, Switzer P, Journel A (2006) Filter-based classification of training image patterns for
spatial simulation. Math Geol 38:63–80

Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.