Best Paper 2001

Petroleum Geoscience 2000

Gradual deformation and Iterative Calibration of Truncated Gaussian Simulations (Suggested running head: Gradual deformation of Truncated Gaussian Simulations)

Lin Y. Hu, Mickaële Le Ravalec and Georges Blanc

ABSTRACT Recently, we proposed the gradual deformation approach for constraining stochastic models to dynamic data (well-tests and production history). In this paper, we review the basic gradual deformation algorithm and extend its application to different types of truncated Gaussian simulations (including non stationary truncated Gaussian and truncated pluri-Gaussian simulations). A case study on the calibration of a reservoir lithofacies model to well-test pressure data illustrates the efficiency of the gradual deformation approach.

INTRODUCTION Since its introduction in the petroleum industry, the truncated Gaussian simulation method (Matheron et al. 1987) has been widely used for modeling reservoir heterogeneity. This method consists in truncating a continuous Gaussian simulation with multiple thresholds to represent the distribution of lithofacies in a heterogeneous reservoir. This method was further extended to non stationary truncated Gaussian simulations (Beucher et al. 1993) and to truncated pluri-Gaussian simulations (Galli et al. 1994). It has been proven to be highly flexible for representing a wide range of geological patterns and shapes. A few methods were developed for conditioning truncated Gaussian simulations to indicator data of lithofacies along wells (Freulon & de Fouquet 1993). There are also methods for incorporating seismic derived information or well-test derived information in truncated Gaussian simulations (Moulière et al. 1997; Hu et al. 1998). However these methods cannot be extended to the integration of production data. Recently, the gradual deformation approach (Hu 2000) was proposed for constraining stochastic models to production history data. This approach is based on the algorithm for gradually deforming a Gaussian-related stochastic model while preserving its spatial variability. This algorithm consists in building a continuous stochastic process, whose state space is the

1

ensemble of the realizations of a spatial stochastic model. The gradual deformation of a realization induces generally smooth variations in the objective function, which measures the mismatch between the production data observed from the oil field and the corresponding data calculated in the realization. Therefore this objective function can be minimized using an efficient optimization algorithm (for instance a gradient based algorithm or the golden section research algorithm). An iterative procedure is used to progressively reach a satisfactory solution. This paper presents the gradual deformation approach applied to truncated Gaussian models. We first describe the gradual deformation of different types of truncated Gaussian simulations. Then we present a case study to illustrate the calibration of a reservoir lithofacies model to welltest pressure data.

GRADUAL DEFORMATION GAUSSIAN SIMULATIONS

OF

TRUNCATED

Gradual deformation of a truncated Gaussian simulation consists in deforming the underlying continuous Gaussian simulation. A gradual deformation of the Gaussian simulation induces a smooth change in the corresponding lithofacies model. The basic algorithm for the gradual deformation of a Gaussian simulation consists in generating two independent realizations y1 and y 2 of the standardized Gaussian random function Y , and then building a realization chain by linear combination: (1) y (t ) = y1 cos t + y 2 sin t It can be shown that, for any t , y (t ) is a realization of Y . Hence, by truncating y (t ) , we obtain a continuous realization chain of the indicator random function. More details and several extensions of this algorithm (e.g., local and structural deformations) are presented in (Hu 2000; Le Ravalec et al. 2000). Here, we will provide examples of the gradual deformation of different variants of truncated Gaussian simulations.

Stationary truncated Gaussian simulations

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 : 1 Y ( x) ≥ s (2) I ( x) =  0 otherwise 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). Figure 1 depicts a truncated Gaussian simulation with three lithofacies and its gradual deformation. In this example, disjoint truncation intervals are used to define lithofacies. This allows direct transitions between any two lithofacies.

Non stationary truncated Gaussian simulations 2

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 kriging 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 randon function (Adler & Thovert 1998). Figure 2 shows a two lithofacies model obtained by truncating a Gaussian simulation with another Gaussian simulation. The black lithofacies may represent channel belts. The gradual deformation is applied to one of the Gaussian simulation, resulting in a gradual deformation of the lithofacies model.

Truncated pluri-Gaussian simulations Truncated pluri-Gaussian simulations are a natural extension of truncated single Gaussian simulations. They involve the combination of several truncated Gaussian simulations. This makes it possible to generate extremely varied lithofacies models. (Galli et al. 1994; Le Loc'h & Galli 1997). Unlike the single 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, presented in the previous paragraph, is a particular example of pluri-Gaussian simulations. In general, the gradual deformation of a truncated pluriGaussian simulation consists in applying the gradual deformation algorithm to each underlying Gaussian simulation. Figure 3a shows the gradual deformation of a three lithofacies model built by truncating a Gaussian simulation and its derivative along the X axis. We used the same method as Armstrong and Galli (1999) to generate the Gaussian simulation and its derivative, except that the truncation rule (Figure 3b) is slightly modified to allow the grey lithofacies (e.g. representing crevasse splays) to be attached to only one side of the black lithofacies (e.g. representing channels).

Remark on the conditioning to well data 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 & Geman 1984; Freulon & de Fouquet 1993) and then by applying conditioning kriging 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. It can be noticed that in this particular case, the truncated Gaussian method becomes equivalent to the probability field method (Srivastava 1992; Froidevaux 1993) for simulating lithofacies. Indeed, 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. In the following, we keep in mind that realizations of a lithofacies model are always conditioned to lithofacies indicator data along wells, and we focus only on the problem of calibration to production data.

3

ITERATIVE CALIBRATION GAUSSIAN SIMULATIONS

OF

TRUNCATED

The calibration of a stochastic model to production data can be formulated as an optimization problem. This involves the definition of an objective function that measures the mismatch between the data observed on the oil field and the corresponding data calculated on the realizations of the stochastic model. Then this objective function is minimized in the ensemble of realizations of the stochastic model. In general, a mean square formulation is used for the definition of the objective function. The basic idea of the gradual deformation approach is to parameterize the reservoir model using formula (1) and to minimize the objective function with respect to the parameter t . Thus, we solve simultaneously two major problems of the stochastic optimization. First we reduce a high dimensional optimization problem to a low dimensional one, and second we preserve the spatial variability of the stochastic model. The procedure is as follow. Starting with an initial realization y1 of Y and another realization y 2 of Y independent of y1 , we build a continuous chain of realization y (t ) by using the formula (1). We minimize the objective function with respect to parameter t. Thus, the problem of N dimensional optimization problem is reduced to a one-dimensional problem. In addition, the optimized vector y (t opt ) is actually a realization of Y . However, the above optimization might not reduce the objective function to a low-enough level. In such a case, we replace the initial realization y1 by the optimized realization y (t opt ) , and we iterate the above procedure until a satisfactory calibration to production data is obtained. When performing optimization with lithofacies models generated, for instance, by the truncated Gaussian method, the objective function is generally not differentiable with respect to the parameter t . The gradient-based optimization algorithms cannot be used. For the above onedimensional optimization, the golden section search method can be used (e.g., PRESS et al. 1992). In the case of pluri-Gaussian simulations, it may necessary to deform simultaneously several Gaussian simulations. This involves several chains of realizations and therefore several parameters to optimize. In such conditions, we can use the simplex method to perform optimization. (e.g., PRESS et al. 1992).

CASE STUDY: CALIBRATING A LITHOFACIES MODEL TO WELL-TEST PRESSURE DATA To illustrate the application of the above stochastic optimization method, consider the calibration of a stochastic reservoir model to well-test pressure data. The construction of the reservoir model is based on geological knowledge of a real oil field. The field is composed of three lithofacies: two reservoir lithofacies (lithofacies 1 and 2) and one non-reservoir lithofacies (lithofacies 3). Table 1 summarizes the petrophysical properties of the three lithofacies. The pluri-Gaussian method was used to simulate the specific lithofacies distribution of the oil field. We first generated lithofacies 3 by truncating a Gaussian realization. Then in the complementary part of lithofacies 3, we generated lithofacies 1 and 2 by truncating another Gaussian realization independent of the previous one. The field is discretized on a regular grid of 60 X 59 X 15 blocks

4

of size 15m X 15m X 1.5m. A spherical variogram model was used for the generation of the Gaussian realizations. The main anisotropic direction is diagonal with respect to the reservoir grid. The variogram ranges of the Gaussian random function used for generating lithofacies 3 are respectively 300m, 80m, 3m in the three anisotropic directions. The variogram ranges of the Gaussian random function used for generating lithofacies 1 and 2 are respectively 150m, 40m, 1.5m in the three anisotropic directions. The proportions of lithofacies 1, 2 and 3 are 6%, 16% and 78%, respectively. Figure 4a shows the middle layer of a realization used as the synthetic reservoir for this study. A numerical well-test was performed with a finite difference well-test simulator (Blanc et al. 1996). The well travels horizontally through the middle layer of the synthetic reservoir along the x axis. Figure 4a shows the perforated section of the well. The diameter of the well is 7.85cm, the well bore storage is 0 and the skins of lithofacies 1, 2 and 3 are 0., 3. and 50, respectively. The numerical well-test was assumed to last 240 days with a constant rate of 5 m3/day (reservoir condition) so as to investigate almost the entire reservoir field. Figure 4b shows the pressure variation and its logarithmic derivative. Our objective is to build realizations of the lithofacies model constrained by lithofacies data along the well and by the well-test pressure curve. The objective function is defined as the sum of the square differences between the pressure responses of the synthetic reservoir and the pressure responses of a realization. Because the dynamic behavior of the reservoir model is mainly controlled by the contrast between the reservoir and non-reservoir lithofacies, the Gaussian realization used for generating the reservoir lithofacies 1 and 2 was fixed at once and only the Gaussian realization used for generating the non reservoir lithofacies 3 was deformed to match the pressure data. Figures 5a and 5b show two examples of the objective function of two realization chains as defined by formula (1). We observe that although the lithofacies model is gradually deformed, the corresponding objective function profiles present discontinuities caused by the discrete nature of the lithofacies model. Note, in particular, the huge drop of the objective function from realization 86 to realization 87 (Figure 5a). A detailed investigation of the realizations 86 and 87 reveals that this discontinuity is due to a change of lithofacies in a single grid block which provides a shortcut to the fluid flow. Figures 6a and 6b show the middle layer of two initial realizations constrained only by the lithofacies data along the well. The numerical well-tests on these realizations result in pressure responses very different from that of the synthetic reservoir (Figures 7a and 7b). Starting respectively from these two independent realizations, and using the iterative optimization method, we obtained the two calibrated realizations (Figures 8 and 9) after a few iterations. Figures 10a and 10b show their respective objective functions as the number of iterations increases.

CONCLUSIONS This paper extends the application of the gradual deformation method to different types of truncated Gaussian simulations (including non stationary truncated Gaussian and truncated pluriGaussian simulations). The gradual deformation of a truncated Gaussian realization yields in general an objective function with discontinuities. This is due to the discrete nature of lithofacies models. Gradient-based optimization methods cannot be used, hence we apply a non gradientbased method like golden section search for minimizing the objective function. Results from the 5

case study demonstrate the efficiency of our approach to calibrate a lithofacies model to well-test pressure data. A set of calibrated realizations can be obtained by repeating the gradual deformation procedure from a set of independent initial realizations. Our presentation is limited to the global deformation of truncated Gaussian simulations with fixed statistical parameters (lithofacies proportions, variogram ranges etc.). However, the method can be extended to local gradual deformation. Local deformation may significantly improve calibration speed in cases where observations are spread over different areas of a field. The deformation of a realization can also be performed by modifying the statistical (structural) parameters of the stochastic model. In other words, calibration can be performed simultaneously with respect to realizations and statistical parameters.

ACKNOWLEDGMENTS This work is financed by the joint research project in reservoir engineering between the Elf Exploration & Production and the Institut Français du Pétrole. We are grateful to A. Journel, O. Dubrule, A. Hurst and a anonymous reviewer for their helpful comments.

REFERENCES ADLER, P.M. & THOVERT, J.F. 1998. Real porous media: Local geometry and macroscopic properties. Applied Mechanics Reviews, 51, No.9, 537-585. ARMSTRONG, M. & GALLI, A. 1999. Derivative based plurigaussian simulations. In: Proceedings of the 5th Annual Conference of the International Association for Mathematical Geology, Trondheim, Norway, 6-11 August 1999. BLANC, G., GUERILLOT, D., RAHON, D. & ROGGERO, F. 1996. Building geostatistical models constrained by dynamic data - a posteriori constraints. Paper SPE 35478, In: Proceedings of the NPF/SPE European 3D Reservoir Modelling Conference, Stavanger, Norway, 16-17 April 1996. BEUCHER, H., GALLI, A., LE LOC’H, G., RAVENNE, C. & HERESIM GROUP 1993. Including regional trend in reservoir modeling using the truncated Gaussian method. In: Soares, A. (ed.), Geostatistics Troia'92, Kluwer Academic Publishers, Dordrecht, 1, 555-566. FROIDEVAUX, R. 1993. Probability field simulation. In: Soares, A. (ed.), Geostatistics Troia'92, Kluwer Academic Publishers, Dordrecht, 1, 73-83. FREULON, X. & DE FOUQUET, C. 1993. Conditioning a Gaussian model with inequalities. In: Soares, A. (ed.), Geostatistics Troia'92, Kluwer Academic Publishers, Dordrecht, 1, 201-212. GALLI, A., BEUCHER, H., LE LOC’H, G., DOLIGEZ, B. & HERESIM GROUP 1994. The pros and cons of the truncated Gaussian method. In: Armstrong, M. & Dowd, P. (eds), Geostatistical simulations, Kluwer Academic Publishers, Dordrecht 217-233. GEMAN, S. & GEMAN, D. 1984. Stochastic relaxation, Gibbs distribution and Bayesian restoration of images I.E.E.E. transactions: Pattern analysis and machine intelligence, 6, 721-741. HU, L.Y., BLANC, G. & NOETINGER, B. 1998. Estimation of lithofacies proportions by use of well and well-test data. SPE Reservoir Evaluation & Engineering, 1, No.1, 69-74. HU, L.Y. 2000. Gradual deformation and iterative calibration of Gaussian-related stochastic models. Mathematical Geology, 32, No.1, 87-108. JOURNEL, A.G. & HUIJBREGTS, CH.J. 1978. Mining Geostatistics, Academic Press, London. LE LOC'H, G. & GALLI, A. 1997. Truncated plurigaussian method: theoretical and practical points of view. In: Baafi, E.Y. & Schofield, N.A. (eds), Geostatistics Wollongong 96, Kluwer Academic Publishers, Dordrecht, 1, 211-223.

6

LE RAVALEC, M., NOETINGER, B. & HU, L.Y. 2000. The FFT moving average (FFT-MA) generator: An efficient numerical method for generating and conditioning Gaussian simulations. Mathematical Geology, 32, No.6. MATHERON, G., BEUCHER, H., DE FOUQUET, C., GALLI, A. & RAVENNE, C. 1987. Conditional simulation of the geometry of fluvio-deltaic reservoirs. Paper SPE 16753, In: Proceedings of the SPE Annual Technical Conference and Exhibition, Las Vegas, 22-25 September 1987. MOULIERE, D., BEUCHER, H., HU, L.Y., FOURNIER, F., TERDICH, P., MELCHIORI, F. & GRIFFI, G. 1997. Integration of seismic derived information for reservoir stochastic modeling using truncated Gaussian approach. . In: Baafi, E. Y. & Schofield, N.A. (eds), Geostatistics Wollongong 96, Kluwer Academic Publishers, Dordrecht, 1, 374-385. PRESS, W.H., TEUKOLSKY, S.A., VETTERLING, W.T. & FLANNERY, B.P. 1992. Numerical Recipes in Fortran 77: The Art of Scientific Computing. 2nd edition, Cambridge University Press. SRIVASTAVA, R.M. 1992. Reservoir characterization with probability field simulation. Paper SPE 26753, In: Proceedings of the SPE Annual Technical Conference and Exhibition, Washington, D.C., 4-7 October 1992.

7

FIGURES

Figure 1: Gradual deformation of a stationary truncated Gaussian simulation with three lithofacies over a 2 dimensional grid of 100X100 nodes. The Gaussian simulation is generated with an isotropic Gaussian variogram. The black lithofacies corresponds to the simulated Gaussian numbers in the interval [−0.5,0.5) , the white lithofacies to Gaussian numbers in (−∞,−1) or in [0.5,1) , the grey lithofacies to Gaussian numbers in [1,+∞) or in [−1,−0.5) .

Figure 2: Gradual deformation of a non stationary two lithofacies model, over a 2 dimensional grid of 100X100 nodes, constructed by truncating a Gaussian simulation with another Gaussian simulation. The Gaussian simulations are generated with anisotropic Gaussian variograms. The

8

black lithofacies corresponds to the situation where the first Gaussian simulation is greater than the second Gaussian simulation, contrarily to the white lithofacies.

channels).

Figure 3a: Gradual deformation of a three lithofacies model, over a 2 dimensional grid of 100X100 nodes, built by truncating a Gaussian simulation and its derivative. The Gaussian simulation is generated with an isotropic Gaussian variogram.

Gaussian

G W

B

W

G Derivative

Figure 3b: Rule for the truncation of a Gaussian simulation and its derivative (W: White lithofacies, B: Black lithifacies, G: Grey lithofacies).

9

Kx (md) Ky (md) Kz

φ

ct (10-5 bar1

(md)

(%)

10.

10.

17

2.1857

lithofacies 2 1.

1.

1.

14

2.0003

lithofacies 3 0.1

0.1

0.001

9

1.8148

Lithofacies

10.

)

1

Table 1: Petrophysical properties of the three lithofacies. (Kx, Ky, Kz = permeability along x, y, z directions; φ = porosity; ct = total compressibility)

1e+2 (bar)

Pressure

1e+1

Well

Facies 1 Facies 2

Derivative

Facies 3

1e+0 1e+4

1e+5

1e+6

Time (s) 1e+7

1e+8

Figure 4b: Pressure and its logarithmic derivative Figure 4a: Middle layer of the 3D synthetic reservoir versus time (in seconds) from a numerical well-test on generated by the truncated Gaussian method. The field the synthetic reservoir. is 900m, 885m and 22.5m along x, y and z directions.

10

100

100

Objective function

10

10

86

1

1

0.1

0,1

87 -1

-0.75

-0.5

-0.25

0

0.25

0.5

t/π π

t/π π

0.01 0.75

Objective function

0,01 1

-1

a

-0,75

-0,5

-0,25

0

0,25

0,5

0,75

1

b

Figure 5: Two examples of the objective function versus the parameter

t /π

of the continuous chain of realizations.

Well

Well

b

a

Figure 6: Two initial realizations (middle layer) of the lithofacies model conditioned to only lithofacies data along the well.

11

1e+2

1e+2

(bar)

(bar)

Pressure

Pressure

1e+1

1e+1

Derivative 1e+0 1e+4

1e+5

1e+6

Derivative

Time (s) 1e+7

1e+0 1e+4

1e+8

a

1e+5

1e+6

Time (s) 1e+7

1e+8

b

Figure 7: Comparison between the pressure curves of the synthetic reservoir and the pressure curves of the initial models.

Well

Well

a

b

Figure 8: Two realizations (middle layer) of the lithofacies model conditioned to lithofacies data along the well and to the well-test pressure curve.

12

1e+2

1e+2

(bar)

(bar)

Pressure

Pressure

1e+1

1e+1

Derivative 1e+0 1e+4

1e+5

Derivative

Time (s)

1e+6

1e+7

1e+0 1e+4

1e+8

a

1e+5

1e+6

Time (s) 1e+7

1e+8

b

Figure 9: Comparison between the pressure curves of the synthetic reservoir and the pressure curves of the calibrated models.

100

100 Objective function

Objective function

10

10

1

1

0,1

0,1

0,01

0,01 Number of iterations

0,001 0

1

2

3

4

5

6

7

8

Number of iterations

0,001 9

0

a

1

2

3

4

5

6

7

8

9 10 11 12 13 14

b

Figure 10: Objective function versus the number of iterations (number of continuous chains of realizations).

13