Multiple‐point geostatistics for modeling subsurface heterogeneity: A comprehensive review

Most of the world's massive fields produce from naturally fractured vuggy carbonate reservoirs. However, modeling such reservoirs is a challenging task due to the co-existence of porous and free-flow domains over a wide range of scales. Realistic characterization of fractured vuggy reservoirs requires preservation of the geological structure and reproduction of patterns. In this study, we model a discrete fracture-vug network using microtomography and multiple-point geostatistics (MPG) modeling technique. For this purpose, micro-CT scan images of a real core is used to generate two micro-scale geological models (2D and 3D reference models, which are used as original training images) with three structures; matrix, fractures, and vugs. The structures in the reference models are also separated into matrix-fracture and matrix-vug to construct two other training images for each of the 2D and 3D original training images. The three training images are then used in MPG to generate multiple equiprobable models that mimic the original networks of the fractures and vugs. The models created using the original training images are called direct models. The others that are generated using the training images after separating the structures (matrix-fracture and matrix-vug), referred to as separable training images, are combined into fracture-vug-matrix models by superimposing the fracture networks from the matrix-fracture models on their corresponding matrix-vug models. These final three-structure models are called the combined models. The spatial continuity, connectivity, patterns reproduction, and flow performance of the generated models are compared with those of the reference models (original training images). The results show that, the separable modeling method produces better results when used in 2D modeling, however, that is not the case for 3D modeling. All the generated realizations preserved the structures and reproduced the spatial continuity of the reference models for 2 and 3D cases. In addition, the patterns are well reproduced in both cases. However, the connectivity and the flow performance of the traditional (direct) MPG method are found to be more accurate in the 3D case than those of the separable MPG method. These results proof that the separable MPG method is not suitable for real 3D naturally fractured vuggy reservoirs modeling applications. However, more detailed sensitivity analysis of the MPG parameters could provide new insights in utilizing this method for 3D reservoir modeling applications.

This article, written by Special Publications Editor Adam Wilson, contains highlights of paper SPE 190087, “Unsupervised Statistical Learning With Integrated Pattern-Based Geostatistical Simulation,” by Q. Li and R. Aguilera, SPE, University of Calgary, prepared for the 2018 SPE Western Regional Meeting, Garden Grove, California, USA, 22–27 April. The paper has not been peer reviewed. This paper presents a new geostatistics modeling methodology that connects geostatistics and machine-learning methodologies, uses nonlinear topological mapping to reduce the original high-dimensional data space, and uses unsupervised-learning algorithms to bypass problems with supervised-learning algorithms. The algorithm presented is a neural topology-preserving pattern-based geostatistical simulation algorithm that integrates the self-organizing map (SOM) concept and its updated version—growing self-organizing map (GSOM)—with an unsupervised competitive learning structure. Introduction In oil and gas reservoir modeling, any model construction faces challenges of limited data to some extent. The heuristic behind all geostatistical techniques is the implicit existence of statistical relationships among available data. “Data” here is a broad term; it could be discrete points, such as porosity or permeability at certain locations, but it also could be training images (TIs), which are used in this work. Using TIs as input data originated with multiple-point geostatistics. The aim was to overcome the limitations of using traditional two-point statistical variograms to describe geological continuity, especially in the case of curvilinear structures, which are quite common in nature, such as in fracture networks and geological fluvial structures. The authors write that geostatistics could benefit from the machine-learning or statistical-learning areas. Machine-learning tasks can be divided into two protocols, supervised learning and unsupervised learning. The difference depends on whether the input data have correct labels or not. For the investigation considered in this paper, after retrieving image patches from TIs, machine-learning algorithms were used to cluster those image patches into different classifications. If the correct clusters to which the image patches belong is known beforehand, the data can be said to have correct labels. Further, a large amount of these labeled data could be used for a model to learn. Supervised learning occurs when the model is finely tuned by guidance of these correct labels using an error-correction process. Unsupervised learning, on the other hand, occurs when neither how many clusters exist nor the correct clusters to which the image patches belong is known. Thus, without correct labels for guidance, this problem is identified as an unsupervised-learning problem. This is the first major drawback of applying machine-learning algorithms to geostatistical simulations. Two other important issues are encountered when performing pattern-based geostatistical simulations. The first is related to the large number of image patches. This large image-patch database contains high pattern redundancy, which could make pattern similarity comparison inefficient. The second issue is the high dimensionality of each image patch, which could typically be described by low-dimensional nonlinear-structure manifolds. To model and visualize these high-dimensional data, a nonlinearity dimensionality reduction should be sought.

Summary Accurate characterization of sub-surface oil reservoirs is an essential prerequisite to the design and implementation of enhanced oil recovery (EOR) scenarios. Specifically, in reservoir characterization, integrating static and dynamic data into reservoir models to construct accurate and realistic models has received considerable attention. Unlike most of the conventional geostatistical approaches of integrating data into reservoir models that are based on semi-variograms (two-point statistics) as a measure of spatial connectivity, a complete multiple-point (MP) statistic framework is presented in this paper. In contrast to two-point statistic methods, MP statistics-based methods are capable of reproducing curvilinear geological structures. The algorithm starts with extracting MP statistics from training images (TI) using an optimal spatial template. After collecting different patterns and building the MP histogram, the pattern reproduction process commences. This process begins from data locations and then grows to fill the whole reservoir domain. The algorithm accounts for three main practical issues: uncertainty in geological scenarios, scanning template and non-stationarity. The MP statistics algorithm (growthsim) is capable of integrating data from multiple data sources. Among these data types is dynamic data or flow history. The conventional approach to integrate production information into reservoir models is by iterative perturbation of the reservoir model until the production history of the reservoir is matched. Iterative methods have been applied till date to random fields that are completely characterized by a two-point co-variance function. In contrast, this paper presents a forward modelling approach that investigates history matching within a MP modelling framework. A novel technique implemented in this research is based on the merging of MPs inferred from history matched and geological models. Pattern growth is performed subsequently by sampling from the merged MP histograms. History matched models using the presented approach show an excellent agreement with underlying geological descriptions and match production history.

The multiple point geostatistics can not reproduce the continuous shape of the channel due to the limitation of the algorithm with random draw of the simulated value. A new method is proposed which integrates the random walk process with multiple point geostatistics. First the random walk process produces the channel main streamline, then the draw of simulated value in multiple point geostatistics is conducted by the streamline. The method solves the continuities problem in channel. The real model is constructed with both the new method and only the multiple point geostatistics and the results show the new method can reproduce the continuous sinuous channel. A cross-validation is implemented and the result is also well, which proves the new method is robust and can be put into real reservoir modeling.

SummaryHistory matching forms an integral part of the reservoir modeling workflow process. Despite the existence of many history-matching tools, the integration of production data with seismic and geological continuity data remains a challenge. Geostatistical tools now are routinely employed for integrating large-scale seismic and fine-scale well/core data. A general framework for integrating production data with diverse types of geological/structural data is largely lacking. In this paper, we develop a new method for history matching that can account for production data constraint by prior geological data, such as the presence of channels, fractures, or shale lenses. With multiple-point (MP) geostatistics, prior information about geological patterns is carried by training images from which geological structures are borrowed, then anchored to the subsurface data. A simple Markov chain is proposed to iteratively modify the MP geostatistical realizations until history match. The method is simple and general in the sense that the procedure can be applied to a wide variety of geological environments without requiring a modification of the algorithm. The method does not make assumptions on the flow model. Synthetic cases are used to assess the flexibility of the proposed method.

In many reservoir-modeling applications, geological uncertainty is treated by considering multiple realizations generated from a specified geological scenario. In reality, however, the geological scenario itself is uncertain, and the use of qualitative criteria for its specification may lead to inaccuracy in flow predictions. In this work, we introduce a systematic procedure for the determination of the most likely a posteriori geological scenario. As is common in geomodeling applications, the scenario is defined in terms of a training image (TI). We introduce continuous parameterizations for uncertain TI attributes such as channel thickness and orientation, and then determine optimal values for these and other key model quantities. Optimality is defined here in terms of the level of agreement between observed data and flow results for realizations generated from a given geological scenario. The optimum scenario is determined using Particle Swarm Optimization. In the second step of the methodology, a set of specific realizations, which provide closer agreement with observed data, is identified using a rejection-sampling method. The workflow is applied to a synthetic channelized system, and the procedure is repeated using several different `true' reservoir realizations to gauge its performance with data from realizations that are more or less representative of the true scenario. Values for TI parameters found by our procedure are shown to be in reasonable agreement with those of the true scenario in all cases considered. In addition, following the identification of specific realizations using rejection sampling, predicted flow results are shown to be of similar quality to those from the true scenario.

Thin, irregularly shaped surfaces such as clay drapes often have a major control on flow and transport in heterogeneous porous media. Clay drapes are often complex, curvilinear three-dimensional surfaces and display a very complex spatial distribution. Variogram-based stochastic approaches are also often not able to describe the spatial distribution of clay drapes since complex, curvilinear, continuous, and interconnected structures cannot be characterized using only two-point statistics. Multiple-point geostatistics aims to overcome the limitations of the variogram. The premise of multiple-point geostatistics is to move beyond two-point correlations between variables and to obtain (cross) correlation moments at three or more locations at a time using training images to characterize the patterns of geological heterogeneity. Multiple-point geostatistics can reproduce thin irregularly shaped surfaces such as clay drapes, but this is often computationally very intensive. This paper describes and applies a methodology to simulate thin, irregularly shaped surfaces with a smaller CPU and RAM demand than the conventional multiple-point statistical methods. The proposed method uses edge properties for indicating the presence of thin irregularly shaped surfaces. Instead of pixel values, edge properties indicating the presence of irregularly shaped surfaces are simulated using snesim. This method allows direct simulation of edge properties instead of pixel properties to make it possible to perform multiple-point geostatistical simulations with a larger cell size and thus a smaller computation time and memory demand. This method is particularly valuable for three-dimensional applications of multiple-point geostatistics.

Multiple-point geostatistical simulation has recently become popular in stochastic hydrogeology, primarily because of its capability to derive geologically reasonable patterns and multivariate distributions from a training image and conditioning training image data to multiple hard and soft data sources. This article resents, evaluates and contrasts the results of using multiplepoint geostatistical simulation for producing geologically realistic models of a Quaternary inner-alpine aquifer. Borehole data, expert-designed geological profiles and training image data were subject to conditional simulation with the Single Normal Equation Simulation (SNESIM) algorithm, one of the most widely used Matrix Product State (MPS) algorithms. The sensitivity of model predictions to the training image and hard as well as soft data input was evaluated. Modeling results indicate that soft conditioning in MPS is a convenient and efficient way for integrating secondary data such as geological expert drawings. Keywords: Multiple point statistics; Secondary data integration; Aquifer heterogeneity; Quaternary over-deepened alpine basin

Reservoirs result from complex and intricate natural formation processes, and are therefore very heterogeneous at multiple scales. As fluid displacements are very sensitive to heterogeneity, it is crucial to understand and characterize how petrophysical properties vary in space to be able to predict oil production. Multiple point statistics has become a common tool to build models representing geological formations comprising objects characterized by complex geometries such as curvilinear shapes. The idea consists in referring to a training image depicting the expected geological structures, and then in generating images reproducing the multiple point statistics inferred from the training image. We investigate a substitute to the commonly referenced multiple-point methodologies that was initially developed in computer graphics: it is called texture synthesis. We modify the way it works so that it can handle reservoir simulation issues. Additionally, we focus on its multiscale variant to make it more efficient. Within this multiscale context, the basic training image is used to create a set of images with decreasing resolutions that are all viewed as training images. Then, a realization is simulated one scale at a time, starting from the level of lowest resolution. When simulating the realization at a given scale, the selection of the value to be assigned to a grid block involves various levels of comparison. We still compare the neighbors of the grid block under consideration with the events extracted from the training image corresponding to the same level of resolution as done for multipoint simulation. In addition, we compare its ancestors with the events extracted from the coarser training images. The ancestors are the values already simulated to populate the coarser grids. In such conditions, the first realization generated at the coarsest scale is able to capture the largest geological objects or heterogeneities. Then, it is used as a constraint that is propagated scale by scale towards the finest resolution grid. This makes it possible to use small size neighborhoods. A few improvements are then suggested to make the algorithm more efficient without deteriorating the simulated images. Last, we present numerical tests performed with complex training images.

Transient well test data conveys significant information about the subsurface heterogeneities in terms of some variations in the well test pressure response curves. It is therefore important to enhance the use of the well test data for building a validated geological model to include the effective reservoir heterogeneities that are reflected on the well test plots. In this work, we present a novel geoengineering workflow for geologically consistent updating of the geostatistical facies models using pressure transient data. We use Multi-Point Statistical (or Geostatistical) simulations (MPS) with conditioning hard and soft data to generate the geostatistical realizations that can preserve the spatial connectivity of the facies. Static model transient tests are then generated using high resolution numerical simulations. The results are compared with the measured well test data for an inversion. The inversion step involves a geologically consistent Probability Perturbation Method (PPM) for perturbing the geostatistical models which are combined with a Gaussian Process (GP) modeling approach for finding the optimum spatial distribution of the facies and the other unknown model parameters. Conditional two-dimensional models of a low-energy anastomosing channelized model are considered in this study. The results show that using such an approach the spatial variation of the facies is maintained and the transitions across the facies boundaries are consistently preserved. In this paper, the geostatistical models are updated simultaneously with other unknown model parameters, including the PPM's parameter (r), facies permeabilities and the non-Darcy D-factor. This multidimensional inversion is efficiently performed by GP in less than 100 compositional simulations. The novelty of this work is to efficiently use the well test data for updating the static models in a fluvial reservoir using a perturbation of the geological models. Moreover, for the first time, a novel optimization method is combined with PPM to consistently update the model parameters with a limited reservoir simulation budget.

TiConverter: A training image converting tool for multiple-point geostatistics

The construction of training images (TIs) depicting the geological prior is one of the most critical step in multiple-point statistics. Geophysical techniques may be used to reduce the uncertainty in the understanding of prior geological scenarios. We developed a methodology to verify the consistency of geophysical data with independently-built TIs. Synthetic geophysical models built from TI scenarios are compared, using multidimensional scaling, with inverted models from field surveys to check if TIs are consistent with geophysical models. Then, the probability of each TI scenario is computed. A cluster analysis enables to determine which parameters used in building the TIs are most impacting the geophysical response. The methodology is tested using ERT to analyze TI scenarios in the Meuse River alluvial aquifer (Belgium)

Reservoir models must be fully updateable: when new data become available, the optimal model identified from a first matching process must be modified to honor also these newly collected data. Therefore, a new matching process has to be run starting from the previously matched model. A facies realization can be modified by varying either facies proportions or facies spatial distributions. The gradual deformation method was introduced to perturb facies distributions (Hu, 2000). This parameterization technique involves a reduced number of parameters and allows for preserving the spatial variability of facies. The resulting changes can be global or localized over a target sub-domain. Local gradual deformation is possible only when applied to the random numbers used to generate the facies realization. Therefore, a prerequisite is to know the random numbers used to produce the existing matched reservoir model. However, in many cases, especially when the second matching study is performed a significant time after the first one, the random numbers are no longer available. The only known object is the realization describing the spatial distribution of facies. At this point, the problem to be solved consists of identifying random numbers, which allow for generating a stochastic realization close to the one populating the existing geological model. In this paper, we develop different methods for identifying the required random numbers and reconstructing the existing reservoir model. Three cases are considered: the existing model was generated from the Sequential Indicator Simulation algorithm, the truncated Gaussian method, or the pluri-Gaussian method. The potential of our methodology is illustrated from a few realistic examples. A particular application shows how an object-based facies realization can even be modified, while considering a reconstruction process based upon the truncated Gaussian method.

Multiple point statistics (MPS) provides a systematic approach for pattern-based simulation of geologic objects from a conceptual training image (TI). The TI encodes the higher-order spatial statistics of the expected connectivity structures through stationary patterns representing the underlying geologic features. The pattern-imitating nature of MPS simulation implies that the simulated facies inherit the spatial structure of the general features in the TI. This property makes the MPS approach very sensitive to uncertainty in the prior TI. Since TIs are constructed using uncertain data and imperfect assumptions, multiple TIs may be necessary to account for the uncertainty and full range of structural variability in facies descriptions. We present a Bayesian mixture modeling approach for adaptively sampling conditional facies from multiple uncertain TIs using a probability conditioning method (PCM). Using the PCM, we invert the flow data to obtain a facies probability map for drawing conditional facies realizations from each TI. The number of samples drawn from each TI is proportional to the weight assigned to them. The TI weights are assigned based on the predictive performance of its corresponding conditional facies realizations. We demonstrate the suitability of the proposed method using numerical experiments in fluvial formations.

Downscaling in remote sensing