Elements for an Integrated Geostatistical Modeling of Heterogeneous Reservoirs

Résumé — Éléments pour une modélisation géostatistique intégrée de réservoirs hétérogènes — L’approche géostatistique pour la modélisation des réservoirs hétérogènes permet, d’une part, d’intégrer l’ensemble des données de natures diverses et d’échelles différentes, et d’autre part, d’évaluer les incertitudes par génération de multiples scénarios possibles. La conception d’un modèle géostatistique de réservoir doit d’abord dépendre de l’environnement de dépôt géologique, ce qui permet de représenter les hétérogénéités majeures qui contrôlent l’écoulement des fluides. Ensuite, les données quantitatives provenant des carottes, des diagraphies, de la sismique, des tests de puits, etc., doivent être utilisées pour l’inférence des paramètres structuraux du modèle. Enfin, le calage des réalisations du modèle aux données hydrodynamiques permet d’augmenter davantage la fiabilité du modèle dans la prévision en production. Cet article présente les éléments d’une méthodologie intégrée pour la modélisation des réservoirs hétérogènes. Nous introduisons d’abord les modèles géostatistiques de base couramment utilisés pour décrire les réservoirs hétérogènes, suivis par un bref aperçu sur l’inférence des paramètres structuraux des modèles. Nous


INTRODUCTION
A reservoir is a porous and permeable subsurface formation, containing one or several fluids (water, gas and oil) under pressure.No reservoir presents the same petrophysical properties everywhere.When a reservoir containing hydrocarbons is discovered, the knowledge on the spatial distribution of petrophysical properties (porosity and permeability) of the reservoir is of primary importance for optimizing its development project.
First of all, we need to identify the depositional environment of the reservoir.One distinguishes two main types of reservoirs: the clastic ones and the carbonate ones.The clastic reservoirs are made primarily of hard and chemically stable grains of quartz.They come from the destruction of eruptive rocks such as granite, according to the traditional cycle of erosion-transport-deposition.The depositional environment of a clastic reservoir can also vary: from continental aeolian sandstones to turbidites, while passing by the fluviatile sandstones (braided channels or bars of meander), the offshore bars, the beaches and the bars of mouth.The carbonates are limestones and dolomites.They are much less stable than quartz and very sensitive to diagenesis: they can be dissolved and reprecipitated, dolomitized, fractured under tectonic constraints.One knows, for example, the limestones built by corals, in marine environment.Each type of depositional environment presents its own geometrical architecture.Outcrop studies help us to distinguish between the various reservoir architectures and to define geological rules for building reservoir models.
In complement with geological conceptual knowledge, there are many tools of investigation to quantify, directly or indirectly, the petrophysical properties of the reservoir.For instance, we have well logs, rock samples, seismic reflection, well seismic, repeated (4D) seismic, well tests, production history, etc.The investigation domain (the scale), the precision and the nature of the information obtained from these tools are very different.For example, the well logs provide lithological or petrophysical information at a centimetric resolution along the wells, while the seismic reflection provides a very dense cover of the horizontal variation of the acoustic properties of the reservoir with a decametric vertical resolution.The interpretation of a well test can give an idea on the average permeability in the vicinity of a well, whereas production data (pressure, watercut, etc.) are rather related to the large-scale (field-scale) connectivity of the reservoir.There is often complementarity between data obtained by different tools of investigation.It is therefore important to integrate all these data in the reservoir model.
The objective of reservoir modeling is to build digital representations (models) of the reservoir in conformity with the geological rules and constrained to all available quantitative information.Also, the investigation on a reservoir is carried out throughout the life of the reservoir from its discovery to the end of its production.We must be able to progressively update a reservoir model with the acquisition of new data.Updating a reservoir model by the integration of more and more available data should make it closer and closer to the reality.However, reality is never completely accessible and there is always a degree of uncertainty, whatever the development stage of a reservoir field.It is thus necessary to build multiple representative models that cover all possible scenarios of the reservoir.In general, the reservoir modeling must: -take into account the geological environment of deposition; -integrate in a geologically consistent way the entire quantitative information; -allow the progressive updating to account for new data; -represent all possible scenarios of the reservoir.
In such context, the probabilistic approach appears to be the most appropriate choice for reservoir modeling.There are already a large number of geostatistical models and methods that provide the potentiality to represent the various depositional environments and to integrate the various types of data.In addition, the nature of the probabilistic approach makes it possible to build multiple representations of the reservoir for the evaluation of uncertainties on production forecasts.
This paper presents the elements of an integrated methodology for modeling heterogeneous reservoirs.We first introduce the basic geostatistical models used to describe the heterogeneous reservoirs.This is followed by an outline on the inference of the model structural parameters.Then, we present an inverse approach based on the gradual deformation method for history matching.

GEOLOGICAL MODELING
The diversity of the geological environments justifies the utility of various types of stochastic models.According to their methods of construction, we can distinguish three types of models: -pixel-based models (e.g., models generated by multi-Gaussian simulations); -object-based models (e.g., the Boolean model); -process-based (or genetic) models (model based on the sedimentary processes).
The geological modeling of a heterogeneous reservoir is often decomposed into several steps and implies several basic models.

Continuous Gaussian Simulations
Let Y(x) (x ∈ D) be a stationary multi-Gaussian random function of order 2. Given that Y(x) has zero mean and unit variance, we define a random function Z(x) (x ∈ D) by the following nondecreasing transformation of Y(x): (1) where G is the cumulative distribution function of the standard Gaussian variable, and F the cumulative distribution function of Z(x).One often utilizes Z(x) for modeling physical properties in field D. In general, the nondecreasing function F -1 {G[•]} is expanded in terms of Hermitian polynomials and the expansion coefficients are determined from the empirical cumulative distribution function of Z(x) (e.g., Lantuéjoul and Rivoirard, 1984).
For some commonly used random function models, the transformation function F -1 {G[•]} can have a simple form.If for instance, Z(x) follows a lognormal distribution, we have then: (2) where m and σ are respectively the mean and the standard deviation of ln[Z(x)].Another example is the truncated Gaussian random function (Matheron et al., 1987) that will be presented more in detail in Section 1.2.
The covariance of Y(x) is determined by that of Z(x), and the simulation of Z(x) turns to be that of Y(x).Several methods are available for simulating Gaussian random functions on a regular or irregular grid.
The matrix decomposition method is based on the Cholesky decomposition of the covariance matrix of the multi-Gaussian vector.This method is rigorous but limited to small simulation grid (the number of nodes of the grid of simulation should not exceed a thousand to generate a realization in a reasonable time).Conditioning to well data can be performed simultaneously with the generation of realizations.
Rather than a simulation method, the turning band method (Matheron, 1973) is a stereological device designed to reduce a multidimensional simulation into 1D ones.1D simulations are generated on lines uniformly spaced in 2-3D.These 1D simulations are then projected onto the grid nodes and averaged to generate higher dimensional realizations.An attractive feature of the turning band method is that 1D simulations can be projected onto all kinds of (regular or irregular) grids (Blanc et al., 1998).Thus, it is possible to modify the grid system without changing the realization of the random function model.
Spectral methods (e.g., Chilès and Delfiner, 1999) are based on the spectral decomposition of the covariance function.The continuous spectral method does not require the definition of a grid system before generating realizations, while discrete spectral method is operational only on regular grid system.Like the turning band method, spectral methods provide nonconditional simulations.
The FFT moving average method (Le Ravalec et al., 2000) combines the moving average method with the fast Fourier transformation.It consists in performing the convolution product of the moving average method in Fourier space, and this makes the moving averaging computationally easy and fast.This method can generate efficiently large Gaussian realizations of any stationary covariance function.
The sequential Gaussian simulation method (e.g., Deutsch and Journel, 1992) is based on the fact that the conditional law of a component of a Gaussian vector knowing the other components remains Gaussian.It is thus possible to generate a realization of a Gaussian vector component by component one after the other.Conditioning to well data is carried out simply by taking these data as the first generated values.Moreover, although the sequential method is particularly suited for simulating Gaussian vectors, it does not require, in theory, the multi-Gaussien property.For example, we can use the sequential algorithm to generate a lithological distribution directly.We will reconsider this point in Section 1.3.
Strictly speaking, the probability field method (Froidevaux, 1993) is a device for conditioning a uniform simulation to well data.In practice, a uniform simulation is often derived from a Gaussian simulation that can be generated by any Gaussian simulation algorithm.Some simulation algorithms (e.g., matrix decomposition, and sequential simulation) generate directly conditional realizations, contrarily to some other algorithms (e.g., FFT moving average method, turning bands).In the later case, a conditioning step is necessary.This can be realized using conditioning kriging (Journel and Huijbregts, 1978).The principle is as follow.
Let y nc (x) be a nonconditional realization of a multi-Gaussian random function Y(x), y * (x) the kriging of Y(x) by using the real data set (y(x 1 ), y(x 2 ), ..., y(x N )), and y * nc (x) the kriging of Y(x) with the simulated data set (y nc (x 1 ), y nc (x 2 ), ..., y nc (x N )).Then y c (x), defined by: (3) is a conditional realization of Y(x).

Truncated Gaussian Simulations
The truncated Gaussian simulation method (Matheron et al., 1987) is often used for describing the lithofacies distribution in a reservoir field (Ravenne et al., 1987).Let Y(x) (x ∈ D) be a standard Gaussian random function defined over the field D, and I(x) (x ∈ D) an indicator random function defined by truncating Y(x) at threshold s: If Y(x) is stationary, then I(x) is also stationary.It can be used to simulate, for instance, the geometry of the spatial distribution of a lithofacies in an oil reservoir.This method can be easily extended to the case of several thresholds in order to build models with several lithofacies (Matheron et al., 1987).In the common practice of this method, each lithofacies corresponds to an interval of the Gaussian variable.The transitions are possible only between the lithofacies whose intervals are contiguous.However, to generate a lithological model with direct transition between any two lithofacies, one can use disjoint intervals to define the lithofacies.
Using regionalized thresholds (Beucher et al., 1993) makes it possible to generate lithofacies models with regional trend.Regionalized threshold functions are derived from regionalized proportion functions of lithofacies.The proportion of a lithofacies is defined over a given area (support) that can vary from the well data support to the whole reservoir field.Proportion functions are built by using well data and seismic derived information (Moulière et al., 1997).
In the case where well test pressure data can be interpreted in terms of permeability data, these data can also be used to estimate lithofacies proportion functions (Hu et al., 1998).The regionalized thresholds can also be a realization of another Gaussian random function (Adler and Thovert, 1998).Figure 1 shows a nonstationary truncated Gaussian simulation of the lithofacies distribution in a turbidite reservoir (offshore Brazil).
Another possibility of enriching the method of truncated Gaussian simulation consists in simultaneously truncating several Gaussian random functions in the definition of the lithofacies, which makes it possible to build extremely varied lithological models (Galli et al., 1994;Le Loc'h and Figure 1 Truncated Gaussian simulation.Vertical cross section of a reservoir model with four lithofacies represented by different colors, offshore Brazil (after Souza, 1997; Reis, 2001).Galli, 1997).Unlike the truncated mono-Gaussian case, truncated pluri-Gaussian simulations do not necessarily impose a sequence of transitions between different lithofacies.Truncating a Gaussian simulation with another Gaussian simulation, mentioned in the previous paragraph, is a particular example of pluri-Gaussian simulations.
In general, lithofacies indicators along wells can be directly accounted for in truncated Gaussian simulations.This is achieved first by generating truncated Gaussian values along wells using the Gibbs sampler (Geman and Geman, 1984;Freulon and de Fouquet, 1993) and then by applying conditioning on the whole Gaussian field.In the particular case where the proportion support is reduced to the well data support, the conditioning to lithofacies indicator data along wells is directly satisfied without the above step.Note that in this particular case, the truncated Gaussian method becomes equivalent to the probability field method (Froidevaux, 1993) for simulating lithofacies.This is because the local distributions involved in the probability field method are equivalent to the local lithofacies proportions, and in practice the probability field used for sampling local distributions is obtained by transforming a Gaussian field into a uniform field.

General Sequential Simulations
Let Z = (Z 1 , Z 2 , ..., Z N ) be the N-dimensional random vector of interest.Z is neither necessarily multi-Gaussian nor necessarily associated with a stationary random function.The sequential simulation approach reduces the problem of generating an N-dimensional random vector into a series of N univariate generation problems.The sequential simulation of Z involves first the definition of an order according to which the N elements of the vector Z are generated one after another.Without loss of generality, we assume that the N components Z 1 , Z 2 , ..., Z N are already arranged according to the defined order and realizations of Z will be generated sequentially from Z 1 to Z N .For each element Z i = (i =1, 2, ..., N), generating its realization requires the following operations: -build the distribution function of Z i conditioned to (Z 1 , Z 2 , ..., Z i-1 ): (5) -draw a value for Z i from the distribution F c (z i ).
In geostatistical practice, the sequential simulation is frequently used for generating realizations of multi-Gaussian vectors and non-Gaussian indicator vectors.In the later case, the conditional probability of the presence of a lithofacies is estimated roughly by indicator kriging (Deutsch and Journel, 1992).
In general, the major difficulty of a sequential simulation lies in the determination of the conditional distributions For the heterogeneous reserves of complex structures, the approximation of these conditional distributions by kriging based on the variogram is far from being satisfactory.The use of the multi-point statistics (Guardiano and Srivastava, 1993;Strebelle and Journel, 2000) makes it possible to capture more complex geological structures compared to two point statistics (i.e., the variogram).However, the inference of the multi-point statistics is much more difficult and requires resorting to training images.These images can be obtained from the description of analogous outcrops or conceived by an objectbased approach or a genetic (process-based) approach.From this point of view, the object-based approach, the processbased approach and the approach of multi-point statistics are complementary.We can regard the approach of multi-point statistics as a way of conditioning the genetic models or the object-based models.

Boolean Simulations
The Boolean model (Matheron, 1967;Stoyan et al., 1995) belongs to the class of the object-based models.An objectbased model is an arrangement of a population of geometrical objects in space.In reservoir engineering, they are often used to describe meandering systems, fracture networks and porous media at the granulometric scale, etc.
In the stationary Boolean model, the objects are distributed in space according to a Poisson point process of constant density.The forms of the objects and their sizes are independent of their locations.Let A a random elementary object (also called primary grain).The construction of a Boolean model consists in generating independent elementary objects A i (i = 1, 2, ...) having the same statistical characteristics as the object A and then implementing them at the points of a Poisson point pattern.The union of these objects denoted X constitutes a random set called Boolean model: The probability for an objet B being in the complementary of X is given by: (7) where |A ⊕ B ∨ | stands for the volume of A dilated by B (Matheron, 1967) and θ the Poisson density.The Boolean model is completely defined by the Formula (7) called its distribution function.In particular, if object B is reduced to a point, one obtains the "porosity" of the Boolean model: That is the proportion of volume occupied by X -.Using the algorithm based on the Markov iteration (Lantuéjoul, 1997), one can simulate a Boolean model in the field D under the condition that two subsets C 1 and C 0 of D belongs respectively to the union of objects X and to its complement X -.The above basic Boolean model can be generalized by combining the objects of different nature or/and by using a regionalized Poisson density.

Random Genetic Simulations
Genetic modeling is based on the sedimentary process.The parameters of the model of sedimentation are often uncertain and can be represented by random variables or random functions.This way of modeling the geological formations is called the random genetic approach.It has the advantage of being able to better describe the complex geological formations than an approach only based on random functions or random sets.Jacod and Joathon (1971) are among the pioneers of using a random genetic approach to describe the geometry of a sedimentary formation.Since then, a lot of work was devoted to this approach for the modeling of the geological phenomena.Hu et al. (1994) developed a methodology to simulate internal geometrical architecture of deltaic depositional system.Figure 2 shows an example of the random genetic simulation of the Roda sand body X (Spain).More recently, Lopez, Galli and Cojan (2001) proposed a method to build a meandering channel system.The major problem of the genetic approach is the respect of the data observed at wells.That limited the application of this type of method to real cases.A solution to this problem probably lies in the approach of the multi-point statistics.We can use nonconditional random genetic simulations as training images for the inference of the conditional probabilities, then to generate conditional realizations by sequential simulation.
Figure 2 Random genetic simulation.Vertical cross section of the internal geometry of the Roda sand body X, Spain (after Hu et al., 1994).

Model Conception
Until here, we presented some basic stochastic models and the algorithms of their simulation.A geological formation is in general a very complex medium and can be difficultly modeled by only one basic model.The geological modeling of a heterogeneous reservoir often composes several major steps and implies several basic models.We first define geological units between which the depositional environment changes.The change of the depositional environment is often associated with large variations of petrophysical properties.Then, in each geological unit, a lithological modeling is performed.This step implies the categorical models such as the truncated Gaussian model, the model generated by sequential indicator simulation, the Boolean model or the genetic models.The choice of a basic model depends on the geological depositional environment.Finally, inside each lithofacies, one builds a petrophysical model by using continuous simulations.The main reason for separating the geological unit modeling in two steps is that it is in general easier to introduce geological knowledge into a lithological model than in a petrophysical model.
Note that the choice and the reliability of a model also depend on the stage in the life of a reservoir.For example, the genetic models or the object-based models are particularly suited for the description of a reservoir at the stage of its appreciation.At this stage, we have few wells and we need to resort more to the geological concepts to build a realistic model.While the pixel-based models make it possible to model spatial variability on smaller scale when we have much more quantitative information at the stage of its development.
Lastly, a reservoir model must be sufficiently flexible to be able to be progressively updated during the acquisition of new data throughout the life of a reservoir.This imposes the choice of the simulation algorithms allowing the local deformations.

PARAMETER INFERENCE
Before generating realizations of a stochastic model, it is necessary to infer the structural parameters of the model.The means, the variances and the variograms are among the structural parameters usually used in practice.They constitute the two first moments of a random function model.In particular, a multi-Gaussian random function is entirely defined by its two first moments.The choice and the fitting of a variogram (structural analysis) were treated in several reference books of geostatistics (Journel and Huijbregts, 1978;Chilès and Delfiner, 1999).There are also literatures devoted to the estimation of the Poisson density of a Boolean model (Schmitt, 1991;Lantuéjoul and Schmitt, 1991;Schmitt and Beucher, 1997;Benito, 2003).The use of higher order moments remains still uncommon in practice because of the difficulty of their inference, except for work by Guardiano and Srivastava (1993) and by Strebelle and Journel (2000) who use multi-point statistics for modeling complexe geological architectures (see also Section 1.3).The inference of the structural parameters is a crucial step of reservoir modeling because a geologically realistic model with false structural parameters will not be valid for production forecasts.It is thus necessary, as of this step, to integrate a maximum of available information.In this section, we consider the estimation of lithofacies proportions to illustrate the use of well, seismic and well-test information for parameter inference.

Estimation of Lithofacies Proportions
The lithofacies proportions constitute the moments of order 1 of a lithological model.In truncated Gaussian simulations, the lithofacies proportions determine the thresholds of truncation of the Gaussian random function.In sequential indicator simulations, one estimates the probabilities of presence of the lithofacies at each node of the simulation grid.These probabilities of presence can be interpreted as the lithofacies proportions of a set of realizations on the point support.In Boolean simulations, the lithofacies proportions are related to the form and the density of objects.When the object shape and the law of the object size are given, the object density is entirely determined by the lithological proportions.The estimation of the proportions (or probabilities of presence) of the lithofacies is then crucial for the construction of a lithological model.The capacity of prediction of a lithofacies model strongly depends on the accuracy of the estimated proportions (or probabilities of presence) of the lithofacies.
In general, the lithological distribution of an oil reservoir varies in a nonstationary way.For estimating proportion functions, we have lithological data at wells.Description in lithofacies along the wells is discretized according to lithological variability and the resolution of logs (from a few centimeters to a few decimeters in general).At any point x α of the discretization, one defines the indicator function of lithofacies n: The covariance C n (h) of each indicator function is inferred from well data or other sources of information (e.g., analogous outcrop data).In the case of few wells, it is necessary to resort to other sources of information (seismic, well tests, etc.) for the estimation of lithofacies proportions.

Contribution of Seismic Information
Seismic information concerns the acoustic properties of the subsurface formation averaged over a support much larger than that of well data.These acoustic properties are in general correlated with the lithology of the formation.Contrary to well data, seismic data are abundant in terms of the number of seismic traces and thus provide a dense coverage in the horizontal plane.However, the vertical resolution of seismic data is still very limited compared to the thickness of a reservoir.One can only expect to extract, from seismic data, a horizontal trend of the lithological variation of the formation on a thickness much larger than the step of discretization along the wells.
The seismic data can be a set of seismic attributes derived from seismic traces.The seismic attributes, in their initial forms (amplitude, coefficient of reflection, acoustic impedance, etc.), can be difficultly used directly to estimate the proportions of lithofacies because of their different nature from the lithological proportions and their often weak individual correlation with lithological properties.It is thus advisable to first extract, from several seismic attributes, synthetic lithological or petrophysical information (probability of presence of a lithofacies at a seismic pixel, proportion of a lithofacies on a given thickness, average porosity on a given thickness, etc.).Then, this synthetic information is used as constraints or auxiliary data for estimating lithifacies proportions (well data being used as principal information).
The extraction of synthetic geological information (e.g., lithological proportions, petrophysical properties) averaged over a vertical interval is the objective of the quantitative statistical calibration of seismic data (Fournier and Derain, 1995).This calibration consists in establishing a statistical relation between the seismic attributes (s 1 , s 2 , ..., s p ) in the vicinities of the wells and the geological properties g 1 , g 2 , ..., g q known at wells and averaged over the vertical interval defined according to the seismic resolution.
(10) This relation will be used to estimate the geological properties (their estimated values with their confidence intervals) from each seismic trace (for which geological information is not known).The Relation (10) can integrate simple relations between the seismic and geological parameters and their spatial relations (via the cokriging) as well.The reliability of the Relation (10) depends on the number of pairs wells-traces used in calibration.
Most statistical methods of calibration assume that the points on which calibration is based come from a homogeneous population.This assumption must be validated before applying the method.A way of carrying out this validation consists in first applying the qualitative calibration of the seismic data.This is to study the correspondences between seismic facies, defined by the morphology of the traces in the level of the studied reservoir, and the geological lithofacies.That makes it possible to divide the reservoir field into several geologically different subfields.Then, one applies quantitative calibration separately in each subfield.
Whatever the calibration method used, one obtains only an average geological information on the interval defined according to the vertical resolution of the seismic data.However, this information is very dense laterally and should contribute to the estimation of the geological properties between often very distant wells.Several forms of kriging for estimating the lithological proportions under seismic constraint are proposed in Moulière et al., (1997).

Contribution of Well Test Information
Well-tests are commonly used to investigate reservoir dynamic behaviors around wells.A well test consists in imposing some flow rate impulse in the reservoir (e.g.start production, change production rate, stop production, etc.) and measuring the pressure responses.These pressure data inform indirectly average flow behavior around wells and can thus contribute to the estimation of the lithological proportions around the wells.The interpretation of the pressure response as a function of time and its derivative often makes it possible to determine an apparent permeability in the investigation area of a well test.For each well test, one can derive an apparent permeability K wt in a known investigation area V.A power-averaging formula can be adopted to relate K wt to the absolute permeabilities of the lithofacies.(11) where K n stands for the absolute permeability of lithofacies n and P n (V) its proportion in V.The exponent ω is to be calibrated for each well test (Alabert, 1989;Noetinger, 1993).
Like seismic data, well test data are average information on support much larger than that of lithological data at wells.If we adopt the Relation (11), the problem is then to estimate the lithofacies proportions under constraints on their (weighted) average.The iterative method under aggregation constraints was proposed to solve this problem (Hu et al., 1998).This method consists in minimizing the sum of the estimation variances of the additive quantities under constraint acting on their linear combination.Haas and Noetinger (1997) propose another method based on cokriging.

HISTORY MATCHING
Once a stochastic model is defined and the structural parameters determined, one could generate multiple realizations.These realizations need to be calibrated to hydrodynamic data, because a realization, inconsistent with the production history, will certainly not credible for production The calibration of a stochastic model to the hydrodynamic data can be formulated as an optimization problem.First, one defines an objective function that measures the difference between the data observed in the reservoir field and the corresponding responses calculated on a realization of the stochastic model.Then, one modifies this initial realization until the objective function becomes small enough.The efficiency of a calibration procedure depends greatly on the way the realizations of the stochastic model are deformed.Also, experience shows that a calibration to the production, without accounting for geological knowledge, does not lead to a predictive model.Thus, a calibration method should have the following characteristics: -Preservation of the spatial variability/continuity of the stochastic model.-Reduction of the number of parameters to optimize.
-Regularity of the objective function with respect to the parameters to optimize.-Coverage of the whole solution space.

Several calibration methods have been proposed in literature.
Swapping the values at two arbitrary points of a realization or deforming a realization by genetic operations (e.g., Sen et al., 1992) violates the spatial continuity of the stochastic model.The use of these algorithms requires including a covariance term in the objective function, resulting in a very tedious optimization.Instead, changing the value at an arbitrary point of a realization, by drawing from its conditional distribution or by independent drawing followed with a convolution, preserves the covariance of the stochastic model.These algorithms are used in a Markov chain Monte-Carlo approach for conditional simulation (e.g., Hegstad et al., 1994;Oliver et al., 1997).However the objective function is either insensitive or discontinuous to this local discontinuous deformation.Consequently, although the covariance is preserved, optimizing such an objective function can still be extremely slow in the context of constraining large models to dynamic data.The pilot point method proposed by de Marsily et al. (1984) consists in first modifying the values at certain pilot points of a realization by accounting for the sensitivities of the objective function with respect to the values at these points.The modifications at the pilot points are then propagated to the whole realization by using conditioning kriging.The efficiency of this method comes from the reduction of calibration parameters and the use of a gradient based optimization algorithm.Nevertheless, deterministic modifications of the values at the pilot points, dictated by an optimization algorithm, may deteriorate the spatial structure of a realization, although in practice this can be limited by setting a minimal distance between the pilot points and a confidence interval of the values at these points (RamaRao et al., 1995).Note also the method proposed by Srivastava (1994) for "visualizing uncertainty" of a stochastic model.This method consists in first generating a realization in a field larger than the reservoir field.Then by moving the reservoir field in the larger field, one obtains a series of correlated realizations of the stochastic model.If we use this method in an optimization procedure for calibrating a reservoir model to nonlinear data, it is in general necessary to generate a complete realization in a field much larger than the reservoir field so as to obtain a satisfactory solution realization.This limits the applicability of the method when dealing with large reservoir models.
The gradual deformation method makes it possible to modify progressively realizations of multi-Gaussian stochastic models while preserving their space variability/continuity (Hu, 2000).This method consists in building a stochastic process whose state space is the ensemble of realizations of a stochastic model.In particular, we propose to build realization chains by combining independent Gaussian realizations.This method is then coupled with an optimization algorithm for calibrating realizations of a stochastic model to hydrodynamic data.In what follows, we present the gradual deformation method in the multi-Gaussian framework and some of its extensions: structural gradual deformation, local or regionalized gradual deformation.We also introduce the generalization of this method to any type of stochastic model.

Gradual Deformation Method
Let Y(t) be a stochastic process whose state space is the ensemble of realizations of the spatial random function Y.For each time t, Y(t) is a Gaussian random function defined in field D. In general, the axis time of the spatio-temporal random function Y(t) does not have any physical signification.Yet, its introduction in the stochastic model allows us to create dependent transition from one realization to another.
There are many ways of building a spatio-temporal random function Y(t).A convenient way of building a spatiotemporal random function Y(t) consists in combining two independent standard Gaussian random functions Y 1 and Y 2 with identical covariance: It is straightforward to show that, for each t, Y(t) has zero mean because that Y 1 and Y 2 have zero mean and that Y(t) shares the covariance of Y 1 and Y 2 due to their independence.cost t and sin t are simply variance normalization coefficients.Furthermore, for each t, Y(t) is a Gaussian random function due to the fact that it is a linear combination of Gaussian random functions.Given two independent realizations y 1 and y 2 of Y 1 and Y 2 , we get a continuous chain of realizations y(t).Hence by applying the Transformation (1), we obtain a continuous chain of realizations z(t) of the random function Z.The basic idea of the gradual deformation method is to calibrate z(t) to hydrodynamic data by fitting parameter t.
Figure 3 shows a chain of realizations of a Gaussian random function with a Gaussian variogram.
Although, throughout this paper, we assume that the random function Y is stationary of order 2, the deformation algorithm summarized by Equation ( 12) remains valid when Y is only intrinsically stationary (Matheron, 1971).As a matter of fact, if Y 1 and Y 2 are two independent intrinsic random functions with identical variogram, Y(t) defined by Equation ( 12) is an intrinsic random function for each t, and it shares the variogram of Y 1 and Y 2 .However, the extension to non-Gaussian-related models is not straightforward, simply because, in general, linear combination of independent non-Gaussian random functions will not preserve the non-Gaussian distribution, although the variogram is preserved.

Linear data like point values or arithmetic mean values of
Y can be directly accounted for in the spatio-temporal model Y(t) by using conditioning kriging.We insist here that the conditioning kriging is always applied to the stochastic process Y(t) but never to the independent random function Y 1 and Y 2 .In the following, we keep in mind that the linear data are directly integrated in the model using conditioning kriging, and we focus only on the problem of calibration to nonlinear hydrodynamic data.

Iterative Calibration Procedure
Let f obs = (f 1 obs , f 2 obs , ..., f p obs ) be the vector of the nonlinear data observed in the physical field and f = (f 1, f 2, ..., f p ) the corresponding vector of the responses of the stochastic model Z.For a given realization z of Z, the values of f 1, f 2, ..., f p are often obtained through numerical simulations.For instance, if z represents a permeability field, f 1, f 2, ..., f p may represent n pressure data which are generally obtained using numerical fluid flow simulation.The problem of constraining the stochastic model Z by the observations consists in generating a realization of Z which reduces an objective function to a low-enough level.We define the objective function as the sum of the weighted quadratic errors of the model responses with respect to the observations on the physical field: (13) where ω i denotes the weight attributed to response f i .The f i (i = 1, 2, ..., p) are functions of the vector z, so is the objective function J.We have then a multivariate optimization problem.
Denote by N the number of components of the vector z.N is often a huge number (10 4 ∼ 10 7 ).Therefore, it is highly tedious to optimize the objective function directly with respect to the vector z.Moreover, in the frame of stochastic optimization, the vector z cannot be arbitrary and it must be a realization of the stochastic model Z.In order of overcome these difficulties, it is convenient to transform Z to a Gaussian random function Y and then to apply the gradual deformation algorithm on Y. Starting with an initial realization y 0 of Y and another realization u 1 of Y independent of y 0 , we build a continuous chain of realizations y 1 (t) by using Equation ( 12).Then we minimize the objective function J with respect to parameter t.In this way, the N dimensional optimization problem is reduced to a onedimensional one, and moreover the optimized vector y 1 (t opt ) is actually a realization of Y.
However, the above optimization might not reduce the objective function to a low-enough level.It is then necessary to repeat the above procedure with a new realization chain y 2 (t) built by combining y 1 (t opt ) and another independent realization u 2 of Y.More generally, the following iterative optimization procedure can be used.At iteration n, the continuous chain of realizations y n (t) is written: where y n-1 is the optimized realization at iteration n -1, and where u n (n = 1, 2, ...) are a series of independent realizations of Y and they are also independent of the initial realization y 0 .Then by minimizing the objective function with respect to parameter t, we get a new realization y n (t opt ) which improves (or at least maintains) the calibration to the nonlinear data.This iterative procedure is stopped when a satisfactory calibration is reached.We note that, at each iteration, the objective function J depends only on parameter t.Therefore, we deal always with a one-dimensional optimization problem.Note also that, due to the mutual independence of u n (n = 1, 2, ...), the solution space can be reached when n becomes large enough and therefore the procedure converges always to a global minimum (i.e., a calibrated realization).
There are several efficient classic algorithms for minimizing the objective function J (e.g., Press et al., 1992).If J is differentiable with respect to t, then gradient based algorithms can be used.These algorithms require the computation of the derivative of the objective function with respect to t.This can be performed using the finite difference method which requires solving the forward problem twice for each t.If the sensitivity coefficients (i.e., the derivatives of the model responses f i (i = 1, 2, ..., p) with respect to the components z j (j = 1, 2, ..., N) of the model vector z) are available, the derivative of the objective function with respect to t is then derived from: (15) There is an extensive literature devoted to the computation of the sensitivity coefficients ∂f i / ∂z j (e.g., 1994).y j (j = 1, 2, ..., N) stand for the components of the Gaussian vector y.dy j / dt can be easily derived from Equation ( 14), although a conditioning kriging term must be included when linear data are involved.dz j / dy j are simply the derivatives of the Transformation function (1).In some cases, for instance when z is a truncated Gaussian model defined by Equation (4), dz j / dy j does not exist.This makes the objective function nondifferentiable with respect to t.Consequently, gradient based optimization algorithms cannot be directly applied and we must turn to algorithms for nondifferentiable optimization problems.For instance, the golden section search method can be used for the above one-dimensional optimization problem.For multivariate optimization, the simplex method can be used (Press et al., 1992).
Lastly, instead of using two independent realizations, it is possible to create a multidimensional process of realizations by using several independent realizations.Examples of calibration to production data by using the combination of multiple realizations are presented in Roggero and Hu (1998).Note that increasing the number of realizations in each iteration allows reducing the number of iterations in the calibration procedure.However, within each iteration, the search for an optimal combination is more expensive.Thus a compromise that depends on each application is required.The above approach was successfully applied to an mature oil field in the Brazilian offshore (Reis, 2001;Feraille et al., 2003).

Structural Gradual Deformation
In many cases, there is not enough data to accurately infer the structural parameters of a stochastic model (e.g., its mean, variance, covariance function, etc.).These structural parameters are often given in terms of intervals or a priori distributions.If the values of these structural parameters are biased, it will be desperate to find a realization calibrated to nonlinear data.Therefore these structural parameters must also be considered as fitting parameters.This involves the deformation of a stochastic model with respect to its structural parameters.Most stochastic simulation algorithms require at least the specification of the covariance before generating realizations.In order to deform a realization while modifying simultaneously its structural parameters, it is necessary to use a stochastic simulation algorithm that separates the generation of random numbers and the imposition of a structure.
Denote by X a standard Gaussian white noise and L the covariance operator, then: is a Gaussian field with the imposed covariance.In this form, it is clear that we can gradually deform the Gaussian white noise X by using the above method and modify at the same time the covariance function.That is: (17 where U and V are two independent standard Gaussian white noises.Figure 4 shows three realizations with the same Gaussian white noise but with different variogram ranges.
Structural gradual deformation.Three Gaussian realizations with the same random seed but with different variogram ranges.
There exist several simulation algorithms that separate the generation of random numbers and the imposition of a covariance function.The method based on the Cholesky decomposition of the covariance matrix is such a method, but limited to small simulations (say less than 1000 nodes).For large simulations, we can resort to the FFT moving average method (Le Ravalec et al., 2000).
All parameters of the covariance function, including ranges and directions of anisotropy, etc., are contained in the covariance operator L. As for the covariance type, instead of choosing somewhat arbitrarily a conventional covariance (exponential, spherical or Gaussian covariance), it is convenient to use the stable covariance family: (18) that exhibits distinct behaviors at origin depending on parameter α (α ∈ (0.2)).For instance, α = 1 corresponds to the exponential covariance and α = 2 corresponds to the Gaussian covariance.When α varies continuously from 1 to 2, we get a continuous chain of realizations whose covariance evolves continuously from the exponential model to the Gaussian model.Consequently, by coupling with an inverse procedure, one may identify the local smoothness (inaccessible from distant point data) of a physical field.Examples of calibrating structural parameters to production data are presented in Le Ravalec-Dupin et al. (2001).

Local Gradual Deformation
When the observations are scatted in different zones of the studied field, the calibration using global deformation may be inefficient, because improving the fitting in one zone may deteriorate the fitting in another zone.This led us to consider the gradual deformation zone by zone.Consider a partition of the field into n zones.Let X be a Gaussian white noise in the whole field and X 1 , X 2 , ..., X n the partition of X in the n zones.As X 1 , X 2 , ..., X n are mutually independent, it is then possible to perform their gradual deformation individually.By applying the covariance operator L on these independent white noises, we obtain a consistent correlated model Y: (19) where t = (t 1 , t 2 , ..., t n ) and where U i and V i for i = 1, 2, …, n are mutually independent Gaussian white noises.For a given set of realizations of U i and V i , we solve an optimization problem of n parameters t 1 , t 2 , ..., t n to obtain a realization which improves (or at least maintains) the calibration to the data.This procedure can be iterated until a satisfactory calibration is reached.
The idea of local (or regionalized) deformation is inspired by the traditional practice of zonation for history matching in petroleum engineering and for model calibration in subsurface hydrology.The traditional zonation method consists in dividing the whole reservoir field in several zones and performing data calibration separately in each zone.This practice does not guaranty the spatial continuity of the reservoir model between zones.However, instead of directly changing the model property like the traditional method, the local deformation method proposed in this paper intervenes in the underlying Gaussian white noise of the reservoir model.Then by applying the covariance operator on the modified Gaussian white noise, the spatial continuity of the stochastic model is preserved.Consequently, this avoids the drawback of the traditional method.Another attractive feature of the local deformation method is its applicability to the updating of an existing reservoir model when new data are available.For instance, when a new well is drilled, only a limited zone around the well need to be updated to fit the new data, while preserving the other parts of the existing model.Figure 5 shows an example of local deformation of a Gaussian simulation.Examples of local calibration to production data are presented in Reis (2001) and in Le Ravalec-Dupin et al. (2001).
Local gradual deformation.Two Gaussian simulations that are different only around the small window.

Generalization to Non-Gaussian Simulations
The gradual deformation algorithm characterized by Equation ( 12) cannot be directly extended to non-Gaussian-related models, simply because, in general, linear combination of independent non-Gaussian random functions will not preserve the non-Gaussian distribution, although the variogram is preserved.For instance, a linear combination of two independent uniform variables does not follow the uniform distribution but a triangular distribution.
However, from the computing point of view, any stochastic simulation procedure can be regarded as an operation that transforms a series of pseudo uniform numbers (generated by a pseudo random number generator) into a set of numbers of interest that are distributed according to the required distribution.On the other hand, uniform numbers can be transformed to Gaussian numbers by the relation: (20) where G is the standard Gaussian cumulative distribution function, and where U is uniform variable between 0 and 1. Therefore in principle, the gradual deformation algorithm for Gaussian models can be used to deform gradually any kind of stochastic models.In particular, we have applied the above principle to general sequential simulations (Hu et al., 2001) and to Boolean simulations (Hu, 2003).Figure 6 shows an example of gradual migration of objects of a Boolean simulation.We observe in this example the progressive destruction of a flow path (in red color) in the upper part of the field and the progressive construction of a flow path in the lower part of the field.

CONCLUSIONS
Modeling heterogeneous reservoirs requires resorting to a probabilistic approach, which allows, on one hand, to integrate data of various natures and scales and on the other hand, to evaluate uncertainties by generating multiple possible scenarios of the reservoir heterogeneity.Building a stochastic reservoir model must account for the geological depositional environment of the reservoir, so as to represent the major heterogeneities that control fluid flow.The quantitative information from wells, seismic and well test, etc. must be used for the inference of the model structural parameters.Constraining model realizations to the hydrodynamic data from production allows to further increase the model reliability for production forecasts.
Several geostatistical models are available.The Gaussianrelated models or more generally pixel-based models have the advantage of being more easily constrained by quantitative measurements.But they are often geologically poor.On the contrary, the object-based models and the processbased (or genetic) models allow to better account for geological concepts but their constraint to quantitative measurements is more difficult.
The object-based models can be regarded as a compromise between the models based on the random functions and the models based on the geological depositional processes.In the object-based Boolean model, for instance, the locations of objects are defined according to a Poisson point process that is a purely mathematical concept, while the shape of objects can have a geological sense.The Boolean model is of great interest when geological objects can be clearly defined like in fault and fracture networks, vacuolar media, meandering channel systems, etc.There is an algorithm for conditioning Boolean simulations to lithological data.Their calibration to production data can be performed by the migration and deformation of objects based on the gradual deformation method.
The process-based (genetic) models mimic geological processes of deposition and thus allow better representing reservoirs geological architecture.The real applicability of these models will depend on their flexibility to be constrained not only to hydrodynamic measurements (well tests, production history) but also to static data (lithological description along wells, seismic maps, etc.).One possibility lies in the statistical pattern recognition or more particularly in the multi-point statistics (Strebelle and Journel, 2000).Genetic simulations can provide nonconditional realizations that can be used as training images for the inference of pattern statistics or multi-point statistics.From these statistics, it is possible to regenerate realizations by sequential simulation.Conditioning sequential simulations to lithological data is immediate.For the calibration to hydrodynamic data, one can Gradual deformation of a Boolean simulation.Chain of realizations of a Boolean model with elliptic objects (in white).
use the gradual deformation method that is perfectly compatible with the algorithm of sequential simulation.
The hydrodynamic calibration of a reservoir model is formulated as an optimization problem.Optimizing a stochastic model defined on a large grid requires reparameterizing the model, on the one hand to reduce the number of parameters to be optimized and on the other hand to preserve the spatial structure of the model.The gradual deformation is a parameterization method that satisfies both of these requirements.This method can be used not only for the global deformation but also for the local deformation of a reservoir model.This allows for progressive updating of the model during the life of a field.Because of uncertainties in model structural parameters, there is also the possibility to consider them as fitting parameters.Moreover, the gradual deformation method is not limited to Gaussian-related models.Different forms of gradual deformation can be integrated in an iterative optimization procedure to constrain any type of stochastic simulations to hydrodynamic data.Examples from synthetic and real case studies show the applicability of the approach.Methods for further acceleration of the above calibration procedure are proposed in Hu and Le Ravalec-Dupin (2003) and in Schaaf et al. (2003).
Another way to accelerate the calibration procedure consists in reducing the time of fluid flow simulation for each realization.That can be done either by upscaling the fine geological model or by using a fast alternative flow simulator.Work of Mezghani and Roggero (2001), Schaaf et al. (2002) are devoted to the integration of an upscaling step into the procedure for calibrating fine scale models, while Le Ravalec-Dupin and Fenwick (2002), Caers (2003) explore the use of an alternate simulator based on streamlines in the procedure of calibration.
Updating a stochastic reservoir model by the integration of more and more available data reduces the distance of the model from the reality.But the real reservoir architecture is never accessible completely and there is always a degree of uncertainty.In this paper, we focused on obtaining a constrained model realization from a nonconstrained realization.The iteration of the same procedure using other nonconstrained realizations generates other constrained realizations of the stochastic reservoir model.All these constraint realizations are coherent with our conception of the geological environment of deposition and with all the available quantitative measurements.Our principal concern was to preserve the model spatial structure inferred from physical measurements.These constrained realizations are to be used to optimize the field development and production project while evaluating associated uncertainties.Nevertheless, one can wonder whether these constrained realizations, even in a sufficient number, are representative of all possible constrained realizations (or the a posteriori distribution in the Bayesian terminology (Tarantola, 1987)) of the adopted reservoir model.Le Ravalec-Dupin et al. (2000) carry out a preliminary investigation on this question.Research is in progress to further clarify this question.
Nevertheless, even if one can generate a set of representative realizations of the constrained stochastic model, one does not necessarily control all uncertain aspects in the process of reservoir characterization.The quantification of uncertainties is always within the framework of an adopted stochastic reservoir model.However, the choice of a model is often somewhat subjective and therefore uncertain.We believe the uncertainty on the choice of a model is partially accounted for by the fact that we can calibrate the structural parameters of the adopted model.In addition, we can also consider several types of plausible models.We thus increase the reliability of the production forecast.But from practical as well as conceptual points of view, no one can claim to cover all uncertain aspects.There is thus always a degree of doubt."Doubt is not a pleasant condition, but certainty is a ridiculous one" (Voltaire).
Figure 3Global gradual deformation.Chain of realizations of a Gaussian random function with a Gaussian variogram (afterRoggero and Hu, 1998).