IBP2598_10 GEOSTATISTICAL HISTORY MATCHING WITH DIRECT TRANSFORMATION OF IMAGES. APPLICATION TO A MIDDLE EAST RESERVOIR Maria Helena S. C. Caeiro1, Amílcar O. Soares 2,Jorge Tiago Santos3, Álvaro Carvalho 4, Luís Guerreiro5. Copyright 2010, Instituto Brasileiro de Petróleo, Gás e Biocombustíveis - IBP Este Trabalho Técnico foi preparado para apresentação na Rio Oil & Gas Expo and Conference 2010, realizada no período de 13 a 16 de setembro de 2010, no Rio de Janeiro. Este Trabalho Técnico foi selecionado para apresentação pelo Comitê Técnico do evento, seguindo as informações contidas na sinopse submetida pelo(s) autor(es). O conteúdo do Trabalho Técnico, como apresentado, não foi revisado pelo IBP. Os organizadores não irão traduzir ou corrigir os textos recebidos. O material conforme, apresentado, não necessariamente reflete as opiniões do Instituto Brasileiro de Petróleo, Gás e Biocombustíveis, seus Associados e Representantes. É de conhecimento e aprovação do(s) autor(es) que este Trabalho Técnico seja publicado nos Anais da Rio Oil & Gas Expo and Conference 2010.

Resumo Este artigo apresenta uma nova metodologia de history matching geoestatístico, que visa integrar os dados dinâmicos dados históricos de produção de vários poços produtores - num modelo estocástico de propriedades internas. Este é um algoritmo interativo que pretende resolver o problema típico dos procedimentos inversos de history matching. A ideia baseia-se numa abordagem de perturbação global de um campo de permeabilidade, em que seja encontrada a calibração desejada através de uma função objectivo. Este trabalho está baseado na abordagem de transformação directa de imagens de permeabilidade (Mata-Lima et al, 2007), mas com um novo método de perturbaçção global, que apresenta uma convergência muito mais rápida à função objetivo. A perturbação (co-simulação directa) é realizada localmente, tendo por base as regiões de influência de cada poço e a função objectivo é multi-critério, combinação linear das diferenças entre os dados simulados de water cut, pressão e produção de óleo. O algoritmo proposto foi implementado num reservatório real no Oriente Médio.

Abstract This paper presents a new methodology for geostatistical history matching that aims to integrate the dynamic data historical production data of several producer wells - into the stochastic model of internal properties. This is an interactive algorithm to solve the typical inverse problem of history matching procedure. It is based on a global perturbation approach of, for example, a permeability field until the desired match, with a predefined objective function, is achieved. This paper is based on the direct transformation of images approach (Mata-Lima et al, 2007) but with a new method of global perturbations, with a much faster convergence to the objective function. The proposed algorithm was implemented in a real Middle East Reservoir. The perturbation (direct co-simulation) is performed locally, based on the influence regions of each well. The objective function is a linear combination of differences between simulated and real data of water-cut, pressure and oil production.

1. Introduction One of the main concerns in petroleum reservoir characterization field is to reduce the uncertainty in reservoir internal characteristics which can be achieved by reducing uncertainty in reservoir modeling through an efficient history matching. History matching is normally an inverse procedure in which an objective function (OF) is iteratively minimized as much as possible. The main difficulties of the problem – given the high number of independent variables this is a typical ill posed problem and there is an highly non-linear relation between the solution and the independent variables- still remain as challenge in this field of work. The proposed methodology is based on a global perturbation approach in the iterative path to reach the objective function. The seminal idea of these global perturbation algorithms comes from Hu et al (2001), with gradual deformation of Gaussian fields. The main advantage of this approach over the

______________________________ 1 Geologic and Mining Engineer – INSTITUTO SUPERIOR TÉCNICO 2 Professor – INSTITUTO SUPERIOR TÉCNICO 3 Student of Physics Engineering – INSTITUTO SUPERIOR TÉCNICO 4 Reservoir Engineer – PARTEX OIL & GAS 5 Reservoir Engineer – PARTEX OIL & GAS

Rio Oil & Gas Expo and Conference 2010 traditional methodologies of history matching is the simplicity of implementation. Other algorithms explored different methodologies for the global perturbation step of the inverse approach: a probability perturbation method (Caers, 2003), the use of direct sequential simulation and co-simulation to transform images of permeability (Mata-Lima et al, 2007). All these approaches follow, basically, the sequence of 4 steps: i) Generation of images of internal properties reflecting the spatial patterns of a geological conceptual model (main geological features, variograms, reference images); ii) Calculation of dynamic responses (dynamic simulator) of such set of images, given the existing wells (producers, injectors); iii) Evaluate the objective function based on a mismatch between the simulated responses and the real dynamic data; iv) Perturbation of initial images, based on the convergence to the objective function and return to i) until a desired match of the objective function is reached.

2. Proposed algorithm for Geostatistical History Matching The proposed algorithm aims to characterize a petroleum reservoir by an iterative algorithm for inverse modeling which integrates dynamic (production) data into stochastic models. It uses direct sequential co-simulation (coDSS) as a global perturbation method of permeability fields, in order to reach the convergence to the objective function. The flowchart of the iterative algorithm is presented in Figure 1 and can be summarized in the following sequence: i) Simulation of a set of realizations of permeability, honoring the well data, seismic information, geological constraints and main spatial patterns as revealed by the variograms; ii) Evaluation of dynamic responses of producing wells for each image/realization; iii) Calculation of the match between simulated dynamic responses and the real production data, through an objective function. If the mismatch is sufficiently low, the process stops. Otherwise a new generation of images is created (step iv) ; iv) Generation of a new set of images. In this step it is proposed a new method of transforming the images of permeability, such as they match of the dynamic response with the real production data but maintaining the main spatial patterns. Firstly at each iteration step, one “best” image is “composed” by a linear combination of the images corresponding to the closest dynamic responses. For this purpose the distances between the responses and from the responses to the real dynamic data (target) are calculated with Modified Hausdorff Distances (MHD), defining a 2-D referential where the simulated images with the closest responses to the target can be represented. This “composed” image does not have the spatial continuity of original images and it is just used as secondary information for the next images to be generated. This is performed with direct sequential co-simulation algorithm that accounts for the conditioning data – the existing permeability well data and the secondary information, of previous iteration, which influence is determined by a correlation measure of the match obtained in the previous step. The process returns to step ii) for the evaluation of the dynamic responses of the new images, and continues until the desired match is achieved through the minimization of the objective function.

2

Rio Oil & Gas Expo and Conference 2010 1st iteration Soft data (permeability model)

Co-DSS

Hard data (wells)

DSS

Hard data (wells)

Global correlation=0.6

Stochastic realizations Composite correlation coefficient cube (correlation data)

Stochastic realizations

Dynamic Simulator ynamic Simulator (Eclipse) Composite permeability image

Equivalent dynamic realizations

Objective function

No

Equivalent dynamic realizations

Real production data

Best responses to OF Calculation of the local correlation coefficients

Eclipse

Objective function

Desirable minimum?

Yes Composite “best” image of permabilities

Best responses to OF

Stop

Figure 1: Procedure of the proposed algorithm

2.1. Global transformation of Images Let us consider Z1(x) the initial image that one wish to transform in another Z2(x) with the same spatial pattern statistics, e.g. variogram and global histogram: C1(h) , γ1(h) , Fz1(z). Applying the Direct Sequential Co-Simulation algorithm, in each location x0 of a regular grid node, the local mean and variance are calculated by simple collocated co-kriging and estimation variance (equation 1).

Z 2 ( x0 ) * − m2 ( x0 ) = ∑ λα ( x0 )[Z 2 ( xα ) − m2 ( xα ) ] + λ1 ( x0 )[Z1 ( x0 ) − m1 ( x0 ) ]

(1)

α

Since the models γ1(h) and γ2(h) are the same, Markov approximation is quite acceptable to this case, i.e., the corregionalization model is fully defined with the correlation coefficient (correlogram at h=0) ρ1,2(0) between Z1(x) and Z2(x), such is in equation 2. ρ1,2(h) = ρ1,2(0) ρ2,2(h)

(2) 3

Rio Oil & Gas Expo and Conference 2010

The collocated cokriging system is written with the univariate correlogram of variable Z2(x), ρ2,2(h) and the coefficient correlation ρ1,2(0) (equation 3).

⎧ Nα ⎪∑ λβ 2 ( x 0 ).ρ 2 , 2 (xα , x β ) + λ1 ( x0 ).ρ 2 , 2 ( xα , x0 ).ρ1, 2 (0 ) = ρ 2 , 2 ( xα , x0 ) ⎪ β =1 ⎨ Nα ⎪ λ ( x 0 ).ρ (x , x )ρ (0 ) + λ ( x ).ρ (0 ) = ρ (0 ) 2, 2 ~0 1, 2 1 0 2, 2 1, 2 β2 β ⎪⎩∑ β =1

(3)

From the initial image Z1(x) one can obtain N images of Z2(x) with the same statistics of spatial continuity of Z1(x). Identically, from Z2(x), a set of Z3(x) could be obtained and so on. The degree of transformation – semblance of both images - is only function of the correlation coefficient ρ1,2(0)=[-1,+1] : the transformed image Z2(x) is as much identical to original Z1(x) as ρ1,2(0) is close to 1; An inverse image of original Z1(x) is obtained as ρ1,2(0) tends to –1. 2.2. Transformation with multiple images Let us consider that one wish to obtain a transformed image Zt(x) based on a multiple set of Ni multiple images Z1(x), Z2(x), …, ZNi(x). The formalism of equation 1 can now be generalized for the direct co-simulation of Zt(x) having Z1(x), Z2(x),…ZNi(x) as auxiliary variables. The collocated cokriging estimator of Zt(x) is presented on equation 4. Ni

Z t ( x0 ) * − mt ( x0 ) = ∑ λα ( x0 )[Z t ( xα ) − mt ( xα ) ] + ∑ λi ( x0 )[Z i ( x0 ) − mi ( x0 ) ] α

(4)

i =1

Since the models γi(h), i=1, Ni, and γt(h) are the same, the application of Markov approximation is also in this case quite appropriated. i.e., the corregionalization models are totally defined with the correlation coefficients ρt,i(0) between Zt(x) and Zi(x). The affinity of the transformed image Zt(x) with the multiple images Zi(x) are determined by the correlation coefficients ρt,i(0). Hence, one can select the images which characteristics we wish to “preserve” in the transformed image Zt(x). 2.3 Local Screening Effect Approximation Let us assume that to estimate Zt(x0) (equation 4) the collocated value Zi(x0) of a specific image Zi(x), with the highest correlation coefficient ρt,i(0), screens out the influence of the effect of remaining collocated values Zj(x0), j ≠ i. Hence, equation (4) can be written with just one auxiliary variable: the “best” at location x0 (equation 5).

Z t ( x0 ) * − mt ( x0 ) = ∑ λα ( x0 )[Z t ( xα ) − mt ( xα ) ] + λi ( x0 )[Z i ( x0 ) − mi ( x0 ) ]

(5)

α

With the local screening effect, Ni images Zi(x) give rise to just one auxiliary variable. The Ni images are merged in on single image based on the local correlation coefficient criterion. 2.4. Stochastic Inversion: a new approach This new approach of stochastic inversion is based on two key ideas: i) Use the sequential direct co-simulation with the screening effect approximation as the method of transforming images which converge to an objective function. ii) To follow the sequential procedure of genetic algorithm optimization to converge the transformed images towards an objective function. 2.5. Objective function In the next step one intend to evaluate the dynamic response of each permeability field realization obtained in a given iteration. This is accomplished by running a deterministic flow simulator (Eclipse 100 Black Oil, Schlumberger) 4

Rio Oil & Gas Expo and Conference 2010 for each permeability model as input, as well as the locations of different producers and injectors wells. The match evaluation is done by comparing the simulated and historic values of well water cut total (WWCT) and well bottom hole pressure (WBHP). After each response is obtained, a “distance” to the real production data is evaluated. In this study one use the Modified Hausdorff Distance for that purpose. The Hausdorff distance is often used to solve problems of "shape matching". Can be applied to any object in 2D or 3D if it can be represented in point set. The Modified Hausdorff distance (MHD) based on the classical Hausdorff distance, by Dubuisson and Jain (1994), is a measure of similarity used for the purpose of object matching, and has better performance because it is easy to compute and does not require the solution of the corresponding problem. The MHD as high discriminatory power - it is robust to outliers that may have resulted from targeting errors and his value should increase with the increasing of the difference between two objects. Suzuki and Caers (2008), proposed the use of MHD to measure the distance of two points in a referential characterized by the dynamic responses and the similarity of the sand channels patterns of each image. Suppose two sets A and B, and ai and bi are any points that belong to the sets A and B, respectively. Being uai ubi and the locations of ai and bi, then the spatial distance between ai and bi can be defined as the norm of the vector ||uai-ubi|| or Euclidean distance between the locations of the points ai and bi. Using this definition, the spatial distance measured between ai and the set B, d(ai, B) is determined as shown on equation 6.

d (a i , B ) = min bi ∈B u ai − u bi

(6)

Thus, the distance between A and B is defined as:

d ( A, B ) =

1 Na

∑

ai ∈A

d (u ai , B )

(7)

f (d ( A, B ), d ( B, A) = max(d ( A, B ), d ( B, A))) A = a1 ,..., a N a

where:

{ } B = {b ,..., b } 1

Nb

Now let us consider the response values of WWCTli (Well Water Cut Total) of simulated realization l along the time periods ti=1, N and the equivalent real data production WWCT at a given well. Modified Hausdorff Distance between the two responses is shown on equation 8.

MHD 1 wwct =

1 N

N

1

∑ d (WWCT i − WWCT )

(8)

i =1

Identical distance is obtained with WBHP (Well Bottom Hole Pressure). The objective function (MHD)total normalizes the distances of the two variables, which are in completely different scales (equation 9). MHD is the value of the modified Hausdorff distance calculated for each simulation of a certain well MHD is the average value of the modified Hausdorff distance and σ is the standard deviation of the modified Hausdorff distance calculated.

MHDtotal =

⎛ MHDwwct − MHDwwct ⎜ ⎜ σ wwct ⎝

⎞ ⎛ MHDwbhp − MHDwbhp ⎞ ⎟ ⎟−⎜ ⎟ ⎜ ⎟ σ wbhp ⎠ ⎝ ⎠ 2

(9)

In order to perform the co-simulation (with the collocated co-kriging approach) one need to convert the MHD into a correlation coefficient value. Hence, to the highest distance value we assume the correlation coefficient CC = 0.

5

Rio Oil & Gas Expo and Conference 2010 The equation 10 is used to access to the equivalent measurent of correlation. Dhi is the smallest MHD of the considered group of simulations and Dhn is the MHD for the current simulation. The ratio R will be variable because the number of considered simulations increases as the number of iterations performed. n

R =

CC= 1-R

Dh i Dh

(10)

In case of considering a two variables objective function (WBHP_ pressure and WWCT_ water cut), the final correlation coefficient is given by the mean of the correlation values corresponding to the simulation that have a smaller objective function. 2.6. Composing “secondary” permeability image and evaluation of the correlation coefficient´s cube The final step of each iteration consists on the composing of the secondary image for the next set of stochastic realizations of the iterative process - co-DSS. The image corresponding to the best response of each well is selected to compose a permeability image. The “best” image is composed of the parts of each simulation corresponding to the best match, inside a grid of polygonal areas of influence of each well (Figure 2).

Layer 55

Figure 2: Representation of the polygonal regions and the composite image of permeabilities

The correlation coefficients are extrapolated inside each region, representing the correlation between the real production data and the simulated ones.

3. Reservoir Description The reservoir used to the validation of the proposed algorithm is a complex and very heterogeneous carbonate reservoir. The field is currently producing under natural depletion with aquifer support. The static model of permeability used in DSS and Co-DSS is a regular grid composed by 4 162 308 cells, with 334x134x93 grid blocks of a single unit. In Figure 3 it is represented one simulated realization of the permeability field, the respective histogram and spatial variograms. The main direction is the (90,0) with 20 of amplitude and it´s isotropic. The spatial variability of permeability,– variograms, are fitted by an exponential model.

6

Rio Oil & Gas Expo and Conference 2010

Data count: 3630681 Mean: 42.6435 Variance: 7991.67 Maximum: 20620.4

Upper quartile: 30.476 Median: 5.40043 Lower quartile: 0.512283 Minimum: 0.01

Figure 3: Characterization of the static model of the permeability field used as reference for the study

4. Results and Discussion The stochastic model (example of Figure 3) reproduces the histogram, variograms and honors the experimental values of porosity and permeability at well data location. The entire set of producer wells defined 47 regions of perturbation (Figure 2), to compose the image to be used as secondary information in co-simulation procedure in each iteration. The geostatistical history matching described in 2 was applied and, in Figure 4, it is represented the evolution of the coefficient correlation (objective function) for each iteration of one well. In general, one obtained satisfactory results with the first or second iteration. The process was stopped in the fourth iteration as it have achieved a global good performance, i.e., approximately 51% of the wells have a coefficient of correlation bigger than 0.5.

Figure 4: Evolution of the correlation coefficient (from objective function) 7

Rio Oil & Gas Expo and Conference 2010

In Figure 5 the curves of well water cut of each iteration are represented. In the 4th iteration a good match is achieved.

WWCTH (well water cut total history) WWCT (well water cut total)_1st iteration

WWCT_2nd iteration WWCT_3rd iteration

WWCT_4th iteration

Figure 5 – Well Water Cut curves

The composite permeability cube, used as secondary information for the next generation of images, for each iteration is show in Figure 6.

Iteration1

Iteration2

Iteration3

Iteration4

Figure 6: Composed permeability models

Figure 7: One example simulation of the composed permeability model 8

Rio Oil & Gas Expo and Conference 2010

The average cube of the 10 realizations of the composed permeability model is shown on Figure 8, the variance cube on Figure 9 and the correlation coefficient cube oh the 4th iteration on Figure 10.

Figure 8: Average cube of the realizations of the composed permeability model

Figure 9: Variance cube of the realizations of the composed permeability model

Figure 10: Correlation coefficient cube of the 4th iteration

Analysing the images of figures 7, 8, 9 and 10 one can conclude that areas of greater variance are areas of lower correlation, i.e., more uncertainty, and smaller areas of well´s influence are generally areas of higher correlation, i.e., less variance (less uncertainty). Generally, areas where the correlations are higher, (variance is smaller) are areas where the permeability also is lower. There are areas too large that may cause an increase in the variability of the popular in this area.

5. Final Remarks An extended version of the geostatistical history matching (Mata Lima et al, 2007) was presented with another alternative for multi-criteria objective function representation, by using a Modified Hausdorff Distance. This showed very promising results in a very complex case study of a real reservoir. The presented methodology showed to have high potential for the history matching of great number of producer wells in very complex and heterogeneous 9

Rio Oil & Gas Expo and Conference 2010 reservoirs. Other improvements can be considered in next developments such as the local regions of influence of each well can be defined from a prior geologic concept (Hoffman and Caers, 2004), drainage area or expert opinion; the influence of each history matching factor (water-cut, pressure) in the objective function can be local, i.e., characterized by different local weights.

7. Acknowledgments Partex Oil and Gas is greatly acknowledged for supporting this research

8. References AANONSEN, S.I., EYDINOV, D. A multiscale method for distributed parameter estimation with application to reservoir history matching. Computational Geosciences, vol. 10, pp. 97-117, 2006. ABACIOGLU, Y., OLIVER, D. & REYNOLDS, A. Efficient reservoir history matching using subspace vectors. Computational Geosciences, vol. 5, pp. 151-172, 2005. CAERS, J. History matching under training-image-based geological model constrains. SPEJ, 2003. CHRISTIE, M., DEMYANOV, V. & ERBAS, D. Uncertainty quantification for porous media flows. Journal of Computational Physics, vol. 217, pp 143-158, 2006. DUBUISSON, M.P & JAIN, A.K. A modified Hausdorff distance for object matching, Proc. International conference on Pattern Recognition, Israel, 1994. HOFFMAN, B.T. & CAERS, J. History Matching with the Probability Perturbation Method: Application to a North Sea Reservoir. 9th European Conference on the Mathematics of Oil Recovery, Cannes, France, 2004. HOFFMAN, B.T. & CAERS, J. Regional Probability Perturbations for History Matching. Journal of Petroleum Science and Engineering, vol. 46, No. 1-2, pp.53-71, 2005. HU, L., Blanc, G. & NOETINGER, B. Gradual deformation and iterative calibration of sequential stochastic Simulations, Mathematical Geology, vol 33, no 4, pp.475-489, 2001. ISAAKS, Edward & SRIVASTAVA, R.Mohan. An Introduction to Applied Geostatistics, Oxford University Press, New York, Oxford, 1989. MATA-LIMA, H. Geostatistic in Reservoir Characterization: from estimation to Simulation methods. Estudios Geológicos-Madrid, vol. 61, No. 3-4, pp. 135-146, 2005. MATA-LIMA, H. Reducing Uncertainty in Reservoir Moddeling through efficient History Matching. Oil Gas European Magazine, 3/2008. MATA-LIMA, H.& SOARES, A. (2007). Geostatistical History Matching with Direct Sequential Transformation of Images, Petroleum Geostatistics, Cascais-Portugal, 9/2007 MATA-LIMA, H.& SOARES, A. Geostatistical History Matching with Direct Sequential Transformation of Images. Petroleum Geostatistics, Cascais-Portugal, 9/2007. MATA-LIMA, J. Reservoir Characterization with iterative direct sequential co-simulation: Integrating fluid dynamic data into stochastic model. Journal of Petroleum Science and Engineering, 2008. SCHEIDT, C. & CAERS, J. Representing spatial uncertainty using distances and kernels. Mathematical Geosciences, vol. 41, pp 397-419, 2009. SOARES, A. Direct Sequential Simulation and Co-simulation. Mathematical Geology, vol. 33, No. 8, pp.911-926, 2001. SOARES, A., ALMEIDA, J. & GUERREIRO, L. Incorporating Secondary Information using Direct Sequential CoSimulation in Stochastic Modelling II, American Association of Petroleum Geologists Pub, 2005. SOARES, A.O. Direct Sequential Simulation and Co-simulation. Mathematical Geology, vol. 33, No. 8, pp.911-926, 2001. SOARES, A.O. Geoestatística para as Ciências da Terra e do Ambiente. 2ª Edição, Colecção Ensino da Ciência e da Tecnologia, IST Press, Lisboa, 2006. SUZUKI, S. & Caers, J. (2008)-A Distance-based prior model parameterization for constraining solutions of spatial inverse problems. Mathematical Geoscience, Springer, 2008. TAVASSOLI, Z. CARTER, J.N. & KING, P.R. An analysis of history matching errors. Computational Geosciences, vol. 9, pp. 99-123, 2005. WANG, Y. & KOVSCEK, A.R. A streamline approach for History-Matching Production Data. Society of Petroleum Engineers, Oklahoma, 2000.

10

Resumo Este artigo apresenta uma nova metodologia de history matching geoestatístico, que visa integrar os dados dinâmicos dados históricos de produção de vários poços produtores - num modelo estocástico de propriedades internas. Este é um algoritmo interativo que pretende resolver o problema típico dos procedimentos inversos de history matching. A ideia baseia-se numa abordagem de perturbação global de um campo de permeabilidade, em que seja encontrada a calibração desejada através de uma função objectivo. Este trabalho está baseado na abordagem de transformação directa de imagens de permeabilidade (Mata-Lima et al, 2007), mas com um novo método de perturbaçção global, que apresenta uma convergência muito mais rápida à função objetivo. A perturbação (co-simulação directa) é realizada localmente, tendo por base as regiões de influência de cada poço e a função objectivo é multi-critério, combinação linear das diferenças entre os dados simulados de water cut, pressão e produção de óleo. O algoritmo proposto foi implementado num reservatório real no Oriente Médio.

Abstract This paper presents a new methodology for geostatistical history matching that aims to integrate the dynamic data historical production data of several producer wells - into the stochastic model of internal properties. This is an interactive algorithm to solve the typical inverse problem of history matching procedure. It is based on a global perturbation approach of, for example, a permeability field until the desired match, with a predefined objective function, is achieved. This paper is based on the direct transformation of images approach (Mata-Lima et al, 2007) but with a new method of global perturbations, with a much faster convergence to the objective function. The proposed algorithm was implemented in a real Middle East Reservoir. The perturbation (direct co-simulation) is performed locally, based on the influence regions of each well. The objective function is a linear combination of differences between simulated and real data of water-cut, pressure and oil production.

1. Introduction One of the main concerns in petroleum reservoir characterization field is to reduce the uncertainty in reservoir internal characteristics which can be achieved by reducing uncertainty in reservoir modeling through an efficient history matching. History matching is normally an inverse procedure in which an objective function (OF) is iteratively minimized as much as possible. The main difficulties of the problem – given the high number of independent variables this is a typical ill posed problem and there is an highly non-linear relation between the solution and the independent variables- still remain as challenge in this field of work. The proposed methodology is based on a global perturbation approach in the iterative path to reach the objective function. The seminal idea of these global perturbation algorithms comes from Hu et al (2001), with gradual deformation of Gaussian fields. The main advantage of this approach over the

______________________________ 1 Geologic and Mining Engineer – INSTITUTO SUPERIOR TÉCNICO 2 Professor – INSTITUTO SUPERIOR TÉCNICO 3 Student of Physics Engineering – INSTITUTO SUPERIOR TÉCNICO 4 Reservoir Engineer – PARTEX OIL & GAS 5 Reservoir Engineer – PARTEX OIL & GAS

Rio Oil & Gas Expo and Conference 2010 traditional methodologies of history matching is the simplicity of implementation. Other algorithms explored different methodologies for the global perturbation step of the inverse approach: a probability perturbation method (Caers, 2003), the use of direct sequential simulation and co-simulation to transform images of permeability (Mata-Lima et al, 2007). All these approaches follow, basically, the sequence of 4 steps: i) Generation of images of internal properties reflecting the spatial patterns of a geological conceptual model (main geological features, variograms, reference images); ii) Calculation of dynamic responses (dynamic simulator) of such set of images, given the existing wells (producers, injectors); iii) Evaluate the objective function based on a mismatch between the simulated responses and the real dynamic data; iv) Perturbation of initial images, based on the convergence to the objective function and return to i) until a desired match of the objective function is reached.

2. Proposed algorithm for Geostatistical History Matching The proposed algorithm aims to characterize a petroleum reservoir by an iterative algorithm for inverse modeling which integrates dynamic (production) data into stochastic models. It uses direct sequential co-simulation (coDSS) as a global perturbation method of permeability fields, in order to reach the convergence to the objective function. The flowchart of the iterative algorithm is presented in Figure 1 and can be summarized in the following sequence: i) Simulation of a set of realizations of permeability, honoring the well data, seismic information, geological constraints and main spatial patterns as revealed by the variograms; ii) Evaluation of dynamic responses of producing wells for each image/realization; iii) Calculation of the match between simulated dynamic responses and the real production data, through an objective function. If the mismatch is sufficiently low, the process stops. Otherwise a new generation of images is created (step iv) ; iv) Generation of a new set of images. In this step it is proposed a new method of transforming the images of permeability, such as they match of the dynamic response with the real production data but maintaining the main spatial patterns. Firstly at each iteration step, one “best” image is “composed” by a linear combination of the images corresponding to the closest dynamic responses. For this purpose the distances between the responses and from the responses to the real dynamic data (target) are calculated with Modified Hausdorff Distances (MHD), defining a 2-D referential where the simulated images with the closest responses to the target can be represented. This “composed” image does not have the spatial continuity of original images and it is just used as secondary information for the next images to be generated. This is performed with direct sequential co-simulation algorithm that accounts for the conditioning data – the existing permeability well data and the secondary information, of previous iteration, which influence is determined by a correlation measure of the match obtained in the previous step. The process returns to step ii) for the evaluation of the dynamic responses of the new images, and continues until the desired match is achieved through the minimization of the objective function.

2

Rio Oil & Gas Expo and Conference 2010 1st iteration Soft data (permeability model)

Co-DSS

Hard data (wells)

DSS

Hard data (wells)

Global correlation=0.6

Stochastic realizations Composite correlation coefficient cube (correlation data)

Stochastic realizations

Dynamic Simulator ynamic Simulator (Eclipse) Composite permeability image

Equivalent dynamic realizations

Objective function

No

Equivalent dynamic realizations

Real production data

Best responses to OF Calculation of the local correlation coefficients

Eclipse

Objective function

Desirable minimum?

Yes Composite “best” image of permabilities

Best responses to OF

Stop

Figure 1: Procedure of the proposed algorithm

2.1. Global transformation of Images Let us consider Z1(x) the initial image that one wish to transform in another Z2(x) with the same spatial pattern statistics, e.g. variogram and global histogram: C1(h) , γ1(h) , Fz1(z). Applying the Direct Sequential Co-Simulation algorithm, in each location x0 of a regular grid node, the local mean and variance are calculated by simple collocated co-kriging and estimation variance (equation 1).

Z 2 ( x0 ) * − m2 ( x0 ) = ∑ λα ( x0 )[Z 2 ( xα ) − m2 ( xα ) ] + λ1 ( x0 )[Z1 ( x0 ) − m1 ( x0 ) ]

(1)

α

Since the models γ1(h) and γ2(h) are the same, Markov approximation is quite acceptable to this case, i.e., the corregionalization model is fully defined with the correlation coefficient (correlogram at h=0) ρ1,2(0) between Z1(x) and Z2(x), such is in equation 2. ρ1,2(h) = ρ1,2(0) ρ2,2(h)

(2) 3

Rio Oil & Gas Expo and Conference 2010

The collocated cokriging system is written with the univariate correlogram of variable Z2(x), ρ2,2(h) and the coefficient correlation ρ1,2(0) (equation 3).

⎧ Nα ⎪∑ λβ 2 ( x 0 ).ρ 2 , 2 (xα , x β ) + λ1 ( x0 ).ρ 2 , 2 ( xα , x0 ).ρ1, 2 (0 ) = ρ 2 , 2 ( xα , x0 ) ⎪ β =1 ⎨ Nα ⎪ λ ( x 0 ).ρ (x , x )ρ (0 ) + λ ( x ).ρ (0 ) = ρ (0 ) 2, 2 ~0 1, 2 1 0 2, 2 1, 2 β2 β ⎪⎩∑ β =1

(3)

From the initial image Z1(x) one can obtain N images of Z2(x) with the same statistics of spatial continuity of Z1(x). Identically, from Z2(x), a set of Z3(x) could be obtained and so on. The degree of transformation – semblance of both images - is only function of the correlation coefficient ρ1,2(0)=[-1,+1] : the transformed image Z2(x) is as much identical to original Z1(x) as ρ1,2(0) is close to 1; An inverse image of original Z1(x) is obtained as ρ1,2(0) tends to –1. 2.2. Transformation with multiple images Let us consider that one wish to obtain a transformed image Zt(x) based on a multiple set of Ni multiple images Z1(x), Z2(x), …, ZNi(x). The formalism of equation 1 can now be generalized for the direct co-simulation of Zt(x) having Z1(x), Z2(x),…ZNi(x) as auxiliary variables. The collocated cokriging estimator of Zt(x) is presented on equation 4. Ni

Z t ( x0 ) * − mt ( x0 ) = ∑ λα ( x0 )[Z t ( xα ) − mt ( xα ) ] + ∑ λi ( x0 )[Z i ( x0 ) − mi ( x0 ) ] α

(4)

i =1

Since the models γi(h), i=1, Ni, and γt(h) are the same, the application of Markov approximation is also in this case quite appropriated. i.e., the corregionalization models are totally defined with the correlation coefficients ρt,i(0) between Zt(x) and Zi(x). The affinity of the transformed image Zt(x) with the multiple images Zi(x) are determined by the correlation coefficients ρt,i(0). Hence, one can select the images which characteristics we wish to “preserve” in the transformed image Zt(x). 2.3 Local Screening Effect Approximation Let us assume that to estimate Zt(x0) (equation 4) the collocated value Zi(x0) of a specific image Zi(x), with the highest correlation coefficient ρt,i(0), screens out the influence of the effect of remaining collocated values Zj(x0), j ≠ i. Hence, equation (4) can be written with just one auxiliary variable: the “best” at location x0 (equation 5).

Z t ( x0 ) * − mt ( x0 ) = ∑ λα ( x0 )[Z t ( xα ) − mt ( xα ) ] + λi ( x0 )[Z i ( x0 ) − mi ( x0 ) ]

(5)

α

With the local screening effect, Ni images Zi(x) give rise to just one auxiliary variable. The Ni images are merged in on single image based on the local correlation coefficient criterion. 2.4. Stochastic Inversion: a new approach This new approach of stochastic inversion is based on two key ideas: i) Use the sequential direct co-simulation with the screening effect approximation as the method of transforming images which converge to an objective function. ii) To follow the sequential procedure of genetic algorithm optimization to converge the transformed images towards an objective function. 2.5. Objective function In the next step one intend to evaluate the dynamic response of each permeability field realization obtained in a given iteration. This is accomplished by running a deterministic flow simulator (Eclipse 100 Black Oil, Schlumberger) 4

Rio Oil & Gas Expo and Conference 2010 for each permeability model as input, as well as the locations of different producers and injectors wells. The match evaluation is done by comparing the simulated and historic values of well water cut total (WWCT) and well bottom hole pressure (WBHP). After each response is obtained, a “distance” to the real production data is evaluated. In this study one use the Modified Hausdorff Distance for that purpose. The Hausdorff distance is often used to solve problems of "shape matching". Can be applied to any object in 2D or 3D if it can be represented in point set. The Modified Hausdorff distance (MHD) based on the classical Hausdorff distance, by Dubuisson and Jain (1994), is a measure of similarity used for the purpose of object matching, and has better performance because it is easy to compute and does not require the solution of the corresponding problem. The MHD as high discriminatory power - it is robust to outliers that may have resulted from targeting errors and his value should increase with the increasing of the difference between two objects. Suzuki and Caers (2008), proposed the use of MHD to measure the distance of two points in a referential characterized by the dynamic responses and the similarity of the sand channels patterns of each image. Suppose two sets A and B, and ai and bi are any points that belong to the sets A and B, respectively. Being uai ubi and the locations of ai and bi, then the spatial distance between ai and bi can be defined as the norm of the vector ||uai-ubi|| or Euclidean distance between the locations of the points ai and bi. Using this definition, the spatial distance measured between ai and the set B, d(ai, B) is determined as shown on equation 6.

d (a i , B ) = min bi ∈B u ai − u bi

(6)

Thus, the distance between A and B is defined as:

d ( A, B ) =

1 Na

∑

ai ∈A

d (u ai , B )

(7)

f (d ( A, B ), d ( B, A) = max(d ( A, B ), d ( B, A))) A = a1 ,..., a N a

where:

{ } B = {b ,..., b } 1

Nb

Now let us consider the response values of WWCTli (Well Water Cut Total) of simulated realization l along the time periods ti=1, N and the equivalent real data production WWCT at a given well. Modified Hausdorff Distance between the two responses is shown on equation 8.

MHD 1 wwct =

1 N

N

1

∑ d (WWCT i − WWCT )

(8)

i =1

Identical distance is obtained with WBHP (Well Bottom Hole Pressure). The objective function (MHD)total normalizes the distances of the two variables, which are in completely different scales (equation 9). MHD is the value of the modified Hausdorff distance calculated for each simulation of a certain well MHD is the average value of the modified Hausdorff distance and σ is the standard deviation of the modified Hausdorff distance calculated.

MHDtotal =

⎛ MHDwwct − MHDwwct ⎜ ⎜ σ wwct ⎝

⎞ ⎛ MHDwbhp − MHDwbhp ⎞ ⎟ ⎟−⎜ ⎟ ⎜ ⎟ σ wbhp ⎠ ⎝ ⎠ 2

(9)

In order to perform the co-simulation (with the collocated co-kriging approach) one need to convert the MHD into a correlation coefficient value. Hence, to the highest distance value we assume the correlation coefficient CC = 0.

5

Rio Oil & Gas Expo and Conference 2010 The equation 10 is used to access to the equivalent measurent of correlation. Dhi is the smallest MHD of the considered group of simulations and Dhn is the MHD for the current simulation. The ratio R will be variable because the number of considered simulations increases as the number of iterations performed. n

R =

CC= 1-R

Dh i Dh

(10)

In case of considering a two variables objective function (WBHP_ pressure and WWCT_ water cut), the final correlation coefficient is given by the mean of the correlation values corresponding to the simulation that have a smaller objective function. 2.6. Composing “secondary” permeability image and evaluation of the correlation coefficient´s cube The final step of each iteration consists on the composing of the secondary image for the next set of stochastic realizations of the iterative process - co-DSS. The image corresponding to the best response of each well is selected to compose a permeability image. The “best” image is composed of the parts of each simulation corresponding to the best match, inside a grid of polygonal areas of influence of each well (Figure 2).

Layer 55

Figure 2: Representation of the polygonal regions and the composite image of permeabilities

The correlation coefficients are extrapolated inside each region, representing the correlation between the real production data and the simulated ones.

3. Reservoir Description The reservoir used to the validation of the proposed algorithm is a complex and very heterogeneous carbonate reservoir. The field is currently producing under natural depletion with aquifer support. The static model of permeability used in DSS and Co-DSS is a regular grid composed by 4 162 308 cells, with 334x134x93 grid blocks of a single unit. In Figure 3 it is represented one simulated realization of the permeability field, the respective histogram and spatial variograms. The main direction is the (90,0) with 20 of amplitude and it´s isotropic. The spatial variability of permeability,– variograms, are fitted by an exponential model.

6

Rio Oil & Gas Expo and Conference 2010

Data count: 3630681 Mean: 42.6435 Variance: 7991.67 Maximum: 20620.4

Upper quartile: 30.476 Median: 5.40043 Lower quartile: 0.512283 Minimum: 0.01

Figure 3: Characterization of the static model of the permeability field used as reference for the study

4. Results and Discussion The stochastic model (example of Figure 3) reproduces the histogram, variograms and honors the experimental values of porosity and permeability at well data location. The entire set of producer wells defined 47 regions of perturbation (Figure 2), to compose the image to be used as secondary information in co-simulation procedure in each iteration. The geostatistical history matching described in 2 was applied and, in Figure 4, it is represented the evolution of the coefficient correlation (objective function) for each iteration of one well. In general, one obtained satisfactory results with the first or second iteration. The process was stopped in the fourth iteration as it have achieved a global good performance, i.e., approximately 51% of the wells have a coefficient of correlation bigger than 0.5.

Figure 4: Evolution of the correlation coefficient (from objective function) 7

Rio Oil & Gas Expo and Conference 2010

In Figure 5 the curves of well water cut of each iteration are represented. In the 4th iteration a good match is achieved.

WWCTH (well water cut total history) WWCT (well water cut total)_1st iteration

WWCT_2nd iteration WWCT_3rd iteration

WWCT_4th iteration

Figure 5 – Well Water Cut curves

The composite permeability cube, used as secondary information for the next generation of images, for each iteration is show in Figure 6.

Iteration1

Iteration2

Iteration3

Iteration4

Figure 6: Composed permeability models

Figure 7: One example simulation of the composed permeability model 8

Rio Oil & Gas Expo and Conference 2010

The average cube of the 10 realizations of the composed permeability model is shown on Figure 8, the variance cube on Figure 9 and the correlation coefficient cube oh the 4th iteration on Figure 10.

Figure 8: Average cube of the realizations of the composed permeability model

Figure 9: Variance cube of the realizations of the composed permeability model

Figure 10: Correlation coefficient cube of the 4th iteration

Analysing the images of figures 7, 8, 9 and 10 one can conclude that areas of greater variance are areas of lower correlation, i.e., more uncertainty, and smaller areas of well´s influence are generally areas of higher correlation, i.e., less variance (less uncertainty). Generally, areas where the correlations are higher, (variance is smaller) are areas where the permeability also is lower. There are areas too large that may cause an increase in the variability of the popular in this area.

5. Final Remarks An extended version of the geostatistical history matching (Mata Lima et al, 2007) was presented with another alternative for multi-criteria objective function representation, by using a Modified Hausdorff Distance. This showed very promising results in a very complex case study of a real reservoir. The presented methodology showed to have high potential for the history matching of great number of producer wells in very complex and heterogeneous 9

Rio Oil & Gas Expo and Conference 2010 reservoirs. Other improvements can be considered in next developments such as the local regions of influence of each well can be defined from a prior geologic concept (Hoffman and Caers, 2004), drainage area or expert opinion; the influence of each history matching factor (water-cut, pressure) in the objective function can be local, i.e., characterized by different local weights.

7. Acknowledgments Partex Oil and Gas is greatly acknowledged for supporting this research

8. References AANONSEN, S.I., EYDINOV, D. A multiscale method for distributed parameter estimation with application to reservoir history matching. Computational Geosciences, vol. 10, pp. 97-117, 2006. ABACIOGLU, Y., OLIVER, D. & REYNOLDS, A. Efficient reservoir history matching using subspace vectors. Computational Geosciences, vol. 5, pp. 151-172, 2005. CAERS, J. History matching under training-image-based geological model constrains. SPEJ, 2003. CHRISTIE, M., DEMYANOV, V. & ERBAS, D. Uncertainty quantification for porous media flows. Journal of Computational Physics, vol. 217, pp 143-158, 2006. DUBUISSON, M.P & JAIN, A.K. A modified Hausdorff distance for object matching, Proc. International conference on Pattern Recognition, Israel, 1994. HOFFMAN, B.T. & CAERS, J. History Matching with the Probability Perturbation Method: Application to a North Sea Reservoir. 9th European Conference on the Mathematics of Oil Recovery, Cannes, France, 2004. HOFFMAN, B.T. & CAERS, J. Regional Probability Perturbations for History Matching. Journal of Petroleum Science and Engineering, vol. 46, No. 1-2, pp.53-71, 2005. HU, L., Blanc, G. & NOETINGER, B. Gradual deformation and iterative calibration of sequential stochastic Simulations, Mathematical Geology, vol 33, no 4, pp.475-489, 2001. ISAAKS, Edward & SRIVASTAVA, R.Mohan. An Introduction to Applied Geostatistics, Oxford University Press, New York, Oxford, 1989. MATA-LIMA, H. Geostatistic in Reservoir Characterization: from estimation to Simulation methods. Estudios Geológicos-Madrid, vol. 61, No. 3-4, pp. 135-146, 2005. MATA-LIMA, H. Reducing Uncertainty in Reservoir Moddeling through efficient History Matching. Oil Gas European Magazine, 3/2008. MATA-LIMA, H.& SOARES, A. (2007). Geostatistical History Matching with Direct Sequential Transformation of Images, Petroleum Geostatistics, Cascais-Portugal, 9/2007 MATA-LIMA, H.& SOARES, A. Geostatistical History Matching with Direct Sequential Transformation of Images. Petroleum Geostatistics, Cascais-Portugal, 9/2007. MATA-LIMA, J. Reservoir Characterization with iterative direct sequential co-simulation: Integrating fluid dynamic data into stochastic model. Journal of Petroleum Science and Engineering, 2008. SCHEIDT, C. & CAERS, J. Representing spatial uncertainty using distances and kernels. Mathematical Geosciences, vol. 41, pp 397-419, 2009. SOARES, A. Direct Sequential Simulation and Co-simulation. Mathematical Geology, vol. 33, No. 8, pp.911-926, 2001. SOARES, A., ALMEIDA, J. & GUERREIRO, L. Incorporating Secondary Information using Direct Sequential CoSimulation in Stochastic Modelling II, American Association of Petroleum Geologists Pub, 2005. SOARES, A.O. Direct Sequential Simulation and Co-simulation. Mathematical Geology, vol. 33, No. 8, pp.911-926, 2001. SOARES, A.O. Geoestatística para as Ciências da Terra e do Ambiente. 2ª Edição, Colecção Ensino da Ciência e da Tecnologia, IST Press, Lisboa, 2006. SUZUKI, S. & Caers, J. (2008)-A Distance-based prior model parameterization for constraining solutions of spatial inverse problems. Mathematical Geoscience, Springer, 2008. TAVASSOLI, Z. CARTER, J.N. & KING, P.R. An analysis of history matching errors. Computational Geosciences, vol. 9, pp. 99-123, 2005. WANG, Y. & KOVSCEK, A.R. A streamline approach for History-Matching Production Data. Society of Petroleum Engineers, Oklahoma, 2000.

10