Reservoir Modeling Using Multiple-Point Statistics

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.

<p><span>In the beginning of the 2000's [1], multiple-point statistics (MPS) was introduced as a novel geostatistical approach to explore the variability of natural phenomena in a realistic way by observing and simulating data patterns, sensibly improving the preservation of connectivity and shape of the modeled structures.</span></p><p><span>A usual requirement for MPS is the presence of complete and representative training images (TI), showing clear and possibly redundant examples of the studied structures. But in the everyday practice, this information is often partially or scarcely available, strongly limiting the use of MPS.</span></p><p><span>In this presentation we start with an overview of MPS strategies proposed to overcome training data limitations. We consider different examples of multisite rain-gauge networks containing sparse data gaps, with the goal of estimating the missing data, using the same incomplete dataset as TI [2]. Another considered study case regards the use of 2D training images of geological outcrops used to reconstruct a 3D volume of fluvioglacial deposits [3]. </span></p><p><span>We then consider a common problem in hydroclimatological studies: the bias correction of weather radar images with ground rainfall measurements. This is a typical no-TI problem where there is no example of unbiased grid image to train MPS. In this case, we propose a novel pattern-to-point approach, where we create a catalog of local grid patterns, each one associated to a rainfall measurement. This way the MPS algorithm 1) selects ungauged locations, 2) searches similar grid patterns in the catalog, and 3) projects the linked historical ground measurements at the ungauged locations.</span></p><p><span>From early results, this technique seems to recover hidden spatial patterns which correct the highly non-linear bias by extracting information from the pattern-to-point catalog. This is a first step for MPS towards the use of TIs integrating variables of different dimensionality, opening a new methodological path for future research.</span></p><p> </p><p><span>BIBLIOGRAPHY</span></p><p><span>[1] Strebelle, S. "Conditional simulation of complex geological structures using multiple-point statistics." Mathematical geology 34.1 (2002): 1-21.</span></p><p><span>[2] Oriani, F. et al. "Missing data imputation for multisite rainfall networks: a comparison between geostatistical interpolation and pattern-based estimation on different terrain types." Journal of Hydrometeorology 21.10 (2020): 2325-2341.</span></p><p><span>[3] Kessler, T. et al. "Modeling fine</span><span>‐</span><span>scale geological heterogeneity—examples of sand lenses in tills." Groundwater 51.5 (2013): 692-705.</span></p>

Extensive outcrop data or photographs of present day depositions or even simple drawings from expert geologists contain precious structural information about spatial continuity that is beyond the present tools of geostatistics essentially limited to two-point statistics (histograms and covariances). A neural net can be learned to collect multiple point statistics from various training images, these statistics are then used to generate stochastic models conditioned to actual data. In petroleum applications, the methodology developed can be a substitute for objects based-algorithms when facies geometry and reservoir continuity are too complex to be modeled by simple object such as channels. The performance of the neural net approach is illustrated using training images of increasing complexity, and attempts at explaining that performance is provided in each case. The neural net builds local probability distributions for the facies types (categorical case) or for petrophysical properties (continuous case). These local probabilities include multiple point statistics learned from training images and are sampled with a Metropolis-Hastings sampler which also ensures reproduction of statistics coming from the actual subsurface data such as locally varying facies proportions. P. 321

Multiple-point statistics are widely used for the simulation of categorical variables because the method allows for integrating a conceptual model via a training image and then simulating complex heterogeneous fields. The multiple-point statistics inferred from the training image can be stored in several ways. The tree structure used in classical implementations has the advantage of being efficient in terms of CPU time, but is very RAM demanding and then implies limitations on the size of the template, which serves to make a proper reproduction of complex structures difficult. Another technique consists in storing the multiple-point statistics in lists. This alternative requires much less memory and allows for a straightforward parallel algorithm. Nevertheless, the list structure does not benefit from the shortcuts given by the branches of the tree for retrieving the multiple-point statistics. Hence, a serial algorithm based on list structure is generally slower than a tree-based algorithm. In this paper, a new approach using both list and tree structures is proposed. The idea is to index the lists by trees of reduced size: the leaves of the tree correspond to distinct sublists that constitute a partition of the entire list. The size of the indexing tree can be controlled, and then the resulting algorithm keeps memory requirements low while efficiency in terms of CPU time is significantly improved. Moreover, this new method benefits from the parallelization of the list approach.

Building a training image with Digital Outcrop Models

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.

We propose a smooth formulation of multiple-point statistics that enables us to solve inverse problems using gradient-based optimization techniques. We introduce a differentiable function that quantifies the mismatch between multiple-point statistics of a training image and of a given model. We show that, by minimizing this function, any continuous image can be gradually transformed into an image that honors the multiple-point statistics of the discrete training image. The solution to an inverse problem is then found by minimizing the sum of two mismatches: the mismatch with data and the mismatch with multiple-point statistics. As a result, in the framework of the Bayesian approach, such a solution belongs to a high posterior region. The methodology, while applicable to any inverse problem with a training-image-based prior, is especially beneficial for problems which require expensive forward simulations, as, for instance, history matching. We demonstrate the applicability of the method on a two-dimensional history matching problem. Starting from different initial models we obtain an ensemble of solutions fitting the data and prior information defined by the training image. At the end we propose a closed form expression for calculating the prior probabilities using the theory of multinomial distributions, that allows us to rank the history-matched models in accordance with their relative posterior probabilities.

We seek to develop a methodology enabling fast geostatistical simulations honoring both geophysical data and a complex prior model. Particularly, we consider a multiple-point statistics (MPS) framework in which a training image (TI) describes the available prior knowledge. Accurate posterior sampling is then possible by using a so-called extended Metropolis algorithm in which proposals are drawn from the prior using sequential geostatistical resampling. Such a Markov chain Monte Carlo (MCMC) algorithm will eventually locate and sample proportionally to the posterior distribution, however, it is often exceedingly slow and typically demands millions of MCMC iterations before the posterior is sampled sufficiently. We are developing a methodology in which the MPS simulation is built up iteratively pixel-by-pixel starting from an empty grid. At each pixel, multiple proposals are generated using an MPS algorithm and the proposals are accepted proportionally to the likelihood considering conditioning data in terms of linear averages (for instance geophysical data). The likelihood function is generally intractable as it depends on the pixels that have not yet been sampled. We approximate the likelihood function using a Gaussian model in which the posterior mean and covariance are updated sequentially as the simulation builds up. The posterior statistics are approximated by running the algorithm multiple times (sequentially or in parallel). Considering crosshole first-arrival ground-penetrating radar data, we assess the accuracy of our methodology both for multi-Gaussian priors for which analytical posteriors are available and for more complex training images against the extended Metropolis method. Our approach is inherently approximate due to the use of a finite training image, a finite number of candidates for each pixel and the need to approximate intractable likelihood functions. Nevertheless, preliminary results are promising as this method allows directly obtaining a reasonable estimation at a reduced computational cost compared to MCMC.

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)

When using coarse-resolution remote sensing images, super-resolution reconstruction is widely desired, and can be realized by reproducing the intrinsic features from a set of coarse-resolution fraction data to fine-resolution remote sensing images that are consistent with the coarse fraction information. Prior models of spatial structures that encode the expected features at the fine (target) resolution are helpful to constrain the spatial patterns of remote sensing images to be generated at that resolution. These prior models can be used properly by multiple-point statistics (MPS), capable of extracting the intrinsic features of patterns from prior models such as training images, and copying them to the simulated regions using hard and soft conditional data, or even without any conditional data. However, because traditional MPS methods based on linear dimensionality reduction are not suitable to deal with nonlinear data, and isometric mapping (ISOMAP) can reduce the dimensionality of nonlinear data effectively, this paper presents a sequential simulation framework for generating super-resolution remote sensing images using ISOMAP and MPS. Using four different examples, it is demonstrated that the structural characteristics of super-resolution reconstruction of remote sensing images using this method, are similar to those of training images.

A nearest neighbor multiple-point statistics method for fast geological modeling

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

Accurate prediction of petroleum reservoir performance requires reliable models of the often complex reservoir heterogeneity. Geostatistical simulation techniques generate multiple realizations of the reservoir model, all equally likely to be drawn. Traditional to geostatistics, geological continuity is represented through the variogram. The variogram is limited in describing complex geological structures as it measures correlation between rock properties at two locations only: it is a two- point statistic. Reservoir analogs such as outcrops can serve as training images depicting the interpreted geological structure. Due to scarcity of well data, the variogram models are often borrowed from such training sets. However, the same training images could be utilized to extract more complex information in the form of multiple point statistics measuring the joint dependency between multiple locations. This paper compares a traditional variogram-based geostatistical model versus a novel geostatistical method utilizing multiple point statistics borrowed from training images. The comparison is made on the basis of flow performance for a typical North Sea reservoir. To obtain such comparison a "true" reference reservoir is generated using object-based simulation that depicts the complex intertwining of fluvial channels. Next, a different but similar reservoir is generated, termed the "training reservoir". The latter is used to extract the necessary structural information, be it variograms or multiple-point statistics, to build multiple geostatistical models of the "true" reservoir conditioned to sparse well data. A water flood flow scenario with inverted 5-spot pattern is simulated using ECLIPSE on the true reference and the various geostatistical models. Water breakthrough characteristics and water saturation distributions are used for comparison.

The reconstruction of porous media is of great importance in predicting fluid transport properties, which are widely used in various fields such as catalysis, oil recovery, medicine and aging of building materials. The real three-dimensional structural data of porous media are helpful to describe the irregular topologic structures of porous media. By using multiple-point statistics (MPS) to extract the characteristics of real porous media acquired from micro computed tomography (micro-CT) scanning, the probabilities of structural characteristics of pore spaces are obtained first, and then reproduced in the reconstructed regions. One solution to overcome the anisotropy of training images is to use real 3D volume data as a training image (TI). The CPU cost and memory burden brought up by 3D simulations can be reduced greatly by selecting the optimal multiple-grid template size that is determined by the entropy of a TI. Moreover, both soft data and hard data are integrated in MPS simulation to improve the accuracy of reconstructed images. The variograms and permeabilities, computed by lattice Boltzmann method, of the reconstructed images and the target image obtained from real volume data are compared, showing that the structural characteristics of reconstructed porous media using our method are similar to those of real volume data.

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.