A differential scheme for elastic properties of rocks with dry or saturated cracks

Summary

Differential effective medium (DEM) theory is applied to the problem of estimating the physical properties of elastic media with penny-shaped cracks, filled either with gas or liquid. These cracks are assumed to be randomly oriented. It is known that such a model captures many of the essential physical features of fluid-saturated or partially saturated rocks. By making an assumption that the changes in certain factors depending only on Poisson's ratio do not strongly affect the results, it is possible to decouple the equations for bulk (K) and shear (G) modulus, and then integrate them analytically. The validity of this assumption is then tested by integrating the full DEM equations numerically. The analytical and numerical curves for both K and G are in very good agreement over the whole porosity range. Justification of the Poisson ratio approximation is also provided directly by the theory, which shows that as porosity tends to unity, Poisson's ratio tends towards small positive values for dry, cracked porous media and tends to one-half for liquid-saturated samples. A rigorous stable fixed-point is obtained for Poisson's ratio, ν_c, of dry porous media, where the location of this fixed-point depends only on the shape of the voids being added. Fixed-points occur at for spheres and ν_c≃πα/18 for cracks, where α is the aspect ratio of penny-shaped cracks. These theoretical results for the elastic constants are then compared and contrasted with results predicted by Gassmann's equations and with results of Mavko and Jizba, for both granite-like and sandstone-like examples. Gassmann's equations do not predict the observed liquid dependence of the shear modulus G at all. Mavko and Jizba predict the observed dependence of the shear modulus on the liquid bulk modulus for a small crack porosity and a very small aspect ratio, but fail to predict the observed behaviour at higher porosities. In contrast, the analytical approximations derived here give very satisfactory agreement in all cases for both K and G. For practical applications of this work, it appears that the ratio of compliance differences is approximately independent of the crack porosity for a given rock, but the constant is usually greater than for granites, while general statements concerning sandstones are more difficult to make.

cracks, differential effective medium theory, granite, liquid saturation, rock mechanics, sandstone

1 Introduction

The elastic moduli of rocks are dependent on: the mineral properties and distribution; porosity type, magnitude and distribution; and the state of saturation. Two major theoretical approaches have been developed to address the problem of estimating elastic moduli from knowledge of the rock constituents and the microstructure: (1) effective medium theory, which assumes separate pores and cracks that may or may not be connected and (2) poroelastic theory, which assumes that significant portions of the pores and cracks are connected. Effective medium theories, which include the classical bounds of Voigt (1928) and Reuss (1929) and Hashin & Shtrikman (1961, 1962) as well as estimates obtained from self-consistent theories (e.g. Budiansky 1965; Hill 1965; Berryman et al. 1980a,b), require parameters characterizing the pore shape and distribution. Alternatively, poroelastic constitutive equations (Biot 1941; Gassmann 1951) are phenomenological and do not require characterization of matrix and pore space geometry. However, they contain the fundamental result, sometimes in disagreement with experiment (Coyner 1984; Coyner & Cheng 1985; Mavko & Jizba 1991; Berryman & Wang 2001; Berryman 2002), that the shear modulus is always independent of the saturation state (Berryman 1999). Although the lack of a shear dependence on the saturating fluid bulk modulus can be correct for special but real microgeometries (such as purely spherical pores) and very low modulation frequencies, this predicted lack of dependence is often contradicted by high-frequency experiments (above ∼1 kHz), and especially so in rocks with crack porosity. As a result, Biot's theory has been modified in various ways. For example, Mavko & Jizba (1991) partition porosity into ‘soft’ and ‘stiff’ porosity fractions to account for the change of both bulk modulus and shear modulus with fluid saturation.

A large literature has developed to address the many issues related to cracked elastic and poroelastic media, in the dry, saturated, and partially saturated cases. Recent comprehensive reviews of the literature on the analysis of cracked elastic materials include Kachanov et al. (1992), Nemat-Nasser et al. (1993) and Kushch & Sangani (2000), as well as the textbook by Nemat-Nasser & Hori (1993). Some of the notable work on dry cracked solids using techniques similar to those that will be employed here includes Zimmerman (1985), Laws & Dvorak (1987), Hashin (1988) and Sayers & Kachanov (1991). Pertinent prior work on both dry and saturated cracked rocks includes Walsh (1969), Kuster & Toksöz (1974), Budiansky & O'Connell (1976), O'Connell & Budiansky (1974, 1977), Walsh & Grosenbaugh (1979), Hudson (1981, 1986, 1990), Henyey & Pomphrey (1982), Mavko & Jizba (1991), Berryman (1992), Hudson et al. (1996) and Pointer et al. (2000). Among the most popular effective medium approaches are the two implicit schemes: (1) the self-consistent (Budiansky 1965; Hill 1965; Korringa et al. 1979; Berryman 1980a, b) and (2) the differential effective medium (Bruggeman 1935; Cleary et al. 1980; Norris 1985; Avellaneda 1987). Other than their relative ease of computation, the main reason for their popularity is that they are the only ones known to be realizable (i.e. they correspond to an actual microgeometry), and therefore have the property that they can never violate rigorous bounds—unlike some of the popular explicit schemes such as that of Kuster and Toksöz, which is known to violate bounds in some cases. Berryman & Berge (1996) summarize the status of the most commonly used explicit schemes, and conclude that while they can be used, they also must be carefully restricted to sufficiently low volume fractions of inclusions. Of the two implicit and realizable methods, the differential effective medium theory has some unique advantages for analysis that have not been stressed previously, and that we will develop more fully here.

The purpose of this paper is to obtain approximate analytical results for the elastic moduli of dry and fully saturated cracked rock based on differential effective medium (DEM) theory (Bruggeman 1935; Cleary 1980; Walsh 1980; Norris 1985; Avellaneda 1987). Penny-shaped cracks have been used extensively to model cracked materials (Walsh 1965, 1969; Willis 1980; Kachanov 1992; Smyshlyaev et al. 1993), but the penny-shaped crack model is itself an approximation to Eshelby's results (Eshelby 1957; Wu 1966; Walsh 1969) for oblate spheroids having a small aspect ratio. In order to obtain some analytical formulae that are then relatively easy to analyse, a further simplifying assumption is made here that certain variations in Poisson's ratio with the change of crack porosity can be neglected to first order. The consequences of this new approximation are checked by comparison with numerical computations for the fully coupled equations of DEM. The agreement between the analytical approximation and the full DEM for cracked rock is found to be quite good over the whole range of computed porosities. Justification for the approximation is provided in part by an analysis of the actual variation of Poisson's ratio, and some further technical justifications are also provided in two of the three appendices.

For simplicity, the main text of the paper treats materials having only crack porosity, and we consider these models to be granite-like (by which term we mean to imply only the presence of a fairly uniform host material with cracks and do not imply anything regarding other aspects of diagenesis) idealizations of rock. A third appendix shows how the results of the main text change if the model is treated instead as a sandstone-like material (again implying only that the composite certainly has some large aspect ratio pores present as well as the small ones that most concern us here, with no other intended implications of the terminology) having finite stiff, nearly spherical porosity in addition to the soft, crack porosity.

2 Differential Effective Medium Theory

Differential effective medium theory (Bruggeman 1935; Cleary et al. 1980; Walsh 1980; Norris 1985; Avellaneda 1987) takes the point of view that a composite material may be constructed by making infinitesimal changes in an already existing composite. As mentioned in the introduction, there are only two effective medium schemes known at present that are realizable, i.e. that have a definite microgeometry associated with the modelling scheme. The differential scheme is one of these (Cleary et al. 1980; Norris 1985; Avellaneda 1987)—and one version of the self-consistent approach (Korringa et al. 1979; Berryman 1980a, b; Milton 1985) is the other. This fact, together with the associated analytical capabilities (including ease of computation and flexibility of application), provide a strong motivation to study the predictions of both of these schemes and the differential scheme in particular. We can have confidence that the results will always satisfy physical and mathematical constraints, such as the Hashin–Shtrikman bounds (Hashin & Shtrikman 1961, 1962).

When the inclusions are sufficiently sparse that they do not form a single connected network throughout the composite, it is appropriate to use the DEM to model their elastic behaviour (Berge et al. 1993). We assume that the host material has moduli K_m and G_m. The inclusion material has moduli K_i and G_i. Then, the effective bulk and shear moduli of the composite are parametrized by K*(y) and G*(y) when the volume fraction of the inclusion phase is y. The equations governing the changes in these constants are then well known to be  

(1)

and  

(2)

where the scalar factors P*ⁱ and Q*ⁱ, will be explained in the following paragraph, y is porosity, which equals the inclusion volume fraction here, and the subscript i again denotes the inclusion phase. We assume that the reader is somewhat familiar with this approach, and will therefore not dwell on its derivation, which is easily found in many places including, for example, Berryman (1992). 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=φ, or possibly over the whole range to y= 1. When integrating in this way, we imagine the result is simulating cracks being introduced slowly into a granite-like solid. The same procedure can be used for a sandstone-like material, assuming this medium has a starting porosity y=φ₀ with K*(φ₀) =K_s and G*(φ₀) =G_s. Integration then proceeds from y=φ₀ to 1. This introduction of crack (or soft) porosity into a material containing spherical (or stiff) porosity is conceptually equivalent to the porosity distribution model of Mavko & Jizba (1991). For simplicity, we will treat the granite-like case here, but the changes needed for other applications are not difficult to implement, and are treated in Appendix A.

The factors P*ⁱ and Q*ⁱ appearing in (1) and (2) are the so-called polarization factors for bulk and shear modulus (Eshelby 1957; Wu 1966). These depend, in general, on the bulk and shear moduli of both the inclusion, the host medium (assumed to be the existing composite medium * in DEM), and on the shapes of the inclusions. The polarization factors usually have been computed from Eshelby's well-known results (Eshelby 1957) for ellipsoids, and Wu's work (Wu 1966) on identifying the isotropically averaged tensor based on Eshelby's formulae. These results can be found in many places including Berryman (1980b, 1995) and Mavko et al. (1998).

The special case of most interest to us here is that of penny-shaped cracks, where  

(3)

and  

(4)

where α(0 < α < 1) is the crack (oblate spheroidal) aspect ratio,  

(5)

and where the superscript * identifies constants of the matrix material when the inclusion volume fraction is y. This formula is a special limit of Eshelby's results not included in Wu's paper, but apparently first obtained by Walsh (1969). Walsh's derivation assumes G_i/G_m≪ 1 and allows K_i/K_m≪ 1, with these approximations being made before any assumptions concerning the smallness of the aspect ratio α. By taking these approximations in the opposite order, i.e. letting the aspect ratio be small first and then making assumptions concerning smallness of the inclusion constants, we would instead obtain the commonly used approximation for discs. However, this latter approximation is actually quite inappropriate for the bulk modulus when the inclusion phase is a gas such as air (since then the ratio K_i/K_m≪ 1) or for the shear modulus when the inclusion phase is any fluid (since then G_i≡ 0), as the formulae become singular in these limits. This is why the penny-shaped crack model is commonly used instead for cracked rocks.

In general the DEM eqs (1) and (2) are coupled, as both equations depend on both the bulk and shear modulus of the composite. This coupling is not a serious problem for numerical integration. Later in the paper, we will show results obtained from integrating the DEM equations numerically. However, the coupling is a problem in some cases if we want analytical results to aid our intuition. We will now present several analytical results for first the bulk and then the shear moduli, and then we will compare these results with the fully integrated DEM results later on.

2.1 Some analytical results for K*

We now assume that the inclusion phase is a fluid so that K_i=K_f and G_i=G_f= 0. The fluid can be either a liquid or a gas. We consider three cases: (1) liquid inclusions, K_f≫παγ_m; (2) gas inclusions: K_f≪παγ_m; and (3) general inclusions, K_f≃παγ_m. Case 1 corresponds to liquid inclusions, case 2 to gas inclusions and case 3 to a circumstance in which we do not want to limit ourselves to the assumptions of either of the previous two cases, or in which the crack aspect ratio is tuned to the fluid modulus.

2.1.1 Liquid inclusions: K_f≫παγ_m

In this limit, it is somewhat more convenient to rewrite the DEM equations in terms of compliances, rather than stiffnesses, so we have  

(6)

The only terms that couple the equations for bulk to shear modulus have been readily neglected in this case, since P*ⁱ≃K*/K_f. Thus, we expect little if any deviation between the analytical results and the full DEM for the liquid-saturated case. We are treating here the granite-like case such that the limit of zero inclusion volume fraction corresponds to K*(0) =K_m, i.e. the bulk modulus of the pure solid. Then, integrating eq. (6) from y= 0 to φ (φ is the resulting porosity in the composite medium) gives directly  

(7)

which may be rearranged as  

(8)

Eq. (8) can also be obtained as the small-φ limit of Gassmann's equation when the saturating fluid is a liquid. Gassmann's result for the bulk modulus (Gassmann 1951) is expressible as  

(9)

where θ= 1 −K_dry/K_m is the Biot–Willis parameter (Biot & Willis 1957) and B is Skempton's coefficient (Skempton 1954)  

(10)

Expanding eq. (9) for small φ gives eq. (8) to first order in φ. Note, however, that Gassmann's full eq. (9) has the further advantage that it is valid for all values of K_f (right down to zero), not just for values in the liquid range.

Eq. (8) is also the result of Mavko & Jizba (1991) for a granite-like material under a high confining pressure so that the crack-like pores are closed. Their result is stated for a sandstone-like material including both crack-like pores and other pores. However, since we have not considered the presence of any other pores except the crack-like pores in this argument, the correct comparison material is just the mineral matrix.

Appendix A shows how to obtain the result of Mavko & Jizba (1991) from a modified DEM scheme.

2.1.2 Gas inclusions: K_f≪παγ_m

For this limit, the stiffness form and the compliance form of the DEM equations are of equal difficulty to integrate, and a further complication arises owing to the presence of shear modulus dependence in the term γ_m in P. We are going to make an approximation (only for analytical calculations) that γ*≃K*[3(1 − 2ν_m)/4(1 −ν²_m)], so the effect of variations in Poisson's ratio away from ν_m for the matrix material is assumed not to affect the results significantly (i.e. to first order) over the range of integration. Without this assumption, the DEM equations for bulk and shear moduli are coupled and must be solved simultaneously (and therefore numerically in most cases).

With this approximation, the equation to be integrated then becomes  

(11)

where  

(12)

(Note that b≃O(α) since ν_m will be in the range 0.05–0.4 for most minerals.) The result of the integration is  

(13)

This result seems to show a very strong dependence of K* on the aspect ratio and Poisson's ratio through the product α (1 − 2ν). However, we show in Appendix B that there exists a special (or critical) value of Poisson's ratio ν_c that serves as a point of attraction during the integration so that ν→ν_c≃πα/18. This result implies that only the dependence on α is significant for gas inclusions.

It seems that this decoupling approximation might have a large effect for a dry system, but an exact decoupling can be achieved in this case (see Appendix B). The result shows that the only significant approximation we have made in eq. (12) is one of order 2 (ν_c−ν_m) and this term is of the order of 20 per cent of b, and usually much less, for all the cases considered here.

2.1.3 General inclusions: K_f and παγ_m arbitrary

Making the same approximations as in the previous case for γ_m, but making no assumption concerning the relative size of K_f and the aspect ratio, we find that DEM gives  

(14)

which can be rewritten in the form  

(15)

It is now easy to check that eq. (15) reduces to eq. (7) when b→ 0 and that eq. (15) reduces to eq. (13) when K_f→ 0.

2.2 Analytical results for G*

We now consider the same three cases for application of DEM to estimating the shear modulus G*.

2.2.1 Liquid inclusions: K_f≫παγ_m

In this limit, the polarization factor for shear is given by  

(16)

where  

(17)

In this case, we have approximated γ*≃G*/2(1 −ν_m) in order to decouple the G* equation from that for the bulk modulus K*. Note that as y→ 1, we anticipate ν*→ 0.5, so for ν* in the usual range from 0 to 0.5 the factor (1 −ν) varies by at most a factor of 2. Therefore, the condition on K_f is not affected. The parameter c depends on a factor (1 −ν)/(2 −ν), which changes at most by a factor of . Thus, we expect some small deviations between the analytical formula and the full DEM for G* in the liquid-saturated case.

Also note that we could argue, in this limit, that the first term on the right-hand side of eq. (16) is dominant (since c∼α≪ 1), and therefore the second term should be neglected. However, for purposes of comparison with Mavko & Jizba (1991), it will prove helpful to retain the second term.

Integrating the DEM equation, we have  

(18)

In the limit of small c (i.e. small α) and φ→ 0, we have  

(19)

which should be contrasted with the result of Mavko & Jizba (1991) for the same problem  

(20)

Because we need some other results to permit the analysis to proceed, a thorough comparison of the present results with the Mavko and Jizba formula will be postponed to Section 4 on the ratio of compliance differences.

2.2.2 Gas inclusions: K_f≪παγ_m

In this second limit, the equation for G* is especially simple, since  

(21)

is a constant (under our constant Poisson's ratio approximation). Again note that d∼O(α). The DEM equation is then integrated to obtain  

(22)

which should be compared with eq. (13). Within the analytical approximation, we will use eq. (22) as our defining equation for G_dry, and note that we can then replace the volume fraction factor 1 −φ by  

(23)

whenever it is convenient to do so.

Our decoupling approximation for the shear modulus in this case turns out to be somewhat better than the corresponding one for the bulk modulus. The result in Appendix B shows that the only significant approximation we have made in eq. (21) is one of order 0.7 (ν_c−ν_m) and this term is of the order of 7 per cent of d or less for all the cases considered here. The relative error is therefore about one-third of that made in the case of the bulk modulus.

2.2.3 General inclusions: K_f and παγ_m arbitrary

In this more general case, we have  

(24)

where  

(25)

Again, the DEM equations can be easily integrated and yield  

(26)

where c is defined in eq. (17) and d is defined in eq. (21). Then, it is easy to check that the two previous cases are obtained when α→ 0 and K_f→ 0, respectively.

3 Theoretical Examples

We now consider some applications of these formulae. We take quartz as the host medium, having K_m= 37.0 GPa and G_m= 44.0 GPa. Poisson's ratio is then found to be ν_m= 0.074.

For liquid saturation, the shear modulus goes to zero as the crack volume fraction increases, while the bulk modulus approaches the bulk modulus of the saturating liquid, which we take as water here (K_f= 2.2 GPa). This means that the effective value of Poisson's ratio increases towards ν*= 0.5 as the crack volume fraction increases, and thus the approximation that ν* is constant clearly does not hold for this case. We therefore expect that the greatest deviations of the analytical approximation should occur for the case of liquid saturation.

In contrast, for the dry case, both the shear modulus and the bulk modulus tend towards zero as the crack volume fraction increases. Thus, since the trends for both moduli are similar, the approximation of constant Poisson ratio might hold in some cases, depending on whether bulk and shear moduli go to zero at similar or very different rates with increasing crack volume fraction.

Hadley (1976) found that Westerly granite has crack aspect ratios ranging from about 0.0001 almost up to 1.0, with a mode around α= 0.001. We will therefore restrict our choice of examples to a subset of this range, picking discrete values of α= 0.1, 0.01, 0.001 and, when we want to show overall trends, we consider 0.001 ≤α≤ 1 for oblate spheroids.

We show three cases in Figs 1–6: (1) α= 0.1 for Figs 1 and 2. (2) α= 0.01 for Figs 3 and 4. (3) α= 0.001 for Figs 5 and 6. The first two cases are easily integrated for DEM. We use two Runge–Kutta schemes from Hildebrand (1956): eqs (6.13.15) and (6.14.5). When these two schemes give similar results to graphical accuracy, we can be confident that the step size used is small enough. If they differ or if either of them does not converge over the range of crack volume fractions of interest, then it is necessary to choose a smaller step size for integration steps. We found that a step size of Δy= 0.01 was sufficiently small for both α= 0.1 and 0.01, while it was necessary to decrease this step size to Δy= 0.001 for the third case, α= 0.001. (Still smaller steps were used in some of the calculations to be described later.)

Figure 1.

The bulk modulus for dry and liquid-saturated cracked porous media with α= 0.1. A full DEM calculation is shown as a solid line for the saturated case and as a dot-dashed line for the dry case. The analytical approximations in the text are displayed as a dashed line for both dry and saturated cases. Agreement between the full DEM calculations and the analytical approximation is excellent in both cases. Gassmann's prediction is shown by the dotted line.

Figure 2.

The shear modulus for dry and liquid-saturated cracked porous media with α= 0.1. A full DEM calculation is shown as a solid line for the saturated case and as a dot-dashed line for the dry case. The analytical approximations in the text are displayed as a dashed line for both dry and saturated cases. Agreement between full DEM calculations and the analytical approximation is again excellent in both cases. The prediction of Mavko & Jizba (1991) is shown by the dotted line. Gassmann (1951) predicts G_dry and G_sat are the same at all porosities for very low frequency responses.

Figure 3.

Same as in Fig. 1 for α= 0.01.

Figure 4.

Same as in Fig. 2 for α= 0.01. Note that the Mavko–Jizba agreement is poor except at low porosities (≲ per cent).

Figure 5.

Same as in Fig. 1 for α= 0.001. The results of Gassmann (1951) are in very good agreement with DEM for this case.

Figure 6.

Same as in Fig. 2 for α= 0.001. Again, note that the Mavko–Jizba prediction is in poor agreement except at very low porosities (≲0.2 per cent).

The results show that our expectations for the agreement between the analytical and numerical results are in concert with the results actually obtained in all cases. The analytical approximation gives a remarkably good estimate of the numerical results in nearly all cases, with the largest deviations occurring—as anticipated—for the intermediate values of crack volume fraction in the cases of liquid saturation for the bulk modulus estimates. We consider that the results of Figs 1–6 are in sufficiently good agreement that they provide cross-validation of both the numerical and the analytical methods.

For the saturated case, we anticipated little if any deviation for the bulk modulus between the analytical results and the full DEM, as is observed for α= 0.1 and 0.01. Larger deviations are found for α= 0.001. We also observed the anticipated small deviations for the shear modulus between the analytical formula and the full DEM.

Note that Gassmann's predictions for bulk modulus are in very good agreement with the numerical DEM results for saturated cracks and α= 0.001. However, the predictions of Gassmann for shear modulus (i.e. that the shear modulus does not depend on the fluid bulk modulus) are clearly violated in all cases.

For the dry case, we anticipated that the analytical shear modulus formula would be a somewhat better approximation of the full DEM than that for the bulk modulus. Both approximations were expected to be quite good. These results are also observed in the figures.

3.1 On improvements

The analytical results obtained here for the dry case could be improved somewhat in several different ways. Instead of replacing ν* by ν_m, we could have replaced it by the fixed-point value ν_c obtained in Appendix B. Since the fixed-point is an attractor and the values rapidly approach ν_c for small but finite volume fractions, this approximation would guarantee an improved approximation over most of the range of crack volume fraction. However, it will make the approximation a little worse in the very small-volume-fraction region. It has been and will continue to be a significant advantage for our analysis to have formulae valid in the small-φ limit, so we have chosen not to do this here. Alternatively, instead of choosing either of the extreme values of ν*, we could use their mean, their harmonic mean or their geometric mean, etc., with similar benefits and drawbacks. Or, we could make direct use of the results from Appendix B for the decoupled equation for Poisson's ratio. This approach will improve the results over the whole range of volume fractions, but will complicate the formulae considerably. We want to emphasize, however, that our goal here has not been to achieve high accuracy in the analytical approximation, but rather to gain an insight into what the equations were computing and why. Having accomplished this even with the simplest approximation ν*≃ν_m, we do not think it fruitful to dwell on this issue and we will therefore leave this part of the subject for now. For the interested reader, some additional technical justifications of the analytical approximation are provided in Appendix C.

Next, we want to make more detailed comparisons between these results and those of Gassmann (1951) and of Mavko & Jizba (1991) in the remainder of the paper.

4 Ratio of Compliance Differences

We have already seen that there are several advantages of the differential scheme presented here for the purposes of analysis. Another advantage will soon become apparent when we analyse the ratio of the compliance differences  

(27)

This ratio is of both theoretical and practical interest. It is of practical interest because it is often easier to measure bulk moduli, and it would therefore be possible to estimate the shear behaviour from the bulk behaviour if the ratio R were known to be either a universal constant, or a predictable parameter. Mavko & Jizba (1991) show that this ratio is given by when the differences between the dry and the starred quantities are caused by a small amount of soft (crack-like) porosity that is liquid filled for the starred moduli. The derivation of this ratio makes it clear that the value of is actually an upper bound on R(0), i.e. a value that cannot be exceeded for such systems, but also a value that clearly is not achieved for many systems lacking such soft porosity. In particular, it was already known by Mavko & Jizba (1991) that R(0) ≃ 0 when the microgeometry of all the porosity is spherical. The crack-like porosity in Mavko and Jizba's model has finite compressibility normal to its plane and is incompressible in the plane of the crack. Thus, their soft porosity can be thought of as cracks for which the aspect ratios approach zero. Goertz & Knight (1998) have also made a parameter study, showing that a related ratio (RG_m/K_m) is generally less than for oblate spheroids and it tends to zero as the aspect ratios of the oblate spheroids approach unity. It would be helpful to see this behaviour directly in the equations.

The purpose of this section is twofold: (1) to derive the Mavko–Jizba result for R(0) analytically and (2) to show, furthermore, that something definite can be said concerning how R(y) changes for small but finite values of y > 0. The second goal is achieved by considering the Taylor series expansion R(y) ⋍R(0) +y dR}(0)/dy, for small values of y.

4.1 Derivation of R(0)

Each of the four material constants appearing in eq. (27) can be computed/estimated using the DEM. However, R is normally defined only in the limit of very small values of soft porosity, in which case both the numerator and the denominator tend to zero. This type of limit is well known in elementary calculus, and the result is given by L'Hôpital's rule:  

(28)

From this form of R, it is now quite easy to relate the ratio to P and Q discussed earlier. In particular, we find that  

(29)

and  

(30)

and therefore that  

(31)

(For sandstones, we could instead evaluate eq. (31) at y=φ₀ and ν(φ₀) =ν_s. It is only the soft, crack-like porosity that needs to be very small for eq. (31) to be applicable.) Eq. (31) is an exact expression for the ratios of these two slopes when the calculation starts at y= 0 and ν(0) =ν_m. It depends only on the aspect ratio α and Poisson's ratio ν_m of the mineral. It shows a sublinear decrease of R with increasing α, and the value of R reaches zero when α_c= 4(1 −ν_m)/3π. Because the formulae used for the penny-shaped crack model are valid only for very low aspect ratios, this latter behaviour should not be taken literally. We do expect R to decrease as the aspect ratio increases, and the trend should be to zero, but this null value should only be achieved at α= 1. This is the type of behaviour observed, for example, by Goertz & Knight (1998). We will check the quantitative predictions by making a numerical study here for oblate spheroids as a function of the aspect ratio. The results will be similar to those obtained by Goertz & Knight (1998), but not identical for several reasons: (1) Goertz and Knight plot RG_m/K_m (instead of R) for the Mori–Tanaka method (Benveniste 1987); (2) the R values presented here are for an infinitesimal change in soft porosity; and (3) the present calculation is (therefore) actually not dependent on the type of effective medium approximation used, only on the Eshelby (1957) and Wu (1966) polarization factors P and Q.

The appropriate expressions for P and Q for oblate spheroids can be found in Berryman (1980b). We repeat the analysis given above in eqs (28)–(31) step by step for oblate spheroids. The results are shown in Fig. 7, together with the results obtained using the penny-shaped cracks as presented already in eq. (31). We see that the results agree completely for α smaller than about 0.001, and are in qualitative agreement over most of the rest of the range. As already discussed, the penny-shaped crack model is a limiting approximation for the oblate spheroids, and deviations from the curve for oblate spheroids do not have physical significance; they merely indicate the degree of error inherent in this choice of approximation scheme. The results for oblate spheroids should be considered rigorous.

Figure 7.

The small porosity limit R(0) of the ratio of compliance differences as a function of the aspect ratio for oblate spheroids and for the penny-shaped crack approximation to oblate spheroids. Note that the asymptotic value for small α is in both cases, in agreement with Mavko & Jizba (1991).

4.2 Derivation of dR/dy|_y=0

The result so far is quite limited because it tells us about the value of the ratio in eq. (27) only for extremely small values of soft porosity. This result would be of more practical value to us if we also knew something concerning the general behaviour of R(y) for finite values of y. The desired information is contained in the first derivative dR/dy|_y=0, which can also be computed analytically, as we now show.

The correct expression for dR/dy at y= 0 must be obtained carefully, since the numerator and denominator are again both vanishing at the same rate with y→ 0. Then, by using either L'Hôpital's rule as before, or equivalently by taking a Taylor series expansion of the numerator and denominator around y= 0, we find that  

(32)

at y= 0.

A complete calculation of all the terms in the numerator of eq. (32) is not necessary if we choose to restrict our attention to the leading-order terms. It is clear that dR/dy∼O(1/α), and—since the denominator is itself of O(1/α)—we need to track only those terms in the numerator that are of O(1/α²).

We find  

(33)

and  

(34)

Because of the difference in the numerator of eq. (32) and the fact that to the order (in powers of α) at which we are now working , we find that the expression for the second derivative of the bulk compliance difference in eq. (33)exactly cancels the first group of terms in eq. (31). So it is only the second group of terms in eq. (34) that survives (to this order) in the numerator of dR/dy.

The computation has therefore been reduced to finding a convenient expression for  

(35)

valid at y= 0. To simplify matters further, we note that (in this limit) because of eq. (28) we can replace the derivative of the shear compliance difference by that of the bulk compliance difference times , and it is most convenient to do so because of the same factor in the denominator of eq. (32). With this goal in mind, the only remaining calculation we need now is an expression for d(ν_dry−ν*)/dy in terms of derivatives of compliance differences. It follows from formulae such as eq. (B1) in Appendix B that  

(36)

Combining eq. (36) with the previous expressions, we have  

(37)

Thus, the final expression for the desired quantity is  

(38)

The result eq. (38) is compared with our numerical results for the ratio of the compliance differences with α= 0.001 in Fig. 8. Agreement is excellent for the slope in the vicinity of y= 0.

Figure 8.

The derivative dR/dy evaluated at y= 0 gives the slope of the ratio of compliance differences R(y) for small values of y with aspect ratio α= 0.001. See eq. (38).

Our discussion of these results is presented at the end of the final section. The next section will provide some preliminary comparisons with experiment.

5 Some Comparisons with Experiment

Although the main thrust of this paper is to present some new theoretical results, we will nevertheless make some preliminary comparisons with known experimental results on rocks in order to clarify the status of the theory and how it is to be used.

5.1 Qualitative analysis

The data that are most commonly used in this context are the data of Coyner (1984) and Coyner & Cheng (1985). For example, these data are the main ones used in the comparisons between theory and experiment by Mavko & Jizba (1991). Qualitatively it is straightforward to see that predictions of Mavko and Jizba for shear wave velocity in Westerly granite (data from Coyner (1984)), although significantly better than the Biot–Gassmann predictions, are also consistently below the data. This trend is entirely consistent with the results determined here. Assuming that the main effects of interest here in the Westerly granite are a result of the presence of small aspect ratio cracks, the corrections to coming from eq. (31) can be neglected, and then eq. (38) shows that the leading corrections to the Mavko and Jizba results are positive corrections to the ratio . This fact translates into an anticipated underprediction of the shear wave velocity; it is not difficult to verify from the equations that if the value of is used instead of the correct, higher value, then the results are always underpredicted. On the other hand, for sandstones, Mavko and Jizba find that the predicted values of shear wave velocity vary from below to above the measured data. This variation is also expected from the present results since, for sandstones, we expect a larger range of aspect ratios to contribute to the results. In this case, formula (31) is important because it implies that the larger aspect ratios will tend to decrease the value of R, whereas eq. (38) implies that the smaller aspect ratios will tend to increase the R value. The tradeoff between these two effects will be complicated in sandstones, and therefore it is not surprising that the observed results are mixed.

5.2 Quantitative analysis

To provide a more quantitative comparison with data, we now use Coyner's data (Coyner 1984; Coyner & Cheng 1985) to compute values of the ratio of compliance differences directly. From such an analysis, we can then see by how much the measured R differs from the nominal value and whether the trend is in agreement with the more detailed predictions given here. We have performed this analysis for both Westerly granite and for Chelmsford granite, as the data for these two cases are expected to be most consistent with the theory presented in this paper, and in particular we expect the formula (38) to be relevant to both these cases.

To make use of eq. (38) for a real granite, we need to consider how it should be modified when there is a distribution of aspect ratios present in the rock. Recall now that the derivation makes explicit use of an assumption that 1/α is a large number (such as 100 or higher) and that 1/α² is a very large number (such as 10⁴ or higher). So we can consider a distribution of aspect ratios having as its largest member α≃ 0.01. For the sake of argument, we will assume that these aspect ratios constitute the majority of the porosity in the system. If this is not the case, the results we show here do not really change, but the notation becomes cumbersome as we then have to keep track of the difference between total porosity and the porosity associated only with these small aspect ratio cracks. With this limitation in mind, our result is  

(39)

where ƒ(ν_m) = (1 + 2ν_m+ 4ν²_m)/(2 −ν_m)², and the sum is only over small aspect ratio cracks. As stated, we are assuming here that φ=∑_iφ_i. Recall that for penny-shaped cracks, φ_i=N_i(4πL³_i/3V)α_i, where L_i and N_i are, respectively, the radius of the penny-shaped crack and the number of cracks having that aspect ratio and radius, while V is the total volume. We see that the ratio φ_i/α_i=N_i(4πL³_i/3V), is independent of the aspect ratio itself, but not of the number and size of the cracks. (This fact is very useful if we have data collected over a range of external pressures, because then according to Walsh's analysis (Walsh 1965), the aspect ratio should be a linear function of the applied pressure, so φ_i and α_i both decrease linearly and their ratio is virtually constant.) Another way of looking at this is to treat the sum over ratios φ_i/α_i as a measure of an effective aspect ratio, being given specifically by  

(40)

Thus, the effective aspect ratio α_eff is the harmonic mean of the pertinent individual aspect ratios, when weighted by porosity in this way. Thus, by measuring R using Coyner's data, we can invert the data for such an effective aspect ratio, without needing to know any of the φ_i and α_i values explicitly, but we would need to know φ.

In Fig. 9, we display a range of the measured data from Coyner (1984) and Coyner & Cheng (1985) for R(y) of Westerly granite and Chelmsford granite, using the theory to determine the corresponding value of φ/α_eff. This step is required since other data are not currently available to determine this ratio independently. The important points to notice are: (1) that most of these data points lie significantly above , as is anticipated by the present analysis and (2) that, for the differential pressure range displayed here (5–50 MPa), the points for each rock tend to cluster in a fairly narrow range of values for both R and φ/α_eff, especially so considering that these data do not always have the precision we would like for the present calculations.

Figure 9.

A plot showing the range of R values found in Coyner's data for Westerly granite and Chelmsford granite. Since the appropriate values of the ratio φ/α are not known, the theoretical expressions (39) and (40) are used to determine the effective ratio φ/α for each measurement, which is also the reason for the apparently perfect fit to the theory. The data presented lie in the differential pressure range 5–50 MPa.

To complete our data comparisons, next we show a cross-plot of the ratio of compliance differences R versus Poisson's ratio for the dry samples of Westerly granite and Chelmsford granite from Coyner's data again in Fig. 10. As before, the data presented lie in the differential pressure range 5–50 MPa. While Poisson's ratio is a fairly smooth function of differential pressure, R is more sensitive since it is computed from the ratios of the differences of measured values. There does not appear to be any significant correlation between Poisson's ratio and R, judging by this figure. The main conclusion seems to be that Chelmsford granite has a relatively constant R≃ 0.6, while Westerly has the value R≃ 0.3, both of which are above the predicted value of Mavko & Jizba (1991).

Figure 10.

A plot showing the correlation of R values with those of Poisson's ratio ν for dry samples in Coyner's data for Westerly granite and Chelmsford granite. The data presented lie in the differential pressure range 5–50 MPa. While Poisson's ratio is a fairly smooth function of differential pressure, R is more sensitive since it is computed from the ratios of differences of measured values. The main conclusion from this figure seems to be that Chelmsford granite has a relatively constant R≃ 0.6 while Westerly has the value R≃ 0.3, both of which are above the predicted value of of Mavko & Jizba (1991).

We will now return to our main line of argument in the final section. Except for results on Poisson's ratio mentioned earlier and presented in Appendix B, we must now leave further analysis of the data to future work.

6 Discussion and Conclusions

6.1 Discussion

We began the paper by pursuing the differential effective medium predictions for the bulk and shear moduli in a cracked material in which the cracks can be either gas-saturated (dry) or liquid-saturated. The DEM equations can be integrated numerically without serious difficulty for the exact model of oblate spheroids of arbitrary aspect ratio, but the full formulae for oblate spheroids are rather involved. In order to make progress on analytical expressions, part of the effort was directed towards studying the penny-shaped crack model of Walsh (1969). This model is not too difficult to analyse if an additional approximation is entertained. The problem for analysis is that the ordinary differential equations for bulk and shear moduli are coupled. If they can be decoupled either rigorously or approximately, then they can be integrated analytically. We accomplished the decoupling for the penny crack model by assuming that changes in Poisson's ratio occurring in those terms proportional to the aspect ratio are negligible to first order. This permits the decoupling to occur and the integration to proceed. We subsequently checked the analytical results against the full DEM integration for penny-shaped cracks, which showed that the analytical results were in quite good agreement with the numerical results.

Then, trying to understand why the analytical results worked so well, we studied the behaviour of Poisson's ratio for the same system, and found that, as the porosity increases, for the dry systems Poisson's ratio tends to a small positive value of the order of πα/18, where α is the aspect ratio, and for liquid-saturated systems it tends towards ½ in all cases. These results permit more detailed error estimates for the analytical formulae showing that errors will always be less than about 5–20 per cent, depending on the aspect ratio and the porosity value.

We have also shown that the Mavko & Jizba (1991) proportionality factor of relating the differences in shear compliances to the differences in bulk compliances for cracked systems with small amounts of soft porosity is indeed an upper bound on R(0) and that this upper bound is approximately achieved for α≤ 0.001. This proportionality factor decreases monotonically with increasing aspect ratio for oblate spheroids, and vanishes identically for spheres at α= 1.

On the other hand, when the amount of soft porosity is not vanishingly small, our results show that dR/dy is always positive and proportional to 1/α. Fig. 8 shows that R(y) starts at and increases by about 50 per cent when the soft porosity is φ≃ 5πα/4 ≃ 0.0039. So 0.4 per cent soft porosity makes a substantial difference to this proportionality factor. Furthermore, the ratio is clearly not bounded above universally by the value , as one might have supposed prior to the present work.

6.2 Conclusions

The analytical approximation made in this paper seems to be very effective at capturing the first-order behaviour of the bulk and shear moduli for the differential effective medium approximation to cracked porous media in both the dry and saturated cases. The resulting formulae may not be rigorous in all cases, unlike Gassmann's formulae for low-frequency behaviour, but these formulae nevertheless have a wider range of approximate validity (considering both porosity and frequency ranges) than either Gassmann's (1951) or Mavko & Jizba's (1991) results. We believe these results will at the very least provide some helpful insight into the behaviour of these complex systems, and may also provide a stepping stone towards more general formulae in the future.

Some of the results, including those of Section 4 and Appendix B, are essentially exact (for α≤ 0.001, which is required so that P and Q for penny-shaped cracks are in accurate agreement with Eshelby's results for oblate spheroids) and independent of the DEM scheme. Any result obtained here at a single value of porosity and for which derivation makes direct use of Eshelby's formulae for ellipsoidal inclusions is not influenced at all by the implicit microgeometry associated with DEM, and therefore has general validity.

A major result of the paper is contained in Section 4, where we show that the ratio of compliance differences R(y), defined in eq. (27), can be expanded about y= 0. Combining the results for the values of R(0) and dR(0)/dy when α is small, we have obtained the two-term Taylor expansion  

(41)

where ƒ(ν_m) = (1 + 2ν_m+ 4ν²_m)/(2 −ν_m)². In order for the truncated expansion to be valid, we must have the condition that y/α∼O(1) or equivalently y∼O(α). This result is true for any effective medium theory based on penny-shaped cracks, not just the DEM theory presented here. Although R(0) [defined as the limit of R(y) for very small amounts of soft porosity] does indeed take the value of found by Mavko & Jizba (1991), the slope dR/dy|_y=0 is, however, positive and inversely proportional to the aspect ratio α. Thus, for small α, the first correction term in the Taylor series expansion can be very significant, and furthermore, when α is very small, these corrections are significant for similarly very small soft porosity values. The Mavko & Jizba (1991) result therefore appears not to be a bound at all, but rather a mid-range estimate, being too large for hard and/or spherical porosity (see Fig. 7) and too small when significant amounts (say φ≃ 0.004) of soft porosity are present (see Fig. 8).

For practical applications of this work, the first-order result appears to be that the ratio of compliance differences R is approximately constant for a given rock, but the constant is usually greater than for granites. We will need to do a separate, and more detailed, data analysis for sandstones since their results appear to be mixed, both in theory and in practice.

Acknowledgments

We thank Gary Mavko for a very helpful discussion of his own results regarding fluid effects on the shear modulus. We thank Jack Dvorkin for a discussion of Poisson's ratio in the weak frame limit. We thank Brian Bonner and Steven Carlson for discussions concerning experimental data and crack models for rocks. We thank Mike Batzle for constructive criticism. We thank Brad Artman and Patricia A. Berge for helpful comments that improved the manuscript. Work performed under the auspices of the US Department of Energy by the University of California Lawrence Livermore National Laboratory under contract no W-7405-ENG-48 and supported specifically by the Geosciences Research Programme of the DOE Office of Energy Research within the Office of Basic Energy Sciences, Division of Engineering and Geosciences. SRP was supported in part by the sponsors of the Stanford Rock Physics and Borehole consortium. HFW was also supported in part by OBES grant DE-FG02-98ER14852.

Appendix

Appendix A: Sandstone-like sample calculations

The main focus of the paper is on the effects of the addition of cracks to pre-existing materials. When the cracks are added to a homogeneous background, we think of this as being a granite-like material. This case has been treated in the main text. To show the generality of the method, we want to give a brief treatment in this appendix of the sandstone-like situation of a material having a pre-existing porosity φ₀. This porosity may itself be either liquid-saturated or gas-saturated (dry). For simplicity, we will assume here that φ₀ is liquid-saturated when liquid-saturated cracks are to be added and dry when dry cracks are to be added. A key assumption for the liquid-saturated, sandstone-like material is that there is no local fluid flow between the pre-existing porosity and the newly introduced cracks by the DEM procedure.

There are two further alternatives to be considered. First, cracks may be added randomly to the pre-existing material. Secondly, cracks may be added preferentially to the porosity-free host material. We discuss these two cases in turn.

Random addition of cracks

Random addition of inclusions is the case considered in the main text, so the DEM equations themselves do not change. We use (1) and (2) as before, but the range of integration changes to the interval starting at y=φ₀ up to y=φ=φ₀+φ_crack. The only differences in the resulting formulae are that: (a) everywhere the factor (1 −φ) appeared before it is now replaced by the ratio (1 −φ)/(1 −φ₀) and (b) everywhere K_m and G_m appeared these material constants are replaced by K(φ₀) and G(φ₀). So, for example, eq. (7) becomes  

(A1)

When φ_crack≪ (1 −φ₀), we can rewrite this expression as  

(A2)

This formula still differs from the result of Mavko & Jizba (1991), but the difference is nevertheless expected because their derivation does not assume random placement of cracks. We can resolve this discrepancy when we make use of an assumption of preferential addition of cracks in the next subsection.

Preferential addition of cracks

The factor of (1 −y) on the left-hand sides of both eqs (1) and (2) arises from the need to account for the fact that, when an inclusion is placed in a composite, the volume of the inclusion replaces not only host material, but also some of the other inclusion material previously placed in the composite. When y is the inclusion volume fraction, the remaining host volume fraction is (1 −y). So random replacement of dy of the composite medium only replaces (1 −y) dy of the host material. Replacing instead dy/(1 −y) of the composite then gives the correct factor of dy host replacement; thus, the factor of (1 −y) is required in (1) and (2) for random inclusion placement at finite values of y.

If we now assume instead that the inclusions are placed preferentially in pure host material (and this becomes progressively harder to do in practice for larger integrated overall inclusion fractions y), then the DEM equations must be modified to account for this situation.

For example, with preferential addition of inclusions, it is clear from the preceding considerations that DEM eq. (6) is replaced by  

(A3)

Integrating eq. (A3) gives  

(A4)

The validity of this result clearly depends on φ₀ being sufficiently small so that it is possible to find enough pure host material to which cracks can be added ‘randomly’. Taking φ₀→ 0 guarantees satisfaction of the requirement, but the approximation must eventually breakdown as φ₀→ 1.

Eq. (A4) is almost the corresponding result of Mavko & Jizba (1991). Mavko & Jizba use as their comparison state the dry porous material, assuming that no cracks are present or that, when present, they are closed owing to applied external pressure. We can also obtain the same result using eq. (A3), but now K*(φ₀) =K_dry(φ₀), so the integration has a different starting value from that in the previous paragraph. Then, we find  

(A5)

Eq. (A5) is exactly the corresponding result of Mavko & Jizba (1991). Although the right-hand sides of eqs (A4) and (A5) are identical, the results differ, i.e. K*≠K*_MJ, since the assumed host material is fluid-saturated in the first case and dry in the second case.

Appendix B: Poisson's ratio for dry cracks

When the cracks are taken to be dry, so that K_i=G_i= 0 in (1) and (2), it turns out that an elegant decoupling of the DEM equations is possible (also see Zimmerman 1985; Hashin 1988). If we consider the parameter ratio G*/K*= 3(1 − 2ν*)/2(1 +ν*), we find that it satisfies the equation  

(B1)

Furthermore, it is generally true for dry inclusions (not just for penny-shaped cracks) that both P*ⁱ and Q*ⁱ are functions only of the same ratio G*/K*, or equivalently of Poisson's ratio ν*. Thus, we can solve eq. (B1) for either ν* or the ratio of moduli, without considering any other equation.

We notice that the dimensionless polarization factors P and Q are both often close to unity, and furthermore it may happen that, for special values of Poisson's ratio, these terms are equal, so P*ⁱ=Q*ⁱ. If this occurs for some critical value ν*=ν_c, then eq. (B1) guarantees that this value of Poisson's ratio will be preserved for all higher values of porosity, since the right-hand side vanishes initially, and therefore always thereafter in the integration scheme. Such a critical value is usually called a fixed-point of the equations, and such fixed-points can be either stable or unstable. If they are unstable, then a small deviation from the critical point causes a rapid divergence of Poisson's ratio from the fixed-point. If they are stable, then a small deviation produces a situation in which the value of Poisson's ratio gradually (asymptotically) approaches the critical value. When this happens, we say the fixed-point is an attractor. For the DEM eq. (B1), a fixed-point that is an attractor will only be reached in the limit φ→ 1, but the value of Poisson's ratio will change fairly rapidly in the direction of the attractor when the first cracks are added to the system. Such behaviour of Poisson's ratio has been noted before by Zimmerman (1994) and by Dunn & Ledbetter (1995), among others.

For dry penny-shaped cracks, we have  

(B2)

which has a fixed-point approximately (using one step of a Newton–Raphson iteration scheme) at  

(B3)

This shows that, when α is very small, Poisson's ratio for the dry cracked material tends toward small positive values. For somewhat larger values of α, Poisson's ratio approaches a value proportional to α and of the order of πα/18.

For comparison, consider spherical void inclusions (see Berryman (1980) for the general expressions for P and Q). Then, we have  

(B4)

which has a fixed-point at  

(B5)

This result has been remarked upon previously by Zimmerman (1994). Similarly, considering needle-shaped void inclusions, we have  

(B6)

which has a fixed-point at  

(B7)

Dunn & Ledbetter (1995) have shown that all the prolate spheroids have critical Poisson's ratios close to that for spheres. Since needles are the extreme limit of prolate spheroids, we see that eqs (B5) and (B7) are in agreement with their results.

Dunn & Ledbetter (1995) have shown that disc-shaped inclusions (which are achieved by taking oblate spheroids to the α= 0 limit) have a critical Poisson's ratio of ν_c= 0. This result and the others obtained above are collected for comparison in Table B1.

Table B1.

Fixed-points of equation (B1) for commonly considered inclusion shapes. Also listed for comparison are the ratios of compressional (v_p) and shear (v_s) wave velocity, where .

To clarify the behaviour of the solution of eq. (B1), we will do an approximate analysis by expanding the right-hand side around ν= 0 and also note that for small ν, G*/K*≃ 3(1 − 3ν*)/2. Then, eq. (B1) becomes  

(B8)

which can be integrated easily to yield  

(B9)

where ν_m is the starting, or in our case the mineral, value of Poisson's ratio, and  

(B10)

A more precise, and therefore more tedious, analysis of the right-hand side of eq. (B2) gives the improved approximation eq. (B3) for the asymptotic value of ν_c.

In Fig. B1, we show the actual results for Poisson's ratio from the full DEM in the same three examples shown in Figs 1–6. The starting value of Poisson's ratio for quartz is ν_m= 0.0742. For comparison, Table B2 contains a listing of various Poisson's ratios for minerals that could be important in rocks in order to show the range of behaviour observed in nature. Except for different starting locations, we expect the qualitative behaviour of the curves for Poisson's ratio to follow closely that observed in Fig. B1 and eq. (B9) for all cases.

Figure B1.

Asymptotic behaviour of Poisson's ratio as a function of the crack volume fraction for three values of α: 0.1, 0.01, 0.001. The asymptotic value for saturated samples is always . For dry samples, the asymptotic value depends on the geometry of the inclusion, and therefore on α for cracks. The limiting value ν_c≃πα/18 is a stable attractor of the DEM equations, as is observed in this figure.

Table B2.

Typical values of Poisson's ratio for various solid materials contained in rocks. See, for example, Mavko et al. (1998).

(Technical notes concerning the dry case. (1) For α= 0.1, the Runge–Kutta scheme used to solve the coupled DEM equations for K* and G* was sufficiently accurate that ν* could be computed from these values. However, for α= 0.01 and 0.001, the accuracy obtained was not sufficient, so we instead used the same Runge–Kutta scheme but applied it directly to eq. (B1). This approach gave very stable results. (2) The approximations used here assume that the resulting bulk and shear moduli are always much larger than those of gas or air, so that taking K_air≃ 0 is a sensible approximation. When this is not true, i.e. when the effective moduli are so weak that they become comparable in magnitude to K_air, then the results revert back to those for saturated media and the asymptotic result approaches ν_c= 0.5.)

Fig. B2 compares the results for oblate spheroids with those for penny-shaped cracks; both curves were obtained by finding the zeros of P−Q numerically. To provide additional insight, the curve (which was obtained using the functional form of eq. (B3) and fitting the coefficient in the denominator at α= 1) is also shown. We see that the results for penny-shaped cracks deviate substantially from those of oblate spheroids as α→ 1, but they come into agreement at lower values of α≤ 0.001. The deviations from the results for oblate spheroids, again, are not physical and should simply be viewed as artefacts introduced by the very low aspect ratio limiting procedure used to obtain the approximate formulae for penny-shaped cracks.

Figure B2.

Poisson's ratio fixed-point ν_c as a function of α found numerically for oblate spheroids and penny-shaped cracks, and also for penny-shaped cracks using the analytical expression ν_c= 2πα/(36.0 + 2.245πα). The two curves for penny-shaped cracks are nearly indistinguishable on the scale of this plot. The correct fixed-point for spheres (α= 1) is , and this value is attained in the α→ 1 limit by the curve for oblate spheroids.

Some data analysis

It is commonly stated that it is not possible to determine actual values of aspect ratio distributions from compressibility data alone (Walsh 1965; Mavko & Nur 1978). This is probably true in general, owing to the observed (here as well as elsewhere) relative insensitivity of K to α. However, the present results show a strong sensitivity of G to α, and therefore of ν to α. Eqs (B1) and (B2) can thus be used to show directly just how strongly ν* should depend on α.

Assuming as in the main text that 1/α is a large number so terms of O(1) (independent of 1/α) can be neglected, we find that eqs (B1) and (B2) imply  

(B11)

(Compare eq. (B11) with eq. (B8).) Although this equation has a form that can be integrated explicitly analytically, it will be more useful to us to make a small porosity φ approximation and then easily arrive at the result  

(B12)

where g(ν_m) = (1 −ν²_m) (3 −ν_m)/(2 −ν_m) and the ratio φ/α has the same significance as in eqs (40) and (39) of the main text.

The useful properties of eq. (B12) are that the right-hand side depends only on the pure host material Poisson ratio ν_m, and the left-hand side is determined by measuring v_p and v_s for dry samples of the rock. This means that there are fewer variables than there are in the determination of R and so we imagine that it should be possible to find estimates of the ratio φ/α from data by making use of eq. (B12).

Again considering Coyner's data (Coyner 1984; Coyner & Cheng 1985), we find ν results for dry samples of Westerly granite and Chelmsford granite in the range 5–50 MPa. These results are shown in Fig. B3. It is not clear why the predicted values of φ/α found in this way differ from those found in Fig. 9 using R as the estimator. Nor is it clear which estimates are likely to be more correct. Since the R values are found by taking the ratio of two differences of measured quantities and since the R values are, in fact, observed to be a bit erratic (whereas the ν values are generally smoothly varying as the differential pressure changes), it seems likely that the values found here using the ν measurements are more trustworthy, but clearly more work is needed to confirm this.

Figure B3.

Plot showing the range of ν values found in Coyner's data for dry samples of Westerly granite and Chelmsford granite. Since the appropriate values of the ratio φ/α are not known independently, the theoretical expression (B12) is used to determine the effective ratio φ/α for each measurement, which is also the reason for the apparently perfect fit to the theory. The data presented lie in the differential pressure range 5–50 MPa.

Appendix C: Technical justification of the approximation for γ

It is inherent in the mathematical form of all DEM schemes that they always give correct values and slopes of the curves for small values of the inclusion volume fraction, and that they always give the right values (but not necessarily correct slopes) at high volume fractions. We see that these expectations are fulfilled in all the examples shown here.

The approximations made in the text to arrive at analytical results were chosen as a convenient means to decouple the equations for bulk and shear moduli, which are normally coupled in the DEM scheme. For the liquid-saturated case, the approximations for bulk modulus are very good for all values of aspect ratio, but for shear modulus the exponent determined by eq. (17) can deviate as much as a factor of . The value chosen is the maximum value possible, guaranteeing that the analytical approximation will always be a lower bound for this case.

In contrast, for the case of dry cracks, the approximations for the shear modulus are expected to be somewhat better than those for the bulk modulus. The analytical approximation is again expected to be a lower bound for the full DEM result for the shear modulus. Analysis for the bulk modulus is more difficult in this limit as it requires checking that the ratio G*/K* remains finite as the porosity φ→ 1, and this would be difficult to establish if Poisson's ratio were going to , as it does for the liquid-saturated case. However, Appendix B shows that Poisson's ratio actually tends to a value of about ν_c≃πα/18, so there is no singularity in the K* behaviour for this case. This feature is also confirmed by the numerical results.