1 Introduction

Various problems emerge in the process of reservoir modeling and simulation basically due to complex geological heterogeneity. Accurate modeling of heterogeneous nature of reservoir permeability, porosity, etc. will cause better prediction of hydrocarbon production from oil/gas fields (Ringrose and Bentley 2016; Kargozarfard et al. 2019). New insights into many geoscience problems, and especially those in the petroleum industry, are brought by accounting for the spatial relationships of the reservoir property data (Rezaee et al. 2013; Koneshloo et al. 2018; Yang et al. 2016). This is routinely done via kriging-based methods including the sequential Gaussian simulation. (Dubrule 1989; Yong et al. 2010; Yang et al. 2016). The major disadvantage of kriging-based methods is the excessive smoothness of the resultant property for the reservoir (Sahimi 2011). Another important aspect to be considered is the scarcity of hard data of the reservoir that postulates the application of multiple point statistics (MPS) alongside kriging-based methods (Caers 2001; Mariethoz 2018). “Bull eye” is the structure in which the property increases or decreases radially around the hard data point location. It is occasionally produced in property modeling around the hard data points and is categorized as an artifact (Ayzel et al. 2017). “Artifacts” are all unrealistic structures or patterns that may be produced because of the failure of image processing methods. The best method is the one which produces least possible artifacts within its calculations.

Recently, new methods called multiple point statistics (MPS) have been applied for this purpose as a cutting-edge category of methods for reservoir property modeling. MPS is a category of texture synthesis algorithms and is borrowed from image processing techniques (Mariethoz and Lefebvre 2014). In other words, MPS is actually the application of texture synthesis algorithms in earth sciences. We can also call them adapted versions of computer graphics methods for geoscience. By this usage, computer graphics are not only concerned with the visual appearances. This interesting usage of them is an open research area in reservoir engineering. In the MPS calculation, modeling the spatial patterns of natural properties, for example reservoir permeability, is based on training images (TIs), which are the explicit examples of the heterogeneity (Mariethoz and Lefebvre 2014). Generally speaking, training images are an ensemble of all spatial variations that must be regenerated by MPS methods (Mariethoz and Caers 2014). MPS methods statistically reproduce the values occurred in the TI with their spatial arrangement. In other words, the high-order statistics inferred from a TI, which is a spatial and conceptual representation of the structures expected in the geological formation, is used for the property modeling (Gardet et al. 2016). An MPS enables the generation of complex, interconnected or continuous geological structures within the stochastic framework (Le Coz et al. 2016).

1.1 MPS methods

The kriging-based methods calculate the reservoir properties by minimizing the variance of residuals while the MPS methods calculate the properties based on the notion of expectation (Mariethoz and Caers 2014). In an overall point of view, all MPS methods conduct their calculations based on expectation of the similar patterns and variation between TI and the reservoir. The essential idea of MPS is to find the most similar pattern in the TI and place it in the desired location. The most similar pattern(s) is/are found based on the calculation of similarity of the previously determined properties between TI and the simulation grid, the reservoir model in this case. The spatial arrangement of previously determined property creates a concept named a “data event.” In other words, by statistical calculation and finding the best match between a TI data event and the reservoir data event, the MPS method fills the reservoir by copying the best match data event into the reservoir at the desired location. In kriging-based methods, this filling is done only by measuring the distances of the desired location from known locations. The involvement of the data event enables the MPS method to create similar patterns to the TI. However, this concept and consequent directional spatial variation in property are ignored in kriging-based methods and only the distance from the known points determines the resulting values of desirable property. Thus, MPS considers the directions in spatial variation during the property modeling while kriging-based methods ignore it. Therefore, various types of property spatial variation could not be modeled by kriging-based methods.

Guardiano and Srivastava (1993) introduced multiple point geostatistical simulation by considering the borrowing concept from TI spatial variation (Guardiano and Srivastava 1993). Although this strategy was appealing, it was inefficient in computational cost aspects of calculations since the algorithm needs to scan the entire TI for each cell, the synonym of pixel in computer graphics, in the reservoir. Strebelle (2002) developed the first efficient MPS algorithm by using tree structures for searching through TI and introduced it as SNESIM (Strebelle 2002). This algorithm first stores all patterns of TI in the storage step. Then, the algorithm fills the reservoir based on the stored patterns and their frequencies. The most acute difficulty for this MPS method is in simulating the continuous property. Afterward, many other modifications were introduced by other researchers (Hoffman and Caers 2007; Caers et al. 2003; Huysmans and Dassargues 2011). Mariethoz et al. (2010) developed direct sampling (DS) method for eliminating the storage step of previous MPS methods (Mariethoz et al. 2010). A number of modifications have been applied for accelerating and augmenting this method (Rezaee et al. 2013; Huang et al. 2013). Among all modifications, advanced-DS changes the size of a data event during property calculation and produces more realistic patterns of property variation (Meerschman et al. 2013). Rezaee et al. (2013) modified the algorithm in order to assign a bunch of values instead of a single value for a single cell (Rezaee et al. 2013). This method called Bunch-DS achieved significant improvement in the cpu time of calculation. Tahmasebi et al. (2012a) engaged the correlation operator into MPS calculations and introduced a method named cross-correlation simulation, CCSIM (Tahmasebi et al. 2012a). CCSIM played a fundamental role for several enhancements of computation cost and time (Tahmasebi and Sahimi 2016b; Tahmasebi et al. 2012b, 2014). The available one is MSCCSIM which calculates the patterns of TI in a “multi-scale search” manner while using the convolution operator with applicable runtime (Tahmasebi et al. 2014). The image quilting (IQ) technique was first introduced by Efros and Freeman (Efros and Freeman 2001). Mahmud et al. (2014) implemented the IQ algorithm in 3D cases for property modeling applications (Mahmud et al. 2014). IQ results in appropriate modeled properties in quite a short cpu time. The Graph-cut method developed by Kwatra et al. supplies another improvement for regenerating more realistic outputs (Kwatra et al. 2003). The basic idea of Graph-cut is to consider the images as a graph and then solve a “max-flow/min-cut” problem for synthesizing the textures within the TI. The advantage of this method is to store the quality of cutting at each step which makes the method more powerful in producing seamless results. More details are addressed in the literature (Kwatra et al. 2003). Li et al. conveyed the message of the Graph-cut concept into geological settings (Li et al. 2016).

The values in the simulation grid, the reservoir in this case, are the selected values of TI. Selected values are chosen based on the consideration of neighboring values in both TI and the simulation grid. Actually, the most similar pattern is selected in the TI for imbedding in the simulation grid. The most similar pattern is found by using a template scanning over the TI. The template is a rigid set of values and locations that can move on the TI and take various spatial distributions of the property values based on the shape of the template. The similarity of the TI values and values in the simulation grid is defined by the distance concept, which is obtained by the difference of values in the TI and those of simulation grids by an explicit formula. For more details, see (Mariethoz and Caers 2014). After calculating the distance value, the MPS method chooses that part of the TI with the lowest distance and pastes it (or maybe a masked version of it) to the simulation grid. All mentioned MPS methods cause the texture of permeability to become better and better, but still they possess a tile-based nature for the property calculations. Tile-based nature is started with an empty reservoir and populating it piece by piece until all the cells of the reservoir attain the property. Figure 1 shows the basic concept of the MPS approach for calculating the reservoir property by finding the most similar pattern in the TI in each iteration of the calculation. The next generation of MPS methods is the ones borrowing the texture optimization technique from Kwatra et al. (2005). Pourfard et al. (2016) developed this generation of methods for MPS application with parallelism and ended in introducing PCTOSIM (Pourfard et al. 2016). In this new generation of MPS methods, the method starts with random values of property, in contrast to the empty reservoir at the beginning of previous generation of MPS. The algorithm makes these random values step by step similar to the existing patterns and values of TI. Figure 2 illustrates the texture optimization process beginning from random values and ending in the textures similar to TI.

Fig. 1
figure 1

General approach of MPS methods is depicted thorough the CCSIM method for finding the most similar pattern of the training image (TI). Section (e) of the figure is the training image, or TI. The first pattern is selected randomly (Step a), and consequent patterns are selected one after another based on previous one(s) and their overlap region with previously pasted pattern, step b–d (Tahmasebi et al. 2014)

Full size image
Fig. 2
figure 2

General workflow of a new generation of MPS methods via texture optimization. Having a TI (in computer graphic this equals to exemplar) and randomly distributed values, named noise in this figure, the method tries to optimize the random variable to be similar to the existing texture of the training image (TI) or exemplar. Started with noise and ended in level 3 (Kopf et al. 2007)

Full size image

1.2 MPS applications

MPS is used in different disciplines ranging from soil science (Meerschman et al. 2014), prediction of the occurrence of rainfall (Oriani et al. 2014), water resources modeling, i.e., remote sensing, (Mariethoz et al. 2012), hydrogeology (Huysmans and Dassargues 2012), medical imaging (Pham 2012) to the fluid flow through underground formations (Okabe and Blunt 2004; Huysmans and Dassargues 2009; Xu et al. 2012; Tamayo-Mas et al. 2016; Lee et al. 2019; Mosser et al. 2018).

MPS methods could also be implemented for inversion problems (Lee et al. 2019). Caers and Hoffman introduced the probability perturbation method (PPM) in this regard (Caers and Hoffman 2006). This method is devoted to combine the soft and hard data for the reservoir characterization. This concept was used for further modifications (Khani et al. 2017; Oliveira et al. 2017; Castro et al. 2009; Li and Caers 2011). Hamdi et al. used MPS in facies modeling to match the well test data (Hamdi et al. 2014).

1.3 Permeability prediction by MPS methods

A lot of research is accomplished by using MPS as the categorical property modeler (Caers and Zhang 2004; Gardet et al. 2016; Comunian et al. 2012). Although early MPS methods are limited to the categorical parameters like SNESIM (Strebelle 2002) and FILTERSIM (Zhang et al. 2006) at early stages of development, nowadays, more rigorous methods are available to use for both categorical (like facies) and also continuous variables (like permeability) of the reservoir. In addition, many papers are published on pore-scale modeling via MPS methods (Hajizadeh et al. 2011; Okabe and Blunt 2007; Naraghi et al. 2017; Ji et al. 2018).

Park et al. investigated challenges for applying TIs in MPS. They showed that the TI must be consistent with the production data (Park et al. 2013). They discussed the qualitative fashion of geoscientists who ignore great detail or anything quantitative in production history in interpreting the reservoir facies heterogeneity. Many considerations should be applied while selecting a TI for reservoir modeling.

In this study, we investigate the effect of horizontal rock heterogeneity, in other words, rock texture, to scrutinize the effect of its texture on the performance of hydrocarbon production. To this end, a modified version of SPE case 10 (Christie and Blunt 2001) is used for the simulation purposes. Having this modified model proposed for fluid flow between production and injection wells, we inspect the front of the injected fluid through the reservoir during oil production. In fact, only reservoir dimension and well rates are modified to inspect the flood front movement regarding rigorous modeling of reservoir permeability. Actually, MPS methods enable us to see what exactly happens regarding the horizontal heterogeneity of permeability. However, this application of MPS has not been investigated while assuming all real reservoir dimensions and reservoir fluid velocities. Previous surveys have mainly focused on the pore-scale models or introduction of a MPS method and have not engaged permeability modeling as a continuous variable at a real reservoir scale with real injection and production conditions which means they have applied facies modeling instead of permeability modeling. Hence, a main challenge is to investigate the effect of rock texture on the reservoir performance in real cases and the efficiency of the provided MPS methods for this purpose.

2 Methodology

2.1 Reservoir model

To investigate the effect of rock permeability texture, the modified version of the SPE 10 comparative solution project is employed by only considering the horizontal heterogeneity of permeability regardless its depth variation. This eliminates the impact of gravity on fluid flow through reservoir. This modification contains identical specifications to the original SPE 10 comparative solution project (Christie and Blunt 2001) including relative permeability, rock properties (except permeability), well constraints for production, simulation block size in any direction, etc. Modifications are applied only to the number of simulation cells in each direction and well injection fluid and injection/production rate. In this study, the number of grids is changed from \(100 \times 1 \times 20\) to \(210 \times 210 \times 1\) in the x-, y- and z-directions, respectively, and also the injection rate is changed from 34.61 MSCF/day gas injection to 50 STB/day water injection. Modification in the number of reservoir grids leads to 5250 ft (210 × 25 ft) horizontal distance in both x- and y-directions of the reservoir between the corners of the reservoir. Two wells are placed at the opposite corners of this square reservoir, one for water injection and one for oil production. Figure 3 shows the well configuration and shape of the studied reservoir. By this new reservoir dimension, the water injection front through the reservoir could be easily observed during reservoir production regarding the horizontal heterogeneity of permeability.

Fig. 3
figure 3

Position of production and injection wells for the studied reservoir model

Full size image

In this study, a training image and predetermined (reference) permeability for the reservoir are selected to investigate the effect of reservoir permeability texture on its production performance. Only 0.015% of the total number of permeability values (44100) are assumed to be known for this reservoir, i.e., values at seven places are known. This extremely low percentage of known values (“hard data” in texture synthesis terminology) is due to sparsity of wells in the reservoir. The places of these known permeability values are selected randomly through this reservoir. Figure 4 shows the training image (TI) for this study, reference permeability and places having known values of permeability for the reservoir.

Fig. 4
figure 4

a Training image (TI) used for this study. b Predetermined, in other words, “reference,” reservoir permeability. c Places of known values of permeability

Full size image

2.2 Permeability texture synthesis

The aim of this study is to investigate the effect of horizontal heterogeneity of permeability which represents the effect of permeability texture on the reservoir performance. To do this, several available methods of permeability modeling including a kriging-based (SGSIM method) and six MPS methods are employed to construct the permeability of the reservoir. The six MPS methods are Bunch-DS (Rezaee et al. 2013), image quilting (Mahmud et al. 2014), advanced-DS (Mariethoz et al. 2010), MSCCSIM (Tahmasebi et al. 2014), PCTOSIM (Pourfard et al. 2016) and Graph-cut (Li et al. 2016).

These seven methods are applied to create the reservoir permeability for 44,093 places with seven known permeabilities of reservoir as hard data. The methods borrow all patterns from TI (Fig. 4a training image (TI) used for this study. b Predetermined, in other words, “reference,” reservoir permeability. c Places of known values of permeability), while the SGSIM method uses the information of this training image for variogram model calculations. In other words, the variogram model for permeability modeling is calculated based on this TI. A Gaussian variogram with a range of 13 simulation cells in each grid direction is matched based on the experimental variogram model. By using this representative variogram of the TI, SGSIM will be able to create the realization. Anisotropic variograms do not make sense because the TI has resulted in similar ranges for different directions. Anisotropic variograms may only elongate the shapes in the cases with sparse hard data which is not acceptable having such non-trending TI in this study. The advantage of MPS is to consider large- and small-scale variations in permeability simultaneously which cannot be emulated by using an anisotropic variogram or even a truncated Gaussian simulation. The truncated Gaussian simulation uses lithotypes to deliver a facies map of the reservoir instead of a permeability map by truncating a Gaussian realization into different lithotypes (Hu et al. 2001). Lithotypes are created based on the opinion of the modeler, and there exist no eligible criteria behind the boundary selection of lithofacies. In the advanced usage of truncated methods, the pluri-Gaussian method is implemented for this purpose. The pluri-Gaussian is able to recreate complex arrangements of a number of lithotypes (Oliver et al. 2008). The two last methods are not applicable for creating a realization of permeability since they end up in facies as a categorical output and not permeability as a continuous one. The aim of the present study is to investigate the fluid flow through a continuous spectrum of permeability within the reservoir. Therefore, if permeability values are assigned on the basis of the modeled facies, it will be a self-contradictory calculation. This (self-contradiction) will arise from missing a part of the TI information by boundary selection of lithotypes as a preprocessing for pluri-Gaussian and then reassigning the permeability, perhaps by random values in the specific range of permeability values, to post-process it. The results of permeability modeling for all seven methods in terms of one realization are illustrated in Fig. 5.

Fig. 5
figure 5

Result of permeability modeling in one set of 25 realizations by different methods for the reservoir. a SGSIM method. b Bunch-DS method. c Image quilting method. d Advanced-DS method. e MSCCSIM method. f PCTOSIM method. g Graph-cut method. Comparison is not the objective of this study

Full size image

To evaluate permeability modeling methods, all these methods are used to generate an ensemble of realizations which consist of twenty-five 210 × 210 images (as the reservoir) from a 70 × 70 training image (as an outcrop of the reservoir or derived from seismic calculation). A synthesized texture image with bigger size than that of TI has not been studied before considering all real reservoir conditions. All seven methods employed in this study model the permeability as a continuous variable and are not similar to the previous studies where facies were modeled as categorical variables, e.g., study of Ren et al. (2019). Difficulties in real reservoir simulation still remain. This study avoids these issues by focusing on characterization of the permeability from a fluid flow point of view. Figure 6 shows the flowchart of this paper.

Fig. 6
figure 6

Flowchart of this paper. The objective of this study is not to compare these methods for acceptance or rejection of each of them. For comparing these methods, interested readers should contact the developer of each method. The objective is to reveal new usage of them for continuous parameters (in this case, permeability) by introduction of the permeability texture concept in real reservoir simulation cases. MPS methods are famous for facies modeling which is not relevant to this study. In fact, this study proposes the usage of multiple point statistics alongside the kriging-based methods in a complementary manner for modeling the permeability

Full size image

Tan et al. provided a measure for finding the similarity of textures (Tan et al. 2014). This comparison method calculates distances of cluster-based histograms of patterns by a statistical measure of distance termed the Jensen–Shannon divergence. More details are available in the literature (Tan et al. 2014). After calculating the distances, multidimensional scaling (MDS) is subsequently applicable for visualizing purposes (Honarkhah and Caers 2010). Using the Tan et al. method and MDS calculation, one could plot the similarity of the central part of the obtained realizations compared to the TI. This method of visualization for showing the similarity of patterns needs images in the same size for training and also realization. Thus, the center part of each realization is selected. The central part of each reservoir, which starts with cell number 71 and ends in cell number 140 in both x- and y-directions, is selected to find the similarities with the training image.

In this part, first the output of six MPS methods and the SGSIM method for modeling the permeability of the reservoir is analyzed, and second, the effect of permeability texture is discussed by running the reservoir simulation for all seven permeability modeling methods and predetermined permeability of the reservoir.

3 Permeability modeling results

Patterns of the TI are obvious for all methods for permeability modeling except the SGSIM method as shown in Fig. 5. The SGSIM method leads to circular regions of low (or high) permeability in some areas with smooth variation in permeability, but other methods create the patterns of the TI in the reservoir. This smoothness is reported in various studies of MPS method developments (Tahmasebi et al. 2012a; Mahmud et al. 2014; Mariethoz and Caers 2014). Patterns similar to that of TI are observed in Fig. 5b-g, which belongs to MPS methods, while Fig. 5a, which belongs to the kriging-based method, lacks TI patterns and has been constructed with the routine method used in the oil and gas industry. This pattern preservation is studied for simulation purposes in the next part of this study and will be plotted alongside the routine kriging-based SGSIM method. Unlike the other five MPS methods, the PCTOSIM method results in blurred permeability because of its optimization paradigm (Pourfard et al. 2016). In the constructed permeabilities of the other five methods, some artifacts including a noisy value in smooth areas or visible shapes of patchiness are also observed. For more information about patchiness, refer to the literature (Tahmasebi and Sahimi 2016a).

The synthetic reservoir and TI both have a bimodal histogram of permeability. Figure 7 shows the histogram of one set (out of 25) of the methods for predicting the reservoir permeability. Examining the histograms of the figure reveals the similarity of synthesized images to each other, while the PCTOSIM method shows a quite smooth histogram compared to the other ones. This is due to the different nature of this method to the other five MPS methods. The PCTOSIM method calculates permeability values by optimizing an initially random texture which should be in accordance with TI patterns, but other MPS methods fill the initial empty reservoir with values and patterns stored in the TI. In this aspect, five MPS methods conduct the permeability modeling like SGSIM by starting with an empty reservoir and ending in completed filled reservoir with permeability values drawn from the TI, whereas PCTOSIM starts with random values of reservoir permeability at the beginning and ends with patterns existing in the TI. Thus, sharp modes are not seen in the PCTOSIM histogram. MPS methods should be able to synthesize large, i.e., “dense” in texture synthesis terminology, grids of various reservoirs for reservoir simulation purposes. This study shows that MPS methods are suitable for situations where grid numbers are higher than TI for the permeability modeling. Note that in almost all previous studies, the grid numbers of TI and synthesized image were the same (Mahmud et al. 2014; Tahmasebi et al. 2014). By visual judgment of all histograms in Fig. 7, it could be said that all predicted synthesized images of MPS methods have similar histograms compared to TI and the predetermined permeability values except the SGSIM method and the more smoothed one for the PCTOSIM method. Thus, all methods could be employed for further simulation purposes. Figure 8 shows the variogram of all methods in addition to that of the TI and the reference case. According to this figure, the general shape of the variograms is similar for all synthesized models and the reference case. Correlation length of the training image and the reference case, around 16 simulation cells, is comparable to all MPS-related permeability fields as well as the SGSIM value. It could be concluded that no difference is observed from a variogram point of view.

Fig. 7
figure 7

Histogram of all methods for modeling the reservoir permeability in one set of 25 realizations, training image and predetermined reservoir permeability. SGSIM method and all MPS methods have 44,100 values (210 × 210 image size is the size of studied reservoir) for this histogram. The training image (TI) histogram has different numbers of bars because this image has 4900 values (70 × 70 image size of TI). The training image has the same trend of histogram compared to all synthesized permeabilities and the reference case. Comparison is not the objective of this study

Full size image
Fig. 8
figure 8

Variograms of all synthesized reservoirs as well as the training image and the reference case. The horizontal black dotted line shows the value of the variance for each reservoir permeability. A red vertical line shows the value of the correlation length. The unit for the horizontal axis is the number of simulation cells. Comparison is not the objective of this study

Full size image

4 Reservoir simulation results

Synthesized permeability of six MPS methods along with the SGSIM method results allows us to examine the reservoir performance regarding different permeability modeling (in other words, synthesis) approaches. Having the output of six MPS and SGSIM methods, the reservoir simulator calculates the water flooding front of this reservoir for 46,200 days. This covers the whole life of reservoir.

Permeability modeling was conducted several times in the reservoir studies for considering the uncertainty of reservoir prediction. In this study, we scrutinize the permeability modeling with a novel implementation of MPS in order to parallel the traditional kriging-based method output with MPS method output. Careful observation of the fluid flow is obtained during the reservoir simulation by understanding its dependence on horizontal heterogeneity of permeability. By using MPS methods, more accurate variation in permeability is modeled for the simulation purposes. This is shown in Fig. 5 by visualizing the patterns of TI and synthesized permeabilities.

Simulation results (top views) of oil production at various times in terms of water saturation are shown in Fig. 9 for one set of 25 realizations calculated for each method. As depicted in Fig. 9, differences in water front movement between SGSIM and MPS methods are obvious. Based on this figure, fluid flow through porous media is very similar to the reference case of permeability (predetermined permeability) and MPS methods, especially the Graph-cut method. The edge of the flood front is the same for the reference and MPS methods but differs from that of SGSIM and PCTOSIM methods. Employed MPS methods are able to do this because they preserve the patterns of TI to some extent. However, the kriging-based method, which only considers distances for permeability calculation regardless of direction, is unable to fabricate the variation patterns of permeability. PCTOSIM is not able to mimic the edge of front because of the blurred effect in its output. Note that all input parameters of this reservoir simulation are identical for the reference case and all seven methods of permeability modeling. Basically, only the texture of permeability (high order of spatial variation in permeability) differs from case to case.

Fig. 9
figure 9

Top view for water saturation in the studied reservoir simulation for SGSIM and all six MPS methods and the reference permeability of the reservoir at three different times of simulation. Each row belongs to each method labeled on the left. The left column, middle and the right one are the early, middle and the end of simulation time, respectively. Comparison is not the objective of this study

Full size image

It is seen in Fig. 9 that trapping and bypassing patterns of oil during water flooding are not similar to the reference case for the SGSIM method. However, other MPS methods except PCTOSIM show the same volumetric amount and behavior of trapping as the reference case during the reservoir performance in visual judgment. Trapping depends not only on the horizontal heterogeneity of permeability but also on the velocity of displacing and displaced fluids and their mobility which is identical for all eight cases of simulation. Thus, the only reason for trapping is the heterogeneity of permeability, i.e., permeability texture, through the reservoir.

Figure 9 shows rather similar flood front advancements for all saturations including low, average and high water saturations, especially for the Graph-cut method. High water saturations, which occurred around the water injection well, were shown by red and dark red colors. SGSIM does not show the same transition from high to low values of water saturation from the injection well toward the edges of the front. It was an arc shape transition for SGSIM. Usually, saturations at quite high or very low values of water saturations do not cause any trouble for history matching purposes. On the contrary, blocks with saturations in the middle (not too much oil or water in the simulation cell) usually cause problematic history matching calculations. History matching and model calibration are, in essence, identical calculations applied, respectively, in reservoir engineering and hydrology disciplines to match the output of models with real data from earth layers. This figure shows that in middle saturations, permeability texture causes the front to pass through the reservoir which is identical to the reference case having the same permeability texture. In other words, the valuable information of permeability texture definitely emulates the flood front in real cases.

By observing water saturations at the middle and end of simulation carefully, it can be realized that dot-like zones of oil trapping are identical for the reference case and MPS methods except PCTOSIM. This dot-like trapping is because of permeability texture of the reservoir which could only be reproduced by the MPS methods. PCTOSIM creates blurred permeability which will not lead to the correct simulation of oil trapping in an objective point of view in terms of the reservoir performance. Having artifacts of MPS methods like the ones seen in the resulted permeability of five MPS methods except PCTOSIM is rather acceptable for this case from this objective point of view.

Figure 10 shows oil production for 25 realizations alongside their P10 and P90 values as well as the reference curve. P10 and P90 are percentiles of 25 realizations in each time of simulation for ten and ninety percent of data. It shows that the reference case lies most of the time between P10 and P90 values of the ensemble of realizations which suggests the correctness of these methods by avoiding deviation from the reference curve trend. Meanwhile, SGSIM, advanced-DS and PCTOSIM reveal less uncertainty in their results by having rather similar values of P10 and P90. Figure 11 shows the aforementioned explanation for the field water production rate (FWPR). The reference curve for water production lies again between P10 and P90 curves. However, the diversity of realization results depends on the selected method.

Fig. 10
figure 10

Field oil production rate (FOPR) for all tested methods. For all models, the oil rate starts with approximately 50 STB/day, and afterward, it declines. The realization values are drawn in background. The light blue line is the P90 values, and the red line is the P10 values for each ensemble of realizations. The reference value is shown by black dots with yellow shading. Comparison is not the objective of this study. Noted that the horizontal axis for this figure covers the whole life of the reservoir

Full size image
Fig. 11
figure 11

Field water production rate for all tested methods. For all models, the water production starts approximately 20,000 days after the start of simulation. The realization values are drawn in the background. The blue line is the P90 values, and the red line is P10 values for each ensemble of realizations. The reference value is shown by black dots with yellow shading. Again, comparison is not the objective of this study. Noted that the horizontal axis for this figure covers the whole life of the reservoir

Full size image

Figure 12 shows the difference between P10 and P90. The beginning of water production (water breakthrough) causes a sharp rise in this difference. This difference could be regarded as an indication of uncertainty. More uncertainty is expected as this difference rises. It is noted that a sharp decrease occurs after start of water production, especially for the PCTOSIM method (blue arrow) and the advanced-DS method (black arrow). This behavior of rapid reduction in uncertainty has great importance in reservoir engineering problems which is not possible with the SGSIM method. A decreasing trend in difference between P90 and P10 shows the validity of this study since, as the time passed, less oil remains in the reservoir to be produced and consequently less uncertainty would occur.

Fig. 12
figure 12

Difference between P90 and P10 values of field oil production rate for all methods. P90 and P10 are shown in Fig. 10 with light blue and red curves. Noted that the horizontal axis for this figure covers the whole life of the reservoir

Full size image

Figure 13 presents the same values for the water production rate. PCTOSIM shows the most certain behavior since it has the sharpest decrease, i.e., least uncertainty (blue arrow). advanced-DS also shows promising results in this perspective (black arrow). It should be noted that the horizontal axis for this figure covers the whole life of the reservoir. In other words, this representation reveals the enduring effect of permeability texture on reservoir performance and related uncertainties.

Fig. 13
figure 13

Difference between P90 and P10 values of field water production rate for all methods. P90 and P10 are shown in Fig. 11 with light blue and red curves. Noted that the horizontal axis for this figure covers the whole life of the reservoir

Full size image

Uncertainty expected in terms of the total oil production is shown in Fig. 14 for this reservoir. PCTOSIM has less uncertainty than SGSIM for all times of simulation. Although other MPS methods preserve the permeability texture, their prediction of total oil production is not certain. While SGSIM does not preserve the permeability texture, it supplies acceptable results. Advanced-DS is the next certain one showing a fairly comparable uncertainty.

Fig. 14
figure 14

Difference between P90 and P10 values of field total oil production for all methods. Noted that the horizontal axis for this figure covers the whole life of the reservoir

Full size image

Figure 15 shows the P50 values of ensembles with the reference curve. This figure shows quite a good match of P50 with the reference curve, especially for SGSIM, image quilting, advanced-DS and PCTOSIM methods. Figure 16 shows the same curves for the water production rate. These two figures suggest the accuracy of results regarding the reference case. P50 is defined such that half of realizations are higher while the other half are lower than it for any time step of simulation. Therefore, P50 can act as the measure of accuracy and reliability of the resulted ensemble of realizations. According to these two figures, accurate and reliable match is obtained by both SGSIM and MPS methods at the first glance for the whole life of the reservoir.

Fig. 15
figure 15

Evaluation of P50 and the reference curve for all methods. The red curve is the reference reservoir oil production rate, and the blue one is P50 of the ensemble of realizations for each method

Full size image
Fig. 16
figure 16

Evaluation of P50 and the reference curve for all methods. The red curve is the reference reservoir water production rate, and the blue curve is P50 of the ensemble of realizations for each method

Full size image

Figure 17 shows the absolute error. In this semilog figure, the vertical axis is the difference between the two curves of Fig. 15 which is proportional to the inverse of accuracy. All methods show the same trend in their accuracies. There is a decrease in accuracy for all methods which happens because of water breakthrough at around 20,000 days of simulation.

Fig. 17
figure 17

Absolute error values of all methods for the field oil production rate between the reference case and P50 of the ensemble of realizations

Full size image

Figure 18 illustrates the aforementioned errors for the total oil production of the reservoir. Advanced-DS and Graph-cut are the most accurate methods among others with an error of 7000 and 6000 STB at most at approximately 37,500 and 22,500 days, respectively. Next is MSCCSIM which is stable between 10,000 and 15,000 STB error and then increases. Advanced-DS, MSCCSIM, image quilting and PCTOSIM are all more accurate than the SGSIM method since their curves are under SGSIM. Finally, Bunch-DS is less accurate than SGSIM.

Fig. 18
figure 18

Absolute error values of all methods for the field total oil production between the reference case and P50 of the ensemble of realizations

Full size image

To visualize the similarities of central part of permeability modeling outputs with patterns of the training image, the Jensen–Shannon divergence calculations and multidimensional scaling are carried out for this purpose. Consequently, distances between the point of the training image and realization points show the similarities of output for each method regarding the patterns inside the training image. Closer nodes to the training image node are more similar to patterns of the TI than others. Figure 19 shows the result of this kind of visualization by considering image quilting and MSCCSIM methods. A different scattering in visualization is observed which means these two methods deliver different levels of similarity to the training image in their results.

Fig. 19
figure 19

2D plot of the first and the second reduced dimensions resulted from multidimensional scaling. This is a representation of central part of the realizations, by considering the training image as the center of coordinates, for image quilting and MSCCSIM methods

Full size image

Figure 20 shows the same visualization for advanced-DS and Bunch-DS methods as the family of direct sampling methods. Both methods show the same trend in their scattering around the training image representative node. In both Figs. 19 and 20, nodes with different distances to the training image node have different degrees of similarity to the training image. In addition, the arrangement of realizations in different directions regarding the training image point (central point in these figures) reveals the ability of MPS methods to generate reliable permeability values which are quite useful for the subsequent uncertainty quantification calculations of the reservoir.

Fig. 20
figure 20

2D plot of the first and the second reduced dimensions resulting from multidimensional scaling. This is a representation of central part of realizations, by considering the training image as the center of coordinates, for advanced-DS and Bunch-DS methods

Full size image

Figure 21 shows this way of visualization for SGSIM, PCTOSIM and graph-cut methods. More similar nodes are shown in this figure for PCTOSIM and graph-cut methods compared to SGSIM method. This is also recognizable from Fig. 5 as one set of 25 realizations. The ability of ensembles of these two methods in predicting accurate total oil production in long-term reservoir simulations (shown as the two best ones in Fig. 18) is due to more realistic results which is also visualized by Tan et al. method for the central part of the realizations.

Fig. 21
figure 21

2D plot of the first and the second reduced dimensions resulted from multidimensional scaling. This is a representation of realizations with respect to the training image for SGSIM, PCTOSIM and graph-cut methods

Full size image

5 Conclusions

This study provides an implementation of the permeability texture effect on the reservoir performance using a modified version of the SPE 10 solution project. All reservoir specifications are identical to the original SPE 10 except the injected fluid, injection rate and the reservoir dimensions. Reservoir dimensions are set such that the effect of horizontal permeability heterogeneity is clearly seen.

Below are important findings of this study:

  • Visual judgment of synthesized permeability reveals fair efficiency of mentioned MPS methods for this reservoir as a pioneer case.

  • From a “histogram matching” point of view, all methods predicted permeability for the studied reservoir with an acceptable general shape of histograms except SGSIM. Meanwhile, PCTOSIM shows a more smoothed histogram than the others.

  • All MPS methods and SGSIM have similar output variograms with close correlation length compared to the reference and training image. This doesn’t make the variogram an informative mean for evaluation of methods.

  • By an objective single realization point of view of this reservoir, MPS methods resulted in a very good match for simulating the trapping phenomena and flood front movement for the reservoir (especially graph-cut method) compared to the reference case except PCTOSIM because of blurred outputs.

  • The amount and volume of trapping in a reservoir are solely due to the reservoir permeability texture under the same conditions of injection/production and fluid velocities. The aforementioned conclusion is not possible by kriging-based methods and PCTOSIM method.

  • The difference between P90 and P10 as a measure of uncertainty shows lower uncertainty for some MPS methods than for the SGSIM method for water and oil production rates.

  • In terms of the total oil production, PCTOSIM and SGSIM have the least uncertainty among others.

  • The difference between P50 and the reference reveals the accuracy of all methods for the total oil production. The value shows that advanced-DS and graph-cut are the most accurate methods and MSCCSIM, PCTOSIM and image quilting are more accurate than the SGSIM method.