Carbonate Pore Shape Evaluation Using Digital Image Analysis, Tomography, and Effective Medium Theory

Abstract

Carbonate rocks have a wide variety of pore shapes and different types of grains, which greatly affect the elastic properties and characteristics of the reservoir. This causes certain difficulties in petroelastic modeling. One of the problems is the scale of the input data, which is then used to build the rock physics model. The paper presents the results of studying three core samples of carbonate rocks of the Upper Devonian and Lower Carboniferous age, which are located in the South Tatar arch (Volga-Ural oil and gas basin (Russia)). To evaluate the structural characteristics of the pore space, the effective medium theory is used. The input data are the results of laboratory studies that include measurements of the velocities of longitudinal and transverse waves, porosity, and thin section and computed tomography analysis. When using the computed tomography, the core samples are analyzed at different resolution (12–37 µm/voxel). The tomography studies of pore space at different scales provide rather different values of porosity and pore aspect ratio. The tomography-based porosity estimations also differ from the experimentally measured porosity (up to 10%). The pore space characteristics provided by different datasets are used to build a rock physics model for the studied rocks that helps to estimate the elastic wave velocities with three different methods of effective medium theory (self-consistent approximation, differential effective medium (DEM), and the Kuster–Toksöz method). A comparison of the velocity estimations with their experimental analogs for dry rocks may indicate the presence of microcracks whose size is beyond the tomography resolution. Improved rock physics models incorporating both pores and microcracks are then used to predict the elastic wave velocities of fluid-saturated rock in a wide porosity range. It is demonstrated that the predicted values significantly differ (up to 30%) from those provided by the rock physics (RP) models constructed without the support of the tomography results. Moreover, other types of models are considered in which the difference in experimental and theoretical velocities is attributed to changes in the host matrix properties as compared to the calcite polycrystal, which are caused by various reasons.

1. Introduction

According to many studies, up to half of the world’s hydrocarbon reserves are concentrated in carbonate reservoirs [1,2,3]. Carbonate reservoirs, as a rule, are very complex in structure, and their behavior is difficult to predict due to the high heterogeneity. In the Republic of Tatarstan (Russia), a high density of deposits concentrated in carbonate reservoirs is observed on the South Tatar arch and the Melekes depression. In these structures, in carbonate rocks, oil reserves amount to 35–40% of all proven reserves. However, recoverable reserves at the level of modern technologies are estimated at only 10–15% [4]. Carbonates are highly heterogeneous and show different pore systems (moldic, interpar-ticle, intraframe, and microcracks), and oil in the pores is quite viscous [5]. Traditional oil production technologies used in siliciclastic reservoirs show low efficiency in carbonate rocks. At present, the oil recovery factor in the fields of the Republic of Tatarstan, localized in carbonate reservoir rocks, is very low and amounts to 0.15–0.25 [6]. One of the main reasons for the low efficiency is the insufficient knowledge of the pore space. This may lead to the incorrect accounting of these characteristics in the technologies of oil development from such formations. The type of pore space in a rock determines not only the ability to store and transport hydrocarbons, but also characterizes the elastic properties of the rock [7,8]. Elastic properties in sedimentary rocks are affected not only by porosity but, more importantly, by the type and shape of the pores. Elastic wave velocities are a measure of the compressibility and stiffness of a rock. This, in turn, depends on the shape and number of pores [9,10,11,12].

On the plots of porosity versus velocity, a clear inverse relationship is revealed: the higher the porosity, the lower the elastic properties. However, quite often, a large scatter of values is observed on the velocity–porosity plot [13,14,15]. An example of high dispersion of P-wave velocity vs. porosity is shown in Figure 1 for carbonates of Lower Carboniferous age, which are located in the South Tatar arch. This leads to large uncertainties in porosity calculations. The weak relationship between porosity and velocity is due to the ability of carbonates to contain a variety of shapes, sizes, and types of pores, which can both increase and decrease the elastic properties [14,16].

Due to the wide variety of pore types, the pore space must be classified taking into account the structural features of the rocks and their petrophysical properties. One of the first who related structural characteristics and petrophysical properties in carbonates was Gustav Archi [17,18]. His classification focused on the evaluation of porosity. A detailed classification of the pore space was proposed by P. Choquette and L. Pray [19]. It is based on the genesis of the pore space and is not related to the petrophysical properties of the rocks. The basic principle of this classification is to divide pores into two classes: pores associated with the rock structure and pores not associated with it. The first group incorporates intergranular and intercrystalline pores, which include caverns that were formed as a result of the dissolution of rock components (moldic) and intraform caverns.

The pore space classification of carbonate rocks, which takes into account petrophysical properties, was proposed by Jerry Lucia [20]. Considering the results of laboratory studies of porosity, permeability, and capillary properties, as well as the description of the rock matrix, Lucia came to the conclusion that the most optimal is the division of void space into two main types: interparticle and vuggy porosity. Vuggy porosity, in turn, is divided into two classes depending on the connectivity of cavities with each other: (1) isolated cavities connected only by interparticle pores and (2) interconnected cavities connected by channels and voids. According to Lucia, vuggy porosity describes pores comparable to or larger than grains and crystals; accordingly, porosity is not interparticle. Caverns are most often formed as a result of the dissolution of fossil shells, particles, grains, or large cavities. Fissures are also included in this group despite the fact that they may occur not only as a result of sedimentation and diagenesis.

The correlation of different petrophysical characteristics of carbonate reservoir rocks with structural and textural features was reported by Kosarev et al. [21]. These characteristics include the bulk density, dynamic and static Young’s moduli and Poisson’s ratios, elastic wave velocities, cohesion, and internal friction angle.

The relationship between elastic wave velocity and porosity in carbonates may be good if the rocks are classified considering the individual characteristics of the deposit and sedimentation. Thus, based on the dependence of the P-wave velocity on porosity and considering various diagenetic processes, four types of rocks were distinguished in [22]: purely microporous limestones, samples with preserved intergranular and moldic pores, and vuggy limestones. In work [23], the authors analyzed the results of laboratory studies on 214 samples of carbonates from the Lower Cretaceous Provence platform (Urgony platform, Southeast France) in order to quantify the effect of pore type and diagenetic changes on acoustic properties. Even though many difficulties arise in the study of carbonates (laboratory errors, scaling issue, etc.), the authors were able to classify the rocks.

In this work, we studied three samples of carbonate rocks. The pore space of the samples was modeled using the effective media theory. The input data were laboratory studies of the internal structures of rocks. For this, thin sections were prepared and experimentally analyzed with computed tomography.

2. Materials and Methods

2.1. Geological Setting

Three core samples were studied, taken from two oil fields. Sample 1 (sampling depth of 952.7 m) is from the field A, which is located on the southwestern slope of the so-called South Tatar arch of the Volga-Ural basin (Figure 2). The geological section of the region is represented by deposits of the Devonian, Carboniferous, Permian, Neogene, and Quaternary systems. The sedimentary cover rests on the rocks of the crystalline basement. The studied sample is confined to the Tournaisian stage of the Lower Carboniferous, which is promising for hydrocarbons and is being developed. The absence of reliable bridges (impermeable layers) allows us to consider the entire Tournaisian stage as a single, hydrodynamically connected reservoir. Porous–permeable layers are represented by oil-saturated limestones, which alternate with compacted, clayey, often fractured limestones.

Samples 2 (sampling depth of 975.3 m) and 3 (sampling depth of 1126.3 m) were taken from field B (Figure 2), which is located on the southeastern slope of the South Tatar arch. Accumulations of hydrocarbons have been established in carbonate deposits of the Upper Devonian and Lower Carboniferous. Sample 2 is confined to Upper Devonian deposits. Sample 3 was taken from the deposits of the Tournaisian stage of the Lower Carboniferous. All samples were drilled in the horizontal direction.

2.2. Petrographic Description

According to the study of the thin section, Sample 1 is represented by a crinoidal wackestone–packstone with peloids (white arrows) and rare foraminifera (black arrows) (Figure 3a). The grains are 40% of the rock volume, and the matrix is 60%. Crinoids are often represented by fragmented grains, and rarely well-preserved grains; peloids are rare, and individual forms and fragments of algae are even rarer. The distribution of grains is uniform. The grains are not supported. The grain size of crinoids is 0.2–0.4 mm, that of peloids is 0.1 mm, and fragments of foraminifera also reach 0.1 mm. The rock matrix is represented by fine-grained calcite, which is uniformly distributed between grains and fills the entire intergrain space. The pore space is represented by intragrain (red arrows) and intergrain (blue arrows) types of porosity. The visible pores are more often isolated and less connected. The pore size varies from 0.005 mm to 0.1 mm. It has a prevailing size of 0.03–0.05 mm.

According to the lithological description of the thin section, Sample 2 is represented by a crinoidal packstone, with fragments of ostracoda and peloids (Figure 3b). Grains make up to 65%, and the matrix makes up 35% of the volume of the rock. Grains are supported and mostly represented by fragments of crinoids and their fragments are up to 1 mm in size. Significantly fewer grains are represented by fragments of ostracod shells, up to 1.2 mm, and peloids, up to 0.3 mm. The matrix is represented by sparry calcite, which occupies the entire space between grains. The pore space is represented by two types of porosity. Separate-vug pore space predominates (green arrows), up to 0.7 mm in size, and there are significantly fewer intergrain pores (blue arrows), up to 0.05 mm. According to the results of the microscopic description of the thin section, Sample 3 is represented by a crinoidal wackstone (Figure 3c). The grains are 30% of the rock volume; the matrix is 70%. Crinoids are represented by fragmented grains (yellow arrows), and rarely well-preserved grains. The distribution of grains is quite uniform. The grains are not supported. Crinoids are often represented by intensely fragmented grains, and rarely well-preserved grains. The grain size of crinoids is 0.1–0.5 mm. The rock matrix is represented by fine-grained calcite, which is uniformly distributed between grains and fills the intergrain space (blue arrows). Pores are rare and isolated, usually unconnected and represented mostly by the secondary dissoluted intergrain type of porosity. The pore size varies from 0.001 mm to 0.05 mm. It has a prevailing size of 0.02 mm.

2.3. Laboratory Measurements

Laboratory petrophysical studies were carried out on three core samples. The samples were cylinders with a diameter of 30 mm. The ends of the specimens were ground to within 0.01 mm. Before measurements, the samples were dried in an oven at 60 °C for at least 72 h. Before drying, the hydrocarbon was extracted from the samples using a Soxhlet extractor. An alcohol–benzene mixture was used as solvents. The bulk density was calculated from the weight of the dry sample and the volume of the cylinder. Open porosity was measured with the use of the PLAST-215 ATM gas porosimeter (Russia). The velocities of longitudinal and transverse waves (P- and S-waves) propagating in rocks under conditions simulating reservoir ones were measured with the use of a PIK UZ-UEP unit (Russia). The measurements were carried out at a confining pressure of 18 MPa and at a temperature of 23 °C. The source emits a signal that enters one of the transducers and excites pulses of longitudinal and transverse waves (a frequency of about 1 MHz). The waves, after passing through the sample, are recorded by a second sensor. The receiver signal is transmitted to the oscilloscope.

Thin sections were analyzed with the help of a translucent polarizing microscope of the research class Axio Imager (Carl Zeiss, Jena, Germany). To determine the mineral phases, observations were used at two positions of the nicols (parallel and across). During the study, lenses with a magnification factor of 5, 10, and 20 were used. For documentation, a regular digital color camera AxioCam ICc 5 from Carl Zeiss and compatible Zen Pro software 3.2 from this company were used.

The process of obtaining a three-dimensional image on X-ray computed tomography (CT) consists of five stages: (1) preparation of an X-ray computed tomography, (2) tomographic reconstruction, (3) evaluation and treatment of X-ray artifacts, (4) segmentation of pore or grain phases, and (5) solving equations based on the targeted properties. First, standard core plugs (30 mm in length and 30 mm in diameter) were scanned with the help of the General Electric v|tome|x S240 micro- and nanofocus X-ray research system for computed tomography. During CT scanning, the following parameters were used: 140 A current, 160 kV voltage, 1200 projections, 2 averages, and 200 ms timing for 1 projection. The obtained voxel resolution was 37.4 µm. After all measurements, from the inner-most central part of all standard samples, cylindrical subsamples with a diameter of 6 mm and a length of 10–14 mm were drilled. Next, all subsamples were rescanned with the following settings: current of 110 µA, voltage of 100 kV, number of projections of 1200, average of 2, and timing for 1 projection of 200 ms. The resolution of computed tomography was 12–16 µm/voxel on small subsamples.

The reconstruction of the 3D models was achieved with the use of phoenix datos|x reconstruction software. The digital analysis was provided in the Avizo 7.1 software. We used a nonlocal means filter (local neighborhood of 3 pixels, search window of 1 pixel, spatial standard deviation of 1, and intensity standard deviation of 0.1) to reduce the noise and improve the image quality. The region of interest (ROI) was extracted from the center of all scanned volumes. Segmentation of porous space was provided with the help of interactive thresholding and based on the difference in grayscale levels for the air and mineral matrix. The porosity coefficient was calculated as the ratio of the segmented pore volume to the total sample volume. The porous space was separated by the throats with the use of a skeleton-aggressive approach [24], and 3D width and 3D length were calculated for each pore. The aspect ratios for each pore were calculated as the ratio of 3D width to 3D length.

2.4. Application of Effective Medium Theory

The effective medium theory gives an opportunity to relate the measured (or effective) properties with the rock composition and characteristics of microstructure. When constructing a rock physics model, first, a studied rock is replaced by a model medium that reflects the principal features of the rock’s microstructure. Then, parameters describing the rock composition and microstructure should be specified. The compositional parameters include the volume concentrations and elastic moduli of components. The microstructural parameters characterize the shape, orientation, and connectivity of components. Finally, a method of the effective medium theory is to be chosen in order to relate the studied physical properties (for example, elastic) and the model parameters. The aim of RP models is to predict how the physical properties vary for different scenarios in changing rock composition and/or microstructure. For example, we can estimate how the elastic wave velocities change if we fill the pores with another fluid, or increase/decrease the porosity, or change the pore shape, or insert other minerals and organic material in the rock.

For the RP modeling (RPM), we applied three methods of effective medium theory: self-consistent approximation (SCA) of Berryman [25], the differential effective medium (DEM) approach [26], and the Kuster–Toksöz method [27]. These methods are widely used in geophysical practice. Each of the methods relates the effective bulk and shear moduli to the rock composition and shape of ellipsoidal inclusions characterized by aspect ratios. In all of the methods, the inclusions are assumed to be ellipsoids of revolution and only one aspect ratio is required to describe their shape. In the case of voids, the aspect ratio characterizes their relative opening. The DEM and Kuster–Toksöz methods assume the existence of a host matrix. The SCA considers all heterogeneities as ellipsoidal inclusions embedded in the matrix with effective properties. The question definitely arises of how the results depend on the EMT method selected for computations. We performed a special analysis that allowed us to give an answer to this question. We compared the results obtained with the methods for two simple basic models of carbonate rock: (1) calcite matrix with quasi-isometric pores and (2) calcite matrix with cracks. The formulas of the self-consistent approximation [25], DEM [26], and Kuster–Toksöz [27] methods are presented in Appendix A. The formulas of these methods can be also found in [28].

When constructing the RP model, all possible information on rock’s composition, microstructure, and measured physical properties should be taken into account. The results of thin section analysis and tomography studies provide valuable data on rock’s microstructure. Thus, the thin sections give information on the heterogeneity’s type and their mutual distribution in the rock’s volume. However, the shape of heterogeneities (specifically, voids) can be more reliably estimated from the tomography data, since the thin sections provide only 2D information.

For rocks with a complicated and hierarchical inner structure, the constructed RP model may have many unknown parameters. At the same time, the number of available experimental data is limited. This leads to the fact that many sets of model parameters can be found which produce values of effective properties rather close to experimental ones. In order to overcome this obstacle, a sensitivity study of the constructed RP model should be performed. Usually, this helps to simplify the model or decrease a number of unknown parameters. To construct RP models for the studied rocks, we selected the model medium such that it reflected the specific features of the rock’s microstructure recognized by the thin-section analysis and tomography. Thus, we should take into account the existence of porous micritic cement confirmed by the thin section analysis. The calcite grains and pores visible with computed tomography (“visible” pores) are placed in this cement. This RP model contains an unknown parameter, namely, the porosity of the micrite matrix. In some works [29], the microporosity is suggested to determine the difference between the measured porosity and porosity provided by the computed tomography. In our modeling, we also followed this line since the analysis of density of the studied samples suggested that the closed porosity was quite small (around 1%).

To calculate the effective elastic moduli of the micrite matrix, we applied the SCA method. To estimate the elastic moduli of the micritic matrix with calcite grains and “visible” pores, we used the Kuster–Toksöz method, since this method assumes the existence of a host matrix.

We did not distinguish in our RP model between “visible” pores of different size, since the porosity was too small to do that. The difference in size of inclusions is taken into account implicitly in hierarchical RP modeling according to the principle “from small heterogeneities to large heterogeneities” [30]. Thus, previous studies performed for carbonate rocks showed that the “size effect” produces a 3% difference in bulk and shear moduli at a total porosity of 20% when vuggy porosity varies from 1 to 10% [31].

The RPM in our study was used to solve both the forward and inverse problems. The forward problem assumes the calculation of the effective elastic moduli, and then elastic wave velocities from the known composition of samples and microstructural parameters provided by different sets of experimental data we have. The theoretical values of longitudinal and transverse wave velocities (Vp and Vs) are compared with the respective experimental values. This comparison gives a possibility to answer the question of whether the RP model can be further used for the predictions mentioned above. In the case of unacceptable misfit between the theoretical and experimental velocities, the RP model should be modified. Thus, smaller voids can exist in rocks that are not recognizable with the used resolution of tomography (for example, microcracks). In this case, an inverse RP-based problem should be solved to estimate the parameters of the hidden voids.

Voids filled with a contrasting substance compared to the host matrix properties (for example, gas or liquid) have a significant effect on the overall elastic moduli. If the voids have nonisomeric shape (microcracks), even their small volume concentration (less than 1%) can lead to a significant drop in the elastic moduli of rocks (see Section 3.2). Microcracks are thin enough to be detected with tomography. Thus, we should not exclude them from consideration, specifically in the case of packstones. The occurrence of microcracks in packstones was confirmed by previous lithological studies in this area [4,5].

It should be noted that many factors except voids can affect the elastic wave velocities. These factors include the state of microfabrics (intact or stressed grains), grain boundaries, and many others. To take into account the overall effect of these factors we can consider that the elastic moduli of the host matrix are also unknown values to be inverted from experimental data. However, for isotropic rocks, only two experimental velocities are available, namely, Vp and Vs. This means that, even assuming only one set of invisible voids, we should invert four unknowns from these measurements. In this case, the domain of possible solutions of the inverse problem can be quite wide. As shown below, even for the known elastic moduli of the host matrix, this domain is not unique. Therefore, considering a model with hidden voids, we assumed that the elastic moduli of the host matrix are known and equal to those of the calcite polycrystal. The difference between the experimental velocities and velocities obtained with RPM, which includes only the voids provided by tomography, was attributed to invisible voids.

We also considered another type of RP model. This RP model explains the discrepancy in experimental and theoretical velocities obtained with RPM, which incorporates only tomography-visible voids, by various factors changing the elastic moduli of the host matrix. Specifically, one of these factors is the effect of the grain boundaries. In the vicinity of the grain boundaries, the grain material can differ from those of the perfect crystal. Below, we call a substance near the grain faces “intergranular material”. The elastic moduli of this material are included in the list of unknown parameters instead of the characteristics of hidden voids. In this case, the unknown elastic moduli of intergranular material were inverted in two steps. First, we considered that the quasi-isometric voids visible with computed tomography were placed in an isotropic host matrix with unknown elastic moduli. These moduli were inverted from the experimental velocities Vp and Vs. Then, we assumed that the host matrix was also a heterogeneous medium composed of particles of calcite polycrystals surrounded by a soft material with unknown properties. This soft material plays the role of intergranular material. The volume concentration of the soft material should be rather small. We assumed that it did not exceed 1–2%.

For the inverse problem solution, we applied the method of N-dimensional meshes. According to this method, first, for each parameter, a range of its possible variation is specified. Then, the whole range of variation in each parameter is divided into intervals. As a result, we obtained a N-dimensional mesh with nodes at the intersection of the interval’s ends. The forward problem was solved for each node. Then, the sets of parameters providing an “acceptable” discrepancy (a discrepancy comparable to experimental error) between the theoretical and experimental values of elastic wave velocities were chosen as possible solutions to the inverse problem. Finally, the most probable solution was selected among the possible solutions. The most probable solution is the solution that provides the minimum value of the chosen measure of discrepancy in both the velocities Vp and Vs. In this work, the discrepancy measure was taken in the form:

$$Discrepancy measure={\sqrt{0.6{\left[Diff(Vp)\right]}^2+0.4{\left[Diff(Vs)\right]}^2}}^{}·100\%$$

(1)

where Diff(Vp) and Diff(Vs) are the relative differences between the theoretical and experimental values of Vp and Vs. velocities, respectively.

In the RP model with microcracks, N = 3 (mineral matrix, pores, and microcracks). Instead of values of microcrack porosity and aspect ratio, we considered their decimal logarithms. Decimal logarithms were used due to the strong nonlinearity of the dependence of elastic wave velocities on the aspect ratio of cracks and their porosity in the domain of their small values. As a result, the two-dimensional region of possible change in the decimal logarithms of the desired parameters was covered with a mesh. In the mesh nodes, the forward problem of determining the velocities of elastic waves was solved.

When solving the two-step inverse problem for a model with unknown properties of the host matrix, at each step, N = 2.

For the RP model with microcracks, we assumed that a solution was acceptable if the relative difference between the theoretical and experimental velocities did not exceed the threshold value of 5% (a value comparable to the velocity measurement error). If this condition was not met for the specified ranges of parameter values, then the threshold was increased by 1%. Moreover, this increase was carried out sequentially—first, the threshold for the S-wave was increased by 1%. If solutions were not determined again, then the threshold was increased by 1% for the P-wave as well. We observed from the modeling carried out for the three samples that, in the “worst” case, the threshold value for velocity did not exceed 7% for P-wave and 7.5% for S-wave. Note that according to our previous studies, the value of 6–7% was established as the degree of heterogeneity of elastic wave velocities in carbonate rocks using a specially developed method of multilevel ultrasonic sounding [32].

When solving the inverse problem for a model with unknown properties of the host matrix, we considered a discrepancy measure similar to (1) but for elastic moduli instead of velocities. The reason for this is that the density of the intergranular material was unknown. The threshold of 10% for this measure was taken for collecting acceptable solutions.

The most probable solution for the model parameters provided by the inverse problem solution were further used for the predictions of elastic wave velocities for different scenarios of changing the saturation type, porosity, void shape, and mineral composition. In this work, we considered RP-based predictions of velocities for water-saturated rocks in a wide porosity range.

3. Results

3.1. Laboratory Studies

According to the results of laboratory studies, the following petrophysical parameters of the carbonate rocks were obtained: longitudinal and transverse (or P- and S-) wave velocities, bulk density, and porosity. Moreover, the porosity was obtained in two different ways: using the gas permeameter (φ1) and with computed tomography. According to the results of tomography, two porosity values were also obtained: on a standard sample with a diameter of 30 mm (φ2) and on a small sample with a diameter of 6 mm (φ3). For two sample sizes (standard and small), the mode of aspect ratio was estimated from probability density functions (PDFs) shown in Figure 4. The results of laboratory studies are presented in

Table 1. Results of laboratory measurement.

Parameter	Sample 1	Sample 2	Sample 3
Bulk density (g/cm3)	2.388 (±0.008)	2.586 (±0.009)	2.463 (±0.009)
S-wave velocity (km/s)	2.4 (±0.102)	2.9 (±0.128)	2.9 (±0.127)
P-wave velocity (km/s)	4.5 (±0.175)	6.1 (±0.279)	4.9 (±0.181)
Porosity (φ1) (%)	11.49 (±0.5)	3.76 (±0.5)	7.88 (±0.5)
Porosity (φ2) (%)	0.43	2.16	0.13
Porosity (φ3) (%)	7.43	4.53	2.19
Aspect ratio mode (standard sample)	0.50	0.50	0.51
Aspect ratio mode (small sample)	0.50	0.55	0.52

.

When comparing the porosity obtained using different methods and for different rock volumes, it follows that the pores in carbonate rocks were distributed unevenly. This can be also observed on volumetric images provided by the computed tomography (Figure 5 and Figure 6). Figure 6 (left parts) shows the porosity distribution obtained with two different resolutions of computed tomography (around 12 and 37 µm/voxel). The dashed contours on the left indicate the area of the small sample. Voxels for the sample with a 6 mm diameter are shown in green. Voxels for the sample with a diameter of 30 mm are indicated in blue. However, as seen from the PDFs shown in Figure 4, the pore shape distribution was rather similar for different resolutions. The representative elementary volume (REV) was estimated as a volume beginning from which the porosity stayed stable. To find the REV, a volume of the cubic shape was considered inside the sample. We started from the cube with an edge length of 1 mm. We gradually increased the edge length up to 17 mm and calculated the porosity for each cubic volume.

The porosity-based REV analysis performed for the 3D images showed that Samples 1 and 3 exhibited porosity curves that converged to a plateau when the cube edge became close to 4 mm. Sample 2 had a higher REV than the sample’s size and was characterized by a high degree of heterogeneity of the porous structure, which is reflected in Figure 7. Meanwhile, the porosity of this sample stayed stable for the sample size varying from 4 to 12 mm.

Due to the presence of voids of different size, a significant difference in porosity values obtained with various methods was observed.

3.2. Rock-Physics-Based Solution of Forward Problem

When constructing the RPM for the studied samples, the following sets of input data were considered for each sample:

Dataset 1. The porosity is measured with the porosimeter. The pore aspect ratio is estimated (a) on the small sample and then (b) on the big sample.

Dataset 2. The porosity and pore aspect ratio are provided by tomography for the small sample.

Dataset 3. The porosity and pore aspect ratio are provided by tomography for the big sample.

Moreover, at the initial stage of RP modeling, we took into account the amount of micrite matrix recognized for the samples via thin-section analysis.

The elastic moduli of components used for the modeling were as follows. The bulk and shear moduli of the calcite polycrystal were calculated with the self-consistent approximation from the stiffness matrix of the calcite monocrystal presented in [33]. Note that the self-consistent approximation is applicable to the effective elastic properties of polycrystals as well. The bulk and shear moduli were, respectively, equal to 75.1 and 30.3 GPa. The bulk modulus of air was 0.0001 GPa [28]. The density of calcite and air was taken as equal to 2.70 g/cm³ and 0.001 g/cm³ [28].

First, we considered a simple model of rocks consisting of the calcite matrix and “visible” dry pores of equal shape. We applied the SCA, DEM, and Kuster–Toksöz methods in a rather wide porosity range—from 0 to 35%. Our goals were: (1) to compare the results provided by different theoretical methods and (2) to see the difference between the theoretical and experimental results for this model. The results obtained with the three methods are shown in Figure 8 together with the experimental data. The vertical and horizontal bars show the errors in velocity and porosity determination. The errors are, respectively, ±4% and ±5% for Vp and Vs velocities and ±0.5% porosity (in the porosity range up to 16%) as declared in unit test documents.

As can be seen, with the exception of the P-wave velocities for Sample 2, the theoretical velocities obtained with the simple model significantly exceeded the corresponding measured values.

Then, we complicated the RP model by bringing it closer to the results of thin-section analysis. Namely, we incorporated the porous micrite matrix in the model. Samples 1 and 3 had large amount of micrite—60% and 70% of the sample volume, respectively. For Sample 2, the micrite cement was not too high—around 35%. As was mentioned above, the micrite porosity was assumed to be a difference between the measured porosity and porosity provided by the computed tomography. We took the porosity given by the computed tomography performed for the small samples. We found that the porosity of the micrite matrix was 6.7% for Sample 1 and 8.1% for Sample 3. The tomography porosity for Sample 2 was smaller than the measured porosity. Therefore, we considered that the micrite porosity was very low or even absent in this sample. We compared the results obtained with the simple model described above in this section with the two-step model of rocks with porous micritic cement presented in Section 2. The comparison showed that the difference in elastic wave velocities provided by these two models did not exceed 1% for all samples. The reason was that the total porosity was not high. In this case, a distribution of porosity between “visible” and micrite pores did not produce a noticeable effect. Therefore, for further computations performed for the studied samples, we applied the model without a micrite matrix. However, in Section 4, we again consider the micrite-based model for the prediction of elastic properties for higher porosity.

We performed a numerical analysis that allowed us to conclude that an RP model incorporating ten systems of voids characterized by histograms (Figure 4) can be replaced by a model with a single system of voids. The aspect ratio of the voids was equal to the mode values of the histogram (0.50–0.55). Thus, the difference between the Vp and Vs velocities for Sample 1, calculated with the two models for porosity provided by the porosimeter (11.49%), was around 1%. This was much smaller than the experimental errors for Vp and Vs (4% and 5%, respectively).

From the analysis of the obtained modeling results, it follows that for the voids’ aspect ratio, provided by the tomography, and in the range of porosity of the samples 0.58–11.5%, all methods showed similar results. Therefore, only the Berryman self-consistent method was chosen for further modeling.

The overestimated values of the theoretical elastic wave velocities compared to their experimental analogs may be due to the presence of smaller voids, the size of which was beyond the resolution of the tomographic study. However, the volume concentration of these voids should not be too high. Otherwise, the porosity of the rock will not correspond to the experimental data. The voids whose small amount may produce a significant drop in the velocities are microcracks. We assumed that the microcracks were randomly oriented. Even a small content of microcracks (within the experimental error of porosity) can cause a significant decrease in the elastic wave velocities. As a rule, it is not possible to measure the microcrack porosity and estimate their aspect ratio in the experiment. In this case, RPM can be used to estimate these parameters, including not only quasi-isometric voids but also microcracks in the model. Such a model will contain two unknown parameters—microcrack porosity and their aspect ratio. These parameters can be determined by minimizing the discrepancy between theoretical and experimental elastic wave velocities.

Note that microcracks not seen by the computed tomography were also incorporated in an RP model of carbonate rocks by Lima Neto et al. [34] in order to explain the experimental P-wave velocities. Fonseca [35] investigated thin-section images of carbonate rocks with confocal laser scanning microscopy. The porosity and permeability were independently measured for the rocks. It was found that numerical simulations based on image analysis captured the porosity but not the permeability. This may also indicate the presence of invisible voids in carbonate rocks.

When incorporating microcracks in the RP model, again, the question arises about the ambiguity of the simulation results due to the choice of the method. To assess this ambiguity, the elastic wave velocities were calculated for a model of a calcite matrix with randomly oriented cracks of the same shape filled with air. As before, Berryman’s self-consistent approximation, differential scheme (DEM), and Kuster–Toksöz methods were applied. The crack porosity varied from 0 to 0.5%. The aspect ratio of cracks was chosen to be 0.001, 0.01, and 0.1. As can be seen (Figure 9), up to a crack porosity of about 0.15%, the behavior of the velocities as a function of the porosity and even their values obtained with these three methods for all aspect ratios were quite similar. As the porosity increased, the Kuster–Toksöz method gave lower Vp and Vs velocities for an aspect ratio of 0.001. Figure 10 shows the same velocities, but depending on the crack density, which is determined by the formula $\mathsf{\epsilon }=\frac{3\mathsf{\phi }}{4\mathsf{\pi }\mathsf{\alpha }}$, where φ is the crack porosity in fractions of a unit and α is the aspect ratio of fractures. Up to a crack density of 0.2, the elastic wave velocities practically did not show a dependence on the aspect ratio for all the methods used. At the same time, the results obtained with these methods were also very close. Furthermore, for the modeling, we used only the Berryman self-consistent method. We applied it to invert the parameters of microcracks, namely, the crack porosity and aspect ratio, from the experimental data.

3.3. Rock-Physics-Based Inversion for Model with Microcracks

To solve the inverse problem of determining the microcrack porosity and aspect ratio for each sample, three sets of data (Dataset 1, Dataset 2, and Dataset 3), described in Section 3.2, were sequentially considered. The method of N-dimensional meshes described in Section 2.4 was applied for solving the inverse problem.

Figure 11 shows the surfaces of the Discrepancy measure calculated by formula (1) versus the microcrack porosity and aspect ratio obtained for Dataset 1 for each sample. The solutions for Datasets 2 and 3 look similar.

As seen from Figure 11, the range of solutions providing almost the same values of Discrepancy measure for combinations (microcrack porosity, aspect ratio) can be quite wide. Figure 12 shows only acceptable solutions, which can also be quite a lot. For each sample, the solutions are shown for three types of input data (Dataset 1, Dataset 2, and Dataset 3). Trying to find the “most probable” solution can also lead to a certain area of ambiguity. This follows from Figure 13 that shows solutions for the “conditional minimum” values of the Discrepancy measure. More than one and even many such solutions can be found, if we consider similar values of Discrepancy measure, which in the vicinity of the absolute minimum of this characteristic, satisfy a certain condition. Such a condition may be, for example, not exceeding a specified difference in the values of Discrepancy measure. Thus, the solutions shown in Figure 13 for some samples correspond to Discrepancy measure values coinciding up to the second decimal place. However, as follows from Figure 14, the solution becomes almost unambiguous when looking for the crack density instead of the pair (microcrack porosity, aspect ratio).

3.4. Rock-Physics-Based Inversion for Model with Intergranular Material

An example of such an inversion was only demonstrated for Sample 1 and Dataset 1, i.e., the porosity was obtained with the porosimeter, and the aspect ratio of isometric voids was 0.50. This sample had the greatest porosity among others. The formulas of the lower Hashin–Shtrikman bound [28,36] were applied to model the effective elastic moduli of the host matrix represented by calcite crystals surrounded by intergranular material. As mentioned above, this material was introduced to take into account the effect of grain boundaries and the possible imperfection of the crystalline structure in their vicinity. The SCA method was used to calculate the effective elastic properties of the host matrix with quasi-isometric pores. Figure 15 demonstrates the inverse problem solutions for the effective elastic moduli of the host matrix. As seen, the domain of solutions for the host matrix moduli was quite wide. Among the solutions, we selected only those that, at the next step of inversion, allowed us to obtain the result satisfying the following conditions: the Poisson’s ratio of the intergranular material should be positive and smaller than 0.5; the elastic moduli of the intergranular material should be smaller than those of the calcite polycrystal; the bulk modulus of the host material should be greater than that of water; and the total volume concentration of the intergranular material should not exceed 2%.

Commonly, these solutions exhibit a bulk modulus varying from 2.2 to 9 GPa and quite a small shear modulus of around 0.16–0.7 GPa. If we increased the possible volume content of intergranular material up to 5%, the range became wider for the bulk modulus, 8–19 GPa, and narrower for the shear modulus, 0.3–0.4 GPa. In both cases, the Poisson’s ratio was quite high, varying from 0.43 to 0.48. The following questions arise: Is this material realistic or not? Or is it only a formal result of inversion?

3.5. Comparison of Predictive Abilities of RP Models

Each of the constructed RP models was used to predict the elastic wave velocities in rocks saturated with formation water in the range of 0–20% quasi-isometric porosity. The bulk modulus and density of formation water were taken as equal to 2.82 GPa and 1.1 g/cm³ [28], respectively. The microcrack parameters were considered constant and corresponded to the most probable solutions found for each dataset. These parameters are shown in

Table 2. Microcrack parameter corresponding to the most probable solutions.

Name	Crack Porosity, %	Aspect Ratio	Crack Density
Sample 1, Dataset 1	0.01	0.0001	0.23
Sample 1, Dataset 2	0.85	0.0078	0.26
Sample 1, Dataset 3	0.02	0.0001	0.39
Sample 2, Dataset 1	0.07	0.0024	0.07
Sample 2, Dataset 2	0.23	0.0100	0.05
Sample 2, Dataset 3	0.25	0.0066	0.09
Sample 3, Dataset 1	0.10	0.0016	0.15
Sample 3, Dataset 2	0.15	0.0016	0.22
Sample 3, Dataset 3	0.03	0.0003	0.26

. For the fluid substitution, the Gassmann method was used [28,37], i.e., a prediction of velocities for the frequencies of the field experiment was assumed. Figure 16 and Figure 17 show the prediction of longitudinal and transverse wave velocities obtained for each sample as a function of the total porosity. The total porosity is equal to the sum of the quasi-isometric and microcrack porosity. The prediction was carried out for the microcrack parameters obtained for the Datasets 1–3 shown in Table 2.

As follows from the analysis of the results, the estimates of elastic wave velocities obtained for the same value of total porosity differed by no more than 14 and 16% for the Vp and Vs velocities, respectively. Moreover, the maximum discrepancy in the predicted velocities was observed for the sample with the highest porosity (Sample 1) and between the datasets exhibiting the maximum difference in the porosity values.

A similar inversion was carried out for an RP model containing only one type of voids, without dividing them into pores and microcracks. In this case, the data of tomography were not taken into account. For all samples, Dataset 1 was selected for the inversion of model parameters. It is interesting to note that for all three samples, the inverted aspect ratio of voids was almost the same and varied around 0.1, which was much less than the aspect ratios provided by the tomographic studies (0.50–0.55). As for the previous case, the inverse problem had more than one solution (on the order of 10–20). However, the values of the aspect ratio in all solutions were obtained in the vicinity of 0.1. For Sample 1, the solutions were obtained at a given error level of 5%. For Sample 2 and Sample 3, the error levels that made it possible to obtain solutions for the voids’ aspect ratio were 7.5 and 8%. These error values were comparable to those used for the model with two types of voids discussed above.

A comparison of the predicted velocities is shown in Figure 17. The single aspect ratio model gave greatly overestimated elastic wave velocities. The discrepancy in the estimates of velocities for the two models could reach 30%. Moreover, the maximum discrepancy was observed for the rock with the highest porosity (Sample 1). The use of such a “rough” model with one type of voids, not verified by tomography data, can lead to significant errors in the interpretation of field experiment data.

A similar prediction was carried out for two RP models in which the elastic moduli of the host matrix were assumed to be unknown and inverted from the elastic wave velocities. These models did not contain microcracks. The models were constructed for Sample 1. Dataset 1 was used for the inversion of the model parameters. For both RP models, the host matrix moduli were selected from the accepted solutions shown in Figure 15. These moduli (bulk and shear) were, respectively, 59.8 and 16.6 GPa for the first RP model (host 1) and 43.4 and 18.5 GPa for the second model (host 2). Below, we call the models host-models. The first model assumes that the intergranular material exists between the calcite grains. The second model does not consider any particular reasons for the decreased moduli of the host matrix compared to the calcite polycrystal. These elastic moduli of the matrix in host 2 model provide the most probable solution to the inverse problem.

The prediction results obtained for these two models are shown in Figure 18 for Sample 1 together with the results found for another two RP models. One of the models is the dual-porosity model with pores and microcracks and the other is an RP model (the “as simple as possible” model) with one system of voids (the same as in the upper-left plot in Figure 17). All the results were found for Dataset 1.

As seen, the choice of host matrix in the models host 1 and host 2 did not considerably affect the prediction results. Moreover, the difference between the velocities provided by the dual-porosity model and the host models did not exceed 10% in the whole porosity range. All the predictions based on the tomography significantly differed from those provided by the “as simple as possible” model (for a single aspect ratio for all voids).

In all the predictions described above, we ignored the existence of the micrite matrix in the samples due to the fact that the porosity was not too high. To analyze when the distribution of porosity between the micrite matrix and “visible” pores became noticeable, we performed the following modeling for Sample 1 and Sample 3. In the first model, we fixed the “visible” porosity and varied the micrite porosity from 0 to 40%. In the second model, we fixed the micrite porosity and varied the “visible” porosity in a wide range. The fixed values of “visible” and micrite porosity corresponded to those found for the sample. This means that the “visible” porosity was 7.4% for Sample 1 and 2.2% for Sample 3. The micrite porosity was, respectively, 6.7% and 8.1% for these samples. The total porosity range was up to 35%. The rocks were water-saturated.

A comparison of elastic wave velocities obtained with two models is shown in Figure 19. The maximum difference in velocities caused by the redistribution of porosity between the micritic matrix and “visible” pores was seen for a porosity of 35%. The difference rapidly increased with porosity and, finally, reached 54% for P-waves and 63% for S-waves, respectively, for Sample 3. This was the sample with the highest micrite content. The difference in the velocities became greater than 5% beginning from a porosity of 17%.

4. Discussion

As mentioned earlier, samples of different stratigraphic age were studied. Two samples (Sample 1 and Sample 3) were Carboniferous and Sample 2 was Devonian. The Lower Carboniferous samples had lower density and higher porosity than the Upper Devonian sample. This was due to the structure of the pore space. The pores in Sample 1 and Sample 3 were mainly represented by the secondary dissolute type of intergranular porosity, and in Sample 2, the separate cavernous pore space prevailed. According to the results of the RPM, the aspect ratio of microcracks in the Carboniferous samples was lower than in the Devonian sample.

The tomography studies carried out for the samples demonstrated that the pore space had a hierarchical structure. For two samples (Sample 2 and Sample 3), the tomography performed at different scales (resolutions of 12 microns and 37 microns) indicated a similar pore shape. For all three samples, significantly differing porosity was observed for these two scales. As a rule, the porosity obtained for big samples is much smaller. Thus, for Samples 1 and 3, it was less than 1% (0.58% and even 0.14%). Additionally, for these samples, the porosity measured with the porosimeter was greater, by 3 and 5%, than the porosity provided by tomography for the small scale. This means that a part of the connected porosity could be beyond the tomography resolution. Another possible reason for this discrepancy was that the pores were unevenly distributed within the studied samples. In contrast, for Sample 2, the tomography for small sample exhibited greater porosity as compared to that given by the porosimeter. This could indicate the presence of closed porosity, since the porosimeter measures connected porosity. Moreover, as follows from Figure 5, it also could be due to the uneven distribution of pores in this sample.

For each sample, the RPM was performed for three datasets depending on the tomography resolution and the method of obtaining porosity (measured with the porosimeter or provided by computed tomography). When constructing the RP models, the results of thin-section analysis and computed tomography were incorporated into the models. Specifically, first, the porous micrite matrix was taken into account. For the Carbonaceous samples (Sample 1 and Sample 3), the micrite matrix occupied a huge amount of the rock volume (60–70%). However, even for these samples, the results of RP modeling demonstrated that the samples’ porosity was too small for the model to be sensitive to the porosity distribution between the micrite matrix and the pores seen with computed tomography. This gave us an opportunity to use a simplified RPM without dividing the quasi-isometric porosity into these two parts.

For the RPM, we applied three different methods of EMT including SCA, DEM, and Kuster–Toksöz. We showed that for the range of quasi-spherical porosity of the studied samples (up to 11.5%) and even for higher porosity (up to 20%), all of these methods produce differing results within the experimental error. However, as seen from Figure 8, for greater porosity, the results may differ more markedly. The following question arises: which method gives a more accurate result? It is not so easy to answer this question since all the methods are approximate. Note that the formulas of these methods were derived under different assumptions of what a heterogeneous medium looks like. Thus, SCA assumes that all heterogeneities have an ellipsoidal shape and are embedded in an effective matrix. The DEM proposes a step-by-step scheme of embedding inclusions in the effective matrix obtained in the previous step. This model can produce good results in the case when inclusions differ in size, and inclusions of each next portion are larger than the previous ones. The Kuster–Toksöz method considers a model of a host matrix with isolated inclusions. When applying these methods in practice, the difference in the models should be taken into account. The calculation results should be verified on experimental data. However, the data may contain experimental errors. An interesting and helpful theoretical verification of these methods was presented by Germanovich and Dyskin [38,39]. The authors of these works compared the quadratic terms with respect to the concentration of inclusions in the virial expansion formulas [40] and SCA and DEM formulas. The results obtained by the authors allowed them to conclude that DEM was the most accurate approach.

The theoretical velocities obtained for dry samples are, in general, significantly higher than their experimental analogs. This result indicates that the carbonate rocks may also contain voids whose size is beyond the tomography resolution. The microcrack parameters (aspect ratio and microcrack porosity) were inverted from the measured Vp and Vs velocities. It was shown that a wide domain of solutions for pairs of these parameters existedfor each sample. However, if the microcracks were characterized by a single parameter, namely, the crack density, the solution was rather stable. Note that the existence of microcracks that were not seen with the tomography studies was debatable. However, their occurrence was quite realistic in Samples 1 and 2 that were characterized as packstones, according to previous lithological studies in this area [4,5].

It is interesting that the crack density found for Samples 1 and 3 had similar values that varied between 0.15 and 0.39 depending on the dataset used for inversion (see Table 2 and Figure 12). These values were much higher compared to Sample 2 for which the crack density varied from 0.05 to 0.09. The rocks represented by Samples 1 and 3 were of the Lower Carboniferous period whereas Sample 2 was rock from the Upper Devonian age. Does a relation exist between the microcrack parameters and the age of carbonate rocks? Or is this, perhaps, a coincidence? This is a debatable question to be answered in future works.

The microcrack parameters obtained for different datasets together with the pores’ aspect ratio provided by computed tomography were used for the RP-based prediction of elastic wave velocities (Vp and Vs) for rocks saturated with formation water. The prediction was performed for a porosity up to 20%. It was demonstrated that the predicted velocity values can differ within 16% depending on the dataset used for the microstructural parameters and porosity evaluation.

In our study, we also considered another RP model of the studied rocks which was called “as simple as possible”. This simple model did not take into account either the data of thin-section analysis or the computed tomography data. This model contained only one type of void with the same aspect ratio. The aspect ratio was inverted only on the basis of elastic wave velocities and porosity measured using a porosimeter. Note, that this type of pore shape evaluation is often met in practice. Our modeling showed that all the samples had voids of a similar aspect ratio of around 0.1. Note that this value has nothing in common with those provided by the computed tomography for pores. Moreover, the elastic wave velocities predicted by this simplest model for water-saturated rocks may differ from the respective values found with the model containing pores and microcracks by 30%. Moreover, the difference is rather high within the whole porosity range. This may cause significant errors in the interpretation of field data.

In addition to the RP model with pores and microcracks, two alternative models were considered. These models assumed that the elastic moduli of the host matrix were different from those of the calcite polycrystal. The reason for the difference can be attributed to various factors including the effect of grain boundaries. These models contained only one set of voids, namely, quasi-spherical pores visible with computed tomography. The models differed from each other only in the elastic moduli of the host matrix, which were inverted from the measured elastic wave velocities. Both the host matrices provided an acceptable solution to the inverse problem. The choice of the host matrix moduli did not have a significant effect on the prediction. It was shown that the velocities predicted by these models for water-saturated rock within the porosity range up to 20% did not significantly differ from those of the dual-porosity RP model. The difference was within 10% for both velocities Vp and Vs.

The following question arises: what is the benefit of incorporating the computed tomography results in RP modeling, irrespective of the fact that a large part of voids and possibly other heterogeneities are beyond the tomography resolution? The answer is as follows. These results give evidence for the existence of quasi-isometric pores and provide information on their shape. This helps to separate their effect from the other factors that control the effective elastic properties. This gives an opportunity to make the predictions based on the constructed RPM more reliable even without direct knowledge of other factors that control the effective elastic properties.

5. Conclusions

The porosity of carbonate rocks has a complicated hierarchical structure. Moreover, the porosity can be unevenly distributed within the rock. This leads to the phenomenon that the porosity and pore shape estimated with the use of different methods and at different scales may be significantly different. The most reliable estimation of the porosity can be obtained with direct measurements of porosity (for example with a porosimeter, as in our study). However, even in this case, the closed porosity may remain out of consideration.

A rock-physics-based estimation of elastic wave velocities from the pore shape provided by the computed tomography at a resolution of around 12–37 microns can lead to greatly overestimated velocities, even for porosity measured with a porosimeter. The difference can be attributed to the existence of other factors affecting the effective elastic properties, which could not be detected with computed tomography. In particular, microcracks can be among such factors.

The thin-section analysis and computed tomography provide significant information on the rock fabric including the shape of visible pores. Combining this information with data on the elastic wave velocities and rock physics modeling allows one to reveal the parameters of voids whose size is beyond the tomography resolution, namely, microcracks. Moreover, if, for some reason, other types of RP models are used (without hidden voids), the available information on visible pores allows one to more reliably invert the elastic moduli of the host matrix affected by various factors which are as-yet unknown.

The information on rock microstructure provided by both thin-section analysis and computed tomography, and then included in rock physics models, may greatly enhance the predictive power and reliability of the models.