The petrophysical properties of rock are closely related to its microstructure. Macroscopic properties or behaviors are essentially manifestations of the microstructure, or macroscopic properties depend on the microstructure. When studying the petrophysical properties of heterogeneous materials, Torquato^[1] showed that the macroscopic properties of materials such as elasticity and permeability are affected by microstructure. In micromechanics, the theoretical methods for estimating the elastic properties of rocks generally include analytical modeling and numerical simulation methods. The numerical methods use field emission scanning electron microscopy (FE-SEM) or X-ray tomography to acquire the microstructural information, then use finite element method to mesh the rock microstructure, simulate the loading process, calculate the stress and strain field, and finally infer the elastic properties. A representative work conducted by Arns et al.^[2] used the finite element method to calculate the elastic properties of the digital core of Fountainebleau sandstone. The results showed that the law of the elastic properties changes accords with the Gassmann theory. Numerical simulation methods can accurately reproduce the micro-pore structure of the rock, but meshing is very difficult to accurately fit the pore structure, and the computation is very time consuming. Analytical modeling methods are generally based on Eshelby's equivalent medium theory, and the pores in rocks are approximated by spherical, spheroidal, cylindrical, flat cracks, and hard hollow spherical shells, then the rock effective elastic properties can be calculated analytically. For example, Zimmerman^[3] assumed that the pores in the rock were spheroids, and can be divided into spherical pores, needle-shaped pores, and coin-shaped pores according to different pore aspect ratios. According to linear elastic Hooke's law and Eshelby's equivalent medium theory, Zimmerman also derived microscale rock deformation theory. Compared with the numerical simulation method, the analytical method can greatly save computation costs. However, due to the lack of detailed description of real porous microgeometry, analytical methods are mostly used for theoretical qualitative research.

In recent years, digital rock physics (DRP) was established and has been extensively used to study reservoir rock microstructure, effective properties and fluid transport mechanism^[4,5,6,7,8]. Thanks to its low requirements on core sample size, and high repeatability of experiment and simulation, DRP has become a popular new means to study the reservoir properties and seepage mechanism of unconventional oil and gas at micro-nano scale^[9,10]. Although DRP technology has been developed to be a mature research method and theoretical system in recent years, compared with real reservoir conditions and macroscopic research methods^[11], DRP technology generally does not consider the influence of in-situ stress on microstructure and percolation, and neglects the stress sensitivity, which makes DRP technology still have certain defects in practical application. This paper attempts to study the relationship between the basic elastic properties of the digital core (pore volume compressibility) and the micro-pore structure to provide theoretical support for the micro-deformation of the digital core. In reservoir engineering, the pore volume compressibility is not only an important input parameter in the material balance calculation^[12], but also an important parameter in the diffusivity equation of pressure in percolation mechanics and well test analysis^[13].

Based on the real porous microgeometry of ten representative digital cores including Bentheimer sandstone and so on, this paper analyzes the characteristics of the pore shape parameters, establishes the abstract spheroidal pore model, combines the analytical method and the equivalent medium theory to study the constitutive relationship of microstructure deformation. The effective elastic modulus and pore volume compressibility are solved and the effects of pore structure on the pore volume compressibility are investigated.

The spheroidal pore modeling method was firstly presented by Sadowsky et al.^[14,15], since most pores can be identically treated as spheroids such as sphere, needle-shape cylinder, and cracks etc. For an ellipse with a semi-major axis length of A and a semi-minor axis length of B (Fig. 1a), rotating about its major axis yields a prolate ellipsoid or named ‘needle-shape pore’ (Fig. 1b), while rotating about its minor axis forms an oblate ellipsoid or named ‘coin-shape pore’ (Fig. 1c).

The spheroid aspect ratio ($\rho $) is defined as the ratio of the polar radius and the equatorial radius. A prolate ellipsoid has an aspect ratio larger than 1 ($\rho >1$), and an oblate ellipsoid has an aspect ratio smaller than 1 ($\rho <1$). By solving the elastic problems in spheroidal coordinate systems, we can derive the equation of pore volume compressibility ^[3].

For a microscale pore itself, Zimmerman believes that the external stress exerted on porous medium can be treated as stress from infinity, thus only normal stress needs to be considered, so the pore volume strain is calculated by ^[3]

For a dry hole condition, pore pressure is constant, the above equation can be simplified as

The pore volume compressibility refers to the fractional change in the pore volume of the rock with a unit change in pressure. The pore volume compressibility C_pp can be used to characterize the deformation scale of the pores in reservoir rocks. In the following section, we first derive the pore volume compressibility under the condition of single-pore deformation including prolate ellipsoidal, oblate ellipsoidal and spherical pore deformation, and then extend it to the multi-pore deformation case.

1.1.1. Prolate and oblate ellipsoidal pore deformation

According to Zimmerman’s single-pore elastic deformation theory, the compressive deformation of a single prolate ellipsoidal pore in infinite, isotropic elastic medium can be solved by coordinate transform. The Cartesian coordinate system can be established first with the center of the prolate ellipsoidal pore as the origin, and the z-axis as the polar axis, then convert it into a prolate ellipsoidal coordinate system (Fig. 2).

where l is the half focal length on the polar axis (z axis), l is taken as unit length in the following derivation. The values of $\alpha $, β, and γ range from 0 to ∞, 0 to π, and 0 to 2π respectively.

When α is a constant, the equation represents the surface of the prolate ellipsoid; when β is a constant, the equation represents the hyperboloid; when γ is a constant, the equation represents a plane passing through the z axis. Therefore, the cavity surface is $\alpha ={{\alpha }_{0}}$, the equatorial radius is $\sinh {{\alpha }_{0}}$, the polar radius is $\cosh {{\alpha }_{0}}$, and the aspect ratio is $\rho =\coth {{\alpha }_{0}}$.

It is now convenient to introduce new variables defined by

and the differential arclength in the orthogonal coordinate system has the form of

where the metric coefficient h_i($i=\alpha ,\beta ,\gamma $) is given by

Since the cavity surface is under uniform pressure of unit magnitude, and all stresses must vanish infinitely far from the cavity, the boundary conditions for the problem are

The Boussinesq three-function approach can be used for solving the equations of elastic equilibrium subject to the above boundary conditions, as presented by Sadowsky et al.^[14] The pore volume change is found by integrating the normal component of the displacement over the entire surface of the cavity, then the pore volume compressibility can be solved.

Since the stress variation on the pore was Δσ, the pore volume compressibility is expressed as

where

Substituting equations (8) and (10) into equation (9) lead to the final expression for C_pp:

where

Since the oblate ellipsoidal coordinate system can be converted from the prolate one, the pore volume compressibility of oblate ellipsoids can be obtained directly from equation (11):

where

The volume strain for prolate and oblate ellipsoids are respectively given by

1.1.2. Spherical pore deformation

The volume of a spherical pore is given by ${{V}_{\text{p}}}=\frac{4}{3}\text{ }\!\!\pi\!\!\text{ }{{r}^{3}}$. When the hydrostatic stress dσ is applied at infinity, the radial strain of the spherical pore is calculated by^[16]

and the volume strain of the spherical pore is

The pore volume compressibility of spherical pore can be solved from equation (18):

The compressibility of a single pore of different shapes was discussed above, which was characterized by the volume change of a single pore under a unit pressure difference. The pore volume change of the entire pore system is the sum of the changes of each independent pore volume. However, in real rocks, the stress field of each pore will be affected by surrounding pores, and the existence of surrounding pores will increase the pore compressibility of individual pore.

In the elastic deformation analysis of real porous materials (such as rocks), considering that the macro-mechanical properties (or transport properties) of rocks are closely related to their microstructures, the equivalent medium theory proposed by Eshelby^[17] is often used to solve the effective properties of macroscopic materials from the microstructure calculation. In this paper, the rock is equivalent to a mixture with inclusions in an infinite matrix, that is, the pores are regarded as inclusions, and the grains are regarded as the matrix. The effective elastic properties of rocks are related with the elastic modulus, volume fraction, microgeometry and spatial distribution of pore and matrix. Before solving the pore volume compressibility under multi-pore elastic deformation conditions, the effective modulus of the rock under the multi-pore conditions needs to be solved first.

1.2.1. Estimation of effective moduli

The effective moduli under the multi-pore elastic deformation conditions can be obtained by the self-consistent method, which was first proposed by Hill^[18] and Budiansky^[19] and is the most commonly used equivalent medium method. The self-consistent method treats the inclusions (pores) in the composite material as spheres, and then embeds them in a medium with unknown effective modulus to calculate the effective properties. Wu^[20] considered inclusions as spheroids and proposed a self-consistent model different from spherical inclusions.

This paper considers inclusions (pores) of arbitrary shapes and properties in an infinitely large matrix. The properties of this matrix containing inclusions are equal to the unknown effective properties of the entire heterogeneous material. Assuming that both the inclusion and matrix phases are isotropic, and considering that the spherical inclusions are randomly distributed, the effective bulk modulus and effective shear modulus are calculated as^[21]:

The above formulae are implicit equations and the calculation process requires numerical iteration. For the case of non-spherical inclusions, using the self-consistent model proposed by Wu^[20] and treating the inclusions as spheroids, the modulus of the two-phase mixture can be calculated by:

1.2.2. Estimation of pore volume compressibility

Use the effective shear modulus $G_{\text{SC}}^{*}$ calculated from equation (23) to replace G and substitute into equations (11) and (13), we obtain the calculation formula for the pore volume compressibility of porous media with multiple prolate and oblate ellipsoidal pores:

Considering that the real reservoir rocks all contain various pore shapes, the final pore volume compressibility can be calculated according to the volume proportions of pores with different aspect ratio distributions. Let the core sample pore aspect ratio distribution range be $\left[ {{\rho }_{1}},{{\rho }_{N}} \right]$, where the pore volume with aspect ratio ${{\rho }_{k}}$ is ${{V}_{k}}$ and the pore volume proportion is ${{x}_{k}}$, then:

The pore volume compressibility for mixed multi-pore case is given by

The digital core samples selected in this paper includes high, medium, and low permeability sandstones and some carbonate rocks. Ten digital core samples are Bentheimer sandstone, Doddington sandstone, Berea sandstone, Fountainebleau sandstone, Wilcox tight sandstone, Estaillades carbonate rock, Daqing sandstone, Tarim sandstone, sandstone of Eastern South China Sea and Xinjiang sandstone. The first six digital core samples are published benchmark digital core data ^[22], and the other four samples are obtained by CT scan using either downhole cores (Daqing and Tarim) or drill cuttings (Eastern South China Sea and Xinjiang). The workflow of establishing spheroidal pore model for Berea sandstone will be introduced here as an example.

The benchmark Berea sandstone digital core used in this study was published by the PERM research group at Imperial College of Technology^[23]. Totally 1024 two-dimensional grayscale images of 1024×1024 pixels were obtained by CT scanning experiment (resolution was set to 2.77 μm). Three- dimensional grayscale image was constructed from two-dimensional slices, then watershed algorithm was applied to extract pore space and convert the data into binary images.

Avizo image processing software was used to separate the pore space and label the individual pores (Fig. 3, different colors are used to distinguish the individual pores). The principle of moment of inertia is applied for parameter analysis on all labeled pores. For each pore, the principal axis length, center of gravity position, and Euler angle of the second standard central moment spheroids are calculated. Then the real pores can be approximated by the spheroids with the same characteristic parameters, and the real pore space is converted into the spheroidal pore model (Fig. 4, different colors are used to distinguish individual pores).

From the spheroidal pore model of Berea Sandstone, the pore aspect ratio distribution was generated (Fig. 5). From Fig. 5, one can observe that the range of the pore aspect ratio in Berea sandstone is 0.2 to 4.9. According to equation (27), the pore volume compressibility of the Berea sandstone can be calculated by using the pore volume proportions. For the rest of 9 digital core samples, the relationships between the pore aspect ratio and the pore volume proportion are also required to calculate the pore volume compressibility.

The spheroidal pore modeling method proposed in this paper is applicable to most sandstones and carbonate rocks except oolitic limestone. The main reason is that the premise of this method is to assume that most pores in the core are spatial convex bodies, so the realistic pores can be approximated by spheroidal pores to form an abstract pore model. However, for Ketton limestone, one of the oolitic limestones, the SEM pictures (Fig. 6a) and the three-dimensional porous network obtained by CT scan (Fig. 6b) show that most diagenetic particles are spheres or spheroids, and the inter-particle pore spaces are obviously concave, so the proposed method is not applicable for Ketton limestone.

Based on the proposed spheroidal pore modeling approach, the effect of microstructure on pore volume compressibility were firstly studied in two synthetic cases: single-type pores and mixed-type pores.

Here we consider that there are only single-type pores in the core. To cater for the reality of multi-pore existence and various pore aspect ratios in the core, it is assumed that the aspect ratio of prolate ellipsoidal pores and oblate ellipsoidal pores conform to the Gaussian distribution. The bulk modulus of the rock particles (quartz) is 37 GPa, the shear modulus is 44 GPa, and the porosity is 20%. On this basis, the effect of pore aspect ratio on pore compressibility was studied.

First consider the case of prolate ellipsoidal pores, assuming that the range of pore aspect ratio mean is 6.00-20.00, the variance is 0.50-3.00, and the porosity is 20%, we could obtain the relationship between pore volume compressibility and the pore aspect ratio whose distribution is featured by different mean and variance (Fig. 7). Fig. 7 indicates that: (a) when pore aspect ratio of prolate ellipsoidal pores conforms to the Gaussian distribution, pore compressibility increases with increasing pore aspect ratio, and increases faster when pore aspect ratio is less than 10.00. (b) the smaller the variance, the larger the pore compressibility. Under the same aspect ratio mean value, a concentrated aspect ratio distribution leads to a larger pore compressibility and a larger deformation under the same stress.

Similarly, assuming that the aspect ratio of oblate ellipsoidal pores conform to the Gaussian distribution with a mean value of 0.10-0.80, a variance of 0.04-0.20, and the porosity of 20%, we could obtain the relationship between pore volume compressibility and the pore aspect ratio whose distribution is featured by different mean and variance (Fig. 8). Fig. 8 indicates that (a) when the aspect ratio of the oblate ellipsoidal pores conforms to the Gaussian distribution, the pore compressibility decreases with increasing aspect ratio, and the smaller the variance, the faster the decrease. (b) same as the prolate ellipsoidal pore case, the smaller the variance, the larger the pore compressibility, i.e. a concentrated aspect ratio distribution leads to a larger pore compressibility and a larger deformation under the same stress. (c) For a mean value of aspect ratio greater than 0.30, the pore compressibility is not influenced by the variance longer.

Single-type pores rarely exist in real rocks in which mixed-type pores are the common case. Equation (27) can be used for calculating pore volume compressibility for the rock with mixed-type pores. The effects of porosity and the volume proportion of different types of pore on the overall pore compressibility are discussed here.

3.2.1. Effect of the volume proportion of different types of pore

Assuming the porosity is 20%, the aspect ratio of prolate ellipsoidal pores is 2.00, and the aspect ratio of oblate ellipsoidal pores is 0.10, we could calculate the dependence of pore compressibility on the volume proportion of different types of pores (prolate ellipsoid, oblate ellipsoid, and sphere) (Fig. 9). From Fig. 9, we can see that (a) when three types of pores are mixed, pore compressibility has a positive correlation with the volume proportion of the oblate ellipsoidal pores, while has a negative correlation with the volume proportion of spherical and prolate ellipsoidal pores. (b) when the prolate ellipsoidal pores are mixed with spherical pores only, pore compressibility increases with increasing spherical pores. The results indicate that under a certain aspect ratio, rocks with more prolate ellipsoidal pores are less likely to deform, while with more oblate ellipsoidal and spherical pores are more likely to deform, and the oblate ellipsoidal pores has a greater effect than spherical pores.

3.2.2. Effect of porosity

To study the effect of porosity, we assume the aspect ratios of prolate and oblate ellipsoidal pores are 2.00 and 0.50 respectively. Consider 4 mixing modes for different types of pores: (a) the mixing ratio of the oblate and prolate ellipsoidal pores is 1:1; (b) the mixing ratio of the oblate and spherical pores is 1:1; (c) the mixing ratio of prolate ellipsoidal and spherical pores is 1:1; (d) the mixing ratio of prolate ellipsoidal pores, oblate ellipsoidal pores, and spherical pores is 1:1:1, the relationship between pore compressibility and porosity can be obtained according to equation (27) (Fig. 10). Fig. 10 indicates that the pore compressibility increases with increasing porosity, but the compressibility value and its growth rate are greatly different under different mixing modes: (a) For three cases containing oblate ellipsoidal pores, under the same porosity, the pore compressibility is much larger than the mixing mode of prolate ellipsoidal and spherical pores, and the larger the porosity, the more obvious the difference, indicating that rocks containing oblate ellipsoidal pores are more likely to deform, thus the pore compressibility is larger. (b) Compared with the mixing modes of prolate-oblate ellipsoidal pores and spherical-oblate ellipsoidal pores, the coexistence of three types of pores has a smaller pore compressibility, indicating that the coexistence of three types of pore will increase the rock compression strength. (c) The mixing mode of spherical-oblate ellipsoidal pores has a slightly larger pore compressibility than that of prolate-oblate ellipsoidal pores, indicating that spherical pores are more likely to deform than prolate ellipsoidal pores.

The pore space of ten real digital core samples except Berea sandstone is shown in Fig. 11 (the colors used to distinguish individual pores). By applying the methodology presented above, spheroidal pore models were established for each digital core sample. The pore volume proportions with varying aspect ratios were calculated (Table 1). From Table 1 we can see that the aspect ratios of prolate and oblate ellipsoidal pores in the ten samples are (1.11-12.75) and (0.11-0.90) respectively. Since the pore volume compressibility is calculated based on the volume proportions of different types of pores, the volume proportions of different types of pores are more important than the numbers of them. According to the theoretical analysis results above, the volume proportions of different types of pores, and the range, mean value and variance of pore aspect ratios have greater impact on pore compressibility.

Pore volume compressibilites of the ten digital core samples were estimated based on the presented theoretical models and the pore characteristic parameters listed in Table 1. Fig. 12 indicates that the pore volume compressibilities are generally increase with increasing porosity. Among the three samples with larger porosity (No. 1-3 in the Table), the pore volume compressibilities of Bentheimer and Doddington sandstone samples are mainly affected by macro porosity, and Doddington sandstone has a larger pore volume compressibility due to larger porosity. The Estaillades sample has the same porosity as the Doddington sandstone, but the Estaillades sample has more spherical and prolate ellipsoidal pores, so the pore compressibility is relatively small. Among the five samples with medium porosity (No. 4, 7-10 in the Table), the Berea sandstone sample has more spherical pores, lower mean and smaller variance of prolate ellipsoidal pore aspect ratio, higher mean of oblate ellipsoidal pore aspect ratio, which leads to a smaller pore volume compressibility. Among the two samples with smaller porosity (No. 5 and 6), although their pore compressibilities are not much different, the changing trend is slightly different from other samples. Compared with the Fountainbleau sandstone, the Wicox sandstone has more prolate ellipsoidal pores and less oblate ellipsoidal pores in volume proportion. If only analyzed from the perspective of volume proportion, the Fountanbleau sandstone should have a larger pore compressibility, but this is different from the calculation results of the model. Through in-depth analysis of the data, it is found that the results are related to the minimum value of the aspect ratio of oblate spheroidal pores in two ultra-low porosity samples. Although the Wicox sandstone has less oblate spheroidal pores, there are more oblate spheroidal pores whose aspect ratios are within the range between 0.1 and 0.2, while the minimum aspect ratio of the oblate spheroidal pores in the Fountainbleau sandstone is 0.22. This difference ultimately resulted in a slightly higher pore compressibility for the Wicox sandstone. This influencing factor will be analyzed more systematically in future research.

The pore volume compressibility of digital core is related to porosity, pore aspect ratio and volume proportion of different types of pores: (1) The pore compressibility of prolate ellipsoidal pores is positively related to the pore aspect ratio, while the oblate ellipsoidal pores have the opposite situation. (2) If the pore aspect ratio meets the Gaussian distribution for prolate and oblate ellipsoidal pores, with the same mean values, the more concentrated the aspect ratio, the larger the pore compressibility, and the larger the deformation under the same stress conditions. (3) The pore compressibility increases with increasing porosity. (4) When mixed-type pores are coexisting in the core, the pore compressibility is related to the volume proportion of different types of pores. For a determined porosity, the more oblate ellipsoidal and spherical pores, the easier the rock is to deform, and the larger the pore compressibility; the more prolate ellipsoidal pores, the more difficult the rock is to deform, and the smaller the pore compressibility.

The spheroidal pore modeling method is used to estimate the pore volume compressibility of 10 typical digital core samples. The obtained results are consistent with the rules described in the theoretical model. This method is validated to estimate the pore volume compressibility of other digital core samples.

A—semi-major axis of the ellipse, m;

A_i—constant term when solving the displacement integration, dimensionless;

B—semi-minor axis of the ellipse, m;

C_pp^*—mixing pore volume compressibility, Pa^-1;

C_pp—pore volume compressibility for single pore, Pa^-1;

C_pp,ob—pore volume compressibility for oblate ellipsoidal pores, Pa^-1;

C_{pp, pro}—pore volume compressibility for prolate ellipsoidal pores, Pa^-1;

C_pp,s—pore volume compressibility for spherical pores, Pa^-1;

$C_{\text{pp, }k}^{*}$—pore volume compressibility for the pores with aspect ratio of $\rho_{k}$, Pa^-1;

$C_{\text{pp, ob}}^{*}$—pore volume compressibility of porous medium containing oblate ellipsoidal pores, Pa^-1;

$C_{\text{pp, pro}}^{*}$—pore volume compressibility of porous medium containing prolate ellipsoidal pores, Pa^-1;

G—shear modulus, Pa;

G₀—matrix shear modulus, Pa;

G₁—pore shear modulus, Pa;

$G_{\text{SC}}^{*}$—effective shear modulus, Pa;

h, ${{h}_{\alpha }}$, ${{h}_{\beta }}$, ${{h}_{\gamma }}$—metric coefficient, m^-1;

i, j—axis serial number in spheroidal coordinate system;

K₀—matrix bulk modulus, Pa;

K₁—pore bulk modulus, Pa;

$K_{\text{SC}}^{*}$—effective bulk modulus, Pa;

l—semi-focal length on the major axis of prolate ellipsoid, m;

N—number of ranges of pore aspect ratio;

p—pressure, Pa;

$P$, $\overline{P}$—intermediate variables in spheroidal coordinate system, m;

$Q$, $\overline{Q}$—intermediate variables in spheroidal coordinate system, m;

${{Q}_{0}}$—semi-major axis length in spheroidal coordinate system, m;

$\overline{{{Q}_{0}}}$—semi-minor axis length in spheroidal coordinate system, m;

r—radius of spherical pore, m;

R_pro, R_ob—intermediate variable, dimensionless;

s—arc length in spheroidal coordinate system, m;

S_G—pore geometry factor for calculating effective shear modulus, dimensionless;

S_K—pore geometry factor for calculating effective bulk modulus, dimensionless;

U—displacement, m;

V_k—pore volume with the aspect ratio of ${{\rho }_{k}}$, m³;

V_p—pore volume, m³;

${{x}_{k}}$—pore volume proportion with the aspect ratio of ${{\rho }_{k}}$, dimensionless;

x, y, z—axis in cartesian coordinate system, m;

$\alpha $, $\beta $, $\gamma $—axis in spheroidal coordinate system, degree;

${{\alpha }_{\text{0}}}$—$\alpha $ values at the cavity surface in spheroidal coordinate system, degree;

${{\alpha }_{\text{B}}}$—Biot coefficient, dimensionless;

${{\text{ }\!\!\delta\!\!\text{ }}_{m,n}}$—Kronecker function, 1 for m=n, 0 for m≠n, dimensionless;

ε_p—pore volume strain, dimensionless;

ε_p,pro, ε_p,ob, ε_{p, s}—pore volume strains of prolate ellipsoid, oblate ellipsoid and sphere, dimensionless;

$\xi $—pore surface area, m²;

$\rho $—aspect ratio, dimensionless;

${{\rho }_{1}}$—minimum aspect ratio, dimensionless;

${{\rho }_{k}}$—an aspect ratio value within the distribution range, dimensionless;

${{\rho }_{N}}$—maximum aspect ratio, dimensionless;

${{\rho }_{\text{ob}}}$—aspect ratio of oblate ellipsoid pores, dimensionless;

${{\rho }_{\text{pro}}}$—aspect ratio of prolate ellipsoid pores, dimensionless;

$\sigma $—normal stress, Pa;

$\Delta \sigma $—stress variation, Pa;

$\upsilon $—Poisson’s ratio, dimensionless;

$\phi $—porosity, %.