next up [*] print clean
Next: About this document ... Up: Table of Contents




Stress-induced transverse isotropy in rocks

Lawrence M. Schwartz,[*] William F. Murphy, III,[*]
and James G. Berryman

Author has no known email address


The application of uniaxial pressure can induce elastic anisotropy in otherwise isotropic rocks. We consider models based on two very different rock classes, granites and weakly consolidated granular systems. We show (1) that these models share common underlying assumptions, (2) that they lead to similar qualitative behavior, and (3) that both provide a microscopic basis for elliptical anisotropy. In contrast, a finely layered transversely isotropic medium always shows anelliptical anisotropy. In the granular case, we make experimentally verifiable predictions regarding the horizontally propagating modes based on the measured behavior of the vertical modes.


It is well known that the elastic properties of the rocks comprising the earth's sub-surface are highly non-linear (Walsh, 1965a,b,c). In terms of seismic and sonic wave propagation, the signature of this non-linearity is the pressure dependence of the P and S sound speeds. An additional complication is the fact that the sub-surface pressure environment is often anisotropic. Accordingly, it is essential for the proper interpretation of both seismic and sonic measurements that we have a semi-quantitative understanding of the effects of stress anisotropy.

In the present paper we examine the uniaxial response of systems whose elastic properties are isotropic under the application of hydrostatic stress. While the P and S velocities of such systems exhibit clear pressure dependence, the VP/VS ratio is often independent of the applied pressure (Domenico, 1977). However, once the applied stress is uniaxial, these systems exhibit transversely isotropic (TI) behavior and the three VP/VS ratios depend on the applied stress (Nur and Simmons, 1969; Murphy, 1982; Zamora and Poirier, 1990; Yin and Nur, 1992).

We study two distinct types of models that predict pressure induced velocity anisotropy in rocks: (1) granular materials under combined hydrostatic plus uniaxial loading and (2) rocks with randomly oriented cracks under uniaxial loading. The first class of models was developed by Schwartz et al. (1984), Schwartz (1984), and Walton (1987) and is appropriate for weakly consolidated granular materials. The second model was developed by Walsh (1965a,b,c) and Nur (1969; 1971) to describe the properties of granites. While these models are directed toward very different rock classes, we will see that the physical bases of the induced anisotropy are quite similar, as are their qualitative predictions.

In the model developed by Walsh (1965a,b,c) and Nur (1969; 1971), the rock is represented by an isotropic array of penny-shaped cracks. Under uniaxial compression, the normal stress acting on each contact is assumed to vary as
\sigma_n = \sigma \: {\cos}^2 \psi
 \end{eqnarray} (1)
where $\psi$ is the angle between the crack normal to the stress axis. Basically, the ${\cos}^2 \psi$ dependence is the simplest variation that is consistent with the symmetry of the problem. As the value of $\sigma$ is increased, cracks oriented with their normals close to the pressure axis close and the elastic properties of the rock become anisotropic. In the granular models developed by Schwartz (1984), adjacent grains are coupled (at their contacts) by effective normal and tangential springs with force constants Dn and Dt. Once again, under uniaxial stress, these quantities are assumed to vary as
D_n = D^{(0)}_n \left[ 1 + \delta \cos^2 \psi \right],\nonumber\\ D_t = D^{(0)}_t \left[ 1 + \delta \cos^2 \psi \right],
Here $\psi$ is now understood to be the angle between the contact normal and the pressure axis and $\delta$ is proportional to the difference between the uniaxial and transverse applied stresses. Note, that for unconsolidated granular media (e.g., sand packs) it is essential that there be some transverse confining pressure to give the material an underlying elastic integrity (i.e., to prevent the vanishing of D(0)n and D(0)t). The same assumption is required in the formalism developed by Walton (1987) and we show here that the application of his formalism leads to results that are identical to those derived by Schwartz (1984).


Let x1, x2, and x3 be spatial coordinates and u1, u2, and u3 be the displacement components of an elastic wave in a rock. If the stress tensor is $\sigma_{ij}$ and the strain tensor is related to the displacements by
\epsilon_{ii} = {{\partial u_i}\over{\partial x_i}},
 \end{eqnarray} (3)
\epsilon_{ij} = {{\partial u_i}\over{\partial x_j}} +
 ...}\over{\partial x_i}} \quad\hbox{for}\quad i,j = 1,2,3,\,i\ne j,
 \end{eqnarray} (4)
then the general stress-strain relations in elastic rocks are
\sigma_{ij} = c_{ij}\epsilon_{kl} \quad\hbox{or}\quad \epsilon_{ij} = S_{ijkl}\sigma_{kl},
 \end{eqnarray} (5)
where cij is the fourth rank stiffness tensor and Sijkl is the corresponding compliance tensor. Repeated indices are summed in (5), but not in (3). The tensor notation may be conveniently replaced by vectors and matrices using the Voigt notation, whereby the subscripts of stress and strain are mapped according to the prescription $11 \to 1$, $22 \to 2$, $33 \to 3$, $23 \to 4$, $31 \to 5$, and $12 \to 6$.Symmetries take care of the remaining combinations.

Assuming that the symmetry axis is in direction x3, so the x1x2-plane is isotropic, then the relation between stress and strain for such transversely isotropic media becomes
\pmatrix{\sigma_{11}\cr \sigma_{22}\cr \sigma_{33}\cr
\epsilon_{23}\cr \epsilon_{31}\cr \epsilon_{12}\cr}.
 \end{eqnarray} (6)
A similar equation relates $\epsilon_{ij}$ to $\sigma_{ij}$ through the compliance matrix ${\bf S}$, so the stiffness matrix is just the inverse of the compliance matrix:
\pmatrix{C_{11}& C_{12}& C_{13}& & & \cr
 C_{12}& C_{11}& C_{13...
 ...& S_{55}& & \cr
 & & & & S_{55}& \cr
 & & & & & S_{66}\cr}^{-1}.
 \end{eqnarray} (7)
For transversely isotropic materials, there are only five independent constants, although six constants appear in each of these matrices. The remaining condition on the constants is C11= C12+ 2C66 for the stiffnesses, and $S_{11}= S_{12}+ {1\over2}S_{66}$ for the compliances. These conditions follow easily from the requirement that a rotation of the coordinate system in the x1x2-plane should not change the constants for a transversely isotropic material.

Taking the displacement vector ${\bf u}$ to be a plane wave proportional to $\exp i({\bf k}\cdot{\bf x}- \omega t)$ and letting $\rho$ be the density of the medium, representative characteristic dispersion relations for propagating waves are easily shown (Berryman, 1979) to be
\rho\omega_{\pm}^2 = {\textstyle {1\over2}}\biggl\{(C_{11}+C_{5...
+ 4(C_{13}+C_{55})^2k_1^2k_3^2}\biggr\},
when the polarization is normal to direction x2 (so u2 = 0), and similarly
\rho\omega_{sh}^2 = C_{66}k_1^2 + C_{55}k_3^2,
 \end{eqnarray} (9)
when the polarization is purely normal to the direction of propagation (so u1 = u3 = 0). Equation (8) gives two dispersion relations for waves having mixtures of compressional and shear polarizations, being neither pure compressional nor pure shear except in the x1 and x3 directions. In these special directions, $v_+ = \omega_+/k$ is the velocity of a pure compressional (P) wave and $v_- = \omega_-/k$ is the velocity of a pure shear (SV) wave. For intermediate angles, these two waves are known respectively as quasi-P and quasi-SV waves. Equation (9) is always (for all angles) the dispersion relation for a pure shear (SH) wave.


We use the results of Walton (1987) to obtain predicted stiffnesses for a granular material comprised of a dense random packing of spheres with constant radius under hydrostatic confining strain with an additional uniaxial strain applied in the x3 direction. Walton gives general expressions for elastic stiffnesses as a function of an arbitrary applied strain. We will not quote his general expression here, but merely write down the results we obtain from his formula using an applied macroscopic strain of the form
\epsilon_{ij} = \epsilon\delta_{ij} + \Delta\epsilon_3\delta_{i3}\delta_{j3},
 \end{eqnarray} (10)
where $\epsilon$ is a uniform hydrostatic strain and $\Delta\epsilon_3$ is the additional uniaxial strain in the x3 direction. Then, introducing the direction cosines $\xi_i$ for the unit vector connecting the centers of adjacent spheres, we find that the stiffness coefficients depend on an average applied strain given by
\bigl<(-\epsilon_{ij}\xi_i\xi_j)^{1\over2}F\bigr\gt \simeq
 ...+{{\Delta\epsilon_3}\over{2\epsilon}}\bigl<\cos^2\psi F\bigr\gt)
 \end{eqnarray} (11)
where $\bigl<\bigr\gt$ is the average over all possible orientations of the unit vector, $\psi$ is again the angle with respect to the direction of the applied strain while F is a complicated function of direction cosines [see Walton (1987)].

All the formulas that follow are proportional to a constant factor of the form
\gamma= {{3n(1-\phi)(-\epsilon)^{1\over2}}\over{4\pi^2B(2B+C)}},
 \end{eqnarray} (12)
where n is the average number of contacts per spherical particle, $\phi$is the porosity, and
B = {{1}\over{4\pi}}\Biggl[{{1}\over{\mu}}+{{1}\over{\lambda+\m...
 \end{eqnarray} (13)
where $\lambda$ and $\mu$ are the Lamé constants of the mineral composing the spherical grains.

With these definitions, we find that the elastic stiffnesses are
C_{11}= \gamma\Biggl[{{4B}\over{3}}+{{2C}\over{5}} +
 ...\epsilon}}\biggl({{2B}\over{15}} + {{C}\over{35}}\biggr)\Biggr],
 \end{eqnarray} (14)
C_{13}= \gamma\Biggl[{{2C}\over{15}} 
+ {{\Delta\epsilon_3}\over{\epsilon}}\biggl({{C}\over{35}}\biggr)\Biggr],
 \end{eqnarray} (15)
C_{33}= \gamma\Biggl[{{4B}\over{3}}+{{2C}\over{5}} +
 \end{eqnarray} (16)
C_{55}= \gamma\Biggl[{{2B}\over{3}}+{{2C}\over{15}} +
 \end{eqnarray} (17)
C_{66}= \gamma\Biggl[{{2B}\over{3}}+{{2C}\over{15}} +
 \end{eqnarray} (18)
The remaining constant is given by C12= C11- 2C66. Eqs. (14) through (18) are, in fact, the same as those derived by Schwartz (1984) using the orientation dependent force constants (2). To establish this correspondence, we employ the identities
{ {D^{(0)}_t} \over {D^{(0)}_n} } = 
{ 2(1 - \nu) \over (2 - \nu)} = 
{2 B \over 2B + C}
 \end{eqnarray} (19)
and set $2 \delta = \Delta \epsilon_3 / \epsilon$. The parameter $\nu$appearing in (19) is the Poisson ratio of the sphere material.

Figure 1
The pluses (+) and solid curve represent experimental (Murphy, 1982) and calculated (VP/VS)2 ratios for propagation along the pressure axis. The (essentially) level dashed curve is the corresponding ratio calculated for propagation in the transverse direction with the shear wave polarized in the transverse plane. The decreasing dashed curve is the calculated ratio for transverse propagation with shear polarization in the axial direction.

Equations (14) through (18) can be directly employed to calculate the three independent (VP/VS)2 ratios: C33/C55, C11/C66, and C11/C55. In Figure 1, the results of such a calculation are compared with experimental data on packed and well sorted (grain diameters in the range $106 \rightarrow 125 \: \mu m$) Ottawa sand for the first of these ratios. In these measurements, direct uniaxial pressure p3 was applied to the sand pack which was subject to a nearly zero-strain boundary condition the transverse directions. [Note that, in the isotropic limit, these data approach the results obtained by Domenico (1977).] Because the transverse components of the stress were not measured independently, we cannot be certain of the relation between the applied pressure and the parameter $\Delta \epsilon_3 / \epsilon$. Accordingly, we adopted the empirical fit $\Delta \epsilon_3 / \epsilon = 0.47 p_3/(1 + 0.058 p_3)$.Clearly, the agreement with the measured values is excellent and the predicted results for the horizontally propagating modes can be directly tested by experiment.


Next we consider the work of Walsh (1965a,b,c) and Nur (1969; 1971). Originally their model was formulated in terms of crack closure in granites. We believe, however, that models of this type can provide a phenomemological basis for a wide variety of anisotropic rocks. Accordingly, we have developed a version of their formalism, based on more general arguments about relations among the compliance matrix elements for any TI material. Models of this kind are most conveniently treated in terms of the compliance matrix. It is therefore useful to recall the relations between the Sij and the technical constants (Young's moduli Eii, Poisson's ratios $\nu_{ij}$, and shear moduli $\mu_{ij}$):
\pmatrix{S_{11}& S_{12}& S_{13}& & & \cr
 S_{12}& S_{11}& S_{13...
 ...1} & & \cr
 & & & & 1/\mu_{31} & \cr
 & & & & & 1/\mu_{12} \cr}.

The basic model is that cracks are present in the rock and uniformly distributed over angles before the application of either a hydrostatic confining pressure or a uniaxial pressure. When pressure is applied, cracks close if the pressure normal to the plane of a crack exceeds a threshold, i.e., $p \ge \alpha E_o$, where $\alpha$ is the aspect ratio of the crack and Eo is the Young's modulus of the rock with all cracks ``open.'' If a hydrostatic pressure is applied, then the moduli change but the rock remains isotropic as the cracks close uniformly in all directions. However, if uniaxial pressure is applied, then those cracks oriented normal (or nearly normal) to the symmetry axis will be closed preferentially.

It is assumed that the cracked rock is isotropic before the application of uniaxial pressure field $\Delta p_3$. Then, in terms of the isotropic constants for the mineral composing the rock (without cracks) E, $\nu$, and $\mu= E/2(1+\nu)$, and certain factors (in brackets) dependent on the particular element of the compliance, the compliances are given by
S_{11}= {{1}\over{E}}\biggl[1 + m \int N(\alpha) I_{11}(\alpha)\,d\alpha\biggr],
 \end{eqnarray} (21)
\quad S_{13}= -{{\nu}\over{E}}\biggl[1 + m \int N(\alpha) I_{13}(\alpha)\,d\alpha\biggr],
 \end{eqnarray} (22)
S_{33}= {{1}\over{E}}\biggl[1 + m \int N(\alpha) I_{33}(\alpha)\,d\alpha\biggr],
 \end{eqnarray} (23)
S_{55}= {{1}\over{\mu}}\biggl[1 + m \int N(\alpha) I_{13}(\alpha)\,d\alpha\biggr],
 \end{eqnarray} (24)
S_{66}= {{1}\over{\mu}}\biggl[1 + m \int N(\alpha) I_{12}(\alpha)\,d\alpha\biggr].
 \end{eqnarray} (25)
The remaining constant is given by the rotational invariance condition $S_{12}= S_{11}- {1\over2}S_{66}$. The various new terms appearing in these equations are $N(\alpha) =$ the crack aspect ratio distribution, $m = 16(1-\nu_o^2)(5-2\nu_o^2)a^3/6(2-\nu_o)$ where the individual crack volume is given by $4\pi a^3\alpha/3$ (a being the sphere radius). The other factors in the integrands are
I_{11}(\alpha) = I_{12}(\alpha) = {{2\pi}\over{3}}
 \end{eqnarray} (26)
I_{33}(\alpha) = {{4\pi}\over{3}}\cos^3\psi_c(\alpha),
 \end{eqnarray} (27)
where the critical angle $\psi_c(\alpha) = 0$ for crack closure if $\Delta p_3 < \alpha E_o$, while for larger uniaxial pressures
\cos^2\psi_c(\alpha) = {{\alpha E_o}\over{\Delta p_3}}\quad\hbox{for}\quad
\Delta p_3 \ge \alpha E_o.
 \end{eqnarray} (28)
The condition that $I_{11}(\alpha) = I_{12}(\alpha)$ is required for transversely isotropic media, following again from rotational symmetry of the x1x2-plane. Finally, we have another condition that must be fulfilled by transversely isotropic media relating I13 to I11 and I33, i.e.,
I_{13}(\alpha) = {\textstyle {1\over2}}[I_{11}(\alpha) + I_{33}(\alpha)].
 \end{eqnarray} (29)
For values of uniaxial pressure such that $\Delta p_3 < \alpha_{min} E$, all the factors $I_{ij} = 4\pi/3$, since then the critical angle vanishes. These equations have all the required symmetries to be those of a transversely isotropic medium.

For comparison, we note that in the absence of applied pressure all cracks are open and this model assumes the compliances all depend on the crack aspect ratios through a common factor
\biggl[1 + m \int N(\alpha)\,d\alpha\biggr],
 \end{eqnarray} (30)
i.e., the brackets in (21)-(25) are all replaced by (30) and, for example, $S_{11} = 1/E_o = (1/E)[1+m\int N\,d\alpha]$.Similarly, as hydrostatic pressure is applied, the model assumes that the rock remains isotropic as some of the cracks close. Now the common factor is
\biggl[1 + m \int N(\alpha) \cos\psi_c(\alpha)\,d\alpha\biggr],
 \end{eqnarray} (31)
where the critical closure angle is again given by (28). The factor in (31) is particularly important because identity (29) follows from the condition that an isotropic average of the transversely isotropic results must agree with this hydrostatic result, and so together with (27) determines (26).

Figure 2
Calculated (VP/VS)2 ratios based on the crack closure model. The upper curve represents propagation along the pressure axis. The middle curve is for propagation in the transverse direction with the shear wave polarized in the transverse plane. The lower curve is for transverse propagation with shear polarization in the axial direction. The calculations were based on a reasonable set of assumptions for granites: all the cracks were taken to have aspect ratio 550.0, the porosity was 0.05, and the VP and VS values for the rock with all cracks closed were 5.90 and 3.65 km/sec.

Illustrative calculations based on the above equations are presented in Figure 2. Here, as in Figure 1, the same three (VP/VS)2 ratios are displayed as a function of uniaxial stress. Nur (1969) and Nur and Simmons (1969) have shown that, with suitably chosen distributions $N(\alpha)$, this model provides a reasonable fit to experimental measurements on granites. In the present context, our objective is simply to emphasize the qualitative similarity between the predictions of the crack closure and grain contact models. Accordingly, the calculations shown in Figure 2 were based on a particularly simple distribution in which all cracks were assumed to have the same aspect ratio. The common features of Figures 1 and 2 are quite striking, particularly in light of the very different assumptions underlying the two models.


It has been shown by Berryman (1979) that, if the transverse isotropy is due to fine layering, then (C11-C55)(C33-C55) - (C13+C55)2 is always positive. This combination of stiffnesses will be referred to as the anellipticity parameter ${\cal A}$. Effective medium calculations by Hornby et al. (1993) produce stiffness coefficients that also exhibit a positive anellipticity parameter. The significance of these results is that it allows us to determine, via measurements of elastic velocities in transversely isotropic rocks, whether the anisotropy is due to layering or due to stress. Within the grain contact framework, we find
{\cal A}\equiv (C_{11}-C_{55})(C_{33}-C_{55}) - (C_{13}+C_{55})^2 \simeq
 \end{eqnarray} (32)
showing that this characteristic quantity is small (and may vanish identically) for finite uniaxial strains in this theory. In addition, our numerical calculations indicate that the crack closure model also leads to a vanishing anellipticity parameter for the simple delta function distribution considered here. Thus it may well be true that measurements of velocities are able to distinguish between layered (${\cal A}\gt$)and pressure-induced (${\cal A}\simeq 0$) TI media. Whether the accuracy of measured elastic constants will be sufficient to distinguish these cases is uncertain at present, but this issue may be resolved by a higher order analysis of the anellipticity parameter. In particular, the sign of (32) is uncertain, since various terms of second and higher order were neglected in deriving the result. It is possible that (32) is actually negative for stress-induced anisotropy. However, we should emphasize that at present we do not know of any models that consistently give a negative value for the anellipticity parameter.

Also of interest are Thomsen's anisotropy parameters (Thomsen, 1986) $\varepsilon= (C_{11}-C_{33})/2C_{33}$ and $\delta= [(C_{13}+C_{55})^2-(C_{33}-C_{55})^2]/2C_{33}(C_{33}-C_{55})$.Evaluating these parameters for the grain contact model to lowest order in the applied uniaxial stress, we find that
\varepsilon\simeq \delta\simeq - \biggl({{\Delta\epsilon_3}\ove...
 \end{eqnarray} (33)
neglecting terms of second order in $(\Delta\epsilon_3/\epsilon)$. Similarly, in the crack closure model our numerical results indicate that $\delta
\approx \epsilon$. When these two anisotropy parameters are equal, we have the special case known as ``elliptical anisotropy.'' This fact also follows from (32), since $2C_{33}(C_{33}-C_{55})(\varepsilon-\delta) = (C_{11}-C_{55})(C_{33}-C_{55})-(C_{13}+C_{55})^2$.Elliptical anisotropy occurs when the right hand side of this expression vanishes identically. In contrast, it is known that finely layered materials always satisfy $\delta< \varepsilon$, and so have quasi-P and quasi-SV waves that are never elliptically anisotropic. Helbig (1983) has emphasized that finely layered TI media virtually never display elliptical anisotropy. However, for small uniaxial stresses, the present results show that both granular materials and systems with penny-shaped cracks exhibit elliptical anisotropy in both the quasi-P and quasi-SV waves (Nur and Simmons, 1969). The significance of these results for these models of stress induced TI behavior must be carefully evaluated.


We thank Ken Winkler for suggesting the analysis of the crack closing model. We thank Hezhu Yin and Gary Mavko for helpful conversations about related experiments at Stanford. We thank Brian Bonner and Pat Berge for a number of important discussions.


next up [*] print clean
Next: About this document ... Up: Table of Contents
Stanford Exploration Project