Survey

* Your assessment is very important for improving the workof artificial intelligence, which forms the content of this project

Survey

Transcript

Geophysical Journal International Geophys. J. Int. (2010) 182, 375–388 doi: 10.1111/j.1365-246X.2010.04614.x Seismic velocity anisotropy of phyllosilicate-rich rocks: characteristics inferred from experimental and crack-model studies of biotite-rich schist O. Nishizawa1 and K. Kanagawa2 1 Institute for Advanced Industrial Science and Technology, Central 7, 1-1-1 Higashi, Tsukuba, Ibaraki 305-8567, Japan. E-mail: [email protected] of Earth Sciences, Chiba University, Chiba 263-8522, Japan 2 Department SUMMARY Seismic velocity anisotropy of biotite schist (30 per cent-mode biotite) was measured under confining pressures up to 150 MPa. The rock shows weak orthotropy which was altered from transverse isotropy (TI) generated by biotite-preferred orientation. The orthotropy was caused by microfolding in the rock. The velocity increase under confining pressure indicates that most crack planes are aligned parallel to the cleavage planes (silicate sheet) of the oriented biotite minerals. The anisotropy of the rock is basically TI due to both the aligned biotite minerals and cracks, which have a common symmetry axis. We found that other sheet silicate-rich rocks have a similar anisotropy with the biotite schist, in which the TI-type anisotropy is characterized by the slow P- and S-wave velocities along the symmetry axis. This is caused by the preferred orientation of sheet silicate minerals and the extremely slow P- and S-wave velocities along the axis perpendicular to the silicate sheet compared to the directions along the silicate sheet. When rock contains a large percentage of highly oriented sheet silicates, the fast and slow shear waves exchange their polarities at some off-symmetry axis directions, indicating that the qSwave (quasi-S wave) velocity exceeds the SH-wave velocity. The phase velocity distribution of qS wave shows an asymmetry with respect to the angle from the symmetry axis, which is characterized by a bulge in this distribution located near the symmetric axis. This is inherent to most sheet silicate minerals. When crack density of aligned cracks increases, the P-wave velocity along the symmetry axis decreases considerably. The qS-wave phase velocity in the off-axis directions also decreases, in accordance with the decrease of the P velocity along the symmetry axis. The asymmetry of the qS-wave phase velocity distribution increases as the P-wave velocity decreases along the symmetry axis. This relationship can be well understood by means of Berryman’s extended Thomsen approach. Key words: Microstructures; Seismic anisotropy; Wave propagation; Acoustic properties. 1 I N T RO D U C T I O N Most studies about rock anisotropy have been concerned with Pand S-wave velocities in axial directions, which are defined in association with the rock microstructures (Kern & Wenk 1990; Barruol & Mainprice 1993; Kern 1993; Burlini & Kunze 2000), or the Pwave velocities over whole solid angles (Pros et al. 1998). Recently, studies are more concerned with the P- and S-wave velocities in off-axis directions [Generally, they are not the true P and S waves but the quasi-P (qP) and the quasi-S (qS) waves] (Jones & Wang 1981; Hornby et al. 1994; Johnston & Christensen 1994; Hornby 1998; Jakobsen & Johansen 2000; Takanashi et al. 2001; Wang 2002; Sarout & Guéguen 2008). Some of the sedimentary rocks and metamorphic rocks in the shallow crust generally have planar microstructures. The planar C 2010 The Authors C 2010 RAS Journal compilation structures often include oriented phyllosilicate (sheet silicate) minerals with their silicate sheets mostly aligned parallel to the structural plane. The previous studies indicate that most of the anisotropic rocks in the shallow crust have slower P-wave velocities in the direction normal to the structural plane than in the direction parallel to the plane. This suggests that elastic anisotropy of these rocks is transverse isotropy (TI) or approximated as TI with the symmetry axis perpendicular to the plane. At the same time those rocks have fissility that develops parallel to the structural plane. The fissility suggests that rocks contain thin cracks of which planes are aligned parallel to the isotropic plane of TI symmetry. Aligned cracks bring additional TI-type elastic anisotropy (Anderson et al. 1974; Hudson 1981; Nishizawa 1982; Douma 1988; Nishizawa & Yoshino 2001) over the TI anisotropy produced by aligned sheet silicates. In many cases the symmetry axes of both TI are in the same direction. The 375 GJI Seismology Accepted 2010 April 6. Received 2010 April 6; in original form 2008 September 23 376 O. Nishizawa and K. Kanagawa elastic anisotropy in sheet silicate-rich rocks is thus controlled by a synergistic effect from the preferred orientations of sheet silicates and cracks. One of the interesting characteristics in mica-rich rocks or in some shales are that the slow shear waves along the isotropic plane increase at some area in the off-axis directions. The slow shear wave velocity becomes close to or even exceeds another shear wave velocity in the off-axis directions near the symmetry axis (Hornby 1998; Nishizawa & Yoshino 2001; Takanashi et al. 2001; Kendall et al. 2007; Chesnokov et al. 2009). Crossing over of the two shear wave velocity is called a shear wave singularity (Crampin & Yedlin 1981), and it causes problems for seismic ray tracing. Seismologists pay keen attention to the velocity singularity (Vavryčuk 2003). Nishizawa & Yoshino (2001) suggested that the aligned cracks strongly reduce the P-wave velocity along the symmetry axis and that the P-wave velocity reduction affects the off-axis shear wave velocity distribution. It is important to study how cracks affect seismic velocities in sheet silicate-rich rocks, in terms of numerical modelling and experiments. Here we investigate velocity anisotropy of a cracked biotite-rich rock that has an orthotropic anisotropy altered from the TI-type anisotropy produced by crystal-preferred orientation (CPO) of biotite. Since cracks in rock close under confining pressure, we can evaluate the crack effect on seismic velocity by measuring velocities under confining pressure. We interpret the velocity data measured under confining pressure, for axial and some of the off-axis directions, by comparing model calculations based on an exact formulation originally shown by Mura (1982). Then we try to draw general features of anisotropy in cracked sheet silicate-rich rocks based on a new formulation derived by Berryman (2008). 2 NUMERICAL MODELLING O F A C R AC K E D B I O T I T E - R I C H RO C K 2.1 Inclusion model To estimate the overall elastic moduli of a cracked medium, we use the Eshelby’s (Eshelby 1957) composite medium model consisting of a matrix and inclusions. Let the matrix material have a homogeneous strain, eiAj , or a homogeneous stress field, σiAj . When an inclusion appears inside the matrix, the composite medium is subjected to an internal stress state. Eshelby (1957) presented an idea, which replaces an inclusion with the matrix material and gives a fictitious stress-free strain inside the inclusion, eiTj (referred to as eigenstrain in Lin & Mura 1973). Elastic energy change due to inclusions is expressed by eiTj . There are, however, two extreme cases associated with the far field conditions: constant surface displacements and constant external forces (Eshelby 1957; Yamamoto et al. 1981; Nishizawa 1982). Since we assume a virtual homogeneous medium, the two cases correspond to equivalent homogeneous strain and equivalent homogeneous stress, respectively. The change of the elastic energy corresponding to the equivalent strain eiAj or stress σiAj is given by E = ∓(1/2)σiAj eiTj φ (1) (Eshelby 1957), where the sign ∓ corresponds to the two extreme cases at the far field: constant external forces and constant surface displacements, respectively; φ is the volume fraction of inclusion; σiAj is related to homogeneous strain eiAj through the elastic constants of the matrix (ci0jkl ) by σiAj = ci0jkl eklA . Let ci jkl be the elastic constants of the inclusion. The relationship between the stress-free strain eiTj and strain eiAj in the homogeneous matrix is expressed by Eshelby’s tensor, Si jkl T ci0jkl eklT = ci jkl eklA + ci jkl Sklmn emn , (2) where ci jkl is the difference between the elastic constants of the matrix and inclusion ci jkl = ci0jkl − ci jkl . (3) Eq. (2) shows that the fictitious stress-free strain eiTj can be obtained using eiAj , Eshelby’s tensor Si jkl , and the elastic constants of matrix ci0jkl , and inclusion ci jkl . When the matrix is TI, Eshelby’s tensor Si jkl can be calculated by the integrations given by Lin & Mura (1973), which have been employed by several authors (Nishizawa 1982; Douma 1988; Hornby et al. 1994). The same problem was recently studied by Markov et al. (2005). The Eshelby’s tensor in eq. (2) can be calculated, provided that the strain (or stress) inside the inclusion is homogeneous (Nishizawa & Yoshino 2001). The composite of matrix and inclusions is regarded as an equivalent homogeneous medium and its elastic energy is given by the following equations, corresponding to the conditions in eq. (1): 1 1 1 ∗ −1 A A −1 c σi j σkl = ci0jkl σiAj σklA + σiAj eiTj φ , 2 i jkl 2 2 (4) 1 1 1 ∗ A A c e e = ci0jkl eiAj eklA − σiAj eiTj φ, (5) 2 i jkl i j kl 2 2 where ci∗jkl denotes the elastic constants of the equivalent homo−1 geneous medium, and ci∗jkl −1 and ci0jkl denote the inverse of the matrices ci∗jkl and ci0jkl , respectively. The energy given by eqs (4) or (5) is valid only when φ is very small and the inclusions are non-interacting. To obtain the energy change up to the large φ, we start with the initial inclusion-free matrix and insert inclusions with a small volume fraction φ. As illustrated in Fig. 1, we repeat the process by regarding the composite of the previous step as a new equivalent homogeneous matrix material for the next step. This is called the differential equivalent medium (DEM) method or the numerical self consistent (NSC) approach (Yamamoto et al. 1981; Le Ravalec & Guéguen 1996a,b). The equations are given by 1 ∗(n−1) −1 A(n−1) A(n−1) 1 ∗(n) −1 A(n−1) A(n−1) c σi j σi j = ci jkl σi j σkl 2 i jkl 2 1 A(n−1) T(n−1) + σi j ei j φ, 2 1 ∗(n−1) A(n−1) A(n−1) 1 ∗(n) A(n−1) A(n−1) c e ei j = ci jkl ei j ekl 2 i jkl i j 2 1 A(n−1) T(n−1) − σi j ei j φ, 2 ∗(n) (6) (7) ∗(n) −1 where ci jkl and ci jkl are the elastic constants and their inverse (compliances) obtained at the nth step of calculation, respectively. The superscripts (n − 1) of ci∗jkl , eiAj , σiAj and eiTj denote the values of T(n−1) the (n − 1)th step. ei j is obtained from Eshelby’s tensor for the effective homogeneous medium of TI anisotropy at the (n−1)th step, ∗(n−1) which has the elastic constants ci jkl . Assuming that the composite medium is a new homogeneous matrix material having new elastic constants, we proceed to the next step of the calculation. If both the inclusion and matrix have TI symmetry and the inclusions are aligned to keep their symmetry axes parallel to the symmetry axis of the matrix, the new nominally homogeneous matrix will also be C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation Velocity anisotropy of phyllosilicate-rich rocks 377 Table 1. Elastic constants of biotite (GPa). Biotite C 11 C 33 C 44 C 66 C 13 186.0 54.0 5.8 76.8 11.6 Table 2. Elastic constants of biotite-rich rock model (30 per cent biotite, GPa). Crack-free rock Figure 1. Inclusion models: (a) inclusion model of biotite-rich rock and (b) crack model of biotite-rich rock. (a) A small amount of inclusions with hexagonal symmetry (biotite) are embedded into an isotropic matrix material. → The overall elastic properties of composite medium becomes TI and it is regarded as a new homogeneous TI material. Inclusions are shown by dotted line. → Again, more inclusions are embedded into the TI matrix. → The process is repeated until the biotite volume fraction becomes equal to the expected value. (b) The biotite-rich rock is assumed as a composite material having aligned biotite and isotropic matrix. → The rock is assumed to be an equivalent TI homogeneous material. Aligned fluid-filled cracks are embedded in the matrix. This keeps TI symmetry. → The cracked rock is assumed as a new TI homogeneous material, and more aligned cracks are embedded in to the matrix. ↔ The material is equivalent to the cracked biotite-rich rock. The process of calculation is same as the biotite inclusion model in (a) except that the initial matrix is TI, not isotropic. TI. The above two extreme cases overestimate or underestimate the value of the real energy change, and the elastic constants calculated from eqs (6) and (7) give the lower and upper bounds, respectively (Hill 1952). We use an average of the elastic constants obtained from eqs (6) and (7) for this study as suggested by Yamamoto et al. (1981). 2.2 Calculating elastic constants of biotite-rich rock We consider a biotite-rich rock consisting of an isotropic matrix and 30 per cent volume fraction of biotite. We assume the elastic constants of the isotropic matrix as λ = μ = 35 GPa, which are plausible values of the matrix material forming an aggregate C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation C 11 C 33 C 44 C 66 C 13 126.6 81.9 15.8 47.0 24.4 of quartz and feldspar group minerals that are common in ordinary biotite-rich rocks. Biotite belongs to the monoclinic space group, but elastically it can be regarded as hexagonal for which the symmetry axis is perpendicular to the silicate sheet (Vaughan & Guggenheim 1986). We used the elastic constants given in Table 1 (Aleksandrov & Ryzhova 1961). Hereafter, the reduced Voigt notation Ci j (i, j = 1–6) is used for denoting the measured and calculated elastic constants, instead of ci jkl in the previous theoretical formulations. In TI media, there are five independent values out of 12 non-zero matrix elements, C11 , C33 , C44 , C66 and C 13 , with the relationship C12 = C11 − 2C66 . When the symmetry axes of hexagonal crystals align in the same direction, we can calculate the overall elastic properties by employing the method described in the previous section. Fig. 1(a) illustrates the process for calculating elastic constants of the rock having aligned biotite minerals. Controlled by the cleavage plane parallel to the silicate sheet, biotite minerals usually appear as thin crystals, and mineral shapes are approximated as oblate spheroids described by the major axis, a, and the minor axis, c. Since the present method considers static deformation, the size of inclusion does not affect overall elastic properties. Only the aspect ratio α (= c/a) of inclusions affects the effective elastic properties. Nishizawa & Yoshino (2001) pointed out that the crystal shape of biotite affects the overall S-wave anisotropy of biotite-rich rocks. However the crystal-shape effect on the overall S-wave velocity is significant only when α of biotite crystal is 1/20 ≤ α ≤ 1/10. Taking the biotite shapes in the typical biotite-rich rocks into account (Takanashi et al. 2001), we assume the aspect ratio of biotite crystal 1/20. Less than this value, they do not severely affect the velocity anisotropy in rocks. The values in Table 2 are obtained for a rock containing 30 per cent biotite by using the method described by Nishizawa & Yoshino (2001). These are used as the elastic constants of the initial matrix of aligned crack models. We then also assume that the initial crack-free biotite-rich rock is a homogeneous TI material and inserts aligned microcracks into the matrix. Fig. 1(b) illustrates the process for inserting aligned cracks into the TI matrix. In each step of the calculation, we obtain two elastic constants, corresponding to the eqs (6) and (7). The difference between the two elastic constants becomes large as the difference of elastic constants between inclusion and matrix becomes large. The difference is also controlled by the crack aspect ratio and the step volume fraction φ. However, crack density is the unique parameter that includes the effects of both crack shape and crack volume fraction. Crack density is proportional to the value φ/α. We selected the step volume fraction φ so as to keep small values for φ/α. We calculate averages of the two elastic constants obtained from eqs (6) and (7). 378 O. Nishizawa and K. Kanagawa 2.3 Seismic velocities in TI media Seismic velocities of an equivalent homogeneous medium with elas∗(n) tic constants ci jkl are given by the wave equation ρ 2 ∂ 2ui ∗(n) ∂ u k = ci jkl , 2 ∂t ∂ x j ∂ xl (8) where ρ is the rock density and u i is the ith component of the wave particle motion. The equation can be solved for the plane wave propagating with the wavenumber (q1 , q2 , q3 ), where qi denotes the xi component of the wave vector q. The relationship between ω and q is given by the Christoffel equation (Landau & Lifshitz 1959) 2 ∗(n) (9) ρω δik − ci jkl q j ql = 0. Phase velocities are obtained by ω/|q| for the three independent solutions. The corresponding displacement vector gives the polarization direction of seismic waves (Carcione 2007). In what follows, we deal with only TI media. Among the three phase velocities calculated from eq. (9), the polarization direction of the fastest wave generally deviates from the propagation direction except when the wave propagates along the symmetry axis or on the isotropic plane. This wave is referred to as the quasi-P wave, qP, to distinguish from the true P wave of which polarization is directed to the propagation direction. The polarization of one of the S waves generally deviates from the true S-wave polarization that is in the plane perpendicular to the propagation direction. The wave is referred to as the quasi-S wave, qS. Another S wave is the true S wave with polarization parallel to the isotropic plane and perpendicular to the propagation direction. When waves propagate along the symmetry axis or on the isotropic plane, the two shear waves are the true S waves polarized perpendicular to the propagation direction. Along the symmetry axis their velocities are equal but they are different in the isotropic plane. The phase velocities of qS and SH waves are generally different in diagonal directions. Each phase velocity depends only on the angle measured from the symmetry axis, θ. In TI media there are only five independent elastic constants, and eq. (9) can be described in explicit form by using Ci j . The phase velocities of qP and qS waves, Vq P and Vq S , are given by the following equations (Mavko et al. 1998): 1/2 C11 sin2 θ + C33 cos2 θ + C44 + M(θ ) , (10) Vq P = (2ρ)1/2 Vq S = 1/2 C11 sin2 θ + C33 cos2 θ + C44 − M(θ ) (2ρ)1/2 , (11) formulae that approximate the angular dependence of phase velocities, Vq P (θ), Vq S (θ) and VS H (θ) for weak anisotropy. He derived anisotropic parameters ε, δ and γ , (and a non-independent parameter η) referred to as Thomsen’s anisotropic parameters. ε and γ are obtained from the axial P- and S-wave velocities V P (0), V P (π/2), VS (0) and VS (π/2), δ is obtained from one of the off-axis qP or qS velocity. In Thomsen’s approximation, phase velocities are expressed as Vq P (θ) = V P (0)(1 + δ sin2 θ cos2 θ + ε sin4 θ), (14) Vq S (θ) = VS (0)(1 + η sin2 θ cos2 θ), (15) VS H (θ) = VS (0)(1 + γ sin2 θ), (16) δ becomes precise when using the qP velocity near the diagonal direction (θ = π/4) because sin2 θ cos2 θ takes the maximum value in eq. (14). There are relationships between anisotropic parameters, elastic constants and the velocity values along θ = 0, π/2 and π/4. ε= VP (π/2) − V P (0) C11 − C33 , ≈ 2C33 V P (0) (17) γ = C66 − C44 VS H (π/2) − VS H (0) , ≈ 2C44 VS H (0) (18) (C13 + C44 )2 − (C33 − C44 )2 2C33 (C33 − C44 ) VP (π/2) Vq P (π/4) −1 − −1 , ≈ 4 V P (0) V P (0) δ= η in eq. (15) is given by VP (0) η= VS (0) + 4(C13 + C44 )2 sin2 θ cos2 θ . (12) The phase velocity of SH wave, VS H is given by 1/2 . VS H = (C66 sin2 θ + C44 cos2 θ )/ρ (13) 2 (ε − δ), (20) where η controls the magnitude of Vq S in off-axis directions. Vq S takes an extreme value at θ = π/4 (eq. 15). Recently, Berryman (2008) derived extended Thomsen formulae applicable for strong anisotropy. The major extension is that the new approximation removes an implicit restriction of the diagonal symmetry of the angular distributions of Vq S , which is caused by the extreme values at θ = π/4. The exact velocities of qP and qS are given by eqs (10) and (11). We denote the angle of the extreme value of M(θ) by θm . tan2 θm = where M(θ ) is 2 M(θ ) = (C11 − C44 ) sin2 θ − (C33 − C44 ) cos2 θ (19) C33 − C44 . C11 − C44 (21) Vq P and Vq S are approximated by using θm as 2 sin2 θm sin2 θ cos2 θ Vq P (θ) = V P (0) 1 + ε sin2 θ − ( − δ) 1 − cos 2θm cos 2θ (22) and 2 sin2 θm sin2 θ cos2 θ . = VS (0)(1 + η) 1 − cos 2θm cos 2θ 2.4 Thomsen’s approximation and extended Thomsen formulae for TI media Vq S For interpreting field observations, it will be useful to describe phase velocities as simple functions of θ . Thomsen (1986) derived We now examine the applicability of the Thomsen’s and the extended Thomsen’s approximation for biotite-rich rocks. C (23) 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation Velocity anisotropy of phyllosilicate-rich rocks 2.5 Phase velocities of the biotite-rich rock model For calculating the phase velocities in the biotite-rich rock, we assume the rock density 2.75 × 103 kg m−3 . The value is equal to the density of biotite schist measured in the present experiment (ST-3 of Takanashi et al. 2001). The bulk modulus of the crackfilling fluid controls velocity anisotropy. Fluid bulk modulus differs by orders of magnitude between liquid and gas. We assumed the fluid bulk modulus K f = 0.01 and 1 GPa for gas and liquid, respectively, corresponding to typical cases of gas-filled and liquid-filled cracks in underground conditions. To describe the relationship between seismic velocities and cracks, crack density is conventionally used (O’Connell & Budiansky 1974). Crack density corresponds to a population of cracks with specific aspect ratio, and shows the combined effect of crack volume and crack shape on seismic velocities. The crack density is given by each crack porosity φi (the volume fraction of crack) and its crack aspect ratio αi as = N γi φi /αi , (24) i=1 where γi is the geometrical shape factor for each open crack, and N is the total number of cracks. For oblate spheroidal cracks γi = 3/4π. When crack closes under external forces, φi becomes zero. Crack density decreases as cracks close, then velocities increase. Velocity change can be described as a function of crack density (O’Connell & Budiansky 1974; Soga et al. 1978). It is enough to calculate wave velocities for a representative crack shape as a function of crack density because crack density is the unique parameter that controls the velocity in a cracked medium. In the phase-velocity calculations of Fig. 2, we adopted a single value of 0.0025 (1/400) for crack aspect ratio and show results for different crack densities and different fluid bulk moduli. Figs 2(a)–(d) show angular distributions of phase velocity plotted on one of the quadrants formed by the symmetry axis and the isotropic plane of TI symmetry (hereafter simply referred to as the quadrant). The fluid bulk modulus K f is assumed as 0.01 GPa. The thick solid curves correspond to the crack densities 0.02 and 0.06, respectively. The thin solid curves show velocity distributions in the crack-free state to show the velocity change due to cracks. Figs 2(a) and (b) show Vq P , and Figs 2(c) and (d) show both Vq S and VS H . The rightand left-hand panels correspond to crack densities 0.02 and 0.06, respectively. Figs 2(e)–(h) show phase velocities for another typical value of fluid bulk modulus K f = 1 GPa. The panels in the upper and lower rows and the left- and right-hand columns correspond to the same conditions as in Figs 2(a)–(d). The angular dependence of phase velocity for the model matrix (thin solid curves) are characterized by inflected curves in the Vq P distribution, and the diagonal bulge in the Vq S distribution. The diagonal bulge shrinks as the crack density increases (cf. Figs 2c and d), or the fluid bulk modulus becomes smaller (cf. Figs 2d and h). When fluid bulk modulus is small, cracks have more effects on the P-wave velocity anisotropy, that is the velocity difference between the symmetry axis and the isotropic plane, compared to the S-wave anisotropy. Increase of crack density reduces the off-axis qS-wave velocities. The reduction is more pronounced when compared to the reduction of the qS-wave velocity along the symmetry axis and the isotropic plane for the case K f = 0.01 GPa. This moves the crossover point of the Vq S and VS H distribution curves closer to the symmetric axis (typically in Fig. 2d). The difference between VS H and Vq S in the same direction is referred to as the shear wave C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation 379 splitting. Since the velocity of the SH wave is insensitive to crack density, the shear wave splitting is mainly affected by the angular distribution of Vq S . The bulge in the Vq S distribution shrinks with increasing crack density. The shrinkage affects the shear wave splitting. Thus shear wave splitting is sensitive to crack density at some area of the off-axis directions. With increasing crack density, a reversal of the fast and slow shear wave polarization may occur in a certain range of θ when fluid bulk modulus (K f ) is small. Thomsen’s and the extended Thomsen’s approximations are also plotted for all the cases. There are considerable differences in Thomsen’s parameters defined by elastic constants and those calculated from velocities by using eqs (17)–(19). We use axial velocities for calculating ε and γ . δ shows considerable difference between the two values: from elastic constants and from Vq P (π/4). We use δM that is determined from the elastic constants (moduli). In Figs 2(a)–(h), Thomsen’s approximation and the extended Thomsen’s approximation are plotted as the dotted curves and the break curves, respectively. For Vq P , the extended Thomsen’s approximation shows better fits to the velocity distribution curves than Thomsen’s approximation, even for higher crack density and small fluid bulk modulus. For VS H , Thomsen’s approximation and the extended Thomsen’s approximation are the same. The approximation fits well to the SH velocity distribution curves, and the velocities are unaffected by the crack density and the fluid bulk modulus. Thomsen’s approximation by using δM overestimates the actual q S-wave velocity. The extended Thomsen’s approximation always underestimates the maximum values of Vq S . 3 L A B O R AT O RY V E L O C I T Y MEASUREMENTS OF BIOTITE SCHIST 3.1 Velocities under 150 MPa confining pressure To examine crack effects on seismic velocities in the biotite-rich rock model presented above, we made laboratory velocity measurements of Hidaka biotite schist under confining pressures up to 150 MPa. The rock was sampled from Hidaka metamorphic belt, Hokkaido, Japan (the sample ST-3 in Takanashi et al. 2001). The rock shows well-developed foliation plane and lineations in the foliation plane. The cleavage planes of biotite crystals are mostly aligned parallel to the foliation plane (fig. 2 and fig. 4 of Takanashi et al. 2001). Rock density, CPO data and seismic velocity anisotropy at the confining pressure 150 MPa are presented in Takanashi et al. (2001). Anisotropy of the rock is described with respect to the foliation plane. We select x3 -axis perpendicular to the foliation plane, and x1 - and x2 -axes parallel and perpendicular to the lineation, respectively. P- and two polarized S-wave velocities were measured along the three axes and in diagonal directions in the x1 x3 - and x2 x3 -planes. The polarization directions of S waves are parallel and perpendicular to propagation directions. Figs 3(a-1) and (a-2) show the microstructure and the biotite c-axis distribution of the rock. Velocities were measured by the pulse transmission method. We used piezo-electric transducers of compression and shear modes with 5-mm-diameter disc and 6 mm × 3 mm rectangular plate, for measuring P- (qP) and S-wave (qS) velocities, respectively. The travel distances of elastic waves are 50–80 mm. In this distance range, elastic waves are approximated as plane waves and measured velocities are almost equal to the phase velocities (Dellinger & Vernik 1994). Velocities at the confining pressure 150 MPa are shown in Fig. 3(b). The two panels in Fig. 3(c) show angular velocity distributions in the two quadrant sections, x1 x3 and x2 x3 . The distributions 380 O. Nishizawa and K. Kanagawa Figure 2. Angular phase-velocity distributions for qP, qS and SH waves. Values are plotted on one of the quadrants formed by the symmetry axis (vertical direction) and isotropic plane (horizontal direction). The rock contains 30 per cent volume fraction of biotite. Crack aspect ratio is 0.0025 (1/400) for all calculations. The fluid bulk modulus for (a)–(d) is 0.01 GPa, and for (e)–(h), 1 GPa. The left and right columns correspond to low (0.02) and high (0.06) crack densities (CD). The phase velocities for the aligned crack cases are shown by thick curves. Crack-free states and approximations by Thomsen and the extended Thomsen are shown by thin, dotted and dashed curves, respectively. are calculated from the velocity data in each plane at 150 MPa. Those velocity distributions will be discussed in the succeeding subsections. Axial velocities under 150 MPa confining pressures and the biotite CPO data indicate that the velocity anisotropy comes from preferred orientation of biotite minerals. Biotite minerals basically align their cleavage plane parallel to the foliation plane. This produces TI-type seismic velocity anisotropy. In the x2 x3 section, the cleavage planes of biotite show slight misalignments from the foliation plane (Takanashi et al. 2001), forming microscale folding accompanied by crenulation cleavage of the rock. The misalignments of biotite cleavage planes and crack planes from the foliation C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation Velocity anisotropy of phyllosilicate-rich rocks 381 Figure 3. (a-1) A sketch of the microstructure. (a-2) The biotite c-axis distribution of biotite schist. The contour is an arbitrary scale modified from Takanashi et al. (2001). (b) 150-MPa confining-pressure velocities along the axial and the diagonal directions in the x2 x3 - and x1 x3 -planes. The axes are selected in reference to the foliation plane and the direction of lineation. x3 -axis is perpendicular to the foliation plane, and x1 - and x2 -axes are parallel and perpendicular to the lineation, respectively. Two orthogonally polarized shear waves show different velocities in the diagonal directions of the x1 x3 -plane, whereas the two shear wave velocities are equal in the diagonal directions of the x2 x3 -plane. Approximation is made by assuming the x3 -axis as the symmetry axis and the velocities in the x1 x3 -planes represents the velocities in TI anisotropy. plane alter the velocity anisotropy from TI to orthotropic symmetry. To see the effect of cracks on seismic velocity, we approximate the rock anisotropy by a combination of the two TI anisotropies seen in the x1 x3 -plane and the x2 x3 -plane. This kind of simplification was employed in the field seismic data analysis for orthotropic anisotropy (Tsvankin 1997). The P-wave velocity in the x2 direction is slower than that in the x1 direction. This is due to the misalignments of biotite sheets and cracks from the foliation plane. In the x1 x3 -plane, the fast shear wave is polarized in the plane including the symmetry axis in the direction θ = π/4. The detected fast shear wave is considered as the qS wave. A crossover between the VS H and Vq S distribution appears in the x1 x3 -plane because VS H exceeds Vq S in the x1 -axis. The Vq S distribution forms a bulge in the x1 x3 -plane (Fig. 14 of Takanashi et al. 2001), similar to the phase velocity distribution of Fig. 2. On the other hand, the two polarized shear wave velocities along the direction θ = π/4 are almost equal. This suggests that the bulge of the qS phase-velocity distribution in the x2 x3 -plane is very small. The velocities in the x1 x3 -plane under 150 MPa show striking similarities to the velocities in the biotite-rich rock model presented in the previous section. 3.2 Velocity increase from atmospheric pressure to 150 MPa Fig. 4(a) shows the pressure dependence of the axial velocities and the velocities along the θ = π/4 direction in the x1 x3 - and x2 x3 -planes. The P-wave velocity perpendicular to the foliation plane, V P3 , increases by about 1.5 km s−1 (3.78–5.27 km s−1 ) as the confining pressure increases from atmospheric pressure to 150 MPa, whereas the P-wave velocity in the x1 (lineation) direction, V P1 , increases only by 0.5 km s−1 (5.74–6.27 km s−1 ) in the same pressure range. The velocity in the x2 direction increases by C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation 0.7 km s−1 in this pressure range. Both in the x1 x3 - and x2 x3 -planes, the qP-wave velocities along θ = π/4, Vq P (π/4), increases with increasing the confining pressure. Vq P (π/4) becomes close to V P3 at 150 MPa. The axial S-wave velocities are denoted by VSi j where i and j correspond to propagation and polarization directions, xi - and x j axis directions, respectively. S-wave velocities show only a small increase compared to P-wave velocities except for Vq S (π/4) in the x1 x3 -plane. The differences between VS12 and VS13 , and VS21 and VS23 show the intensity of shear wave splitting in the x1 and x2 axial directions, respectively. Those differences do not show remarkable changes with increasing confining pressure. However, in the direction θ = π/4 in the x1 x3 -plane, the difference between Vq S (π/4) and VS H (π/4), does change significantly with increasing confining pressure. In contrast, in the direction θ = π/4 in the x2 x3 plane, Vq S (π/4) and VS H (π/4) are almost equal. 4 DISCUSSION 4.1 Closure of cracks under confining pressure We interpret the pressure dependence of velocity in Hidaka biotite schist by assuming aligned cracks in a TI matrix. Because more cracks close as the confining pressure increases, the increase of seismic velocity under confining pressure results from the decrease of crack density. In isotropic materials, the following equation gives the crack-closing pressure ( pc ) for the cracks having the aspect ratio α (Walsh 1965; Simmons et al. 1974) pc = π Eα , 4(1 − ν 2 ) (25) where E and ν are the Young’s modulus and the Poisson’s ratio of the matrix material, respectively. Since Poisson’s ratio takes the value 382 O. Nishizawa and K. Kanagawa Figure 4. Seismic velocities in the Hidaka biotite schist as functions of confining pressure. (a) Velocities in the x1 x3 -plane. (b) Velocities in the x2 x3 -plane. between 0.2 and 0.3 for the most of crystalline rocks, eq. (25) gives an approximation that pc is nearly equal to Eα. The relationship will hold for the present anisotropic rock matrix, giving only a small change for pc since the effective E is almost the same for the present aligned crack case. If we tentatively assume that E is 60 GPa (a plausible value for crystalline rocks), the cracks with aspect ratio less than about 1/400 will close at 150 MPa. Cracks with the larger aspect ratios still exist above 150 MPa. Technically, the rock is not a crack-free material under 150 MPa, but we regard the rock at 150 MPa as an equivalent homogeneous material even if the rock actually includes cracks and voids. Then we can calculate elastic constants of the rock below this pressure by inserting cracks into this nominally homogeneous matrix. Table 3. Ci j for Hidaka biotite schist (GPa). C 11 C 33 C 44 C 66 C 13 x1 x3 -plane 108.1 76.35 22.33 36.22 19.35 x2 x3 -plane 95.70 76.35 22.33 33.68 29.35 of Hidaka biotite schist. We adopt the measured density, 2.75 × 103 kg m−3 (ST-3 in Takanashi et al. 2001). C 13 was estimated from the measured Vq P (π/4) and other four elastic constants. Table 3 shows the estimated Ci j of Hidaka biotite schist in the x1 x3 - and x2 x3 -planes at 150 MPa where the rock is nominally crack-free. We calculate velocities in each direction as functions of crack density by using those elastic constants. 4.3 Phase velocity as a function of crack density 4.2 TI velocity anisotropy in the x1 x3 - and x2 x3 -planes Schist often tends to split along their foliation plane. This indicates that microcracks mostly distribute along the foliation plane with crack planes almost parallel to the foliation plane. To interpret elastic wave velocities under confining pressure, it will be reasonable to apply the aligned crack model as the first approximation. The experimental velocity data show different velocity increases between x1 x3 - and x2 x3 -planes. The differences are induced by microfolding seen in the x2 x3 -plane, which induce an orthotropy in an original TI-type anisotropy due to CPO of biotite. The velocities under 150 MPa confining pressure indicate two TI velocity anisotropies that are similar to the biotite-rich rock. Although the schist shows an orthotropic velocity anisotropy, we consider two different TI anisotropies in the x1 x3 - and x2 x3 -planes. The difference is induced by the differences in biotite-preferred orientations. We try to interpret velocities under confining pressure in the two planes independently. We obtain the following elastic constants for C11 , C33 , C44 and C 66 for the x1 x3 - and x2 x3 -planes We calculate velocities by using the same method as for the previous biotite-rich rock model by considering the schist as a TI matrix. Calculated velocities in each measured direction are shown in Figs 5(a)–(c) as functions of crack density. Figs 5(a) and (b) show the velocities in the assumed TI media for the x1 x3 - and x2 x3 -planes, respectively. Fig. 5(c) shows the same relationship for the biotite-rich rock model I in Table 2. The horizontal axis scales are reversed in accordance with the inverse proportionality between crack-closing pressure and the crack density as indicated by eqs (24) and (25). The observed velocities increase with increasing confining pressure. Because velocities under confining pressure do not show the quantitative relationship between velocity and crack density, the two horizontal axes in Figs 4 and 5 cannot be compared directly. The total crack density is given by the summation of the crack densities of different aspect ratio cracks as given by eq. (24). Under low confining pressures, cracks with smaller aspect ratios are preferentially closed. Crack density decreases because the smaller α C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation Velocity anisotropy of phyllosilicate-rich rocks 383 Figure 5. Seismic velocities in the three models as functions of pressure: (a) aligned crack model estimated from the velocities at 150 MPa in the x1 x3 -plane; (b) aligned crack model estimated from the velocities at 150 MPa in the x2 x3 -plane; (c) biotite-rich rock model I. values in the denominator in eq. (24) are dropped out. To obtain the quantitative relationships, we need to know the distribution of the crack aspect-ratio population (Simmons et al. 1974). We do not have the data about crack aspect-ratio distribution. But if we assume that the crack porosity is of the same order for different aspect ratios, velocities are more affected in the low confining pressure range because of closure of low aspect-ratio cracks. If the crack density rapidly decreases with increasing confining pressure at the low confining pressures, the rapid velocity increases in the low confining pressure range in Fig. 4 will correspond to rapid decrease of crack density due to closing of thin cracks. Because P-wave velocity is more sensitive to the crack density, increase of P- (or qP-) wave velocity is more pronounced than S-wave velocities. Even considering closure of the aligned thin cracks, there still remain discrepancies in the velocity increase at the low confining pressure range. V P1 and V P2 increase more than the expected values in Fig. 5. We have to consider other crack orientations. Nur & Simmons (1969) showed that closure of randomly oriented thin cracks brings considerable velocity increase in crystalline rocks at the low confining pressure range. Increase in V P1 and V P2 at the low pressure range suggests that a certain population of randomly oriented thin cracks could be superposed to the population of aligned cracks. In addition to the aligned and random crack orientations, another population of crack orientations will be associated with the microfolding, where crack planes are slightly slanted from the x1 x2 plane. Those slanted cracks would affect the axial velocities in a different manner as shown in Fig. 6. When P and SH waves are propagating parallel to crack planes, their velocities are little affected by cracks. For the slanted cracks in Fig. 6, V P1 propagates parallel to the crack plane. In this case, the P wave does not pass through the crack plane. On the other hand, V P2 is affected strongly by the slanted cracks because it passes through the crack planes. V P3 propagates normal to the crack planes of aligned cracks. Slanted cracks slightly change the incident angle of V P3 , but this will have little effect on velocity. Slanted biotite sheets will have slight effect on V P3 because V P3 increases C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation in the off-axial directions around the symmetry axis as shown in Fig. 2. VS12 and VS21 are slightly affected by the slanted cracks. VS12 propagates parallel to the slanted crack plane but its polarization is slightly inclined with respect to the crack plane. On the other hand, the S21 wave propagates through the slanted crack plane with a small angle from the crack plane. The polarization direction of the S12 wave shows slight deviations from the slanted crack plane and the major vibration direction is still parallel to crack plane. S21 wave propagates through the slanted crack plane and is more affected by slanted cracks than S12 wave. VS21 is therefore more affected by cracks than VS12 . 4.4 Angular velocity distributions Ignoring the small differences caused by the cracks with normals directed off the x3 -axis, the relationship between calculated velocities and crack density will correlate to the pressure dependence of the measured velocity. For example, in Fig. 4(a) V P3 becomes close to Vq P (π/4) as confining pressure increases. This is similar to the large increase of V P3 and rather small increase in Vq P (π/4) as shown in Fig. 2(b) with decreasing crack density. A crossover of the two shear wave velocities in the diagonal direction, VS H (π/4) and Vq S (π/4), appears in Fig. 4(a). Vq S (π/4) becomes larger than VS H (π/4) as crack density decreases. This is similar to the change in the qSand SH-wave velocity in off-axis directions with decreasing crack density, as shown in the biotite-rich rock model in Fig. 2(d). We conclude that the experimental results can be basically explained by crack models having aligned thin cracks in the biotite-rich rock matrix having a TI-type anisotropy. We examine velocity anisotropy of cracked biotite-rich rocks in more detail, by plotting phase velocity in all directions. Figs 7(a)–(d) show phase velocity of the biotite schist calculated from the elastic constants in Table 3. As crack density increases, the anisotropy in P wave is more enhanced; the phase velocity surface of qP wave shows a flat-bottomed pan shape, indicating that the qP-wave 384 O. Nishizawa and K. Kanagawa Figure 6. Effect of microfolding on velocities. Microfolding causes tilts of the normals of crack planes and biotite cleavage planes. The tilts affect each velocity in a different way. The most affected velocity is V P2 of which ray path is first parallel to crack planes or biotite cleavage planes but is crossing those planes when the tilts appear. V P2 then becomes smaller than V P1 . velocity is considerably affected by cracks at the range 0 ≤ θ ≤ π/4, but it is little affected at π/4 ≤ θ ≤ π/2. This is the reason why Vq P (π/4) show the smaller velocity increase compared to V P3 with decreasing of crack density as shown in Fig. 4. The experimental results on V P1 , V P3 and Vq P (π/4) support this, because V P3 increases rapidly compared to V P1 and Vq P (π/4) which show almost same velocity increase with respect the confining pressure. VS H (π/4) and Vq S (π/4) come closer to each other as crack density increases, as shown by the thick curves (crack density 0.1) and thin curves (crack-free) in Fig. 7(d). The crossover point of phase velocities is close to the direction π/4 off from the x3 -axis. 4.5 Crack effects on phase velocities in shale Shale is one of the important rocks for reservoir assessment and its anisotropy is a much-discussed subject. Shales generally show TItype anisotropy with strong CPO of clay minerals (Jones & Wang 1981; Hornby et al. 1994) accompanying dominant aligned thin microcracks. Clay minerals have strong velocity anisotropy similar to biotite, with large differences between C 11 and C 33 , and between C 44 and C 66 (Hornby et al. 1994). A number of investigations have been done for interpreting velocity anisotropy in shale (Sayers 2002; Hall et al. 2008; Verdon et al. 2008). Here we apply our crack model to understand velocity anisotropy of shale under confining pressures or under deviatric stresses. First, we have to notice the difference between the present metamorphic rock and shale. The matrices of shales are elastically more compliant than the matrices in metamorphic rocks, and elastic compliances of shale matrices change with the confining pressure (Jones & Wang 1981; Hornby et al. 1994; Hornby 1998; Sayers 2002; Sarout & Guéguen 2008; Voltolini et al. 2009). The velocity increase in shale under confining pressure is not only due to closure of cracks but also due to changes of elastic constants of matrix material by increasing of the contact stiffness between particles under confining pressure. Therefore, it is not easy to extract crack effects only from the velocity data under confining pressures. To interpret the velocity changes in shale under confining pressure, we assume that elastic constants of the matrix is constant and the velocities under confining pressure are only due to closure of cracks. This simple assumption may be useful as the first approximation to understand the wave velocities under confining pressure. We take the two typical cases of velocities in shale based on the models of Hornby et al. (1994): the aligned and the distributed clay-platelet orientations. We assume the bulk modulus of the crackfilled fluid as 0.01 GPa, the case close to gas-filled cracks, which emphasizes crack effects. Figs 8(a-1), (a-2) and Figs 8(b-1), (b-2) show phase velocities for the qP, qS and SH waves. The aligned model shows a remarkable inflection of the Vq P distribution curve and a large bulge of the Vq S distribution and its diagonal asymmetry in the quadrant. On the other hand, the distributed model shows unclear inflection of the Vq P distribution and the small bulge of the qS velocity distribution. In both cases, the axial velocity anisotropy defined by Thomsen’s parameters ε and γ increases with increasing crack density for both P and S waves. With increasing crack density, the off-axis bulge in the Vq S distribution shrinks. In the quadrant velocity distributions, this moves the crossover point between the SH and qS waves (singularity) and the maximum value of the qSwave velocity closer to the symmetry axis. Those affect the off-axis shear wave splitting. The polarity of the fast shear wave reverses in some area, and the magnitudes of shear wave splitting increases or decreases depending on the angle θ. For the distributed orientation model, the inflection in the Vq P distribution and the singularity of shear waves are not clear. The magnitude of shear wave splitting increases with increasing crack density in the propagation directions π/4 < θ ≤ π/2. 4.6 General characteristics in the TI anisotropy of biotite-rich rocks and shales As has been shown here, the strong TI anisotropy in biotite-rich rock is characterized by an increase of Vq S in the off-axis directions, which is manifested as a bulge and a diagonal asymmetry of the angular distribution of Vq S . The magnitude of the off-axis bulge in the Vq S distribution is roughly estimated by Thomsen’s anisotropic parameter η although it does not reproduce the exact distributions. The angle of the maximum Vq S is controlled by θm in eq. (23). Table 4 shows a comparison of anisotropy in several TI-type rocks and crack effects on anisotropy. ε, η and θm control the characteristics of TI anisotropy such as the shape of the Vq P distribution, the magnitude of the bulge in the Vq S distribution and the maximum phase-velocity angle in the velocity distributions of qS wave. The upper segment of the table shows the effect of cracks and the lower segment shows other shales that have weak anisotropy. In the upper segment, we assumed K f = 0.01 GPa and crack density (CD) = 0.1 or 0.06. In all the cases, ε increases with increasing crack density. The two typical cases of totally aligned mineral models, biotite-rich rock model I and Cretaceous shale with aligned clay platelet model (Hornby et al. 1994), have large η values. The large η indicates increase of Vq S in the off-axis directions. With C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation Velocity anisotropy of phyllosilicate-rich rocks 385 Figure 7. Calculated phase velocity distributions in the x1 x3 - plane of biotite schist (Table 3) for crack-free state and cracked state: (a) crack-free case qP wave, (b) cracked case qP wave, (c) crack-free case SH and qS waves and (d) cracked case SH and qS waves. Crack density is 0.1. Cracks are filled with gas (K f = 0.1 GPa) and their normals align parallel to the TI symmetry axis of biotite schist. Figure 8. Calculated phase velocity distributions for clay platelets aligned and distributed orientation models for crack density 0.1 and fluid bulk modulus 0.01 GPa. (a-1), (a-2) aligned model. (b-1), (b-2) distributed model. C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation 386 O. Nishizawa and K. Kanagawa Table 4. Anisotropic parameters and related parameters of biotite-rich rocks and shale rocks. ε − δM V P (0) VS (0) 2 ε δM Biotite-rich rock model I (matrix) Biotite-rich rock model I (CD=0.06)a Hidaka biotite schist Hidaka biotite schist (CD=0.1)a Creta. shale aligned(1) Creta. shale aligned (CD=0.1)a Creta. shale distributed(1) Creta. shale distributed (CD=0.1)a Kaolinite-illite mixture(4) Kaolinite-illite mixture (CD=0.1)a 0.273 0.586 0.208 0.648 0.648 1.364 0.282 0.690 0.766 1.359 −0.254 −0.166 −0.143 0.108 −0.131 0.111 0.063 0.391 0.413 1.204 0.527 0.752 0.351 0.540 0.779 1.253 0.219 0.299 0.346 0.155 5.18 3.91 3.42 2.43 5.63 3.48 3.68 2.58 1.89 1.49 Cotton valley shale(2) Wills point shale(2) Kimmimerridge shale(3) 0.135 0.215 0.344 0.205 0.315 0.177 −0.070 0.100 0.167 2.67 7.47 3.47 Rock θ m (◦ ) ζm 2.732 2.937 1.200 1.313 4.384 4.356 0.807 0.773 0.653 0.231 37.68 31.93 38.43 29.20 31.92 24.47 36.89 29.02 25.85 18.23 0.779 0.785 0.625 0.573 0.735 0.728 0.339 0.299 0.346 0.101 −0.187 −0.747 0.579 36.89 39.27 35.49 −0.156 −0.108 −0.154 η a K = 0.01 GPa, (1) Hornby et al. (1994), (2) Berryman (2008), (3) Hornby (1998), Voltolini et al. (2009), (4) Voltolini et al. (2009) 100 per cent clay model f 50 MPa compression. increasing the crack density, θm decreases, indicating the stronger diagonal asymmetry of the Vq S distribution. With increasing crack density, both Hidaka biotite schist and Cretaceous shale with distributed orientation of clay platelet (Hornby et al. 1994) show the same changes in ε and θm . If we carefully look into the published experimental data, we find that the offaxis qS-wave velocity increases with increasing confining pressure. Sometimes the velocity increase is very small, but it exceeds increases in the SH-wave and the axial S-wave velocities (Jones & Wang 1981; Hornby 1998; Sarout & Guéguen 2008). This is similar to the increases in Vq S in Figs 4(a)–(c). Increase of Vq S under confining pressure can be interpreted as an increase of θm or extension of the bulge of the Vq S distribution with decreasing crack density, as shown in Figs 2, 7 and 8. It should be noted that θm is less than 45◦ in all cases. Because θm is defined by tan2 θm = (C33 − C44 )/(C11 − C44 ), it is less than one because the condition C33 < C11 always holds in anisotropic sheet silicate-rich rocks. Since cracks reduce P-wave velocity along the symmetry axis, C 33 decreases, tan θm decreases with increasing crack density and θm also decreases. By comparing the values of η, θm and ζm in the lower segment with those values of the upper segment, we can roughly estimate the behaviour of the qP- and qS-wave velocity distributions. 4.7 Fitting of the extended Thomsen’s approximation The main point of Berryman’s extended Thomsen method is the following approximation for M(θ ) in eq. (12): M(θ) = (C11 − C44 ) sin2 θ + (C33 − C44 ) cos2 θ [1 − ζ (θ)]1/2 ≈ (C11 − C44 ) sin2 θ + (C33 − C44 ) cos2 θ [1 − ζ (θ)/2] , (26) where ζ (θ ) is given by ζ (θ ) = 4 [(C11 − C44 )(C33 − C44 ) − (C13 + C44 )2 ] sin2 θ cos2 θ . [(C11 − C44 ) sin2 θ + (C33 − C44 ) cos2 θ ]2 (27) The approximation in eq. (26) is valid for small ζ (θ ). ζ (θ ) becomes small when θ is close to zero or π/2. The fitness of the approximation by eq. (26) around the extreme value of ζ (θ ) is controlled by the magnitude of the extreme value, ζm . Values of ζm for biotiterich rock model, Hidaka biotite schist and some shales are listed in Table 4. ζm becomes relatively large for strong preferred orientations of biotite or clay platelets. The difference between the approximation and the exact calculation becomes large in the Vq S distribution when η is large. η and ζm are generally large for a strong anisotropic material like the biotite-rich rock model I and the aligned clay platelet model of the Cretaceous shale. This is the reason why the extended Thomsen’s approximation gives poor fits to the exact velocities in those cases. Hidaka biotite schist shows better fits compared to the above rock models for both matrix and cracked cases because of smaller η and ζm . For most of the distributed orientation models and natural shale rocks, ζm and η are small. In such cases, Berryman’s approximation fits well to the exact calculations. 4.8 Shear wave splitting in the off-axis directions The magnitude of shear wave splitting is measured from the velocity differences between the phase velocities of S H and q S waves. Shear wave splitting has been used for probing subsurface conditions (Crampin 1987; Crampin, et al. 1990). Interpretations based on phase velocities can be applied to teleseismic waves because plane wave assumption is implicitly employed. Here, we tentatively adopt the plane wave assumption and discuss shear wave splitting of phase velocities. Fig. 9 shows the relative differences between VS H and Vq S in per cent, 100 × (VS H − Vq S )/[(VS H + Vq S )/2], as functions of θ for Hidaka biotite schist and Cretaceous shale with distributed orientation of clay platelets. The positive area is the normal splitting zone, VS H > Vq S , whereas the negative area is the reverse, Vq S > VS H , due to the increase of Vq S near the symmetry axis. The boundary of the normal and the reverse zone (zero-crossing) corresponds to the singularity angle. With increasing crack density, the normal zone extends due to the decrease of the off-axis Vq S values. Also θm decreases. This means that the direction of the extreme value of M(θ ) gets closer to the symmetry axis. In Hidaka biotite schist, increasing of crack density causes an exchange of the fast and slow shear wave polarization near θ 55◦ , where the initial slow-wave polarization direction turns to be the fast-wave polarization direction. The distributed orientation model of Cretaceous shale (Hornby et al. 1994) shows a similar trend but the magnitude of shear splitting is smaller. The relative changes of the shear wave splitting between the crack-free state and the cracked state (crack density 0.1) are C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation Velocity anisotropy of phyllosilicate-rich rocks Figure 9. Relative phase velocity differences between two polarized shear waves, VS H and Vq S , for Hidaka biotite schist (solid curve) and the distributed model of cretaceous shale (dashed curve) by Hornby et al. (1994). The thick and thin curves correspond to the cases of crack density 0.1 and crack-free state (matrix), respectively. Fluid bulk modulus is 0.01 GPa for both cases. 10 per cent for the Hidaka biotite schist and 4 per cent for the distributed orientation model of Cretaceous shale, in the directions θ ≥ π/4. In the directions θ < π/4, the relative changes become smaller and come to zero along the symmetry axis. Observed magnitudes of shear wave splitting depend on the path distances of anisotropic media through which seismic waves propagate. Most of the anisotropic layers will not be thick enough to produce a large traveltime retardation of the two shear waves. However, a sophisticated method presented by Miyazawa et al. (2008) observed a 2 per cent shear wave splitting, and suggested that the splitting as small as ∼0.2 per cent can be detected. Considering this, shear wave splitting could be used for monitoring reservoir conditions by observing anisotropy of cap rocks. 5 C O N C LU S I O N S Rocks containing large amounts of sheet silicates or clay platelets often show TI-type anisotropy characterized by slow P- and S-wave velocities along the symmetry axis. This anisotropy is produced by the preferred orientation of sheet silicates or clay platelets. Shear wave anisotropy of such TI-type rocks is mostly characterized by increase of the qS-wave velocity (the slow velocity in axial directions) in the off-axis directions with the direction of the maximum value biased to the symmetry axis. Cracks tend to be formed with their planes parallel to the planes of silicate sheet in sheet silicate minerals. The cracks enhance the existing axial anisotropy caused by CPO of sheet silicate minerals, by reducing the velocities of the P and S waves along the symmetry axis. Reduction of the P-wave velocity along the symmetry axis is closely linked to the change of the off-axis qS-wave velocities. An increase of crack density enhances the diagonal asymmetry of the angular distribution of the qS-wave velocity, and gives changes in the axial dependence of shear wave splitting between SH and qS waves. Under certain conditions, a decrease of the off-axis qS velocity caused by cracks can appear as exchange of polarization direction between the fast and slow shear waves. Although biotite schist is not a major constituent of the Earth’s crustal rocks, the present results represent a typical case for the C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation 387 most common TI-type rocks, such as shales, in the shallower part of the Earth’s crust. Most of the previous measurements on velocity under confining pressure in sheet silicate-bearing rocks could be explained by the aligned crack model. Berryman’s formulation for characterizing seismic anisotropy in TI-type rocks can be successfully applied. Velocity changes due to aligned cracks in TI-type rocks are enhanced when bulk modulus of the crack-filling fluid decreases. For example, gas phase enhances crack effects more than liquid because of the more than three orders of magnitude smaller bulk modulus of gas. The velocity differences between SH and qS waves are more enhanced in the off-axis directions than in the axial directions in sheet silicate-rich rocks when sheet silicate minerals aligned in the same direction. Anisotropy may be useful for monitoring reservoir cap rocks. Thus this study would present important clues to interpret the mechanism of anisotropy in common TI-type rocks in the shallower part of the Earth. AC K N OW L E D G M E N T S We wish to thank M. Takanashi and K. Yasunaga for the velocity measurements under confining pressure. R. Kranz kindly reviewed the manuscript and gave the authors some suggestions about crack density and its effect on velocity. The comments from the reviewers were helpful to improve the manuscript. REFERENCES Aleksandrov, K.S. & Ryzhova, T.V., 1961. The elastic properties of rockforming minerals, II: layered silicates, Izv. Acad. Sci. USSR, Geophys. Ser., 12, 186–189. Anderson, D., Minster, B. & Cole, D., 1974. The effect of oriented cracks on seismic velocities, J. geophys. Res., 79, 4011–4015. Berryman, J.G., 2008. Exact seismic velocities for transversely isotropic media and extended Thomsen formulas for stronger anisotropies, Geophysics, 73, D1–D10. Barruol, G. & Mainprice, D., 1993. A quantitative evaluation of the contribution of crustal rocks to the shear-wave splitting of teleseismic SKS waves, Phys. Earth planet. Inter., 78, 281–300. Burlini, L. & Kunze, K., 2000. Fabric and seismic properties of Carrara marble mylonite, Phys. Chem. Earth (A), 25, 133–139. Carcione, J.M., 2007. Wave fields in real media: wave propagation in anisotropic, anelastic, porous and electromagnetic media, in Handbook of Geophysical Exploration I: Seismic Exploration, Vol. 38, eds Helbig, K. & Treitel, S., Elsevier, Amsterdam. Chesnokov, E.M., Tiwary, D.T., Bayuk, I.O., Sparkman, M.A. & Brouwn, R.L., 2009. Mathematical modelling of anisotropy of illite-rich shale, Geophys. J. Int., 178, 1625–1648. Crampin, S., 1987. Geological and industrial implications of extensive dilatancy anisotropy, Nature, 328, 491–496. Crampin, S. & Yedlin, M., 1981. Shear-wave singularities of wave propagation in anisotropic media, J. Geophys., 49, 43–46. Crampin, S., Boothe, D.C., Evans, R., Peacock, S. & Fletcher, J.B., 1990. Changes in shear wave splitting at Anza near the time of the North Palm Springs earthquake, J. geophys. Res., 95, 11 197–11 212. Dellinger, J. & Vernik, L., 1994. Do travel times in pulse-transmission experiments yield anisotropic group or phase velocities? Geophysics, 59, 1774–1779. Douma, J., 1988. The effect of the aspect ratio on crack-induced anisotropy, Geophys. Prospect., 36, 614–632. Eshelby, J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. Roy. Soc., A241, 376– 396. 388 O. Nishizawa and K. Kanagawa Hall S.A., Kendall, J.M., Maddock, J. & Fisher, Q., 2008. Crack density tensor inversion for analysis of damages in rock frame structure, Geophys. J. Int., 173, 577–592. Hill, R., 1952. The elastic behaviour of a crystalline aggregate, Proc. Phys. Soc., Lond. A65, 349–354. Hornby, B.E., 1998. Experimental laboratory determination of the dynamic elastic properties of wet, dry shales, J. geophys. Res., 103, 29 945–29 964. Hornby, B.E., Schawartz, L.M. & Hudson, J.A., 1994. Anisotropic effectivemedium modelling of the elastic properties of shales, Geophysics, 59, 1570–1583. Hudson, J.A., 1981. Wave speeds and attenuation of elastic waves in material containing cracks, Geophys. J. R. astr. Soc., 64, 133–150. Jakobsen, M. & Johansen, T.A., 2000. Anisotropic approximations for mudrocks: a seismic laboratory study, Geophysics, 65, 1711–1725. Johnston, J.E. & Christensen, N.I., 1994. Elastic constants and velocity surfaces of indurated anisotorpic shales, Surv. Geophys., 15, 481–494. Jones, L.E.A. & Wang, H.F., 1981. Ultrasonic velocities in Creataceous shales from the Williston basin, Geophysics, 46, 288–297. Kendall, J.-M. et al., 2007. Seismic anisotropy as an indicator of reservoir quality in siliciclastic rocks, in Structurally Complex Reservoirs, Vol. 292, pp. 123–136, eds Jolley, S., Barr, D., Walsh, J. & Knipe, R.J., Geological Society, Special Publications, London. Kern, H., 1993. P- and S-wave anisotropy and shear-wave splitting at pressure and temperature in possible mantle rocks and their relation to the rock fabric, Phys. Earth planet. Inter., 78, 245–256. Kern, H. & Wenk, H.-R., 1990. Fabric-related velocity anisotropy and shear wave splitting in rocks from the Santa Rosa mylonite zone, California, J. geophys. Res., 95, 11 213–11 223. Landau, L.D. & Lifshitz, E.M., 1959. Theory of Elasticity, Pergamon, Oxford, (3rd edition in 1985, Butterworth-Heinemann). Le Ravalec, M. & Guéguen, Y., 1996a. High and low frequency elastic moduli for a saturated porous/cracked rock (differential self consistent and poroelastic theories), Geophysics, 61, 1080–1094. Le Ravalec, M. & Guéguen, Y., 1996b. Comments on “The elastic modulus of media containing strongly interacting cracks” by Paul M. Davis and Leon Knopoff, J. geophys. Res., 101, 25 373–25 375. Lin, S. & Mura, T., 1973. Elastic fields of inclusions in anisotropic media (II), Phys. Status Solidi, A, 15, 281–285. Markov, M., Levine, V., Mousatov, A. & Kazatchenko, E., 2005. Elastic porperties of double-porosity rocks using the differential effective medium model, Geophys. Prospect., 53, 733–754. Mavko, G., Mukerji, T. & Dvorkin, J., 1998. The Rock Physics Handbook: Tools for Seismic Analysis in Porous Media, Cambridge University Press, Cambridge. Miyazawa, M., Snieder, R. & Venkataraman, A., 2008. Application of seismic interferometry to extract P- and S-wave propagation and observation of shear-wave splitting from noise data at Cold Lake, Alberta, Canada, Geophysics, 73, D35–D40. Mura, T., 1982. Micromechanics of Defects in Solids, Martinus Nijhoff Publishers, Hague. Nishizawa, O., 1982. Seismic velocity anisotropy in a medium containing oriented cracks-transversely isotropic case, J. Phys. Earth, 30, 331–347. Nishizawa, O. & Yoshino, T., 2001. Seismic velocity anisotropy in mica-rich rocks: an inclusion model, Geophys. J. Int., 145, 19–32. Nur, A. & Simmons, N., 1969. The effect of saturation on velocity in low porosity rocks, Earth planet. Sci. Lett., 7, 183–193. O’Connell, R.J. & Budiansky, B., 1974. Seismic velocities in dry and saturated solids, J. geophys. Res., 79, 5412–5426. Pros, Z, Lokajı́ček, T., Přkryl, R., Špičák, A, Vajdová, V. & Klı́ma, K., 1998. Elastic parameters of West Bohemian granites under hydrostatic pressure, Pure appl. Geophys., 151, 631–646. Sarout, J. & Guéguen, Y., 2008. Anisotropy of elastic wave velocities in deformed shales. Part 1—experimental results, Geophysics, 73, D75–D89. Sayers, C.M., 2002. Stress-dependent elastic anisotropy of sandstones, Geophys. Prospect., 50, 85–95. Simmons, G., Siegfried, R. & Feves, M., 1974. Differential strain analysis: a new method for examining cracks in rocks, J. geophys. Res., 79, 4383–4385. Soga, N., Mizutani, H., Spetzler, H. & Martin, R.J., 1978. The effect of dilatancy on velocity anisotropy in Westerly granite, J. geophys. Res., 83, 4451–4456. Takanashi, M., Nishizawa, O., Kanagawa, K. & Yasunaga, K., 2001. Laboratory measurements of elastic anisotropy parameters for the exposed crustal rocks from Hidaka metamorphic belt, central Hokkaido, Japan, Geophys. J. Int., 145, 33–47. Thomsen, L., 1986. Weak elastic anisotropy, Geophysics, 51, 1954–1966. Tsvankin, I., 1997. Anisotropic parameters and P-wave velocity for orthorhombic media, Geophysics, 62, 1292–1300. Vaughan, M.T. & Guggenheim, S., 1986, Elasticity of muscovite and its relationship to crystal structure, J. geophys. Res., 91, 4657– 4664. Vavryčuk, V., 2003. Parabolic lines and caustics in homogeneous weakly anisotropic solids, Geophys. J. Int., 152, 318–334. Verdon, J.P., Angus, D.A., Kendall, J.M. & Hall, S.A., 2008. The effect of microstrucuter and nonlinear stress on anisotropic seismic velocities, Geophysics, 73, D41–D51. Voltolini, M., Wenk, H.-R., Mondol, N.H., nut Bjorlykke, K. & Jahren, J., 2009. Anisotropy of experimentally compressed kaolinite-illite-quartz mixtures, Geophysics, 74, D13–D23. Walsh, J.B., 1965. The effect of cracks on the compressibility of rock, J. geophys. Res., 9, 381–389. Wang, Z., 2002. Seismic anisotropy in sedimentary rocks, part 2: laboratory data, Geophysics 67, 1423–1440. Yamamoto, K., Kosuga, M. & Hirasawa, T., 1981. A theoretical method for determination of effective elastic constants of isotropic composite, Sci. Rep. Tohoku Univ., Ser. 5, 28, 47–67. C 2010 The Authors, GJI, 182, 375–388 C 2010 RAS Journal compilation