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.
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 . 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 . 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)  and improved version  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 . Second, a calculating method based on Biot theory , 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 , multiple linear regression based on seismic attributes, etc., are practical but unsuitable [25,26], because they do not consider influencing factors like pore geometry . 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.1 Eshelby–Walsh ellipsoidal inclusion theory
The ellipsoidal inclusion theory for dry rocks  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 .
where and are compressibility of dry rock and rock matrix, respectively; is rock porosity; and 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 and 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 .
2.2 Gassmann’s equation and the deduction of porosity inversion method
The most famous Gassmann’s equation  using rock compressibility coefficients is expressed as follows.
where is the rock saturated fluid equivalent compression coefficient, is the rock matrix compression coefficient, is the dry rock compression coefficient, 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 approximate to , which has approved high precision to calculate the rock matrix compression coefficient , then we can get the expression from (3).
Then replacement of with from the expression (1), we get the following expression.
Thus, we get the expression to calculate rock porosity:
where 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 , and C, which can be gained from the log data certainly to calculate the porosity.
, 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 .
A linear equation can be transformed from expression (6) as follows.(7)
Expression (7) is the linear relation between the rock saturated fluid equivalent compression coefficient and the rock porosity. Where(8)
Then 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  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 as the liquid part of Gassmann’s equation
And which can also be expressed as a liquid part of Russell’s expression
where and are p-wave and s-wave impedance, respectively, c is a constant, when , is an initial small value, we can get a string of values.
, the rock pore fluid compression coefficient, can be calculated by the following expression and with the saturated water log data , where H is the depth.(9)
, the rock saturated fluid equivalent compression coefficient, can be calculated by elastic parameters as inversed by prestacked seismic data, and the calculation expression is as described below.(10)
C, the rock pore geometry parameter, is generally calculated by the expression , 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)
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 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 can meet the condition, where makes tiny influence on the result of the expression . Thus, the final porosity calculating expression can be reduced to:(12)
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 by cross-plot fitting (Figure 4a) and using the string-calculating method to calculate a string of (Figure 5) and also calculate the rock pore geometry coefficient (Figure 6) and do cross-plot fittings between and other elastic parameters like bulk modulus, P-wave modulus, Possion’s ratio, and so on (Figure 4b–d); thus, we get the parameter , which is of high correlation to the calculated parameter C and then we can calculate C of the second well H1-2.
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.
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 , , 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%.
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.
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.
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 Dinghui Y, Meixia W, Xiao M. Symplectic stereomodelling method for solving elastic wave equations in porous media. Geophys J Int. 2014;1:560–79.Search in Google Scholar
 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
 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
 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
 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
 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
 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
 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
 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.Search in Google Scholar
 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
 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
 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
 Yingbo Z. Seismic porosity inversion method and its application. Oil Geophys Prospect. 1994;29(3):261–73.Search in Google Scholar
 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
 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
 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.Search in Google Scholar
 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
 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
 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
 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
 Gassmann F. Uber die Elastizitat poroser Medien Vier der Nat. Gesellschaft in Zurich. 1951;96:1–23.Search in Google Scholar
 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
 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
 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
 Liu WL. Reservoir exploitation seismic techniques. Beijing: Petroleum Industrial Press; 1996.Search in Google Scholar
© 2021 Kaiyuan Liu et al., published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.