Research and application of seismic porosity inversion method for carbonate reservoir based on Gassmann’s equation

2.1 Eshelby–Walsh ellipsoidal inclusion theory

The ellipsoidal inclusion theory for dry rocks [28] simulated pore as an elastic ellipsoidal inclusion in rocks and is supposed that a region within an isotropic elastic solid undergoes change of form. In particular, if the region is an ellipsoid, the strain inside it is uniform and may be expressed in terms of tabulated elliptic integrals. An ellipsoidal region in an infinite medium has elastic constants different from those of the rest of the material. And a fundamental form of rock matrix parameters expression based on the Eshelby’s theory is as given below [29].

where β D and β s are compressibility of dry rock and rock matrix, respectively; η is rock porosity; μ D and μ s are shear modulus of rock and rock matrix, respectively; m, n, and α are nondimensional positive constants related with the shape and property of rocks. α is the aspect ratio, the expression for m and n is described below, and K S and μ S are the matrix elastic parameters.

And there is a relation table between Possion’s ratio and m and n as shown in Table 1 when the Possion’s ratio is at the range of 0.15 ≤ σ ≤ 0.40 .

2.2 Gassmann’s equation and the deduction of porosity inversion method

The most famous Gassmann’s equation [31] using rock compressibility coefficients is expressed as follows.

where β E is the rock saturated fluid equivalent compression coefficient, β s is the rock matrix compression coefficient, β D is the dry rock compression coefficient, β p is the rock pore fluid compression coefficient, and η is the rock porosity.

Gassmann’s equation is suitable for processing seismic data because of its low-frequency assumption. It can be used for fluid substitution as well as porosity inversion, but they have different purpose and calculation flow.

In case of diphasic medium, compressibility of pore fluid is much bigger than that of rock matrix, especially in the deep area where the carbonate reservoir has high pressure and small compressibility (Figure 1), and so we can take β p approximate to β p − β s , which has approved high precision to calculate the rock matrix compression coefficient [32], then we can get the expression from (3).

Figure 1

Comparison of the rock compression coefficient in different conditions [30].

Figure 2

Rock slice photo by microscope, with ellipsoidal pore shape.

Figure 3

The predicting porosity working flow.

Then replacement of β s with β D from the expression (1), we get the following expression.

Thus, we get the expression to calculate rock porosity:

where C = α / m is called the rock pore geometry coefficient.

2.3 Calculation of parameters and working flow of the porosity inversion

From the rock porosity expression (6), we can see it needs to calculate some parameters such as β s , β p , β E , and C, which can be gained from the log data certainly to calculate the porosity.

1. β s , the rock matrix compression coefficient, can be calculated in 2 ways that mentioned as below, which means we only calculate one point value to instead the value of the whole well with depth [33].

   A linear equation can be transformed from expression (6) as follows.

   (7) η = β E 1 β p + C β s − β s β p + C = A β E − B

   Expression (7) is the linear relation between the rock saturated fluid equivalent compression coefficient and the rock porosity. Where

   (8) A = 1 β p + C β s B = − β s 1 β p + C β s

   Then β s = − B / A can be deduced from the expression (8).

   However, there is a small problem with the calculating method as upper described. With the depth going down, the lithology is changing in the actual situation, so the linear method of calculating the one-point value instead of a string of values is not suitable for the actual situation. So we employed a method [34] to calculate a string of values corresponding to the logs of wells, which is more suitable to invert porosity. The method is stated as follows.

   We take f 1 as the liquid part of Gassmann’s equation

   f 1 = β 2 M = 1 − K dry K S η K f l + 1 − η K S − K dry K S 2

   And which can also be expressed as a liquid part of Russell’s expression

   f 2 = Z p 2 − c Z s 2 ρ sat

   where Z p and Z s are p-wave and s-wave impedance, respectively, c is a constant, when f 1 − f 2 < B , B is an initial small value, we can get a string of β s values.

2. β p , the rock pore fluid compression coefficient, can be calculated by the following expression and with the saturated water log data [35], where H is the depth.

   (9) β oil = 1 1.19 − 0.362 H + 0.042 H 2 β gas = 1 0.00014 + 0.00946 H + 0.00145 H 2 β w ater = 1 2.02 + 0.304 H − 0.0572 H 2

3. β E , the rock saturated fluid equivalent compression coefficient, can be calculated by elastic parameters as v p , v s , ρ inversed by prestacked seismic data, and the calculation expression is as described below.

   (10) β E = 1 ρ v p 2 − 4 3 v s 2

4. C, the rock pore geometry parameter, is generally calculated by the expression C = α / m , and C is always constant, but in reality C must be verified with the change in the lithology and pore geometry type. We calculate C by expression (6) and the log data without calculating m and α , thus C is more accurate to inverse porosity than when it is a constant. The calculating expression is as described below.

   (11) C = β s ⋅ η β p − β E + β s ( β E − β s ) β p

5. Aspect ratio and the final form of porosity calculating expression.

   In general, we do statistical analysis by observing the zoomed rock slice under microscope to take an average of the aspect ratio α = b / a in different wells, where the pore shape shows as an ellipse (Figure 2), and then calculate m based on the relation between m and Possion’s ratio (see expression 2), finally we get the value of C. From these works, we have known that the range of m is [0.54, 1.2] from the result of rock physics testing analysis. To carbonate reservoir, there is no prospecting meaning as the aspect ratio is smaller than 0.072, which means the pore is approaching a tiny micro-fracture, and this situation is not suitable for reserving oil and gas. But for good reservoirs, the range of the aspect ratio is [0.12, 1], and the expression C β s + 1 β p ≈ C β s can meet the condition, where 1 β p makes tiny influence on the result of the expression C β s + 1 β p . Thus, the final porosity calculating expression can be reduced to:

   (12) η = C β s ( β E − β s )

6. Working flow

We get all of the necessary parameters with which we can predict the porosity of other wells or a certain survey using the calculating expression (13). And the working flow is described in Figure 3.

3.1 Well testing

To test the porosity prediction method, we chose one well whose shear wave slowness was measured and neutron porosity and water saturated log data are to be interpreted yet. Then we calculated the rock pore geometry parameter C and predicted the porosity data of another well nearby. We take the first well named H1-1 to predict the porosity of the second well named H1-2. Following the working flow, we calculate the rock matrix compressibility β s = 0.00928 by cross-plot fitting (Figure 4a) and using the string-calculating method to calculate a string of β s (Figure 5) and also calculate the rock pore geometry coefficient C (Figure 6) and do cross-plot fittings between C and other elastic parameters like bulk modulus, P-wave modulus, Possion’s ratio, and so on (Figure 4b–d); thus, we get the parameter β E , which is of high correlation to the calculated parameter C and then we can calculate C of the second well H1-2.

Figure 4

(a) The cross-plot between β E and the measured porosity of well H1-1, (b) the cross-plot between β E and the pore geometry parameter C. (c) the cross-plot between compressional wave modulus and the pore geometry parameter C.

Figure 5

The calculating result of β s of the well H1-1.

Figure 6

The calculating result of C of the well H1-1.

Thus, we calculate the porosity of the second well H1-2, and the result is shown in Figure 7 as a comparison with the traditional TAR method. We can see that our method has better predicting effect.

Figure 7

Comparison of the predicting porosity result between the new method and the traditional method.

3.2 Field application

Taking the porosity inversion of a survey in China Sea, for example, the target stratum is a carbonate reservoir, which is identified as Neocene stratum formed lately. Following the prediction of the working flow, we reinverted the porosity of this survey using the prestacked inversion result data V p  ,   V s , and ρ of the target reservoir.

The predicting result is shown in Figure 8; from the result of the cross-well (H1) inline section, we can see high correlation to the interpreted porosity log data. And the predicting porosity attribute slice is shown in Figure 9 where a is the inversed porosity slice of carbonate top layer and b is the inversed porosity slice of carbonate bottom layer; from these porosity attribute slice, the porosity intensity distribution of the survey was shown. It is also correlated with the actual situation as the reservoir layer of H1 well and H3 well is filled with water due to its high porosity of about 32%.

Figure 8

Through-well (H1) porosity section correlated with porosity log.

Figure 9

(a) The predicting porosity intensity distribution of the time of the carbonate top layer plus 20 ms, and (b) the predicting porosity intensity distribution of the time of the carbonate bottom layer.