**3. Seismic inversion**

The seismic data is the best source of spatially extensive measurement over the reservoir, but as explained previously, is very coarse vertically. From the wells we can use the acoustic impedance (AI) which is the result of multiplying the lateral interval velocity by the layer density, that later can be transformed in to seismic amplitude, this is possible because the difference between the different geological materials is linked to the response of the seismic signal.

The advantages of working with AI instead of the recorded seismic data are that is layer property rather than an interface property and hence more like geology and therefore it has a more physical meaning and during the inversion process all the well data is tied to the seismic data giving better understanding of the quality of both datasets.

So, the integration of the amplitude seismic data with the acoustic impedance data from wells requires some kind of transformation. There are two methods that transform one in to the other, an inverse method and a direct or forward method.

Both of them use wavelets for the transformation. These can be summarized as onedimensional pulse and the link between seismic data and geology and they are a kind of impulse of energy that is created when a marine air gun or land dynamite source is released during the acquisition of seismic surveys.

The first method, "inverse process", transforms the amplitudes from seismic to acoustic impedance by removing the wavelet from the amplitudes and obtaining a model of the acoustic impedance. The problem of simultaneously inverting reservoir engineering and seismic data to estimate porosity and permeability involves complex processes such as fluid flow through porous media, and acoustic wave propagation and cannot be solved by linear inversion methods as [11]. On the other hand, the process of sending a wavelet into the earth and measure the reflection is called forward process. Here we need to know what it is in the subsurface in terms of geological models divided in layers and each characterized by his acoustic impedance. Forward modeling combines the sequence of differences in the acoustic impedances with a seismic pulse to obtain a synthetic amplitude trace.

In the forward model if we compare those sequences of synthetic amplitude variations with the real ones, we can identify and quantify the differences and try to minimize them.

#### **3.1. Synthetic seismogram**

122 New Technologies in the Oil and Gas Industry

predefined threshold.

At each iterative step one knows how close is one given generated image from the objective, by the global and local correlation coefficients between the transformed traces and the real seismic traces. These correlation coefficients of different simulated images are used as the affinity criterion to create the next generation of images until it converges to a given

In a last step, porosity images can be derived from the seismic impedances obtained by seismic inversion and the uncertainty derived from the seismic quality is assessed based on

For the case of characterization of the reservoir in terms of facies distribution several methods for the integration of the seismic data in facies models have been proposed, several of which rely on the construction of a facies probability cube by calibration of the seismic data with wells. If only post-stack seismic data is considered, is typically inverts the seismic amplitudes into a 3D acoustic impedance cube, and then converts this impedance into a 3D

This facies probability can serve as input of several well-known geostatistical algorithms to create a facies realization, such as the use the cube as locally varying mean on indicator kriging or the use of the tau model in [10], in a multi-point simulation to integrate the facies

While this provides satisfactory results in most cases, the resulting facies realization does not necessarily match the original seismic amplitude from which the acoustic impedance was inverted. Indeed, if one would forward simulate, for example a 1D convolution on a single facies realizations, then this procedure does not guarantee that the forward simulated

Another objective of this work is to present a geostatistical methodology, based on multipoint technique that generates facies realization compatible with the field seismic amplitude data, and therefore this new procedure has two main advantages, matching field seismic amplitude data in a physical sense, not merely in a probabilistic sense and using

The seismic data is the best source of spatially extensive measurement over the reservoir, but as explained previously, is very coarse vertically. From the wells we can use the acoustic impedance (AI) which is the result of multiplying the lateral interval velocity by the layer density, that later can be transformed in to seismic amplitude, this is possible because the difference between the different geological materials is linked to the response of the seismic

The advantages of working with AI instead of the recorded seismic data are that is layer property rather than an interface property and hence more like geology and therefore it has

probability cube with spatial continuity information provided by a training image.

the quality of match between synthetic seismogram and real seismic.

facies probability using a calibration method of choice.

multi-point statistics, not just the two-point statistics (variogram).

seismic matches the field amplitudes.

**3. Seismic inversion** 

signal.

It is known that the propagation of waves or sounds that pass through earth can be explained by an elastic wave equation. One simple, but powerful and commonly used approach for computing the seismic response of a certain earth model is the so-called convolutional model as in [12-14], which can be derived from an acoustic approximation of an elastic equation.

This convolutional model (equation 1) creates a synthetic seismic trace Sy(t), that is the result of convolving a seismic wavelet w(t) with a time series of reflectivity coefficients r(t), with the addition of a noise component n(t), as follows:

$$\mathbf{S}\mathbf{y}(\mathbf{t}) = \mathbf{w}(\mathbf{t}) \,^\*\mathbf{r}(\mathbf{t}) + \mathbf{n}(\mathbf{t}) \tag{1}$$

An even simpler assumption is to consider the noise component to be zero (equation 2), in which case the synthetic seismic trace is simply the convolution of a seismic wavelet with the earth's reflectivity:

$$\mathbf{S}\mathbf{y}(\mathbf{t}) = \mathbf{w}(\mathbf{t}) \,^\*\mathbf{r}(\mathbf{t})\tag{2}$$

The reflection coefficient is one of the fundamental physical concepts in the seismic method. Basically, each reflection coefficient may be thought of as the response of the seismic wavelet to an acoustic impedance change within the earth and represents the percentage of the energy that is emitted compared to the one that is reflected

Mathematically, converting from acoustic impedance to reflectivity involves dividing the differences in the acoustic impedances by the sum of them. This gives the reflection coefficient at the boundary between the two layers.

Each layer can have different rock acoustic impedance, which is a function of two porositydependent rock properties: bulk density and velocity (equation 3). For instance, the P-wave acoustic impedance is given by:

$$\text{AI(t)} \newline \text{= density(t)\* P-Velocity(t)} \tag{3}$$

So the reflectivity coefficient (equation 4) can be computed from:

$$r\_t = \frac{\mathbf{A}\mathbf{I}\_{t+1} - \mathbf{A}\mathbf{I}\_t}{\mathbf{A}\mathbf{I}\_{t+1} + \mathbf{A}\mathbf{I}\_t} \tag{4}$$

where AI represents the acoustic impedance and the subscripts t and t+1 refer to two subsequent layers in the stratigraphic column.

The wavelet is usually extracted from the seismic survey through deconvolution (the methodology of extraction the wavelet from the seismic signal). During the deconvolution, the wavelet extraction does not take in account the low frequency of the earth model, this represents the compaction of the porous media in depth, which in some seismic inversion algorithm and in the inverse modeling, are ignored or only take in account adjacent geology units, not the full vertical analysis.

In this work, the problem is resolved by the algorithm of generation of the acoustic impedance model, consider the trend that are represented in the wells log data, by creating initially a model only conditioned by the wells data (AI).

#### **3.2. Correlation coefficient comparison**

Using the convolution algorithm described previously, one can forward an entire model (or cube if one considers a model as a cartesian grid) of AI created only by using the wells and its spatial continuity to a synthetic seismic amplitude cube that can be compared with the real amplitude seismic cube. This comparison will give the mismatch between those two cubes of amplitudes.

The mismatch between the synthetic and the field amplitudes is calculated as a simple colocated correlation coefficient as in [15-17]. Alternatively, this mismatch could be calculated in the form of a least-square difference between both amplitude cubes.

In a first step, each of the N synthetic seismic amplitude cubes are divided into layers (figure 1).The choice of number of layers is a tuning parameter of the algorithm.

The amount of mismatch between each pair of columns of the synthetic and the real amplitude is calculated and stored in a new cube (figure 1), named the mismatch cube. Essentially, this new cube provides information on which locations in the reservoir are fitting the seismic and which need improvement.

**Figure 1.** Calculus of the mismatch cube.

124 New Technologies in the Oil and Gas Industry

acoustic impedance is given by:

The reflection coefficient is one of the fundamental physical concepts in the seismic method. Basically, each reflection coefficient may be thought of as the response of the seismic wavelet to an acoustic impedance change within the earth and represents the percentage of

Mathematically, converting from acoustic impedance to reflectivity involves dividing the differences in the acoustic impedances by the sum of them. This gives the reflection

Each layer can have different rock acoustic impedance, which is a function of two porositydependent rock properties: bulk density and velocity (equation 3). For instance, the P-wave

> 1 1 AI AI AI AI *t t*

where AI represents the acoustic impedance and the subscripts t and t+1 refer to two

The wavelet is usually extracted from the seismic survey through deconvolution (the methodology of extraction the wavelet from the seismic signal). During the deconvolution, the wavelet extraction does not take in account the low frequency of the earth model, this represents the compaction of the porous media in depth, which in some seismic inversion algorithm and in the inverse modeling, are ignored or only take in account adjacent geology

In this work, the problem is resolved by the algorithm of generation of the acoustic impedance model, consider the trend that are represented in the wells log data, by creating

Using the convolution algorithm described previously, one can forward an entire model (or cube if one considers a model as a cartesian grid) of AI created only by using the wells and its spatial continuity to a synthetic seismic amplitude cube that can be compared with the real amplitude seismic cube. This comparison will give the mismatch between those two

The mismatch between the synthetic and the field amplitudes is calculated as a simple colocated correlation coefficient as in [15-17]. Alternatively, this mismatch could be calculated

in the form of a least-square difference between both amplitude cubes.

*t t*

*t*

*r* 

AI(t) = density(t) \* P-Velocity(t) (3)

(4)

the energy that is emitted compared to the one that is reflected

So the reflectivity coefficient (equation 4) can be computed from:

coefficient at the boundary between the two layers.

subsequent layers in the stratigraphic column.

initially a model only conditioned by the wells data (AI).

**3.2. Correlation coefficient comparison** 

units, not the full vertical analysis.

cubes of amplitudes.

A mismatch cube is calculated for each of the simulated models using their corresponding synthetic seismic models.

#### **3.3. Co-generation of acoustic images**

The generation of an acoustic impedance model is the main objective of this work, but until an optimized one is created, it is submitted to a convergence genetic modification.

These images are generated through stochastic simulation, i.e. generation of different models with the same probabilistic and spatial distribution. This means that if the parameters that are used to the creation of the models are too tight, the produced images are almost identical, if the parameters allow more freedom, these models are totally different.

Using this idea, the proposed methodology uses two approaches for image generation, two different programs are used. The first one is the DSS (Direct Sequential Simulation) for the generation of the first unconditioned images and then for the convergence is used the Co-DSS (Direct Sequential Co-Simulation) that uses a soft or secondary data to simulation conditioning as in [3,18,19].

The second algorithm is the Snesim (Single Normal Equation SIMulator), a multi point simulation algorithm described in [20,21]. This version is used to create the first unconditional facies models. A second version with modification in [2], uses the influence of secondary data to create the models that are conditioned to more information. Both these programs use secondary data for the co-simulations, and this secondary or soft data is the one that makes the algorithm converge.

For the use of Direct Sequential Co-Simulation for global transformation of images as in [17]**,**  let us consider that one wishes to obtain a transformed image Zt(x), based on a set of Ni images Z1(x), Z2(x),…ZNi(x), with the same spatial dispersion statistics (e.g. correlogram or variogram and global histogram): C1(h) , 1(h) , Fz1(z).

Direct co-simulation of Zt(x), having Z1(x), Z2(x),…ZNi(x) as auxiliary variables, can be applied as in [3]. The collocated cokriging estimator (equation 5) of Zt(x) becomes:

$$-Z\_t(\mathbf{x}\_0)^\ast - m\_t(\mathbf{x}\_0) = \sum\_a \lambda\_a \left(\mathbf{x}\_0\right) \left[Z\_t(\mathbf{x}a) - m\_t(\mathbf{x}\_a)\right] + \sum\_{i=1}^{Ni} \lambda\_i \left(\mathbf{x}\_0\right) \left[Z\_i(\mathbf{x}\_0) - m\_i(\mathbf{x}\_0)\right] \tag{5}$$

Since the models i(h), i=1, Ni, and t(h) are the same, the application of Markov approximation is, in this case, quite appropriated, i.e., the co-regionalization 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 with which characteristics we wish to "preserve" in the transformed image Zt(x). So, a local screening effect approximation can be done, by assuming that to estimate Zt(x0) (equation 5) 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 6 can be written with just one auxiliary variable: the "best" at location x0:

$$\sum\_{\alpha} Z\_t(\mathbf{x}\_0)^\ast - m\_t(\mathbf{x}\_0) = \sum\_{\alpha} \lambda\_\alpha(\mathbf{x}\_0) \left[ Z\_t(\mathbf{x}\alpha) - m\_t(\mathbf{x}\_\alpha) \right] + \lambda\_i(\mathbf{x}\_0) \left[ Z\_i(\mathbf{x}\_0) - m\_i(\mathbf{x}\_0) \right] \tag{6}$$

With the local screening effect, Ni images Zi(x) give rise to just one auxiliary variable. The Ni images are merged in on a single image based on the local correlation coefficient criterion.

Basically, the correlation or least mismatch of each of the previously simulated images at each x0 is converted to a single image where the best of each x0 is selected, this way.

In practice, the algorithm chooses in each x0 of the simulated acoustic impedance model the best genes and the correlation that those genes have. Based on figure 2, the process is explained as follows:

The process starts by analyzing each of the mismatch cubes that were created previously, and in each position of the cubes (x0) it compares which has the higher value of correlation or lower mismatch. In this case the example used is the correlation, so, the points 1a and 1b are compared, and since it is the 1a that has the higher value of correlation, that value is copied (1c) to the Best Correlation Cube.

At the same time, since it was from the first simulation that has better coefficient of correlation (CC) with the real seismic, another cube is created with the correspondent acoustic impedance (AI) value (1d). Next, it starts comparing the next set of values from the CC cubes, and the process is the same. With a comparison between the sets 2a and 2b, and since in this case it is the second cube that has the higher value of cc, it is this values that is copied to the Best CC cube (2c). The same is done from the AI cube of this simulation, copying the acoustic impedance values to the Best AI (2d). This process continues until all the N simulated models have been analyzed and to all the values of the cubes. Through this way two new cubes are created, one with the best genes from each acoustic impedance model (Best AI) generated previously and another cube with the confidence factor of each part (Best CC).

126 New Technologies in the Oil and Gas Industry

one that makes the algorithm converge.

variable: the "best" at location x0:

criterion.

explained as follows:

copied (1c) to the Best Correlation Cube.

variogram and global histogram): C1(h) , 1(h) , Fz1(z).

programs use secondary data for the co-simulations, and this secondary or soft data is the

For the use of Direct Sequential Co-Simulation for global transformation of images as in [17]**,**  let us consider that one wishes to obtain a transformed image Zt(x), based on a set of Ni images Z1(x), Z2(x),…ZNi(x), with the same spatial dispersion statistics (e.g. correlogram or

Direct co-simulation of Zt(x), having Z1(x), Z2(x),…ZNi(x) as auxiliary variables, can be

00 0 00 0

 

( )\* ( ) () () () ()

Since the models i(h), i=1, Ni, and t(h) are the same, the application of Markov approximation is, in this case, quite appropriated, i.e., the co-regionalization models are

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 with which characteristics we wish to "preserve" in the transformed image Zt(x). So, a local screening effect approximation can be done, by assuming that to estimate Zt(x0) (equation 5) 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 6 can be written with just one auxiliary

00 0 00 0 ( )\* ( ) ( ) () () () *Zx mx x Zx mx x Zx mx t t*

With the local screening effect, Ni images Zi(x) give rise to just one auxiliary variable. The Ni images are merged in on a single image based on the local correlation coefficient

Basically, the correlation or least mismatch of each of the previously simulated images at

In practice, the algorithm chooses in each x0 of the simulated acoustic impedance model the best genes and the correlation that those genes have. Based on figure 2, the process is

The process starts by analyzing each of the mismatch cubes that were created previously, and in each position of the cubes (x0) it compares which has the higher value of correlation or lower mismatch. In this case the example used is the correlation, so, the points 1a and 1b are compared, and since it is the 1a that has the higher value of correlation, that value is

each x0 is converted to a single image where the best of each x0 is selected, this way.

 

*t t ii i*

 (6)

*t t t t ii i*

*Zx mx x Zx mx x Zx mx*

1

 

*i*

 (5)

*Ni*

applied as in [3]. The collocated cokriging estimator (equation 5) of Zt(x) becomes:

totally defined with the correlation coefficients t,i(0) between Zt(x) and Zi(x).

**Figure 2.** Process of creating "Best Correlation Cube" and "Best Acoustic Impedance cube"

Since the process is an iterative one, in the very first steps of the iterative process, the secondary image (Best AI) do not have the spatial continuity pattern of the primary (simulated AI), as it results from a composition of different parts of a set of simulated images. As the process continues, the secondary image tends to have the same spatial pattern of generated images by co-simulation, because the correlation values are becoming higher and the freedom of the co-simulation is diminishing. Finally these two new data sets are used as soft data for the next iteration.

#### **4. Algorithm description using direct sequential simulation**

The proposed stochastic inversion algorithm is settled with the following key ideas in mind: generation of stochastic images, transformation of the images in synthetic seismograms, chose and keep the best "genes" from each of the images and then use them through the exercise of the genetics algorithm formalism and the stochastic co-simulation to create a new generation of images and the convergence of the process. It can be summarized in the following steps;


At each iterative step one knows how closer is one given generated image from the objective, by the global and local correlation coefficients between the synthetic seismogram and the real seismic. These correlation coefficients of different simulated images are used as the affinity criterion to create the next generation of images until it converges to a given predefined threshold. A simplified diagram is show in figure 3.

**Figure 3.** Diagram of the proposed algorithm

#### **4.1. Case study**

128 New Technologies in the Oil and Gas Industry

function is reached;

**Figure 3.** Diagram of the proposed algorithm

simulation. Instead of individual traces of cells;

impedance in internal reservoir characteristics.

predefined threshold. A simplified diagram is show in figure 3.

from acoustic impedances, with a known wavelet;

seismic by computing, for example local correlation coefficients;

i. Generate a set of initial 3D images of acoustic impedances by using direct sequential

ii. Create the synthetic seismogram of amplitudes, by convolving the reflectivity, derived

iii. Evaluate the match of the synthetic seismograms, of entire 3D image, and the real

iv. Ranking the "best" images based on the match (e.g. the average value or a percentile of correlation coefficients for the entire image). From them, one selects the best parts (the columns or the horizons with the best correlation coefficient) of each image. Compose

vi. The last step of the process is the transformation of the optimized cube of acoustic

At each iterative step one knows how closer is one given generated image from the objective, by the global and local correlation coefficients between the synthetic seismogram and the real seismic. These correlation coefficients of different simulated images are used as the affinity criterion to create the next generation of images until it converges to a given

one auxiliary image with the selected "best" parts, for the next simulation step; v. Generate a new set of images, by direct sequential co-simulation, using the best locals correlation as the local co-regionalization model and return to step ii) starting an iterative process that will end when the match between the synthetic seismogram and the real seismic is considered satisfactory or until a given threshold of the objective A case study of a Middle East reservoir will be presented but only a small part of the full reservoir is studied due to data confidentiality and the coordinates presented are modified.

The field described in [22,23] is a carbonate reservoir with a deposition geology that holds in some zones a strong internal geometry with clinoforms. These are sedimentary deposits that have a sigmoidal or S shape. They can range in size from meter, like sand dunes, to kilometers and can grow horizontally in response to sediment supply and physical limits on sediment accumulation.

This type of geological phenomenon increases the complexity of the internal reservoir, making the geophysics interpretation of seismic a very difficult job, mainly when the resolution of acquisition is very coarse. It also causes a great impact in oil production in wells and reservoir characterization and modeling.

From the calibration of the acoustic impedance data of the wells with the seismic, a wavelet was extracted and there were 19 wells in the study area but only two were used, because only these two have the velocity log, and without the velocity log the acoustic impedance data cannot be calculated. On those wells the acoustic impedance was determined, then calibrated with the seismic and finally upscaled to fit the seismic scale of 4 ms.

The convolution for well 11 is show in figure 4.

**Figure 4.** Calculation of the AI, convolution and match with seismic of well 11.

As one can notice, there is a very good match with the synthetic amplitude (red line on right side) calculated with the acoustic impedance (white line in the middle) from the wells and the real seismic amplitude extracted from the seismic cube (cyan line in the right side).
