ECMOR XI - 11th European Conference on the Mathematics of Oil Recovery

Seismic pressure point sources have been simulated near the surface of sand dunes in Oman. Real topographic elevation data was assembled from recorded field data in the region, in-filled with necessary interpolated data, to form the surface topography implemented. Recent results describing seismic sand curves for the region are employed in the near-surface medium definition, and a constant half space of realistic seismic velocities and densities are used for the simulated sub-sand medium. Topographic effects in the form of scattering and reflections from prominent sand dune topography are confirmed, showing differences between scattering from elongated and normally oriented sand dunes to the excited source orientation. 8th order Finite-Differences (F-Ds) are used to discretize full elastic wave equations in the velocity-stress formulation, reducing gradually the F-D order when approaching all numerical grid boundaries. At the free surface, I employed boundary conditions for arbitrary surface topography derived for the particle velocities, discretized by 2nd order F-Ds, and at the side and bottom boundaries exponential damping is applied for non-reflective boundaries. A rectangular grid is used to model the Cartesian wave equations from a curved grid by adding in the equations extra terms that depend on the arbitrary surface topography and its spatial derivatives. The curved grid is adapted to an arbitrary surface topography with gradually less curvature with depth towards the bottom of the grid, which is planar. The code is parallellized by domain decomposition through MPI (Message Passing Interface) parallellization between processors. This makes possible the time domain modelling of higher frequencies (even relevant for exploration) and larger areas than could be done otherwise. Only in the last 3-5 years has such complex modelling been commercially viable in hydrocarbon exploration.

Anisotropic effective medium theory (for homogenization of fractured/composite porous media under the assumption of complete scale-separation) can be used to map the effects of sub-grid (or sub-seismic) fractures onto the grid-scale. The effective poroelastic stiffness and hydraulic permeability tensors that determines the seismic amplitude versus angle-azimuth (AVAZ) and production data (provided that one has information about the overburden and the porous matrix, as well as suitable tools for seismic forward modelling and fluid flow simulation) can therefore be viewed as functions of a relatively small number of parameters related with the fractures (e.g., fracture density and azimuthal orientation). In this paper, we develop a rock physics-based Bayesian method for joint inversion of seismic AVAZ and production data with respect to these parameters of the fractures. Within this stochastic framework, the expectation value of the effective permeability tensor is given by an integral of an effective medium approximation (derived by using rigorous integral equation or Green’s tensor function methods) weighted by the posterior probability density function over the allowed values of the parameters of the fractures. The present work complements those of Will et al. (2005), Jakobsen et al. (2007) and Jakobsen and Shahraini (2008) in the sense that we are using different inversion methods, seismic attributes and/or (rock physics) models of the fractured reservoir. Our interdisciplinary method for characterization of sub-grid fractures can in principle be applied to rather complicated models of fractured reservoirs (e.g., involving multiple sets of fractures characterized by different fracture densities and orientations, etc.). However, the initial (synthetic) modelling and inversion results presented here are associated with a relatively simple model of a fractured reservoir in which a single set of vertical fractures are located within sub-set of the (piecewise constant) reservoir model (that can in principle be identified using complementary methods). We have analysed synthetic AVAZ and production data contaminated with different amounts of normal distributed noise, and managed to recover the (unknown parameters of the fractures) effective permeability tensors with good accuracy. The results show clearly that the incorporation of seismic anisotropy attributes into the history matching process helps to reduce the uncertainties of the estimated permeability tensors.

In oil industry and subsurface hydrology, geostatistical models are often used to represent the spatial distribution of different lithofacies in the reservoir. Two main model families exist: multipoint and truncated Gaussian models. We focus here on the latter. In history matching of lithofacies reservoir model, we attempt to find multiple realizations of lithofacies configuration, that are conditional to dynamic data and representative of the model uncertainty space. This problem can be formalized in the Bayesian framework. Given a truncated Gaussian model as a prior and the dynamic data with its associated measurement error, we want to sample from the conditional distribution of the facies given the data. A relevant way to generate conditioned realizations is to use Markov Chains Monte Carlo (MCMC). However, the dimension of the model and the computational cost of each iteration are two important pitfalls for the use of MCMC. In practice, we have to stop the chain far before it has scanned the whole support of the posterior. Further more, as the relationship between the data and the random field is non-linear, the posterior can be multimodal. Hence, the chain may stay stuck in one of the modes. In this work, we first show how to reduce drastically the dimension of the problem by using a truncated Karhunen-Loève expansion of the Gaussian random field underlying the lithofacies realization. Then we show how we can improve the mixing properties of classical single MCMC, without increasing the global computational cost, by the use of parallel interacting Markov chains at different temperatures. Applying the dimension reduction and this innovative sampling method lowers drastically the number of iterations needed to sample efficiently from the posterior. We show the encouraging results obtained when applying the methodology to a synthetic history matching case.

Since the first version of Kalman Filter was introduced in 1960 it received a lot of attention in mathematical and engineering world. There are many successful successors like for example Ensemble Kalman Filter (Evensen 1996) which has been applied also for reservoir engineering problems. The method proposed in [Zhou et al. 2007] draws together the ensemble filtering ideas and an efficient covariance representation, and is expected to perform well in history matching for reservoir engineering. It is the Ensemble Multiscale Filter. The EnMSF is a different way to represent the covariance of an ensemble. The computations are done on a tree structure and are based on an ensemble of possible realizations of the states and/or parameters of interest. The ensemble consists of replicates that are the values of states per pixel. The pixels in the grid are partitioned between the nodes of the finest scale in the tree. A construction of the tree is led by the eigenvalue decomposition. Then, the state combinations with the greatest corresponding eigenvalues are kept on the higher scales. The updated states/parameters using the EnMSF are believed to keep geological structure due to localization property. It comes from the filter’s characterization where the pixels from the grid (e.g. permeability field) are distributed (in groups) over the finest scale tree nodes. We present a comparison of covariance matrices obtained with different setups used in the EnMSF. This sensitivity study is necessary since there are many parameters in the algorithm which can be adjusted to the needs of an application; they are connected to the tree construction part. The study gives the idea of how to efficiently use the EnMSF. The localization property is discussed based on the example where the filter is run with a simple simulator (2D, 2 phase) and a binary ensemble is used (the pixels in the replicates of permeability have two values only). Several possible patterns for ordering the pixels are applied.

The ensemble Kalman filter (EnKF) method has proven to be a promising tool for reservoir model updating and history matching. However, because the EnKF requires that the prior distributions for the parameters to be estimated are Gaussian, or approximately Gaussian, the application to facies models, and channel systems in particular, has been a challenge. In this paper we suggest two different approaches for parameterization of the facies models in terms of Gaussian perturbations of an existing "best guess" model. Method 1 is inspired by level set methods, where surfaces (here facies boundaries) are implicitly modelled through a level set function, normally defined as the signed distance from the nearest surface. Model realizations are generated by adding a Gaussian random field to the initial level set function. Method 2 is based on a similar idea, but the realizations are generated by adding a Gaussian random field to a smoothed indicator function. The smoothing is performed with a Shapiro filter, which is fast and simple to use for any number of dimensions. The performance of the methods is illustrated using a 3D model inspired by a real North Sea fluvial reservoir. It is shown that realistic facies model realizations may be generated from realizations of a Gaussian random field. Based on one realization from the prior for each of the methods, two sets of synthetic production data were generated. Then the prior model ensemble generated with each of the methods were conditioned to each of the two data sets using EnKF. Reasonably good matches were obtained in all cases, including those where the true model is not a realization from the same statistical model as is used to generate the prior realizations.

Oil and gas reservoirs are commonly complex heterogeneous natural structures that are described by means of several direct or indirect field measurements involving different physical processes operating at different spatial and temporal scales. Seismic techniques provide a description of the large scale geological structures. In some cases, they can help to characterize the spatial fluid distribution, which knowledges can in turn be used to improve the oil recovery strategy. In practice, these measurements are always expensive and due to their indirect, local and incomplete nature, an exhaustive representation of the entire reservoir cannot be attained. Several uncertainties are always remaining, that must be conveniently propagated in the modeling workflow to deduce in turn uncertainties of the production forecasts. Those uncertainties are essential when setting up a reservoir development scenario. A typical issue is to choose between several oil recovery scenario, or position and trajectory of a new well. Due to the cost of the associated field operations, it is essential to model the risks due to the remaining uncertainties. It is within this framework that devising strategies allowing to set-up optimal data acquisition schemes can have many applications in oil or gas reservoir engineering, or in the recently considered CO2 geological storages involving analogous technologies. We present a method allowing us to quantify the information that is potentially provided by any set of measurements. Applying a Bayesian framework, we quantify the information content of any set of data using the so called Kullback-Leibler divergence between posterior and prior distributions. In the case of a gaussian model where the data depend linearly on the parameters, analytic formulae are given and allow us to define the optimal set of time acquisitions. The redundancy of information can also be quantified, highlighting the role of the correlation structure of the prior model.

The goal of history matching is to construct one or more reservoir models that match observed field behaviour and use those models to forecast field behaviour under different operating conditions. In constructing the model to be matched, there are a number of choices to be made, such as the type of geological model to construct and the number of unknown parameters to use for history matching. We often choose a single geological model and vary a set of parameters to give the best fit to the data (e.g. production history). There are two areas of concern with this approach. The first is that there are usually a number of possible geological models which are all plausible to some degree, and we may be able to make better forecasts by using this information. The second is that increasing the number of unknown parameters may give a better fit to the data, but may not predict so well if the model is over-fitted. The goal of this paper is to examine the application of two techniques to handle these problems. The first technique uses the concept of minimum description length (MDL), which was developed from information theory, and quantifies the trade-off between model complexity and goodness of fit. The second technique is called Bayesian Model Averaging, and was developed to improve the reliability of weather forecasts using codes developed at different centres by constructing a suitable average of the models which takes into account not only the uncertainty forecast by each model but also the between model variance. Both techniques are illustrated with real reservoir examples. The MDL approach is shown on a simple reservoir model based on a Gulf of Mexico Field. The BMA approach is shown on a field with a moderate number of injectors and producers.

It is well-known that history matching of reservoir models with production measurements is an ill-posed problem, e.g. different choices for the history matching parameters may lead to equally good history matches. We analyzed this problem using the system-theoretical concept of structural identifiability. This allows us to analytically calculate a so-called information matrix. From the information matrix we can determine an identifiable parameterization with a significantly reduced number of parameters. We apply structural identifiability analysis to single-phase reservoir simulation models and obtain identifiable parameterizations. Next, we use the parameterization in minimizing an objective function that is defined as the mismatch between pressure measurements and model outputs. We also apply the structural identifiability analysis to an object-based parameterization describing channels and barriers in the reservoir. We use the iterative procedure to determine for reservoir models with 2025 grid block permeability values a structurally identifiable parameterization of only 13 identifiable parameters. Next, we demonstrate that the parameterization leads to perfect history matches without the use of a prior model in the objective function. We also demonstrate the use of the identifiable object-based parameterization, leading to geologically more realistic history matches.

We consider an efficient AVA stochastic inversion algorithm which performs inversion of equal angle offset volumes to volumes of Vp, Vs and density. Exact solution of Zoeppritz equations is used in order to determine reflection coefficients. Using combination of Simulated Annealing and Metropolis-Hasting sampler, we generate a set of equiprobable high-resolution volumes of elastic properties which are consistent with seismic data, well measurements and reproduce detailed geological features. Stochastic realizations generating process is constrained to reproduce seismic data, well data, low-frequency elastic trends, 3D variograms for elastic parameters, estimates of probability distributions of (Vp,Vs,Den), estimates of vertical alteration of elastic properties, etc. The multi-trace approach and extensive use of a priori geological knowledge greatly reduces influence of seismic noise on the inversion results for each of the solutions obtained. Multiple grid approach is used to properly treat large-scale variograms features. For reproduction of vertical alteration of elastic properties we use multi-point statistics. The training dataset is formed from log-derived properties. This eliminates the need for building an artificial training model. Application of multi-point statistics has a greater potential for ensuring more accurate account of vertical variations of acoustic properties, for obtaining geologically valid inversion solutions and for decreasing the uncertainty degree. The seismic inversion algorithm has been efficiently parallelized. It uses a shared memory computer and inverts simultaneously different 1D models on different processors preserving spatial correlation of elastic properties. AVA stochastic inversion results are used for cascaded stochastic simulation of lithology and fluid units, for simulation of reservoir properties, uncertainty analysis. Based on Bayesian approach we can take into account of existing uncertainties, in particular, uncertainties caused by inaccuracy of the mathematical model that relates elastic parameters and reservoir properties, and uncertainties caused by non-uniqueness and inaccuracy of seismic inversion results.

The experimental design methodology is widely used to quantify uncertainty in the oil and gas industry. This technique is adapted for uncertainty quantification on non-spatial parameters, such as bubble point pressure, oil viscosity, and aquifer strength. However, it is not well adapted for the case of geostatistical (spatial) uncertainty, due to the discrete nature of many input parameters as well as the potential nonlinear response with respect to those parameters. One way to handle spatial uncertainty is called the joint modeling method (JMM). This method, originally proposed in a petroleum context by Zabalza (2000), incorporates both non-spatial and spatial parameters within an experimental design framework. The method consists of the construction of two models, a mean model which accounts for the non-spatial parameters and a dispersion model which accounts for the spatial uncertainty. Classical Monte-Carlo simulation is then applied to obtain the probability density and quantiles of the response of interest (for example the cumulative oil production). Another method to quantify spatial uncertainty is the distance kernel method (DKM) proposed recently by Scheidt and Caers (2007), which defines a realization-based model of uncertainty. Based on a distance measure between realizations, the methodology uses kernel methods to select a small subset of representative realizations which have the same characteristics as the entire set. Flow simulations are then run on the subset, allowing for an efficient and accurate quantification of uncertainty. In this work, we extend the DKM to address uncertainty in both spatial and non-spatial parameters, and propose it as an alternative to the joint JMM. Both methods are applied to a synthetic test case which has spatial uncertainty on the channel representation of the facies, and non-spatial uncertainties on the channel permeability, porosity, and connate water saturation. The results show that the DKM provides a more accurate quantification of uncertainty with fewer reservoir simulations. Finally, we propose a third method which combines aspects of the DKM and the JMM. This third method again shows improvement in efficiency compared to the JMM alone.