1 INTRODUCTION

The elastic properties of the rock depend on the microscopic structure of pore system significantly. Most rocks usually have two or even more than two different pore types, such as pore, crack, cavity, etc.; the complex pore system makes the relationship between the velocity and porosity of the rock highly scattered(Baechle et al., 2008; Sayers, 2008). Therefore, it needs to establish more accurate multiple-porosity rock physical models to characterize the variations of elastic moduli with porosity for porous rock.

Effective medium theory, such as Kuster-Toksöz theory and differential effective medium theory, has been used to study the elastic properties of porous rock. Kuster and Toksöz(1974)derived expressions for bulk and shear moduli of multiple-porosity rock by using wave scattering theory, in which the effects of elasticity, volume content and pore shape of inclusions are taken into account. But Kuster-Toksöz theory does not consider the interaction between different pores, and requires that the ratio of porosity to pore aspect ratio has to be much smaller than 1. Kuster-Toksöz theory has been shown to violate the upper and lower Hashin-Shtrikman bounds in some cases(Berryman, 1980). The differential effective medium(DEM)theory models two-phase composites by incrementally adding inclusions to the matrix phase, and is applied to the determination of the effective elastic properties of porous rocks that are dry or saturated by fluid through the numerical solution of differential equations(Berryman, 1980; Norris, 1985; Zimmerman, 1985; Berryman, et al., 2002). The DEM has the property that it can never violate rigorous bounds compared with the Kuster and Toksöz theory(Berryman, 1980; Berryman and Berge, 1996; Berryman et al., 2002). However, since the ordinary differential equations for bulk and shear moduli are coupled, it is more difficult to integrate them to yield accurate analytical formulae for the bulk and shear moduli. Therefore, the behavior of the bulk and shear moduli can only be accurately simulated through the numerical solutions of differential equations. Moreover, the effective elasticities of porous medium depend on the adding order of the inclusions, i.e. the pores or cracks with different aspect ratios.

Xu and White(1995)developed a double-porosity(stiff and soft pores)model(now is called Xu-White model)for shaly s and stones based on Kuster-Toksöz theory, which assumed that the volume contents of different pores depend on the s and and clay content respectively. But the model inherited the problems existing in Kuster-Toksöz theory, and had the computational efficiency problem for its iterative algorithm. Keys and Xu(2002)simplified the Xu-White model and obtained an approximation for dry-rock bulk and shear moduli by assuming a constant dry-rock Poisson’s ratio in coupled equations, which improved the computational efficiency. However the Poisson’s ratio or modulus ratio derived from their dry-rock approximation is not constant as they presupposed. Xu and Payne(2009)promoted the Xu-White model for shaly s and stones to carbonates and established the Xu-Payne model, but the model also had the same problems as the Xu-White model.

To solve the problem that DEM theory lacks analytical formulae, many researchers have made some efforts and achieved a series of achievements. Zimmerman(1991)derived the implicit analytical results for the elastic moduli of dry rock with spherical pores from DEM theory. Berryman et al.(2002)obtained approximate analytical results for the elastic moduli of dry and fully saturated cracked rock based on DEM theory. Li and Zhang(2010, 2012)derived analytical solutions of the bulk and shear moduli for dry rock with the three specific pore shapes: spherical pores, needle-shaped pores and penny-shaped cracks in order to decouple the differential equations by applying an analytical approximation for dry-rock modulus ratio. Similarly, Li and Zhang(2011)obtained analytical formulae for the bulk and shear moduli of dry rocks with ellipsoidal pores based on DEM theory. Those DEM analytical formulae were proposed based on the single-porosity system, and were used for the effective pore aspect ratio inversion and shear wave velocity prediction(Li et al., 2013).

The purpose of this paper is to obtain analytical formulae for the bulk and shear moduli of multiple-porosity dry rocks with ellipsoidal pores based on effective medium theory and use them to characterize the dry-rock moduli dependency on porosity. We first derive the differential schemes for Kuster-Toksöz theory from the Kuster-Toksöz expressions for bulk and shear moduli. Next we derive analytical solutions of the bulk and shear moduli for multiple-porosity dry rock with ellipsoidal pores based on DEM theory by applying an analytical approximation for dry-rock modulus ratio, in order to decouple the equations for bulk and shear moduli. And then we use the theoretical model for checking these approximations by comparison with numerical computations for the differential equations of the DEM, and the agreement between the analytical approximations and the full DEM is discussed. Finally the analytical approximations are applied to predict the variations of the bulk and shear moduli with respect to porosity for dry s and stone data, which were measured by Han et al.(1986)in the laboratory.

2 DIFFERENTIAL SCHEMES FOR KUSTER-TOKSÖZ THEORY AND THEIR NUMERICAL SOLUTION

By using a long-wavelength first-order scattering theory, Kuster and Toksöz(1974)derived the expressions for the effective bulk modulus, K_KT^*, and shear modulus, G_KT^* of rock with N inclusions

where K_m and G_m are the bulk and shear moduli of matrix, respectively; v_i is the volume percentage of the ith inclusion in porosity φ, K_i and G_i are the bulk and shear moduli of the ith inclusion, the scalar factors P^mi and Q^mi are the so-called polarization factors which describe the effect of the ith inclusion in a background medium m, ${\xi _m} = \frac{{{G_m}(9{K_m} + 8{G_m})}}{{8({K_m} + 2{G_m})}}$, and the superscript * denotes the effective medium.

Kuster-Toksöz theory and differential effective medium theory can be combined to calculate the elastic properties of multiple-porosity rock(Xu and White, 1995; Keys and Xu, 2002). Differential effective medium theory takes the point of view that a composite material may be constructed by making infinitesimal changes in an already existing composite(Norris, 1985; Berryman et al., 2002; Mavko et al., 1998). Therefore, when a tiny porosity dy is added into the composite material whose porosity is y, K_m and G_m become K^*(y) and G^*(y) before adding, respectively, K_KT^* and G_KT^* become K^*(y + dy) and G^*(y + dy) after adding, respectively, P_*i and Q_*i become P^*i and Q^*i which describe the effect of the ith inclusion in an effective medium *, in Eq.(1). So Eq.(1)can be rewritten as

Because dy [ \ll \] 1 in Eq.(2), we have $\frac{{3{K^*}(y)+ 4{G^*}(y)}}{{3{K^*}(y + dy)+ 4{G^*}(y)}} \approx 1$, $\frac{{3{G^*}(y)+ 4{\xi ^*}(y)}}{{3{G^*}(y + dy)+ 4{\xi ^*}(y)}} \approx 1$ Letting dK^*(y) = K^*(y + dy)- K^*(y), dG^*(y)= G^*(y + dy)- G^*(y) and dividing both sides of Eq.(2)by the dy, we obtain

Eq.(3)is the differential scheme of Kuster-Toksöz theory. Zimmerman(1984)proposed the similar differential scheme for rock with single spherical pores, which assumes that when a new inclusion is introduced, it always replaces the host matrix material(Zimmerman, 1984; Mavko et al., 1998), and Eq.(3)is called the Zimmermann’s version for DEM theory. While Norris(1985)assumed that it replaces on average either the host matrix material or the inclusion material with probabilities 1 - y and y(Norris, 1985; Mavko et al., 1998), so dy is replaced by dy/(1 - y)in Eq.(3)(Mavko et al., 1998), Eq.(3)becomes The equation above is the Norris’s version of DEM theory for multiple-porosity rock, when N = 1, it is the same as the traditional DEM theory(Norris, 1985; Mavko et al., 1998; Li and Zhang, 2010). These equations are typically integrated starting from porosity y = 0 with values K^*(0)= K_m and G^*(0)= G_m, which are assumed here to be the mineral moduli values for the single homogeneous solid constituent. Integration then proceeds from y = 0 to the desired highest value y = φ.

In general the effective differential Eqs.(4a) and (4b)are also coupled, as these two equations depend on both the bulk and shear moduli of the composite. Therefore, it is difficult to obtain the analytical formulae of the bulk and shear moduli for dry rock by integrating the DEM equations.

First we give a comparison between the analytical formulae and numerical results of the differential Eqs.(3) and (4). We take quartz as the host medium, let K_m = 37 GPa and G_m = 44 GPa. There are three pore shapes in the rock, and their pore aspect ratios α are 1.0, 0.5 and 0.01 respectively, and their volume percentages are 98.99%, 1% and 0.01% respectively, in order to meet the requirement that the ratio of porosity to aspect ratio has to be much smaller than 1 for each inclusion in Kuster and Toksöz theory(Kuster and Toks&vouml;z, 1974). We use the fourth-order Runge-Kutta scheme to get numerical solutions of the differential Eqs.(3) and (4), the step size used in numerical integration is set to △y = 0.1. We assume that there are two different materials in the cavities, one is solid inclusions and the other is dry void. The elastic moduli of solid inclusions are K_i = 14 GPa and G_i = 10 GPa, while those of dry voids are set to 0. We show comparisons of the analytical and numerical results for two versions of differential equations in Fig. 1, and pink points represent the numerical results of Eq.(3), blue points represent the numerical results of Eq.(4). Furthermore, we give the theoretical estimation of Kuster-Toksöz theory which is displayed in black points and the Hashin-Shtrikman bounds which are displayed in red lines. Figs. 1a and 1b show the comparisons of the bulk and shear moduli for solid inclusions, and Figs. 1c and 1d show the comparisons of the bulk and shear moduli for dry voids. The results show that the theoretical estimation of Kuster-Toksöz theory and the upper Hashin-Shtrikman bound are identical, the numerical results of Eq.(3)are close to Kuster-Toksöz theory when the porosity is low, but as the porosity increases the difference between them increases and the numerical results of Eq.(3)violate the Hashin-Shtrikman upper bound and have no physical meaning. Thus Eq.(3)can only be used to simulate the elastic properties of multiple-porosity rock when the porosity is low, while Eq.(4)is always within the Hashin-Shtrikman bounds. As Eq.(3)is in conflict with the Hashin-Shtrikman bounds when the porosity is high, we only use Eq.(4)to calculate the effective moduli of multiple-porosity rock.

3 DEM ANALYTICAL FORMULAE FOR DRY MULTIPLE-POROSITY ROCK

For dry rock, the elastic moduli of the ith inclusion K_i = 0 and G_i = 0(1 ≤ i ≤ N), Eq.(4)becomes

and In general Eqs.(5a) and (5b)are also coupled, as two equations depend on both the bulk and shear moduli of the composite(Berryman, 1980; Li and Zhang, 2011). Therefore, it is difficult to obtain the analytical formulae of the bulk and shear moduli for dry multiple-porosity rock by integrating the DEM equations. Using the similar assumption of Li and Zhang(2010, 2011)that there is an almost linear relationship between the difference P_eff^* -Q_eff^* and the effective modulus ratio K^*(y)/G^*(y), we can obtain the analytical formula of modulus ratio for dry multiple-porosity rock(see Appendix A). where K_dry and G_dry are the bulk and shear moduli for dry multiple-porosity rock respectively, constants a and b are intercept and gradient respectively, which meet P_eff^* - Q_eff^* = a + bK^*(y)/G^*(y)(Li and Zhang, 2010), a =$\sum\limits_{i = 1}^N {{v_i}{a_i}} $, b=$\sum\limits_{i = 1}^N {{v_i}{b_i}} $constants ai and bi are intercept and gradient of P^*i - Q^*i for the ith inclusion respectively(as there are no fluids in inclusions, the inclusions are actually pores). Constant bi is simply set to the first derivative of P^*i - Q^*i where K^*/G^* = K_*/G_* and then constant ai can be calculated from P^*i - Q^*i = a + bK^*(y)/G^*(y)(Li and Zhang, 2011, 2012). The expressions of P^*i - Q^*i for ellipsoidal pores and its first derivative for porous rocks with ellipsoidal pores are given in Appendix A of Li and Zhang’s paper (Li and Zhang, 2011).

In order to decouple the equations for bulk and shear moduli, substituting Eq.(6)into(5a), and then, integrating it from y = 0 to φ gives directly the analytical results for K_dry, and then, dividing the K_dry by Eq.(6)gives the corresponding analytical results for G_dry, we have(see Appendix B)

where ${S_{0i}} = \frac{{2 - 3{g_i} - 3{\theta _i}}}{{4{\theta _i} - 6\theta _i^2 - 4{g_i}}}, {S_{1i}} = \frac{{{\theta _i} - {g_i}}}{{2 - 3{g_i} - 3{\theta _i}}}, {S_{2i}} = \frac{4}{3}, {S_{3i}} = \frac{{2({\theta _i} - {g_i})}}{{3(2{\theta _i} - 3\theta _i^2 - 2{g_i})}}$, αi is the pore aspect ratio of the ith ellipsoidal inclusion.

4 THEORETICAL AND EXPERIMENTAL EXAMPLES

Eqs.(7) and (8)can be used to study effective elastic properties of dry rocks with two or more kinds of pore shapes. Most current studies focus on the double-porosity rock(Xu and White, 1995). For convenience of comparison, we here also conduct the theoretical and laboratory data analysis by using the analytical formulae of bulk and shear moduli for dry double-porosity rock.

First we give a comparison between the analytical formulae and numerical results of the differential equations. We take quartz as the host medium, let K_m = 37 GPa and G_m = 44 GPa respectively which is the same as Fig. 1. There are two kinds of pore types, one represents pore and the other represents crack, and their pore aspect ratios α are 0.15 and 0.005 respectively.

To characterize the effect of different pores or cracks on the elastic properties of rock, we give three different volume percentages of pores and cracks, which are(90%, 10%), (50%, 50%) and (10%, 90%)respectively. We also use the fourth-order Runge-Kutta scheme to get numerical solutions of the differential Eq.(5). Fig. 2 shows the comparisons of the analytical and numerical results, in which points represent the numerical results, solid line represents the analytical results calculated by Eqs.(7) and (8)whose constant b is simply set to the first derivative of the difference between the polarization factors when K^*/G^* = K_m/G_m = 1.1892 and constant a can be then calculated from P^*i - Q^*i = a + bK^*(y)/G^*(y)(see Appendix A). The results show that the analytical results of bulk and shear moduli are in sufficiently good agreement with the numerical results in three different volume percentages of pores and cracks. As the volume percentage of cracks increases, the effective elastic moduli of the rock decrease, and the decrease rate of the effective moduli as the volume percentage of cracks increases from 10% to 50% is substantially larger than that as the volume percentage of cracks increases from 50% to 90%.

We use 69 dry s and stone data of Han(1986)to analyze the predictions of the analytical formulae. The data were obtained by measuring wave velocities in dry-rock samples at confining pressure 40 MPa and ultrasonic frequency 1 MHz and 0.6 MHz respectively. The porosity of these samples ranges from 5% to 30%, and the clay content range in these samples is between 0 and 50%. Fig. 3a shows the relation between porosity and P-wave velocities of Han’s dry s and stone data. The measured results indicate that P-wave velocities of dry s and stone decrease as porosity or clay content increases. The bulk and shear moduli of quartz used in the analysis were taken to be 40 GPa and 43.7 GPa, respectively. The bulk and shear moduli of clay were taken to be 20.9 GPa and 6.85 GPa, respectively. We use the double-porosity DEM analytical formulae to calculate the elastic properties of s and stone data, in which the pore aspect ratios of stiff pores and soft pores are 0.18 and 0.04, respectively. It needs the volume percentages of two different kinds of pores to predict the elastic moduli with porosity using Eqs.(7) and (8). Fig. 3b shows the variations of the bulk moduli of 10 dry clay-free s and stones with porosity in different pore volume percentages. We can see that when the volume percentages of soft pores V_cp = 0 there are only two points on the predicted line. And the point(porosity is 6%)is far away from the predicted line, when the volume percentage of soft pores increases to V_cp = 25%, the predicted line is in agreement with that point. It shows that there is only two points having no soft pores, while the other eight points having both stiff and soft pores to varying growth levels.

We estimate the volume percentage of soft pores by minimizing the error between theoretical predictions by Eqs.(7) and (8) and experimental measurements. Fig. 3c shows the relationship between the inverted volume percentages of soft pores and clay content. In general, the porosity of the rock is larger, the volume percentage of soft pores is smaller, but the volume percentage of soft pores is not significantly correlated with clay content. We obtain the shear wave velocity from the inverted volume percentages of soft pores by Eqs.(7) and (8). Fig. 3c shows the comparison between the predicted and measured shear wave velocities, and the error bars with 95% confidence interval are shown. Observe that for the majority of samples the analytical formulae are in good agreement with the measured data, and it shows that the pore structure of s and stone can be characterized by the double-porosity DEM analytical model.

5 CONCLUSIONS

By applying the Kuster-Toksöz theory this paper established the differential equations of both Zimmermann’s and Norris’s versions for multiple-porosity rock, and derived analytical solutions of the bulk and shear moduli for dry rock from the differential equations of Norris’s version. The analytical formulae are functions of the bulk and shear moduli of the mineral grain, porosity, pore volume percentage and pore geometry. The pore geometry is characterized by a linear relation whose intercept a and gradient b might be approximately calculated from the difference between the effective polarization factors and its first derivative where the effective modulus ratio is equal to the mineral grain modulus ratio.

The differential effective medium model of multiple-porosity rock avoids the problem that the estimated elastic properties depend on the adding order of different pore-type inclusions when the classical DEM equations are used for multiple-porosity rock. When the porosity type is single, the multiple-porosity DEM equations become the classical DEM equations.

Numerical results show that the elastic moduli predicted by the DEM equations of Norris’s version never violate Hashin-Shtrikman bounds while those predicted by the DEM equations of Zimmermann’s version violate bounds in some cases, so the latter can only be used to predict the elastic moduli for the porous rock with very low porosity.

Theoretical models and laboratory data tests show that the predictions of the analytical formulae are in good agreement with the numerical solutions of the full DEM equation. The DEM analytical formulae for multiple-porosity rock can accurately predict the variations of elastic moduli with porosity. The dual-porosity DEM analytical model can be used to invert for the pore volume percentages of both cracks and pores by minimizing the error between theoretical predictions and experimental measurements. The inversion results for the s and stone experimental data show that there is no direct correlation between the crack volume percentage and clay volume.

6 Appendix A Derivation of the Dry Rock Modulus Ratio of Multiple-Porosity Rock

Here we first define the effective polarization factor, which is the sum of the polarization factors weighted by their volumes of the pores contained in the medium, i.e.

and Therefore, Eq.(5)becomes and Thus, Eq.(A3)–(A4)are completely consistent with the classical DEM differential equations of single-porosity rock(Norris, 1985; Mavko, et al., 1998; Li and Zhang, 2010), their analytical approximations can be obtained by using the similar methods as Li and Zhang(2011, 2012)did.

Subtracting Eq.(A4)from Eq.(A3)yields

Li and Zhang(2010)assumed that there is an almost linear relationship between the difference P^*i-Q^*i and the effective modulus ratio K^*(y)/G^*(y) for dry single-porosity rock, i.e. for the ith pore in the multiple-porosity rock, we have where constants ai and bi relate to pore shapes and a_i < 0, b_i ≥ 0. This assumption is also suitable to multiple-pore system, because where $a = \sum\limits_{i = 1}^N {{v_i}{a_i}}, b = \sum\limits_{i = 1}^N {{v_i}{b_i}} $

Substituting Eq.(A7)into Eq.(A5)yields

Furthermore, Eq.(A8)is decomposed into

Integrating Eq.(A9)fromy = 0 to φ, and applying the initial condition we have where 0≤ φ < 1. Eq.(A11)can be easily rearranged to give

7 Appendix B Derivation of the Dry Bulk and Shear Moduli of Multiple-Porosity Rock

According to the appendix A shown in Li and Zhang’spaper(2011), for dry spheroidal void the polarization factor of the ith pore of multiple-porosity rock is

where α_i is the pore aspect ratio of the ith spheroidal void.

From Eqs.(5a) and (B1)we obtain

The right-h and side of Eq.(B2)is furthermore decomposed into where Exp and ing the last term around y = 0 on the right h and side of Eq.(B3)using first-order Taylor series, we have For dry rock, K^*(y)/G^*(y) in Eq.(B8)is equal to K_dry(y)/G_dry(y). Substituting Eq.(6)into Eq.(B8), then integrating it from y = 0 toφ, and since 0 ≤ y ≤ φ and 0 ≤ φ ≤ 1, omitting the quadratic term with respect to (i.e., omitting the second term within the square bracket in Eq.(B8)), we have Eq.(B9)can be rearranged to give Dividing Eq.(B10)by Eq.(6)gives

8 ACKNOWLEDGMENTS

Thank the editors and anonymous reviewers for many insightful comments and suggestions that improve this paper. This work was supported by the National Natural Science Foundation of China(41404088) and the Important National Science and Technology Specific Projects of China(2011ZX05001).