Multiple Point Statistics: A Review

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.


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 problems. 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 (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 heterogeneity and non-smooth phenomena. One consequence is the underestimation and overestimation for low and high values, respectively. This problem becomes evident 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 simulations 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). 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 simulations. 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 generated models using these techniques cannot represent an accurate representation of any physics-based simulations (e.g. flow, grade distribution, contaminate forecasting and etc.).
In the MPS methods, the spatial statistics are not either extracted using variogram, 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 demonstrated. Next, different MPS methods are explained and the advantages and disadvantages associated with each method are demonstrated. Finally, some avenues for future research are discussed.

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 simulation is to produce realizations z u α ð Þ, α = 1, . . . , N f g , conditioned to n available data and reproducing a given multivariate distribution. For this aim, the multivariate distribution is decomposed into a set of N univariate conditional cumulative distribution functions (ccdfs): ð Þ conditioned to a set of n original data and N − 1 ð Þ previously simulated values.

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 previously 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 transforms or histogram anamorphous through Hermite polynomials.

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 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 I u ð Þ f g= p is the stationary proportion of a given category. The indicator variogram can be expressed as: ð30:2Þ Usually, the categorical variables expressed as a set of K discrete categories that z u ð Þ 0, . . . , k − 1 f g . Therefore, the indicator value for each of the defined classes can be expressed as: 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: We can rewrite the above equation for categorical variables by using simple Kriging as: where E I u, k ð Þ f gis 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 random 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).

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.

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 traditional 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 sciences 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 , g a 2D model generated using the process-based techniques , h a 3D model generated by the process-based methods  algorithm to provide any further alterations. The results of object-based simulation methods are one of the best and most accessible sources for TIs. • Process-based Methods: Process-based methods (Biswal et al. 2007(Biswal et al. , 1999Bryant 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 methods are, however, computationally expensive and require considerable calibrations. 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.

Simulation Path
Geostatistical techniques are conducted on a simulation grid G, which is constructed 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).

Random Path
Random path is one of the most commonly used visiting path in sequential simulation 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 unbiasedness around the conditioning point data.

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. (2014Tahmasebi et al. ( , 2012a implemented a raster path along a fast similarity computation and achieved high-quality realizations. Such paths usually produced high quality realizations that can barely be produced using the random path algorithms. 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 algorithms, 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.

Some Other Definitions
Simulation Grid (G): a 2D/3D computation grid on which the geostatistical modeling 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.

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 through the hybrid frameworks will result in an efficient algorithm. Such a combination 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).

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: where O j and S j 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 where E i 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, 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 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) 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 function. 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 statistical methods have also used the effect of iterations on removing the artifacts (Pourfard et al. 2017;. It should be noted that the sequential algorithms can be parallel using different strategies which are not discussed 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 algorithm. 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: 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 s k , k = 1, . . . , K f g , which usually represents facies unit. This algorithm, like any other conditional techniques, calculates the joint probability over n discrete points using: 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 ðk 1 , . . . , k n Þ at the locations s u + h 1 ð Þ, . . . , s u + h n ð Þ. The above probability is replaced with the following equation in SNESIM: where N n and c d n ð Þ denote the total number of patterns in the TD and number of replicates for the data event d n = s u + h n ð Þ= s k α , α = 1, . . . , n f g . 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 . 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  Wu et al. (2008) 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 predefined 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 . 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.   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)

Pattern-Based Algorithms
Pixel-based algorithms can have problems to preserve the continuity of the geological structures. To palliate this, some pattern-based methods have been developed 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 distance for finding most similar pattern. The algorithm can be summarized as follows. 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) 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 considerable 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 TI Realization Fig. 30.10 The results of cumulants for modeling of two complex channelized systems. The results are taken from Dimitrakopoulos (2010, 2011) 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 calculating 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 TI Realization TI Realization Fig. 30.11 The results of SIMPAT. These results are taken from Arpat (2005) 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 Fig. 30.12 The results of FILTERSIM. These results are taken from Zhang et al. (2006) 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  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 neighborhood is used, which leads to poorly connected features in the outcome realizations.

TI Realization
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. Secondly, 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: with i ∈ 0 ½ T x + D x − 1Þ and j ∈ ½0 T y + D y − 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 L x × L y , with x ∈ 0, . . . , D x − 1 f gand y ∈ 0, . . . , D y − 1 È É . An OL region of size D x × D y and a data event D T 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 cutting 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 modeling in digital rock physics (Karimpouli et al. 2017;Karimpouli and Tahmasebi 2015;Tahmasebi et al. , b, c, 2017a.

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.

Realization #1
Realization #2 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.

TI Realization
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: 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 distributions defined earlier (see the SNESIM algorithm). Thus, a normalization terms, namely n x and n p , should be included such that n x + n p = 1. Note that n x and n p 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: 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  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 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). 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

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 conditioning 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 functions 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 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. 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.