Misalignment of the orientation of fractures and the principal axes for P and S waves in rocks containing multiple non-orthogonal fracture sets

Abstract

Azimuthal anisotropy in rocks can result from the presence of one or more sets of partially aligned fractures with orientations determined by the stress history of the rock. A shear wave propagating in an azimuthally anisotropic medium splits into two components with different polarizations if the source polarization is not aligned with the principal axes of the medium. In the presence of two or more non-orthogonal sets of vertical fractures, the symmetry of the medium may be approximated as monoclinic with a horizontal plane of mirror symmetry if, in the absence of fractures, the rock is transversely isotropic (TI) with the symmetry axis perpendicular to the bedding plane. For such a medium, the fast and slow polarization directions for vertically propagating shear waves are not parallel or perpendicular to any of the fracture planes but lie in directions given by the principal axes of a second-rank fracture compliance tensor. This tensor is independent of the normal compliance of the fractures. It follows that for vertical propagation of shear waves in a vertically fractured TI medium, any number of arbitrarily oriented vertical fracture sets is equivalent to two mutually perpendicular fracture sets, provided that the seismic wavelength is large compared to the size of the fractures. For offsets typical of surface seismic acquisition, the P-wave velocity at fixed offset varies with azimuth as cos2(φ − φ₀ ), where φ is the azimuth measured with respect to the fast polarization direction for a vertically polarized shear wave; φ₀ depends on both the normal and the shear compliance of the fractures and may differ from zero if the ratio of the normal to shear compliance of the fractures varies significantly between fracture sets. If this ratio is similar for all fractures, φ₀ ≈ 0 and the principal axes of the variation in P-wave velocity with azimuth for fixed offset are determined by the principal axes of the second-rank fracture compliance tensor.

anisotropy, cracked media, fractures, shear-wave splitting.

Introduction

The purpose of this paper is to examine the misalignment of the orientation of fractures and the principal axes for P- and S-wave propagation which may occur in sedimentary rocks containing two or more non-orthogonal fracture sets. An understanding of this relation is important for determining fracture orientation in fractured reservoirs using seismic data. Natural fractures in reservoirs play an important role in determining fluid flow during production, and hence the density and orientation of fractures are of great interest (Reiss 1980; Nelson 1985). The fractures with the largest apertures at depth tend to be oriented normal to the direction of the minimum in situ stress and may therefore lead to significant permeability anisotropy in the reservoir (Sayers 1990).

In the absence of fractures, sedimentary rocks may often be described, to a good approximation, as being transversely isotropic (TI) with a symmetry axis oriented perpendicular to the bedding plane (Jones & Wang 1981; Tosaya 1982; Hornby 1994). However, many formations contain fractures with orientations determined by the stress history of the rock rather than by the orientation of the bedding plane. A rock possessing both an anisotropic fabric and a preferred orientation of vertical fractures will display azimuthal anisotropy, the seismic velocities varying with azimuth about the symmetry axis of the unfractured sediment.

A vertically propagating shear wave in an azimuthally anisotropic medium will split into two components with different polarizations if the source polarization is not aligned with the principal axes of the medium (Crampin 1985). In the case of a single set of aligned, vertical fractures, the fast shear wave (S1) for vertical propagation has particle motion parallel to the fracture planes, whilst the slow shear wave (S2) has particle motion perpendicular to the fracture planes. A measurement of the fast polarization direction can therefore be used to determine the orientation of the fractures. The presence of shear-wave splitting in sedimentary basins is now well established (Lefeuvre 1994). Examples cited by Lefeuvre (1994) include the use of converted P-to-S waves in surface-to-surface measurements (Harrison & Stewart 1993), converted P-to-S waves in surface-to-borehole measurements (Garotta & Granger 1988; Lefeuvre & Queen 1992), pure S waves in surface-to-surface measurements (Alford 1986; Lynn & Thomsen 1990; Davis & Lewis 1990, Mueller 1991) and pure S waves in surface-to-borehole measurements (Cliet et al. 1991; Winterstein & Meadows 1991; MacBeth & Crampin 1991; Lefeuvre et al. 1992; Lefeuvre et al. 1993; Queen et al. 1992).

Although shear waves are considered to be more reliable indicators of fracture orientation than P waves, considerable interest remains in the use of P waves for determining fracture orientation, since these form the basis of most commercial seismic surveys (Lefeuvre 1994; Lynn et al. 1994). Chang & Gardner (1992) suggested that the fracture orientation of a subsurface fracture zone may be determined by analysing P-wave interval velocities, and a method for inverting the results for fracture orientation has been proposed (Sayers 1995). Grechka & Tsvankin (1996) showed that the azimuthal variation in P-wave normal moveout (NMO) velocity, which describes the variation in P-wave velocity with small offset, is elliptical for arbitrary anisotropy. For a single set of vertical fractures, the principal axes of this ellipse are aligned parallel and perpendicular to the fractures (Sayers 1995; Grechka & Tsvankin 1996; Tsvankin 1996). A measurement of the azimuthal variation in P-wave NMO velocity may therefore be used to determine the orientation of a single set of vertical fractures.

Although the orientation of a single set of vertical fractures can be determined by measuring the fast and slow shear-wave polarization directions, or the azimuthal variation in P-wave NMO velocity, many reservoirs contain two or more sets of fractures with different orientation (Reiss 1980; Nelson 1985). Frequently, the different fracture sets are not orthogonal, but have some angle other than 0° or 90° between their strikes (Reiss 1980; Nelson 1985). For two or more non-orthogonal sets of vertical parallel fractures within a transversely isotropic rock having a vertical axis of symmetry, the symmetry of the rock can be approximated as monoclinic (Winterstein 1990). For such a medium the fast and slow polarizations for a vertically propagating shear wave will not coincide with any of the fracture planes but will lie at some intermediate angle (Liu et al. 1993; MacBeth 1996). The principal axes of the P-wave NMO ellipse will also not coincide with any of the fracture planes. In this paper, the relation between the orientation of the fractures and the fast and slow shear-wave polarization directions in a medium containing two or more sets of non-orthogonal vertical fractures is investigated, and the deviation between these directions and the principal axes describing the variation of P-wave velocity with azimuth is examined.

Vertical Propagation In Monoclinic Media

It is assumed that the medium is monoclinic with a horizontal plane of mirror symmetry. A sedimentary rock containing several sets of fractures with normals lying in the bedding plane is an example of a medium with a single plane of mirror symmetry if, in the absence of fractures, the rock is transversely isotropic with the symmetry axis perpendicular to the bedding plane (Winterstein 1990). It is convenient to choose the coordinate axes x₁ ,x₂ ,x₃ with the normal to the mirror plane parallel to x₃ . A monoclinic medium with a horizontal plane of reflection symmetry has the elastic stiffness tensor

 

(1)

in the conventional (two-subscript) condensed 6 × 6 matrix notation (Nye 1985). For an arbitrary choice of axes having x₃ perpendicular to the mirror plane, the elastic behaviour of the material is apparently described by the 13 independent c_ij . However, a choice of the axes x₁ and x₂ exists such that c₄₅ = 0 and therefore the number of independent elastic constants in this coordinate frame is reduced to 12 (Norris 1989; Helbig1993). This choice corresponds to a rotation by α about the x₃ axis with

 

(2)

(Helbig1993).

For vertical propagation in a monoclinic medium, the Christoffel equation gives

 

(3)

(Fryer & Frazer 1987), where q is the vertical slowness. The root q² = ρ/c₃₃ represents a pure P mode. The shear slownesses are roots of the equation

 

(4)

For the choice of axes for which c₄₅ = 0, eq.(4) factorizes completely, and the two shear waves S1 and S2 have vertical slownesses given by

 

(5)

 

(6)

Seismic Anisotropy Of Fractured Reservoirs

In this section, expressions for the elastic stiffnesses of a transversely isotropic medium containing fractures of arbitrary orientations are given following the approach of Sayers & Kachanov (1991, 1995). This approach has been shown to agree with numerical simulations (Kachanov 1992). Other approaches may also be used. For example, Hudson (1990, 1991) and Peacock & Hudson (1990) have extended the theory of Hudson (1981) for aligned circular cracks to arbitrary crack orientation distributions.

The effective elastic compliance tensor s_ijkl of a rock containing fractures relates the average strain ε_ij over a representative volume V to the average stress components σ_ij :

 

(7)

For fractures, ε_ij may be written in the form

 

(8)

where sijkl_b is the compliance tensor of the unfractured background rock, which may be of arbitrary anisotropy; S_r is the surface of the rth fracture lying within V; n_i are the components of the local unit normal to the fracture surface, which may in general be curved; and square brackets [] denote jump discontinuities in the displacement; see, for example, Sayers & Kachanov (1991, 1995). Note that eq.(8) is applicable to finite, non-planar fractures in the long-wavelength limit, i.e. the applied stress is assumed to be constant over the representative volume V .

When the fractures are approximately planar, and their unit normal is denoted by n, with components n_i , a linearity assumption is conveniently introduced through a ‘fracture compliance tensor’B (Kachanov 1992) with components B_ij such that, for the rth fracture,

 

(9)

Eq.(8) then gives

 

(10)

where Δs_ijkl is the change in compliance due to the presence of fractures. Note that this equation is exact for a parallel set of infinite fractures (Schoenberg & Sayers 1995).

The simplest assumption concerning the properties of the fractures is to let the normal compliance of the fractures be given by B_N and the tangential compliance by B_T . With this assumption,

 

(11)

so that the fracture behaviour is invariant with respect to rotation about an axis normal to the fracture. This behaviour is also assumed in the penny-shaped crack model of Hudson (1981). In general, the tangential compliance of the fracture will be a function of direction in the plane of the fracture (Schoenberg & Douma 1988). In this paper it is assumed that this variation can be neglected. It follows that (Sayers & Kachanov 1995)

 

(12)

where

 

(13)

and

 

(14)

Note that α_ij and β_ijkl are symmetric with respect to all rearrangements of the indices, so that, for example, β₁₁₂₂=β₁₂₁₂ , β₁₁₃₃ = β₁₃₁₃ , etc.

If B_N = B_T for all fractures, Δs_ijkl is completely determined by the second-rank tensor α_ij . Δs_ijkl is then orthotropic, with principal axes coinciding with the principal axes of α_ij (Vakulenko & Kachanov 1971; Kachanov 1980; Sayers & Kachanov 1991). If B_N ≠ B_T the fourth-rank tensor β_ijkl is also required and may cause deviations from orthotropy.

Multiple Sets Of Vertical Fractures

For vertical fractures having normals lying within the x₁ x₂ plane, the non-vanishing α_ij and β_ijkl are α₁₁ , α₁₂ , α₂₂ , β₁₁₁₁ , β₁₁₁₂ and β₁₁₂₂ = β₁₂₁₂ , β₁₂₂₂ and β₂₂₂₂ .

Since α_ij is a second-rank tensor it may be diagonalized with an appropriate choice of axes. With this choice, α₁₂ = 0, and the non-vanishing independent Δs_ijkl follow from eq.(12):

 

(15)

 

(16)

 

(17)

 

(18)

 

(19)

 

(20)

 

(21)

 

(22)

To determine the c_ijkl it is convenient to transform the s_ijkl to the conventional (two-subscript) condensed 6 × 6 matrix notation, 11→1, 22→2, 33→3, 23→4, 13→5, 12→6, with factors 2 and 4 introduced as follows (Nye 1985):

 

 

 

The compliance tensor of the fractured medium resulting from eqs (15)–(22) has the form

 

(23)

for a transversely isotropic background medium.

The stiffness tensor follows upon matrix inversion and is given by eq.(1), with components

 

(24)

 

(25)

 

(26)

 

(27)

 

(28)

 

(29)

 

(30)

 

(31)

 

(32)

 

(33)

 

(34)

 

(35)

 

(36)

where

 

(37)

Shear-wave polarization

If the background medium is transversely isotropic with a vertical axis of symmetry, it follows from eq.(34) that the choice of axes x₁ and x₂ which diagonalizes α_ij also gives c₄₅ = 0. It follows that the fast and slow polarization directions for a vertically propagating shear wave are aligned with x₁ and x₂ chosen in this way. The fast and slow directions for vertically propagating shear waves therefore coincide with the principal directions of α_ij . The azimuth of the fast direction φ_S1 is therefore given by

 

(38)

Thus for vertical propagation of shear waves in a vertically fractured TI medium, any number of arbitrarily oriented vertical fracture sets is equivalent to two mutually perpendicular fracture sets, provided that the seismic wavelength is large compared to the size of the fractures.

It can be seen from eqs (33), (35), (19) and (20) that, for this choice of axes, the shear-wave slownesses for vertical propagation in a vertically fractured medium depend only on the components α₁₁ and α₂₂ of the second-rank tensor α_ij . The shear-wave anisotropy may be defined in terms of the velocities v_S1 and v_S1 of the fast and slow vertically propagating shear waves as follows:

 

(39)

where μ = c_44b = c_55b is the shear modulus governing vertical shear-wave propagation in the transversely isotropic background medium.

For the case of a single set of aligned, vertical fractures, the fast shear wave (S1) for vertical propagation has particle motion parallel to the fracture planes, whilst the slow shear wave (S2) has particle motion perpendicular to the fracture planes. A measurement of the fast polarization direction can therefore be used to determine the orientation of the fractures. In the presence of two or more sets of non-orthogonal vertical fractures, the fast and slow polarization directions for vertically propagating shear waves are not parallel or perpendicular to any of the fracture planes but lie in directions given by the principal axes of the second-rank fracture compliance tensor α_ij . A measurement of the fast and slow polarization directions therefore defines a convenient choice of axes for which c₄₅=0. The deviation of the fast and slow shear-wave polarization directions from the strike of the fractures will be illustrated below for the case of two non-orthogonal vertical fracture sets.

P-wave anisotropy

The variation of the qP phase velocity with angle in a weakly anisotropic medium of arbitrary symmetry has been studied recently by Sayers (1994), Gajewski & Psencik (1996) and Mensch & Rasolofosaon (1997), and, for orthotropic media, by Tsvankin (1996). For the present case, the equation of Gajewski & Psencik (1996) reduces to

 

(40)

Here n_i are the components of the unit vector in the slowness direction, and should not be confused with the components of the fracture normal used elsewhere in this paper; v_P (0) is the P-wave speed for propagation along Ox₃ ; and n₁ = cosφsinθ, n₂ = sinφsinθ and n₃ = cosθ, where θ is the polar angle measured from the x₃ direction and φ is the azimuthal angle measured from the x₁ direction.

For small polar angle θ, eq.(40) may be written in the form (Sayers 1994)

 

(41)

where

 

(42)

Here

 

(43)

 

(44)

 

(45)

Eq.(42) may be written as

 

(46)

where

 

(47)

The principal axes for P- and S-wave propagation will only coincide if φ₀ is zero. This occurs only if c₃₆ = 0 (see eq.45).

It is interesting to compare eq.(41) with Thomsen's result for near-vertical propagation in a transversely isotropic medium with a vertical axis of symmetry (Thomsen 1986):

 

(48)

where

 

(49)

δ may be written in the form (Sayers 1994)

 

(50)

where

 

(51)

The quantity χ(φ) is therefore a generalization of χ to the case of azimuthal anisotropy resulting from multiple vertical fracture sets. For small anisotropy, these results are in agreement with the results of Grechka & Tsvankin (1996), who showed that the azimuthal variation in P-wave NMO velocity is elliptical for arbitrary anisotropy. If the background medium is isotropic, with Lamé elastic parameters λ and μ and Poisson's ratio ν, and the fracture-induced anisotropy is small, eqs (24)–(38) simplify considerably to give

 

(52)

 

(53)

 

(54)

The principal axes are then given by

 

(55)

In the presence of a single set of vertical fractures, the principal axes for P waves coincide with the crack strike and normal directions and can be used to predict these directions. In the presence of multiple vertical fracture sets, the principal axes for P and S waves coincide only if β₁₁₁₂ + β₁₂₂₂ = 0. Thus the principal directions determined from the azimuthal variation in P-wave velocity determined from seismic data may be different from those determined from the fast and slow shear-wave polarization directions. This occurs because the normal and shear compliances of the fractures may be different. If they are the same, the principal axes of P and S waves coincide.

Two sets of parallel vertical fractures

To illustrate the above results, consider the important case of two sets of parallel vertical fractures (Kachanov 1992) with azimuths φ₁ and φ₂ in an isotropic medium as shown in Fig.1. The assumption of isotropy of the unfractured rock is expected to be a reasonable approximation for most reservoir sandstones and limestones. It follows from eqs (13) and (14) that

 

(56)

and

 

(57)

Let Δφ = φ₂ − φ₁ be the difference in azimuth of the strikes of fracture sets 1 and 2 (Fig.1) and let A₁ = (1 − η)A, A₂ = ηA. Fig.2 shows the azimuth φ_S1 of the fast polarization direction for vertically propagating shear waves, given by eq.(38), measured with respect to the average azimuth (φ₁ + φ₂ )/2 of the two fracture sets as a function of η. This is independent of the properties of the isotropic background medium. For small Δφ, φ_S1 varies approximately linearly with η, in agreement with the empirical results found by Liu et al. (1993) for a range of intersecting vertical fracture systems, and with the approximation of MacBeth (1996). However, for larger Δφ the fast direction varies non-linearly with fracture density, with the fast polarization direction remaining closer to the dominant fracture strike than is predicted by the linear approximation.

Figure 1.

Two fracture sets having azimuths w1 and w2, with average azimuth (w1+w2)/2 and azimuthal difference Dw=w2-w1. The azimuth of the fast polarization direction for a vertically propagating shear wave is denoted by wS1.

Figure 2.

Azimuth wS1 of the fast polarization direction for vertically propagating shear waves measured with respect to the average azimuth (w1+w2)/2 of the two fracture sets as a function of g.

Figure 3.

Contour plots of w0 given by eq. (55) as a function of B1/A1 and B2/A2 for A1=A2 and Dw=60°.

Figure 4.

Contour plots of w0 given by eq. (55) as a function of B1/A1 and B2/A2 for A1=A2 and Dw=30°.

Consider now the effect of the fourth-rank tensor β_ijkl on the variation in P-wave velocity v_P (θ,φ) with azimuth φ for fixed offset. The values of β₁₁₁₂ and β₁₂₂₂ in eq.(55) are zero if the two sets of fractures are orthogonal, or if A₁ = A₂ and B₁ = B₂ so that both fracture sets have equal weight. Figs 3–5 show contour plots of φ₀ given by eq.(55) as a function of B₁ /A₁ and B₂ /A₂ for the cases (1) A₁ = A₂ , Δφ = 60° (Fig.3), (2) A₁ = A₂ , Δφ = 30° (Fig.4) and (3) A₁ = 2A₂ , Δφ = 60° (Fig.5). The principal axes for small-offset P-wave propagation are seen from eq.(55) to depend only on the Poisson's ratio of the background rock, which is assumed to take the value ν = 1/4, corresponding to a V_P /V_S ratio of 1.732. Brine-saturated sands and shales typically have values of V_P /V_S larger than this, and therefore have values of ν > 1/4. The values of φ₀ given in the figures therefore underestimate the deviations which may occur for brine-saturated rocks, since φ₀ increases with increasing ν (eq.55). However, sandstones saturated with gas or light oil may have values of ν < 1/4. For fractures, B_N /B_T is expected to be less than or equal to one. For gas-filled open fractures B_N /B_T ≈ 1, but a lower ratio of B_N /B_T may result from the presence of cement or clay within the fractures, or from the presence of a fluid with non-zero bulk modulus. For B_N /B_T < 1, it is expected from eqs (13) and (14) that B₁ /A₁ and B₂ /A₂ in eqs (56) and (57) will be negative. It can be seen in Figs 3-5 that φ₀ is small if B₁ /A₁ ≈ B₂ /A₂ . This will occur if the ratio B_N /B_T is the same for all fractures. If this ratio varies between fracture sets, φ₀ may be non-zero and the principal axes of the variation in P-wave velocity with azimuth for fixed offset may depart from the fast and slow polarization directions for a vertically propagating shear wave.

Figure 6.

Conclusions

A shear wave propagating in an azimuthally anisotropic medium splits into two components with different polarizations if the source polarization is not aligned with the principal axes of the medium. In the presence of two or more non-orthogonal sets of vertical fractures, the symmetry of the medium may be approximated as monoclinic with a horizontal plane of mirror symmetry if, in the absence of fractures, the rock is transversely isotropic with the symmetry axis perpendicular to the bedding plane. For TI media containing multiple sets of vertical fractures, the fast and slow polarization directions for a vertically propagating shear wave are given by the principal axes of the second-rank fracture compliance tensor α_ij defined by eq.(13). This tensor is independent of the normal compliance of the cracks. It follows that for vertical propagation of shear waves in a vertically fractured TI medium, any number of arbitrarily oriented vertical fracture sets is equivalent to two mutually perpendicular fracture sets, with strikes parallel to the principal directions of α_ij , provided that the seismic wavelength is large compared to the size of the fractures. For two non-orthogonal sets of parallel fractures, with azimuthal difference Δφ less than about 30°, the fast polarization direction is a weighted average of the strikes of the two fracture sets, in agreement with the conclusions of Liu et al. (1993) and MacBeth (1996). However, for larger Δφ the fast direction varies non-linearly with fracture density, with the fast polarization direction remaining closer to the dominant fracture strike than is predicted by the linear approximation.

For offsets typical of surface seismic acquisition, the P-wave velocity at fixed offset varies with azimuth as cos2(φ − φ₀ ), where φ is the azimuth measured with respect to the fast polarization direction for a vertically polarized shear wave (Sayers 1994; Gajewski & Psencik 1996). The value of φ₀ depends on the ratio of the normal to the shear compliance of the fractures, B_N /B_T . If this ratio is the same for all fractures, φ₀ = 0 and the principal axes of the variation in P-wave velocity with azimuth for fixed offset are determined by the principal axes of the second-rank fracture compliance tensor. If this ratio varies between fracture sets, φ₀ may be non-zero and the principal axes of the variation in P-wave velocity with azimuth for fixed offset may depart from the fast and slow polarization directions for a vertically propagating shear wave. For gas-filled, open fractures B_N /B_T ≈ 1, but a lower ratio of B_N /B_T may result from the presence of cement or clay within the fractures, or from the presence of a fluid with non-zero bulk modulus.