History matching of statistically anisotropic fields using the Karhunen-Loeve expansion-based global parameterization technique

Traditional ensemble-based history matching method, such as the ensemble Kalman filter and iterative ensemble filters, usually update reservoir parameter fields using numerical grid-based parameterization. Although a parameter constraint term in the objective function for deriving these methods exists, it is difficult to preserve the geological continuity of the parameter field in the updating process of these methods; this is especially the case in the estimation of statistically anisotropic fields (such as a statistically anisotropic Gaussian field and facies field with elongated facies) with uncertainties about the anisotropy direction. In this work, we propose a Karhunen-Loeve expansion-based global parameterization technique that is combined with the ensemble-based history matching method for inverse modeling of statistically anisotropic fields. By using the Karhunen-Loeve expansion, a Gaussian random field can be parameterized by a group of independent Gaussian random variables. For a facies field, we combine the Karhunen-Loeve expansion and the level set technique to perform the parameterization; that is, for each facies, we use a Gaussian random field and a level set algorithm to parameterize it, and the Gaussian random field is further parameterized by the Karhunen-Loeve expansion. We treat the independent Gaussian random variables in the Karhunen-Loeve expansion as the model parameters. When the anisotropy direction of the statistically anisotropic field is uncertain, we also treat it as a model parameter for updating. After model parameterization, we use the ensemble randomized maximum likelihood filter to perform history matching. Because of the nature of the Karhunen-Loeve expansion, the geostatistical characteristics of the parameter field can be preserved in the updating process. Synthetic cases are set up to test the performance of the proposed method. Numerical results show that the proposed method is suitable for estimating statistically anisotropic fields.     

Application of the EnKF and Localization to Automatic History Matching of Facies Distribution and Production Data

The performance of the ensemble Kalman filter (EnKF) for continuous updating of facies location and boundaries in a reservoir model based on production and facies data for a 3D synthetic problem is presented. The occurrence of the different facies types is treated as a random process and the initial distribution was obtained by truncating a bi-Gaussian random field. Because facies data are highly non-Gaussian, re-parameterization was necessary in order to use the EnKF algorithm for data assimilation; two Gaussian random fields are updated in lieu of the static facies parameters. The problem of history matching applied to facies is difficult due to (1) constraints to facies observations at wells are occasionally violated when productions data are assimilated; (2) excessive reduction of variance seems to be a bigger problem with facies than with Gaussian random permeability and porosity fields; and (3) the relationship between facies variables and data is so highly non-linear that the final facies field does not always honor early production data well. Consequently three issues are investigated in this work. Is it possible to iteratively enforce facies constraints when updates due to production data have caused them to be violated? Can localization of adjustments be used for facies to prevent collapse of the variance during the data-assimilation period? Is a forecast from the final state better than a forecast from time zero using the final parameter fields?To investigate these issues, a 3D reservoir simulation model is coupled with the EnKF technique for data assimilation. One approach to enforcing the facies constraint is continuous iteration on all available data, which may lead to inconsistent model states, incorrect weighting of the production data and incorrect adjustment of the state vector. A sequential EnKF where the dynamic and static data are assimilated sequentially is presented and this approach seems to have solved the highlighted problems above. When the ensemble size is small compared to the number of independent data, the localized adjustment of the state vector is a very important technique that may be used to mitigate loss of rank in the ensemble. Implementing a distance-based localization of the facies adjustment appears to mitigate the problem of variance deficiency in the ensembles by ensuring that sufficient variability in the ensemble is maintained throughout the data assimilation period. Finally, when data are assimilated without localization, the prediction results appear to be independent of the starting point. When localization is applied, it is better to predict from the start using the final parameter field rather than continue from the final state.     

Continuous Facies Updating Using the Ensemble Kalman Filter and the Level Set Method

We present a methodology based on the ensemble Kalman filter (EnKF) and the level set method for the continuous model updating of geological facies with respect to production data. Geological facies are modeled using an implicit surface representation and conditioned to production data using the ensemble Kalman filter. The methodology is based on Gaussian random fields used to deform the facies boundaries. The Gaussian random fields are used as the model parameter vector to be updated sequentially within the EnKF when new measurements are available. We show the successful application of the methodology to two synthetic reservoir models.     

A dynamo in a torus as an explanation of magnetic fields in the outer rings of galaxies

It is currently generally believed that magnetic fields in the disks of spiral galaxies are generated by the dynamo mechanism, which is based on the joint action of differential rotation and the alpha effect, associated with turbulent motions in the interstellar gas. Together with their disks, outer rings are also encountered in galaxies, where magnetic fields may be present. In earlier studies, the generation of magnetic fields has been described in a planar approximation, whose essence is that the size of rings perpendicular to the plane of the galaxy is much smaller than their size in the radial direction. However, it is plausible that these sizesmay sometimes be comparable, so that it would be more logical to suppose that a ring has a toroidal form. A model for a dynamo in a toroidal ring is constructed in this study. This model describes the magnetic field using two functions, corresponding to the toroidal component of the field and the part of the vector potential characterizing its poloidal component. The possible generation of magnetic field in various cases is shown, with both quadrupolar symmetry (close to the fields obtained in the planar approximation) and dipolar symmetry (when two layers with oppositely directed magnetic fields form in the ring). The parameter values for which the generation of fields with one or the other type of symmetry is possible are estimated. The results can also be used to describe the evolution of the magnetic fields in other toroidal astrophysical objects.     

Application of EM algorithms for seismic facices classification

Identification of the geological facies and their distribution from seismic and other available geological information is important during the early stage of reservoir development (e.g. decision on initial well locations). Traditionally, this is done by manually inspecting the signatures of the seismic attribute maps, which is very time-consuming. This paper proposes an application of the Expectation-Maximization (EM) algorithm to automatically identify geological facies from seismic data. While the properties within a certain geological facies are relatively homogeneous, the properties between geological facies can be rather different. Assuming that noisy seismic data of a geological facies, which reflect rock properties, can be approximated with a Gaussian distribution, the seismic data of a reservoir composed of several geological facies are samples from a Gaussian mixture model. The mean of each Gaussian model represents the average value of the seismic data within each facies while the variance gives the variation of the seismic data within a facies. The proportions in the Gaussian mixture model represent the relative volumes of different facies in the reservoir. In this setting, the facies classification problem becomes a problem of estimating the parameters defining the Gaussian mixture model. The EM algorithm has long been used to estimate Gaussian mixture model parameters. As the standard EM algorithm does not consider spatial relationship among data, it can generate spatially scattered seismic facies which is physically unrealistic. We improve the standard EM algorithm by adding a spatial constraint to enhance spatial continuity of the estimated geological facies. By applying the EM algorithms to acoustic impedance and Poisson’s ratio data for two synthetic examples, we are able to identify the facies distribution.     

A robust cylindrical fitting to point cloud data

Environmental, engineering and industrial modelling of natural features (e.g. trees) and man-made features (e.g. pipelines) requires some form of fitting of geometrical objects such as cylinders, which is commonly undertaken using a least-squares method that—in order to get optimal estimation—assumes normal Gaussian distribution. In the presence of outliers, however, this assumption is violated leading to a Gaussian mixture distribution. This study proposes a robust parameter estimation method, which is an improved and extended form of vector algebraic modelling. The proposed method employs expectation maximisation and maximum likelihood estimation (MLE) to find cylindrical parameters in case of Gaussian mixture distribution. MLE computes the model parameters assuming that the distribution of model errors is a Gaussian mixture corresponding to inlier and outlier points. The parameters of the Gaussian mixture distribution and the membership functions of the inliers and outliers are computed using an expectation maximisation algorithm from the histogram of the model error distribution, and the initial guess values for the model parameters are obtained using total least squares. The method, illustrated by a practical example from a terrestrial laser scanning point cloud, is novel in that it is algebraic (i.e. provides a non-iterative solution to the global maximisation problem of the likelihood function), is practically useful for any type of error distribution model and is capable of separating points of interest and outliers.     

Geostatistical Modelling of Cyclic and Rhythmic Facies Architectures

A Pluri-Gaussian method is developed for facies variables in three dimensions to model vertical cyclicity related to facies ordering and rhythmicity. Cyclicity is generally characterised by shallowing-upward or deepening-upward sequences and rhythmicity by the repetition of facies at constant intervals along sequences. Both of these aspects are commonly observed in shallow-marine carbonate successions, especially in the vertical direction. A grid-free spectral simulation approach is developed, with a separable covariance allowing a dampened hole-effect to capture rhythmicity in the vertical direction and a different covariance in the lateral plane along strata, as in space-time models. In addition, facies ordering is created by using a spatial shift between two latent Gaussian functions in the Pluri-Gaussian approach. Rapid conditioning to data is performed via Gibbs sampling and kriging using the screening properties of separable covariances. The resulting facies transiograms can show complex patterns of cyclicity and rhythmicity. Finally, a three dimensional case study of shallow-marine carbonate deposits at outcrop shows the applicability of the modelling method.     

A probabilistic parametrization for geological uncertainty estimation using the ensemble Kalman filter (EnKF)

In the past years, many applications of history-matching methods in general and ensemble Kalman filter in particular have been proposed, especially in order to estimate fields that provide uncertainty in the stochastic process defined by the dynamical system of hydrocarbon recovery. Such fields can be permeability fields or porosity fields, but can also fields defined by the rock type (facies fields). The estimation of the boundaries of the geologic facies with ensemble Kalman filter (EnKF) was made, in different papers, with the aid of Gaussian random fields, which were truncated using various schemes and introduced in a history-matching process. In this paper, we estimate, in the frame of the EnKF process, the locations of three facies types that occur into a reservoir domain, with the property that each two could have a contact. The geological simulation model is a form of the general truncated plurigaussian method. The difference with other approaches consists in how the truncation scheme is introduced and in the observation operator of the facies types at the well locations. The projection from the continuous space of the Gaussian fields into the discrete space of the facies fields is realized through in an intermediary space (space with probabilities). This space connects the observation operator of the facies types at the well locations with the geological simulation model. We will test the model using a 2D reservoir which is connected with the EnKF method as a data assimilation technique. We will use different geostatistical properties for the Gaussian fields and different levels of the uncertainty introduced in the model parameters and also in the construction of the Gaussian fields.     

Capturing and modelling metre‐scale spatial facies heterogeneity in a Jurassic ramp setting (Central High Atlas,Morocco)

Each simulation algorithm, including Truncated Gaussian Simulation, Sequential Indicator Simulation and Indicator Kriging is characterized by different operating modes, which variably influence the facies proportion, distribution and association of digital outcrop models, as shown in clastic sediments. A detailed study of carbonate heterogeneity is then crucial to understanding these differences and providing rules for carbonate modelling. Through a continuous exposure of Bajocian carbonate strata, a study window (320 m long, 190 m wide and 30 m thick) was investigated and metre‐scale lithofacies heterogeneity was captured and modelled using closely‐spaced sections. Ten lithofacies, deposited in a shallow‐water carbonate‐dominated ramp, were recognized and their dimensions and associations were documented. Field data, including height sections, were georeferenced and input into the model. Four models were built in the present study. Model A used all sections and Truncated Gaussian Simulation during the stochastic simulation. For the three other models, Model B was generated using Truncated Gaussian Simulation as for Model A, Model C was generated using Sequential Indicator Simulation and Model D was generated using Indicator Kriging. These three additional models were built by removing two out of eight sections from data input. The removal of sections allows direct insights on geological uncertainties at inter‐well spacings by comparing modelled and described sections. Other quantitative and qualitative comparisons were carried out between models to understand the advantages/disadvantages of each algorithm. Model A is used as the base case. Indicator Kriging (Model D) simplifies the facies distribution by assigning continuous geological bodies of the most abundant lithofacies to each zone. Sequential Indicator Simulation (Model C) is confident to conserve facies proportion when geological heterogeneity is complex. The use of trend with Truncated Gaussian Simulation is a powerful tool for modelling well‐defined spatial facies relationships. However, in shallow‐water carbonate, facies can coexist and their association can change through time and space. The present study shows that the scale of modelling (depositional environment or lithofacies) involves specific simulation constraints on shallow‐water carbonate modelling methods.     

Gaussian Cosimulation: Modelling of the Cross-Covariance

Whenever two or more random fields are assumed to be correlated in reservoir characterization, it is necessary to generate valid cross-covariance models to describe the relationship. The standard methods for constructing covariance matrices for correlated random fields are not very general. In particular, they do not allow one to specify different auto-covariance models for the two fields. It is not possible, for example, for one field to have a Gaussian auto-covariance and the other an exponential auto-covariance, unless the two fields are uncorrelated. The standard approaches also do not allow for nonsymmetric cross-covariance functions. In this report, I present a straightforward method of cosimulation based on the square root of the auto-covariances. The same approach is used for constructing cross-covariance models for the variables. The approach is quite general and does not require symmetry of the cross-covariance. The modelling of the cross-covariance is illustrated with gamma ray and spontaneous potential logs.     

Isotropic Variogram Matrix Functions on Spheres

This paper is concerned with vector random fields on spheres with second-order increments, which are intrinsically stationary and mean square continuous and have isotropic variogram matrix functions. A characterization of the continuous and isotropic variogram matrix function on a sphere is derived, in terms of an infinite sum of the products of positive definite matrices and ultraspherical polynomials. It is valid for Gaussian or elliptically contoured vector random fields, but may not be valid for other non-Gaussian vector random fields on spheres such as a χ ², log-Gaussian, or skew-Gaussian vector random field. Some parametric variogram matrix models are derived on spheres via different constructional approaches. A simulation study is conducted to illustrate the implementation of the proposed model in estimation and cokriging, whose performance is compared with that using the linear model of coregionalization.     

Multimodal ensemble Kalman filtering using Gaussian mixture models

In this paper we present an extension of the ensemble Kalman filter (EnKF) specifically designed for multimodal systems. EnKF data assimilation scheme is less accurate when it is used to approximate systems with multimodal distribution such as reservoir facies models. The algorithm is based on the assumption that both prior and posterior distribution can be approximated by Gaussian mixture and it is validated by the introduction of the concept of finite ensemble representation. The effectiveness of the approach is shown with two applications. The first example is based on Lorenz model. In the second example, the proposed methodology combined with a localization technique is used to update a 2D reservoir facies models. Both applications give evidence of an improved performance of the proposed method respect to the EnKF.     

Non-stationary Geostatistical Modeling Based on Distance Weighted Statistics and Distributions

A common assumption in geostatistics is that the underlying joint distribution of possible values of a geological attribute at different locations is stationary within a homogeneous domain. This joint distribution is commonly modeled as multi-Gaussian, with correlations defined by a stationary covariance function. This results in attribute maps that fail to reproduce local changes in the mean, in the variance and, particularly, in the spatial continuity. The proposed alternative is to build local distributions, variograms, and correlograms. These are inferred by weighting the samples depending on their distance to selected locations. The local distributions are locally transformed into Gaussian distributions embedding information on the local histogram. The distance weighted experimental variograms and correlograms are able to adapt to local changes in the direction and range of spatial continuity. The automatically fitted local variogram models and the local Gaussian transformation parameters are used in spatial estimation algorithms assuming local stationarity. The resulting maps are rich in nonstationary spatial features. The proposed process implies a higher computational effort than traditional stationary techniques, but if data availability allows for a reliable inference of the local distributions and statistics, a higher accuracy of estimates can be achieved.     

Bayesian Gaussian Mixture Linear Inversion for Geophysical Inverse Problems

A Bayesian linear inversion methodology based on Gaussian mixture models and its application to geophysical inverse problems are presented in this paper. The proposed inverse method is based on a Bayesian approach under the assumptions of a Gaussian mixture random field for the prior model and a Gaussian linear likelihood function. The model for the latent discrete variable is defined to be a stationary first-order Markov chain. In this approach, a recursive exact solution to an approximation of the posterior distribution of the inverse problem is proposed. A Markov chain Monte Carlo algorithm can be used to efficiently simulate realizations from the correct posterior model. Two inversion studies based on real well log data are presented, and the main results are the posterior distributions of the reservoir properties of interest, the corresponding predictions and prediction intervals, and a set of conditional realizations. The first application is a seismic inversion study for the prediction of lithological facies, P- and S-impedance, where an improvement of 30% in the root-mean-square error of the predictions compared to the traditional Gaussian inversion is obtained. The second application is a rock physics inversion study for the prediction of lithological facies, porosity, and clay volume, where predictions slightly improve compared to the Gaussian inversion approach.     

Sedimentological evaluation of general circulation model simulations for the “greenhouse” Earth: Cretaceous and Jurassic case studies

Conceptual climate models, based on the workings of the present-day climate system, provided a first-order approach to ancient climate systems. They are potentially very subjective in character. Their main drawback was that they involved the relocation of continents beneath a stable atmospheric circulation modelled upon that of the present. General circulation models (GCMs) use the laws of physics and an understanding of past geography to simulate climatic responses. They are objective in character. However, they require super computers to handle vast numbers of calculations. Nonetheless it is now possible to compare results from different GCMs for a range of times and over a wide range of parameterisations. GCMs are currently producing simulated climate predictions which compare favourably with the distributions of climatically sensitive facies (e.g. coals, evaporites and palaeosols). They have been used effectively in the prediction of oceanic upwelling sites and the distribution of petroleum source-rocks and phosphorites. Parameterisation is the main weakness in GCMs (e.g. sea-surface temperature, orography, cloud behaviour). Sensitivity experiments can be run on GCMs which simulate the effects of Milankovitch forcing and thus provide insights into possible patterns of climate change both globally and locally (i.e. provide predictions that can be evaluated against the rock record). Future use of GCMs could be in the forward modelling of sequence stratigraphic evolution and in the prediction of the diagenetic characteristics of reservoir units in frontier exploration areas. The sedimentary record provides the only way that GCMs may themselves be evaluated and this is important because these same GCMs are being used currently to predict possible changes in future climate.     

On the local multiscale determination of the Earth’s disturbing potential from discrete deflections of the vertical

As a first approximation, the Earth is a sphere; as a second approximation, it may be considered an ellipsoid of revolution. The deviations of the actual Earth’s gravity field from the ellipsoidal “normal” field are so small that they can be understood to be linear. The splitting of the Earth’s gravity field into a “normal” and a remaining small “disturbing” field considerably simplifies the problem of its determination. Under the assumption of an ellipsoidal Earth model, high observational accuracy is achievable only if the deviation (deflection of the vertical) of the physical plumb line, to which measurements refer, from the ellipsoidal normal is not ignored. Hence, the determination of the disturbing potential from known deflections of the vertical is a central problem of physical geodesy. In this paper, we propose a new, well-promising method for modelling the disturbing potential locally from the deflections of the vertical. Essential tools are integral formulae on the sphere based on Green’s function with respect to the Beltrami operator. The determination of the disturbing potential from deflections of the vertical is formulated as a multiscale procedure involving scale-dependent regularized versions of the surface gradient of the Green function. The modelling process is based on a multiscale framework by use of locally supported surface curl-free vector wavelets.      

重磁场的变差函数特征及应用

用变差函数研究重磁场的区域变化特征.变差函数的变程反映重磁场的相干范围, 块金效应反映随机干扰, 基台值反映变异程度.重磁场的理论模拟说明: 重力场的相干范围大于磁场, 重磁场变程主要取决于场源深度, 浅源重磁场变差函数近似为球状模型或指数模型, 深源重磁场近似为连续性更好的高斯模型.磁场场源深度近似等于变程的一半, 重力场场源深度近似等于变程的四分之一.湖北大冶铁矿垂直分量磁异常具有几何各向异性, 北西-南东走向, 变差函数推测磁铁矿平均深度为250m.磁异常小波多尺度分解细节和逼近部分磁场具有协调几何各向异性, 变差函数的各阶场源深度估计结果与功率谱估计结果吻合.      

A Spatial Correlation-Based Anomaly Detection Method for Subsurface Modeling

Spatial data analytics provides new opportunities for automated detection of anomalous data for data quality control and subsurface segmentation to reduce uncertainty in spatial models. Solely data-driven anomaly detection methods do not fully integrate spatial concepts such as spatial continuity and data sparsity. Also, data-driven anomaly detection methods are challenged in integrating critical geoscience and engineering expertise knowledge. The proposed spatial anomaly detection method is based on the semivariogram spatial continuity model derived from sparsely sampled well data and geological interpretations. The method calculates the lag joint cumulative probability for each matched pair of spatial data, given their lag vector and the semivariogram under the assumption of bivariate Gaussian distribution. For each combination of paired spatial data, the associated head and tail Gaussian standardized values of a pair of spatial data are mapped to the joint probability density function informed from the lag vector and semivariogram. The paired data are classified as anomalous if the associated head and tail Gaussian standardized values fall within a low probability zone. The anomaly decision threshold can be decided based on a loss function quantifying the cost of overestimation or underestimation. The proposed spatial correlation anomaly detection method is able to integrate domain expertise knowledge through trend and correlogram models with sparse spatial data to identify anomalous samples, region, segmentation boundaries, or facies transition zones. This is a useful automation tool for identifying samples in big spatial data on which to focus professional attention.