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

Kaiyuan Liu; Li Qin; Xi Zhang; Liting Liu; Furong Wu; Le Li

doi:10.1515/geo-2020-0147

Article Open Access

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

Kaiyuan Liu , Li Qin , Xi Zhang , Liting Liu , Furong Wu and Le Li

Published/Copyright: January 29, 2021

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Open Geosciences Volume 13 Issue 1

Abstract

Carbonate rocks frequently exhibit less predictable seismic attribute–porosity relationships because of complex and heterogeneous pore geometry. Pore geometry plays an important role in carbonate reservoir interpretation, as it influences acoustic and elastic characters. So in porosity prediction of carbonate reservoirs, pore geometry should be considered as a factor. Thus, based on Gassmann’s equation and Eshelby–Walsh ellipsoidal inclusion theory, we introduced a parameter C to stand by pore geometry and then deduced a porosity calculating expression from compressional expression of Gassmann’s equation. In this article, we present a porosity working flow as well as calculate methods of every parameter needed in the porosity inverting equation. From well testing and field application, it proves that the high-accuracy method is suitable for carbonate reservoirs.

Keywords: porosity inversion; Gassmann’s equation; compression coefficient; carbonate reservoir; pore geometry

1 Introduction

Porosity is one of the most important parameters in predicting reservoirs of oil and gas exploration to describe the capacity of reservoir storage. Thus, inverting porosity from seismic data makes it more critical in oil and gas exploration. This article attempts to invert porosity based on seismic data.

Many parameters, such as pore size, aspect-ratio, tortuosity, and pore complexity, play important roles in the relationship between porosity and p-wave velocity [1]. Acoustic velocity in carbonates has been verified to be dependent on pore geometry [2,3,4,5]. Acoustic wave velocity is not only correlated with porosity but more with pore geometry through experimental measurement of carbonate rock samples [6,7]. Carbonate rocks with different pore geometries have quite different acoustic velocities and permeabilities [8,9]. The exact relationship was proposed between pore geometry and acoustic attributes through digital image analysis using cross-polarized light [10,11]. Recent studies have claimed that the main geometric variable of interest is the pore aspect ratio ( α ), which is statistically computed from images by the ratio between the longest and the shortest semi-axes of the ellipse that encloses the pore [12,13]. Permeability, porosity, and acoustic properties are influenced greatly by pore–space parameters of carbonate rocks from X-ray tomography images and digital image analysis [14,15]. Pore sizes, tortuosity, and edginess of the pore are more important in affecting the rock stiffness, and thus the acoustic behavior of carbonates [16]. In one word, pore geometry is important for carbonate porosity inversion.

Nowadays, there are two main methods to inverse porosity. First, the time–average relationship (TAR) [17] and improved version [18] provide us with linear methods to achieve inverse porosity. Then the numbers of relationships were proposed to inverse porosity in complicated situations [19,20,21,22]. Anyway, in fact we do cross-plot fittings between porosity and other elastic attributes to get the corresponding relationships to a specified survey, as carbonate rocks frequently exhibit less predictable velocity–porosity relationships compared to those found in siliciclastics because of complex and heterogeneous pore structures [16]. Second, a calculating method based on Biot theory [23], which need not calculate well data and 3D seismic data, but we must calculate so many parameters, making it hard to widely use the method. The nonlinear method like neural networks’ ANN group method of data handling [24], multiple linear regression based on seismic attributes, etc., are practical but unsuitable [25,26], because they do not consider influencing factors like pore geometry [27]. In our work, we approach a new less parameter needed method based on Gassmann’s equation with the influence of the pore geometry. Through well-predicting experiments and survey application, we can get good correlation to the actual well data, which proves that our method is suitable and effective.

2 Methods

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].

(1) β D = β s 1 + m η α μ D = μ s 1 1 + n η α

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.

(2) m = K S ( 2 K S + 4 μ S ) / [ π μ S ( 3 K S + μ S ) ] n = π 15 8 ( 3 K S + 4 μ S ) 3 K S + 2 μ S + 4 ( 3 K S + 4 μ S ) 3 K S + μ S

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 .

Table 1

The relationship table between Possion’s ratio and m and n [30]

σ	0.15	0.2	0.25	0.3	0.35	0.4
m	1.20	1.10	0.99	0.86	0.72	0.54
n	0.38	0.36	0.34	0.33	0.31	0.29

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.

(3) 1 β E − β s = 1 β D − β s + 1 ( β p − β s ) η

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).

(4) 1 β E − β s = 1 β D − β s + 1 β p η

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.

(5) 1 β E − β s = 1 β s m α η + 1 β p η

Thus, we get the expression to calculate rock porosity:

(6) η = C β s + 1 β p ( β E − β s )

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.

β 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.
β 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
β 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
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
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 )
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 Results

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 {\beta }_{\text{E}} and the measured porosity of well H1-1, (b) the cross-plot between β E {\beta }_{\text{E}} and the pore geometry parameter C. (c) the cross-plot between compressional wave modulus and the pore geometry parameter C.$

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 {\beta }_{\text{s}} of the well H1-1.$

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.

4 Conclusions

In this work of predicting porosity, we explored a new method based on the rock physics. We have done so much research and analysis of the Gassmann’s equation combined with Eshelby–Walsh ellipsoidal inclusion theory. Thus, we exhibited the porosity inversion method. The well testing and field application proved the accuracy of our new method. From the work done in this work, we got the following conclusions:

The method of Gassmann’s equation combined with Eshelby–Walsh ellipsoidal inclusion theory is suitable for predicting porosity for carbonate reservoirs, especially the burial depth of reservoir is deep enough, where the rock pore fluid compression coefficient is obviously much bigger than the rock matrix compression coefficient.
Every parameter can be calculated from the well logs by statistical analysis, and the correlation between elastic attributes and parameters determined the resolution of inverted porosity.
The porosity predicting working flow is simple and easy to manipulate, and it is necessary for a high-resolution prestacked seismic data to improve the porosity predicting accuracy.

Acknowledgments

This work was supported by the project (NO. 2016ZX05023) of China Major national special projects and the Major special project (NO. 2019A-3307) of CNPC. The authors thank all the assessment experts and editors of the journal.

Conflict of interest: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Data availability statement: Some or all data, models, or code generated or used during the study are proprietary or confidential in nature and may only be provided with restrictions.

References

[1] Castro DDD, Rocha PLFD. Quantitative parameters of pore types in carbonate rocks. Rev Brasilra De Geofísica. 2013;31(1):125–36.10.22564/rbgf.v31i1.251Search in Google Scholar

[2] Anselmetti FS, Eberli GP. Controls on sonic velocity in carbonates. Pure Appl Geophys. 1993;141(2/3/4):287–323.10.1007/978-3-0348-5108-4_6Search in Google Scholar

[3] Anselmetti FS, Luthi S. Quantitative characterization of carbonate pore systems by digital image analysis. AAPG Bull. 1998;82(10):1815–36.Search in Google Scholar

[4] Baechle GT, Weger RJ, Massaferro L, Eberli GP. The role of macroporosity and microporosity in constraining uncertainties and in relating velocity to permeability in carbonate rocks. Soc Explor Geophys Expand Abstr. 2004;23(1):1662–5.10.1190/1.1845149Search in Google Scholar

[5] Hommel J, Coltman E, Holger C. Porosity – permeability relations for evolving pore space: a review with a focus on (bio-)geochemically altered porous media. Transp Porous Med. 2018;124(2):589–629.10.1007/s11242-018-1086-2Search in Google Scholar

[6] Anselmetti FS, Eberli GP. The velocity-deviation log: a tool to predict pore type and permeability in carbonate drill holes from sonic and porosity or density logs. AAPG Bull. 1999;83:450–66.Search in Google Scholar

[7] Lima Neto IA, Misságia RM, Ceia MA, Archilha NL, Oliveira LC. Carbonate pore system evaluation using the velocity-porosity-pressure relationship, digital image analysis, and differential effective medium theory. J Appl Geophys. 2014;110:23–33.10.1016/j.jappgeo.2014.08.013Search in Google Scholar

[8] Sun YF. A two-parameter model of elastic wave velocities in rocks and numerical avo modeling. J Comput Acoust. 2004;12(4):619–30.10.1142/S0218396X04002432Search in Google Scholar

[9] Cheng P, Liang N, Li R, Lan H, Cheng Q. Influence of pore structure and solid bitumen on the development of deep carbonate gas reservoirs: a case study of the longwangmiao reservoir in Gaoshiti–Longnusi area, Sichuan Basin, SW China. Energies. 2020;13:1–20.10.3390/en13010001Search in Google Scholar

[10] Weger RJ, Eberli GP, Baechle GT, Massaferro JL, Sun YF. Quantification of pore structure and its effect on sonic velocity and permeability in carbonates. AAPG Bull. 2009;10(10):1297–317.10.1306/05270909001Search in Google Scholar

[11] Dinghui Y, Meixia W, Xiao M. Symplectic stereomodelling method for solving elastic wave equations in porous media. Geophys J Int. 2014;1:560–79.10.1093/gji/ggt393Search in Google Scholar

[12] Weger RJ. Quantitative pore/rock type parameters in carbonates and their relationship to velocity deviations [PhD Dissertation Thesis]. Florida: University of Miami, Coral Gables; 2006. p. 158–67.Search in Google Scholar

[13] Mavko G, Mukerji T, Dvorkin J. The rock physics handbook: tools for seismic analysis in porous media, 2nd edn. New York: Cambridge University Press; 2009.10.1017/CBO9780511626753Search in Google Scholar

[14] Archilla NL, Missagia RM, Hollis C, de Ceia MAR, McDonald SA, Lima Neto IA, et al. Permeability and acoustic velocity controlling factors determined from X-ray tomography images of carbonate rocks. AAPG Bull. 2016;100(8):1289–309.10.1306/02251615044Search in Google Scholar

[15] Bashah S, Pierson B. The impact of pore geometry and microporosity on velocity-porosity relationship in carbonates of central Luconia, Sarawak. Proceedings of the AAPG international conference and exhibition, Singapore, September 2012; 2012.Search in Google Scholar

[16] Bashah NSI, Pierson B. The impact of pore geometry and microporosity on velocity-porosity relationship in carbonates of central Luconia, Sarawak. AAPG international conference and exhibition; 2012.Search in Google Scholar

[17] Wyllie MRJ, Gregory AR, Gardner LW. Elastic wave velocities in heterogeneous and porous media. Geophysics. 1956;21(1):41–70.10.1190/1.1438217Search in Google Scholar

[18] Raymer DS, Hunt ER, Gardner JS. An improved sonic transit time-to-porosity transform. Proc SPWLA 21st Ann Meeting; 1980. p. 1–13.Search in Google Scholar

[19] Hamilton EL, Bachman RT, Berger WH, Johnson TC, Mayer LA. Acoustic and related properties of calcareous deep-sea sediments. J Sediment Petrol. 1982;52(3):733–53.10.1306/212F803F-2B24-11D7-8648000102C1865DSearch in Google Scholar

[20] Han DA, Nur A, Morgan D. Effect of porosity and clay content on wave velocities in sandstones. Geophysics. 1986;51(11):2093–107.10.1190/1.1893163Search in Google Scholar

[21] Kahraman S, Yeken T. Determination of physical properties of carbonate rocks from P-wave velocity. Bull Eng Geol Environ. 2008;67(2):277–81.10.1007/s10064-008-0139-0Search in Google Scholar

[22] Xie M, Mayer KU, Claret F, Alt-Epping P, Jacques D, Steefel C, et al. Implementation and evaluation of permeability-porosity and tortuosity-porosity relationships linked to mineral dissolution-precipitation. Comput Geo Sci. 2015;19(3):655–71.10.1007/s10596-014-9458-3Search in Google Scholar

[23] Yingbo Z. Seismic porosity inversion method and its application. Oil Geophys Prospect. 1994;29(3):261–73.Search in Google Scholar

[24] Nahser MA, Wang YH. Porosity prediction using the group method of data handling. Geophysics. 2011;76(5):o15–22.10.1190/geo2010-0101.1Search in Google Scholar

[25] Emilson PL, Alexandre CV. 3D porosity prediction from seismic inversion and neural networks. Comput Geosci. 2011;37(8):1174–80.10.1016/j.cageo.2010.08.001Search in Google Scholar

[26] Linqi Z, Zhang C, Guo C. Calculating the total porosity of shale reservoirs by combining conventional logging and elemental logging to eliminate the effects of gas saturation. Petrophysics. 2018;59(2):162–84.10.30632/PJV59N2-2018A4Search in Google Scholar

[27] Zhu L, Zhang C, Zhou X, Wei Y, Huang Y, Zhang C. Take the inversion of permeability of tight gas reservoir with the combination of deep belief kernel extreme learning networks and NMR logging data. Interpretation. 2017;37(7):1–43.Search in Google Scholar

[28] Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc Roy Soc Lond. 1957;A241:376–96.10.1007/1-4020-4499-2_18Search in Google Scholar

[29] Walsh JB. The effect of cracks on the compressibility of rock. J Geophys Res. 1965;70(2):381–9.10.1029/JZ070i002p00381Search in Google Scholar

[30] Chen Y, Huang TF. Seismic wave propagation in fractured media: equivalent medium theory. In: Rock Phys (in Chinese). Hefei: Press of University of Science and Technology of China; 2009. p. 403–50.Search in Google Scholar

[31] Gassmann F. Uber die Elastizitat poroser Medien Vier der Nat. Gesellschaft in Zurich. 1951;96:1–23.Search in Google Scholar

[32] He XL, He ZH, Wang RL, Wang XB, Jiang L. Calculations of rock matrix modulus based on a linear regression relation. Appl Geophys. 2011;8(3):155–62.10.1007/s11770-011-0290-4Search in Google Scholar

[33] Jiang L, Wen XT, HE ZH, Huang DJ. Pore structure model simulation and porosity prediction in reef-flat reservoir. CHJ Geophys (Chinese). 2011;54(6):1624–33.10.1002/cjg2.1623Search in Google Scholar

[34] Lin K, Xiong XJ, Yang X, He ZH, Cao JX, Zhang XH, et al. Self-adapting extraction of matrix mineral bulk modulus and verification of fluid substitution. Appl Geophys. 2011;8(2):110–6.10.1007/s11770-011-0278-0Search in Google Scholar

[35] Liu WL. Reservoir exploitation seismic techniques. Beijing: Petroleum Industrial Press; 1996.Search in Google Scholar

Received: 2020-05-19

Revised: 2020-10-16

Accepted: 2020-10-21

Published Online: 2021-01-29

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/geo-2020-0147

Keywords for this article

porosity inversion; Gassmann’s equation; compression coefficient; carbonate reservoir; pore geometry

Creative Commons

BY 4.0