1 INTRODUCTION

The Gassmann equations are important tools for shear wave velocity prediction,fluid substitution and sonic velocities evaluation^[1]. However,how to determine the elastic moduli of dry-rock frame in Gassmann equations is a problem not completely solved yet. Therefore,geophysicists present many methods for calculating them. There are two kinds of popular models,one is the empirical model,and the other is the effective medium model.

The empirical model assumes simply that the elastic moduli are related to porosity,not taking into account the effect of pore geometry. The critical porosity model and Pride model^{[2, 3]} are two common empirical models,which have their own assumptions on the modulus ratio of dry rock. The critical porosity model assumes a constant modulus ratio of dry rock,independent of porosity,while the Pride model assumes that the modulus ratio of dry rock increases with porosity. However,it has been proven that these assumptions are not fully reasonable; and laboratory data have shown that the modulus ratio of dry rock not only increases with the porosity,but also decreases with the porosity^{[4, 5, 6, 7]}. Therefore,there are limitations while applying these models to estimate elastic properties of the rock.

The effective medium models,such as the Kuster-Toksöz theory^[8] and differential effective medium theory^[9],establish the relationship among the elastic moduli of dry rock and porosity,pore shapes. Rock physics theory indicates that the velocity is strongly affected by pore shapes,which are the main cause for divergence between velocity and porosity^{[10, 11]}. Li and Zhang systematically derived the elastic modulus formulae of dry-rock frame from the differential effective medium theory,which are referred to here as the DEM analytical model. These formulae can characterize the dry-rock moduli dependency on porosity and pore shapes^{[7, 12, 13]}. The effective medium models have been widely used in the shear velocity prediction^{[10, 14]}. These models commonly describe the pore shapes of the rock as the ideal ellipsoid^{[8, 10, 15]}. However,the pore aspect ratio which is used to characterize pore shapes is generally difficult to be obtained from laboratory or field measurements. The popular Xu-White model assumes that the porosity of the rock is divided into two parts,one is the stiff pore and the other is the compliant pore,and their volumes are respectively determined from sand and clay content and their aspect ratios are constant^[10]. These assumptions have negative effects on the accuracy of shear velocity prediction. Therefore,some researchers proposed a method to improve accuracy by applying a depth-dependent variation of pore aspect ratios which is determined by the empirical relationship in the Xu-White model^{[16, 17]},but the empirical relationship depends on the samples used,and cannot be randomly extended.

This paper presents a method for pore structure evaluation and shear wave velocity prediction in terms of the effective pore aspect ratio inverted from the P-wave (and S-wave) velocity based on the DEM analytical model^[13]. First we introduce the Gassmann equations,then introduce the analytical formulae,which are the function of the pore aspect ratio,of the bulk and shear moduli of dry rock derived from the differential effective medium theory. Next we apply the nonlinear global optimization algorithm to invert the effective pore aspect ratio by using the Gassmann equations. Finally we study the pore aspect ratio variations with the effective pressure on the granite and sandstone data from laboratory measurements,and perform the effective pore aspect ratio inversion and S-wave velocity prediction on the laboratory sandstone data^[18] and the logging data.

2 METHODOLOGY 2.1 Gassmann Equations

Upon the low frequency assumption,Gassmann^{[1, 19]} gives the formulae of the bulk and shear moduli of the fluid-saturated rock

and where K_dry,K_m and K_f are respectively the bulk moduli of dry-rock frame,matrix and pore fluid,G_dry is the shear modulus of dry-rock frame,and $\phi $ is the total porosity.

After the bulk and shear moduli of the fluid-saturated rock are determined by the Gassmann equations,the P- and S-wave velocities of the fluid-saturated rock can be expressed as

where V_p and V_s are respectively P- and S-wave velocities,and ρ_sat is the bulk density of the fluid-saturated rock.

2.2 Analytical Formulae of Elastic Moduli of Dry Rock

When Gassmann equations (Eqs.(1) and (2)) are used to study the elastic properties of the fluid-saturated rock,the elastic modulus of dry-rock frame must be given in advance,but there are no unified formulae about the elastic moduli of dry-rock frame. Usually these parameters can be determined by the following ways: (1) Rock samples are directly measured in the laboratory^{[18, 20]}. This method has high accuracy,but is generally difficult to obtain the values expected. (2) Make prediction by the empirical relationship^{[2, 3]} or effective medium theory^{[8, 9, 10, 12, 13]},their prediction accuracy depends on the model assumptions^{[12, 21]}.

The dry-rock model,based on empirical relationships,assumes simply that the elastic moduli are related to porosity,regardless of the pore structure. The critical porosity model which is the most common empirical model assumes that both bulk and shear moduli of dry rock vary linearly with porosity for porosities less than a given critical value. Nur (1998) suggested that the relationship between the elastic moduli of matrix and dry-rock frame can be written as^[2]

where $\phi $_c is the critical porosity,and G_m is the shear modulus of matrix. The critical porosity is taken as $\phi $_c ≈ 0.40 for sandstone.

Pride et al. (2004) suggested that the bulk and shear moduli of dry consolidated sandstone can be expressed as^[3]

where c is a consolidation parameter that represents the degree of consolidation between grains,which is taken as 2 < c < 20 for sandstone.

Both the critical porosity model and the Pride model have limitations: there is an implicit and inappropriate assumption on modulus ratio of dry rock in these models. According to Eqs.(5)-(8),the critical porosity model implicitly assumes a constant modulus ratio of dry rock,irrespective of porosity,i.e.,the ratio of P-wave to S-wave velocity V_p/V_s of dry rock is equal to the velocity ratio of matrix,but this assumption sometimes does not match with laboratory measurements^{[4, 5, 21]}. Whereas the Pride model implicitly assumes that the modulus ratio of dry rock increases monotonically with the porosity because the consolidation coefficient is always positive. But the modulus ratio for sandstone may increase with porosity (such as quartz sandstone) or decrease with porosity (such as feldspar sandstone)^{[6, 7]}.

The inclusion theory,such as differential effective medium theory,suggests that the modulus ratio of dry rock is related not only to the porosity but also to the pore geometry^{[9, 15]}. The numerical solutions of differential effective medium theory can accurately simulate the modulus ratio of dry rock with porosity for ideal ellipsoidal pore and penny-shaped crack. Li and Zhang obtain DEM analytical models of the modulus ratio and elastic moduli of dry rock by decoupling the two coupled ordinary differential equations for bulk and shear moduli^{[7, 12, 13]}. For ellipsoidal pores,their analytical formulae are written as

where ${{S}_{0}}=\frac{2-3g-3\theta }{2(2\theta -3{{\theta }^{2}}-2g)},{{S}_{1}}=\frac{\theta -g}{2-3g-3\theta },{{S}_{2}}=\frac{4}{3},{{S}_{3}}=\frac{2(\theta -g)}{3(2\theta -3{{\theta }^{2}}-2g)},\theta =\left\{ \begin{matrix} \frac{\alpha }{{{(1-{{\alpha }^{2}})}^{3/2}}}[arccos(\alpha )-\alpha \sqrt{1-{{\alpha }^{2}}}] & (\alpha <1),\\ \frac{\alpha }{{{({{\alpha }^{2}}-1)}^{3/2}}}[\alpha \sqrt{{{\alpha }^{2}}-1}-cos{{h}^{-1}}(\alpha )] & (\alpha >1),\\ \end{matrix} \right.g=\frac{{{\alpha }^{2}}}{1-{{\alpha }^{2}}}(3\theta -2)$,and they are functions of the pore aspect ratio α. Constants a and b are intercept and gradient respectively,which meet P^*i - Q^*i = a + bK^*(y)/G^*(y). Constant b is simply set to be the first derivative of P^*i - Q^*i and constant a can be then calculated from P^*i - Q^*i = a + bK^*(y)/G^*(y) where K^*/G^* = K_m/G_m. The expressions of P^*i - Q^*i and their first derivatives for porous rocks with ellipsoidal pores are given in Appendix A of the Ref.^[13]. The P^*i and Q^*i are called the polarization factors of the effective bulk and shear moduli K^* and G^*,and they usually have been computed from Eshelby’s well-known results and Wu’s results^{[1, 15]}.

Eqs.(9) and (10) implicitly define the modulus ratio of the dry rock as

Figure 1a shows the modulus ratio with porosity for dry quartz sandstone with α=1,0.2,0.1,0.05,0.01,0.001 by using Eq.(11). We take quartz as the host medium,having K_m=37 GPa and G_m=44 GPa. As seen in Fig.1a,the modulus ratio of dry rock increases with porosity when α > 0.1; and the modulus ratio of dry rock is almost constant when α = 0.1,independent of the porosity; and the modulus ratio of the dry rock decreases with porosity when α < 0.1. Therefore,assumptions that the modulus ratio of dry rock is constant in the critical porosity model and it increases with porosity in the Pride model are special cases of Eqs.(9)-(11). Fig.1b shows the modulus ratio varying with porosity for dry quartz and feldspar sandstone with α = 1 by using Eq.(11). The elastic moduli of feldspsr sandstone is taken as K_m = 37.5 GPa and G_m = 15 GPa. As seen in Fig.1b,for sandstones with the same pore shapes,the modulus ratio of dry feldspar rock decreases with porosity,significantly different from the quartz sandstone. Knackstedt et al. (2003) got similar results by using the finite element method to estimate the elastic properties of the model system^[6]. Fig.1 shows that Eq.(11) can be used to characterize the dry rock modulus ratio varying with porosity for various rocks; and it overcomes the limitations of the empirical models,and has wider applications than the conventional models,such as the critical porosity model.

2.3 Inversion of Pore Aspect Ratios

Eqs.(9) and (10) show that the elastic moduli of dry rock are not only the function of the elastic moduli of the mineral and porosity,but also related to the pore aspect ratio. Substituting Eqs.(9)-(10) into Eqs.(1)-(4),we establish a nonlinear relationship between the wave velocities and the pore aspect ratio. We not only obtain the forward solutions for P- and S-wave velocities from the elastic moduli of the mineral and the porosity and the saturation as well as the pore aspect ratio,but also invert the P- and S-wave velocities to obtain the pore aspect ratio. In order to apply Eqs.(1)-(4) and (9)-(10) to estimate the pore aspect ratio,here we simply treat the rock pores as an ideal ellipsoid with a single aspect ratio,the objective function is

where K_sat^cal and G_sat^cal are respectively the bulk and shear moduli of the predicted fluid-saturated rock,W_p,W_s are the weighting factors,W_p + W_s = 1. When only the P-wave data is available,W_p = 1,W_s = 0. When both the P- and S-wave data are available,W_p = 0.5,W_s = 0.5. We solve Eq.(12) by a nonlinear global optimization: simulated annealing algorithm. When S-wave log is absent,we can predict the S-wave velocity from the effective pore aspect ratio inverted from only the P-wave log. When both P- and S-wave logs are available,we can evaluate these data by the effective pore aspect ratio inverted from both P- and S-wave logs. Because there exist noises in the measured P- and S-wave velocity logging data,the Poisson’s ratio calculated from these data might be negative. Especially,the abnormal values will be present for acoustic travel-times due to the cycle skip particularly in the enlarged diameter section of the wellbore,and S-wave will be affected more seriously. Therefore these data are often not reliable,and must be corrected.

3 REAL APPLICATIONS 3.1 Laboratory Data

We first study the effective pore aspect ratios inverted from both the P- and S-waves varying with the effective pressure. Fig.2a is the results of the effective pore aspect ratios inversion for Coyner’s Westerly granite data^[22] which were obtained by measuring wave velocities and porosity in dry-rock samples at the effective pressure of 5 to 100 MPa. The rock porosity is 0.8% at 0 MPa,decreases slowly with pressure and reduces to about 0.71% at 100 MPa. The results show that the effective pore aspect ratio increases with the effective pressure. The effective pore aspect ratio of granite sample increases from 0.009 to 0.023 when the effective pressure increases from 5 to 100 MPa. Fig.2b is the results of the effective pore aspect ratios inversion for Coyner’s Navajo quartz sandstones data^[22] which were obtained by measuring wave velocities in dry-rock samples at the effective pressure of 0 to 100 MPa. The rock porosity at the effective pressure of 0 to 100 MPa ranges from 11.5% to 11.8%. The results show that the effective pore aspect ratio of the sandstone sample increases with the effective pressure. The effective pore aspect ratio of the sandstone sample increases from 0.08 to 0.15 when the effective pressure increases from 0 to 100 MPa. Fig.2 shows that our method can accurately obtain the pore aspect ratio not only for fractured reservoirs,such as granite,but also for porous reservoirs,such as sandstone. The estimation of the pore aspect ratios of a rock from its measured velocities at different effective pressures was presented first by Cheng and Toksöz (1979)^[23]. Since their method is based on the Kuster-Toksöz model and iterative mode,their computational process is very complicated,and the inversion of the pore aspect ratios spectrum has a non-uniqueness,dependent on the pre-given aspect ratios,especially the aspect ratios of large pores^[24]. On the other hand,their method is limited to the laboratory measurements,difficult to be effectively applied to logging data because the measurements at different pressures can not be simultaneously obtained for the reservoirs interested. Furthermore,Kuster-Toksöz theory itself is suitable for low-porosity rocks,and its application is limited by its theoretical assumptions^{[8, 25]}.

Secondly,we use 75 water-saturated sandstone data of Han^[18] to analyze S-wave velocity predictions from the inverted effective pore aspect ratios. The data were obtained by measuring wave velocities at confining pressure 40 MPa. The porosities of these samples range from 5% to 30%,and the clay content range in these samples is 0 to 50%. The bulk and shear moduli of quartz used in the analysis were taken as 41.42 GPa and 37.87 GPa,respectively. The bulk and shear moduli of clay were taken as 23.98 GPa and 6.96 GPa,respectively. The wet bulk modulus is calculated from its transit time and density that Xu and White provided^[10]. The bulk and shear moduli of shaly sandstones are calculated by using the Voigt-Reuss-Hill averaging method^[1]. Fig.3a shows the relation between the inverted effective pore aspect ratios and porosity,the solid points stand for the inversion only from the P-wave velocity,and the hollow points stand for the inversion from both P- and S-wave velocities. Only three of 75 points are greater than 0.4 and their pore structures are close to the sphere,and they are not shown in Fig.3a. The results show that every hollow point is close to its corresponding solid point and prove that the inversion of the effective pore aspect ratios only from P-wave data is credible and the method of S-wave prediction based on the effective pore aspect ratios is correct. For sandstone,on the whole,the pore aspect ratios are relatively stable (around 0.15) if porosities are greater than 15%. The effective pore aspect ratios have wider distributions if porosities are less than 15%. The higher the clay content,the lower the pore aspect ratios; but if porosities are less than 0.15,there is no clear correlation between the pore aspect ratio and clay content. The P- and S-wave velocities can then be estimated from the inverted effective pore aspect ratio by using Eqs.(9)-(10) and Eqs.(1)-(4). Fig.3b is a comparison between the predicted and measured shear wave velocities. Observe that for the majority of samples the predicted values are in good agreement with the measured,there are only two samples for which the predictions of the shear velocity are deviated from the measured data. This fully shows that the effective elastic moduli of these sandstones can be reasonably well characterized by spheroidal pores. The assumption on porous rock with a single aspect ratio is different from the Xu and White model in which two aspect ratios are used,one for the sand fraction and the other for the clay fraction,but these two aspect ratios have to be fixed. Fig.3c is a comparison between the measured and predicted shear wave velocities by the Xu-White model for the same sandstone samples in Fig.3b^[10]. We see that the predictions of this paper are more accurate than the predictions of the Xu-White model. On the whole the predicted shear velocities by the Xu-White model are overestimated. This fully shows that our method of the shear-velocity prediction by using the pore aspect ratio is more reasonable and credible.

3.2 Logging Data

We evaluate the effective pore aspect ratio inversion and shear wave velocity prediction by logging data from the Bohai Bay. The sandstone reservoirs in this well are located between 2712 and 2880 m of depth where the gas-bearing reservoirs are from 2712 to 2778 m,and the water-bearing reservoirs are from 2778 to 2880 m.

Figure 4a shows the logging and inverted curves of the target. Column one to column three from left to right respectively show the P-wave velocity,density,clay content,porosity and water saturation curves. The inverted effective pore aspect ratio curves are displayed in the fourth column where the inversion only from P-wave is displayed as a red curve and the inversion from both P- and S-wave is displayed as a black curve. These two inversions are almost identical. Fig.4b is the distribution of the inverted effective pore aspect ratios,and they are mainly distributed in the range 0.18 to 0.32,most of which are around 0.2. It indicates that the reservoirs of this well belong to porous reservoirs. In addition,there are several portions of the reservoirs whose effective pore aspect ratios approach to 1. Their pores are close to the sphere and P-wave velocities are relatively high. The fifth column in Fig.4a is a comparison of the measured and predicted shear velocities where the measured shear wave velocity is displayed as a black curve and the predicted one only from P-wave velocity is displayed as a red curve. Note that the predicted shear wave is in good agreement with the measured one,proving that the S-wave velocity reconstructed from the measured P-wave data is reliable and accurate.

4 CONCLUSIONS

This paper proposed a method for the pore structure evaluation and shear wave velocity prediction from the effective pore aspect ratios inverted from the P-wave (and S-wave) velocity. First,we built a link among Pand S-wave velocities,density,porosity,saturation and mineral compositions through Gassmann equations and the DEM analytical model and applied nonlinear global optimization algorithm to find the best estimate for the effective pore aspect ratios; Then,we estimated the S-wave velocity from the inverted effective pore aspect ratios and evaluated the reservoir pore structure by the pore aspect ratios.

Compared with the existing empirical models,the DEM analytical model derived from the differential effective medium theory can characterize pore geometry by pore aspect ratios,overcome the limitations of the empirical models,such as the critical porosity model that only relates the elastic moduli of dry rock to porosity regardless of the pore shape and assumes a constant modulus ratio irrespective of porosity,and accurately describes elastic moduli and its modulus ratio varying with porosity for any kind of lithology.

Applications to the laboratory and borehole data show that the effective pore aspect ratios inverted from velocities can be used to evaluate the pore structure of reservoirs accurately. For fractured reservoirs such as granite,the pore aspect ratio is typically less than 0.025. For porous reservoirs such as sandstone,the pore aspect ratio is usually greater than 0.08. The effective pore aspect ratio only from P-wave is almost exactly the same as that one from both P- and S-wave. And the predicted S-wave from P-wave is in good agreement with the measured shear wave. This fully shows that our method of the shear-velocity prediction by using the pore aspect ratio inverted from P-wave is effective.

Our method on the pore-type evaluation for the reservoir by assuming the rock pores to be an ideal ellipsoid with a single aspect ratio and inverting the pore aspect ratio from P-wave data is more suitable for the reservoirs whose pore type is single pore or crack. For complex reservoirs with multiple pore types,this method can also be used to obtain an effective pore aspect ratio for the S-wave velocity prediction and evaluation. How to invert their individual pore aspect ratios from P-wave data for the reservoir evaluation is an issue to be addressed in the future work.

5 ACKNOWLEDGMENTS

We thank the editors and anonymous reviewers for many insightful comments and suggestions that improve this paper. This work was supported by the National Major Science and Technology Specific Projects of China (2011ZX05001).