Assessing the Probability of Training Image-Based Geological Scenarios Using Geophysical Data

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.

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.

<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>

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

Two approaches are traditionally used to build numerical models for facies distributions within a reservoir. Pixel-based techniques aim at generating simulated realizations that honor the well data values, and reproduce a given variogram which models two-point spatial correlation. However, because the variogram cannot look at spatial continuity between more than two locations at a time, pixel-based algorithms give poor representations of the actual facies geometries. In contrast, object-based techniques allow reproducing crisp geometries, but the conditioning on well data requires iterative "trial-and-error" corrections, which can be time-consuming, particularly when the data are dense with regard to the average object size. This paper presents a new approach that combines the easy conditioning of pixel-based algorithms with the ability to reproduce "shapes" of object-based techniques, without being too time and memory demanding. In this new approach, the complex geological structures expected to be present in the reservoir are characterized by multiple-point statistics, which express joint variability at many more than two locations at a time. Such multiple-point statistics cannot be inferred from typically sparse well data but could be read from training images depicting the expected subsurface heterogeneities. A training image need not carry any locally accurate information on the reservoir; it need only reflect a prior stationary geological/structural concept. Thus training images can be generated by object-based algorithms freed of the constraint of data conditioning. The multiple-point statistics inferred from the training image(s) are then exported to the reservoir model, where they are anchored to the well data using a pixel-based sequential simulation algorithm. This algorithm is tested for the simulation of a turbidite system where flow is controlled by meandering channels with cross-bedding. The training image reflecting the channel patterns is an unconditional realization generated by an object-based algorithm. The final simulated numerical models reproduce these channel patterns, and honor exactly all well data values at their locations. The methodology proposed appears to be practical, general, and fast.

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.

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.

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.

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

For more than half a century, geostatistics has been developed and has increasingly been used for modeling subsurface heterogeneity. Traditionally, geostatistical simulations are based on a random function model defined according to the specificities of the geological formation under investigation. Unlike traditional geostatistics, multiple‐point (MP) geostatistics avoids the explicit definition of a random function but directly infers the necessary multivariate distributions from training images. This confers on MP geostatistics a potential applicability to any geological environment, provided that there is a training image representative of the geological heterogeneity and that the essential features of this training image can be characterized by statistics defined on a limited point configuration. This paper presents a comprehensive review of MP geostatistics. If the principle of MP geostatistics is straightforward and attractive, its industrial applicability largely depends on the implementation methods. The use of a search tree for storing MP statistics is a great step that made MP geostatistics actually applicable in industry. In the meantime, other methods such as multiple grids and image postprocessing are introduced and allow enhancing the reproduction of patterns observed in training images. Because of the sequential procedure, MP simulations can easily honor local hard data. There are methods available for constraining MP simulations to target global statistics (facies proportions) and also for integrating spatial auxiliary constraints that can tremendously improve the spatial features of MP simulations. Furthermore, both the gradual deformation method and the probability perturbation method are compatible with MP geostatistics and allow integrating hydrodynamic data into MP simulations. More recent developments of MP geostatistics include sequentially simulating patterns instead of points and using different geological scenarios (training images) for dynamic data inversion.

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