Frequency-dependent seismic anisotropy in fractured rock

To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 349. Copyright © 1997 Elsevier Science Ltd Copyright © 1997 Elsevier Science Ltd Int. J. Rock Mech. & Min. Sci. Vol. 34, No. 3-4, 1997 ISSN 0148-9062 To cite this paper: Int. J. RockMech. &Min. Sci. 34:3-4, Paper No. 349 FREQUENCY-DEPENDENT SEISMIC ANISOTROPY IN FRACTURED ROCK W. Y i l ; K.T. Nihei2; J.W. R e c t o r l ; S. N a k a g a w a l ; L.R. Myer2; N.G.W. C o o k 1 1 Department of Materials Science & Mineral Engineering, University of California, Berkeley, CA 94720, USA z Earth Sciences Division, E.O. Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA ABSTRACT To examine the effects of multiple, aligned fractures in rock, we have developed a two-dimensional elastic finite difference code for fractured media. Fractures are incorporated into the model explicitly as displacement-discontinuity boundary conditions. The wavefield is computed using a 4th-order staggered grid scheme. Simulations were performed for a broadband explosion point source (center frequency 374 Hz) located at the center of the model. The model consisted of 90 horizontal fractures spaced approximately 1/8 of a wavelength apart. The normal and shear fracture stiffnesses were selected such that the normal incidence transmission coefficient is 0.6. The simulations show strong scattering attenuation of the P-wave in the vertical direction (± to the fractures) and channeling of guided waves in the horizontal direction (ll to the fractures). The same code was also used to model wave propagation in an anisotropic medium with equivalent effective moduli for the 90-fracture system. Significant differences between the amplitudes, velocities, and frequency content of the waves in the explicit and equivalent medium fracture models were observed. These differences result from frequency-dependent time delays and filtering across each fracture and channeling along fractures that are not included in the zero-frequency effective medium description. These effects are especially interesting because they illustrate that the dynamic properties of fractured rock include significant amplitude anisotropy that may prove useful in the characterization of fractured rock. Copyright © 1997 Elsevier Science Ltd KEYWORDS fracture anisotropy • fracture interface waves • transversely isotropic model • explicit fracture model • displacement-discontinuity boundary condition INTRODUCTION Single fractures in rock can give rise to a variety of interesting seismic wave phenomena, including low-pass filtering of transmitted waves, the generation of reflected and converted waves, and the guiding of fracture interface waves (Schoenberg 1980; Pyrak-Nolte et al. 1990a; Pyrak-Nolte et al. 1992). Although some work has been done on seismic wave propagation in rock with multiple fractures (e.g., Pyrak-Nolte et al. 1990b; Schoenberg, Sayers 1995), a comprehensive picture of the seismic wave phenomena produced by multiple fractures has yet to emerge. This paper uses numerical simulations to ISSN 0148-9062 To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 349. Copyright © 1997 Elsevier Science Ltd investigate the effects of multiple, parallel fractures on seismic waves propagating at normal, parallel and oblique incidence to the fractures and, in particular, examines dynamic seismic anisotropy in fractured rock. Conventional approaches for the seismic characterization of multiple, parallel fractures in rock typically utilize zero-frequency effective medium theories. The additional compliance of the fractures are lumped into normal and shear stiffnesses, and the overall stiffness of the fracture+rock system is captured in effective elastic moduli (e.g., White 1983; Schoenberg, Muir 1989; Schoenberg, Sayers 1995). Reducing the properties of the fractured rock system to its static effective properties results in P- and S-wave velocities that vary with respect to fracture orientation and shear wave splitting. However, since this approach is inherently a static (i.e., zero-frequency) approximation, it does not include frequency-dependent amplitude and velocity variations with respect to fracture orientation due to reflection losses across fractures and wave channeling along fractures. The objectives of this paper are to use numerical finite difference simulations: (1) to investigate the effects of multiple, parallel fractures on the amplitudes and velocities of compressional and shear waves, and (2) to examine differences in the static (zero-frequency) and dynamic amplitude and velocity anisotropy. FRACTURE ANISOTROPY Anisotropy in the elastic properties of rock can result from the combined effects of aligned microcracks, mineral grains, bedding planes, and layering that vary with direction. At the scale of a hydrocarbon reservoir, anisotropy may also be present in the form of aligned fractures. Because fractures can significantly affect the flow characteristics of the reservoir, seismic methods for determining the orientation of the fracture sets are of considerable interest to reservoir geophysicists. Approaches for estimating the static anisotropic elastic moduli for rock with aligned fracture sets have been developed by a number of investigators, including White 1983 and Schoenberg, Muir 1989. Recently, the applicability of these zero-frequency theories has been questioned (Pyrak-Nolte et al. 1990b; Frazer 1995). This has led to the development of numerical techniques for computing the dynamic anisotropic properties of fractured rock (Frazer 1995; Coates, Schoenberg 1995). In the following sections, this paper will describe a numerical scheme for modeling the complete wavefield in a medium with multiple parallel fractures. The velocities and amplitudes of this explicit fracture model are compared with those obtained from a static equivalent transversely isotropic (TI) fracture model. DISPLACEMENT-DISCONTINUITY FRACTURE MODEL Fractures in rock are thin, localized regions of compliance. This compliance results in sharp jumps in the normal and tangential displacements across the fracture, the magnitude of these jumps being proportional to the compliance of the fracture and the stress acting on the fracture. For a planar fracture oriented in the x-y plane, the displacement-discontinuity (Schoenberg 1980; Pyrak-Nolte et al. 1990a) boundary conditions for the fracture are given by u2 - u2 = " c K (1) z ISSN 0148-9062 To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 349. u:-u: = Copyright © 1997 Elsevier Science Ltd 'T "~ K (2) (3) ~x ='t:= (4) Here, the stresses are represented by gj, K i is the fracture stiffnesses in units of [Pa/m], and Ui is the particle displacement, with superscripts + and - referring to opposite sides of the fracture. The displacement-discontinuity boundary condition is a generalized boundary condition in the sense that it degenerates to the boundary conditions for a welded interface as ~:--->ooand to those for a traction-free interface as ~:--->0.This boundary condition is essentially a lumped properties approximation that captures the macromechanical stress-displacement behavior of the fracture. Because it is a linear model, it is valid provided that the magnitude of the seismic displacements are insufficient to mobilize nonlinear asperity deformation and frictional sliding, a condition which is typical in most geophysical investigations. Laboratory acoustic measurements of body wave transmission across single fractures (Pyrak-Nolte et al. 1990a) and fracture interface waves along single fractures (Pyrak-Nolte et al. 1992; Roy, Pyrak-Nolte 1995) have established the validity of the displacement-discontinuity model for describing the dynamic properties of fractures. FINITE DIFFERENCE SCHEME FOR FRACTURED ROCK Elastic wave propagation in a medium with multiple, parallel fractures was investigated by using an elastic, two-dimensional, staggered grid finite difference code based on the approach of Virieux 1986. The code uses fourth-order differencing in space, and second-order differencing in time (Levander 1988), +D:o:: ,] I ' "°" p,,., L D:v:-" .= 1 ",~" J ' [D, "" "" D*." , ,, ÷j4, =t.: llm.÷j~, j uv x *,",t "lr'~ + D:~;~q4.j = t_:,, . uv • II~+~ ,l . o C ._. L " ="*"~÷~ ~ i] +D:C:- -+~4'J "~ + C . ,D*v "-~1,, .+/1~, m+~,j +';M/ "m+l.j ~+J'~ , Zttl+~.]+Y2 +C C (5) D*I;v"-~;,, x. ,j |3m÷~ 'j D-v In these equations, v i is the particle velocity, g). are the components of stress, p is the material density, and t is time. The major difference of the staggered grid formulation from standard finite difference schemes is that the different components of the velocity field are not know at the same grid node. In the staggered grid scheme (Figure 1), the horizontal velocity and density are defined at location (m, j), the ISSN 0148-9062 To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 349. Copyright © 1997 Elsevier Science Ltd vertical velocity and density at (re+l~2, j+l/2), the normal stress and elastic constants Cll, C33 and C13 at (m+1/2, j), and tangential stress and the elastic c o n s t a n t C44 at (in, j+1/2). Only half the unknowns are required at each time step. The grid size of the staggered grid method is half the grid size of the conventional finite difference method. The primary advantages of the 4th order staggered grid scheme for this study are its computational efficiency, accuracy, and flexibility by which fractures can be modeled either explicitly, as displacement-discontinuity boundary conditions or implicitly using effective elastic constants. EQUIVALENT ELASTIC PROPERTIES OF FRACTURED ROCK A two-dimensional elastic medium containing a single set of parallel fractures each with fracture s t i f f n e s s e s Ki can be described by a transversely isotropic medium characterized by five elastic constants. The relationships between these constants and seismic velocities can be expressed by (Mavko et al. 1993), CII -~ P V2pll c , , = c,, - 2 0 % C.,. ~ (6) V 2 p,,, P- and SV-wave motions are mixed if the angle of propagation relative to the plane of symmetry is other than 0 or 90 degrees. Therefore, pure P- and SV- waves exist in a TI medium only in these directions. The velocities of quasi P- and quasi S- waves propagating in intermediate directions depend on C13 (which does not affect the velocities of P- and S- waves propagating parallel and perpendicular to the axis of symmetry), and can even produce cusps and triplication's in the wavefront of the quasi S-wave (Helbig 1966). The static equivalent properties of aligned fractures can be obtained by decomposing the elastic compliance of the fractured rock as the sum of the compliance of the host rock and the compliance for the set of parallel fractures. The resulting effective elastic moduli for the fractured rock mass can be represented by a transversely isotropic medium (Schoenberg, Sayers 1995), C. = rn.(1- 8,,) C. = ~ . ( 1 - 8 , , ) (7) where ISSN 0148-9062 To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 349. rn~ = ~ . + 2At b Copyright © 1997 Elsevier Science Ltd t~ = ZT/z" l+ZT#0 (8) Z,, m~ rb=-- 1 + Z,, m s 1 1 Z,, = ~ ZT = HK r H 1 ¢ ,~ where )~b and ~tb are the Lam6's constants, H is the fracture spacing, and •T and KN are the transverse and normal fracture stiffnesses, respectively. As KT and KN --->~ , the fractures become welded and the elastic constants become equal to those of an isotropic elastic medium defined by )~b and ~tb. COMPARISON OF TI M O D E L W I T H EXPLICIT F R A C T U R E M O D E L Finite difference simulations were performed to examine the differences between a model in which fractures are modeled explicitly as displacement-discontinuity boundary conditions and a model in which the fractures are modeled by their equivalent transversely isotropic (TI) properties described by Eq. (7) and (8). Figure 2 shows the numerical model which consists of a square region of 12 wavelengths in size. Simulations were performed using a broadband explosion source (first-derivative of a Gaussian function, central frequency equals 374 Hz) located at the center of the model. At the external boundaries of the models, hybrid absorbing boundary conditions were imposed. The explicit fracture model consisted of 90 horizontal fractures modeled as displacement-discontinuity boundary conditions described by normal and shear stiffnesses, ~:T = ~:N = 7.25 x 109 Pa/m (normal incidence P-wave transmission coefficient=0.6). The fracture spacing is approximately 1/8 of a wavelength (i.e., 7.5 fractures per wavelength). For the TI fracture model, the five elastic constants were selected to exactly represent the static effective properties of the 90 fractures system based on the relationships given in Eq. (7) and (8). Snapshots of the horizontal and vertical components of particle velocity are shown for the TI model and the explicit model at 16 ms (Figure 3). The explicit model shows strong scattering attenuation of the P-waves in the vertical direction (_1_to the fractures) and channeling of waves in the horizontal direction ({ to the fractures). The waves also exhibit triplication of the wavefronts in the horizontal direction. The vertical component of the particle velocity shows a stronger triplication than the horizontal component. In addition to triplication's, the energy of seismic waves appear focused along the vertical direction. A possible explanation for this phenomena is that seismic waves are multiply reflected, resulting in localized resonances in the vertical direction. In contrast, the TI model shows a strong P-wave arrival with an elliptical wavefront and smaller variability in amplitude with direction of propagation. These characteristics correspond to an elliptical anisotropic medium (Helbig 1983). The anisotropy ratio A- v ,~.~ _ v I V ~rmrnt~ (9) VII for both models were calculated from the peak P-wave amplitudes in the horizontal and vertical ISSN 0148-9062 To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 349. Copyright © 1997 Elsevier Science Ltd directions (Table 1). The explicit fracture model shows a small anisotropy ratio, indicating that the amplitude anisotropy is much stronger for this model. P-wave velocities as a function of the direction of propagation relative to the fractures are shown in Figure 4(a). At 0 degrees (parallel to fractures), the P-wave velocities for both models are identical. As the angle of propagation is increased from 0 to 90 degrees, the TI model shows consistently lower velocities than the explicit fracture model. This discrepancy in the velocities can be attributed to the frequency-dependent group time delay that results when a finite frequency wave propagates across a fracture (Pyrak-Nolte etal. 1990a), _ 2(g/Z) (10) 49C/Z) 2 where ~: is the fracture stiffness, Z is the acoustic impedance, and co is the angular frequency. From this equation, it is evident that the maximum time delay occurs in the static limit (i.e., as co-->0). Thus, the velocities in the TI fracture model will always be lower than the explicit fracture model since it is based on elastic constants derived from the (static) equivalent medium approximation. The vertical component of the P-wave amplitude is displayed as a function of the direction of propagation relative to the fractures in Figure 4(b). For an isotropic elastic medium without fractures, the vertical component of displacement for a P-wave generated by an explosion source should increase from zero at 0 degrees to a maximum at 90 degrees. The P-wave amplitudes for both models increase from zero at 0 degrees. The explicit fracture model, however, shows a marked decrease in P-wave amplitude at approximately 30 degrees. This sharp reduction in the P-wave amplitude at angles 30 to 90 degrees results from strong scattering attenuation of the transmitted wave. At normal incidence (90 degrees), the transmission losses resulting from reflections off each fracture can be estimated from the plane wave transmission coefficient I 21N (Pyrak-Nolte et al. 1990a) ITI ~- {, (I1) where N is the number of fractures. At 90 degrees, the calculated results agree well with the predictions of Eq. (11). The amplitude spectra of the vertical component of the P-wave at 90 degrees is displayed in Figure 5. The spectrum for the explicit fracture model is much lower than the TI model due the reflection losses described above. In addition, a 75 Hz shift in the central frequency of the peak amplitude is present in the explicit fracture model. This low pass filtering that results in the explicit fracture model is consistent with the behavior predicted by Eq. (11). A comparison of the wavefields in the TI and explicit fracture models demonstrate that the frequency-dependent nature of fractures can result in significant differences in the velocities and amplitudes. These results suggest that modeling fractured rock using the TI approximation may significantly underestimate the velocities and the amplitude anisotropy. SUMMARY ISSN 0148-9062 To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 349. Copyright © 1997 Elsevier Science Ltd This study has investigated elastic wave propagation in a medium containing multiple, parallel fractures, and, in particular, has examined frequency-related amplitude and velocity anisotropy. This analysis demonstrates that a fractured medium modeled as a static equivalent transversely isotropic (TI) medium will underestimate amplitude anisotropy and frequency-dependent velocities. Future work will focus on quantifying the effects of fracture spacing (relative to wavelength), fracture stiffness, and frequency on the dynamic anisotropic properties of fractured rock. ACKNOWLEDGMENTS This work was carried out with support in part by the Gas Research Institute under Contract No. 5093-260-2663 and by the Director, Office of Energy Research, Office of Basic Energy Sciences under U.S. Department of Energy Contract No. DE-AC03-76SF00098. FIGURES Paper 349, Figure 1. !:n. i÷1) (m. : 1 ) . [] © [m-l= j-l) 4) I , [] q . . . - . . . . . . [] ~/., 0 V, O E~ (m-1,.'2.i.1'.. V-] () ,z (m j-l) Figure 1. Stencil for staggered grid finite difference scheme. Black symbols are for particle velocities and buoyancy (l/p) at time nA. White symbols are for stresses and Lame's coefficients at time (n + 1)A. ISSN 0148-9062 To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 349. Copyright © 1997 Elsevier Science Ltd Paper 349, Figure 2. Explicit Fracture M~_d._e/ N--g0 TI Mode] .90 L . . . . , " "- - - ' - ' - ~ ' ~--.. "~'.'. r'~l .--. H,' ~" ' [.= - - - - " ":.1 ...... .... ... " ' .'..~._~. "~" ' ' . . . . . . . I -liP'- g _ .-;-~.-'~:. ~ liP- ,-4 C'1 - - .~;c,u.'c ~, , ~. ___..'.'_,.~_~:-._ ' ~ ,~ " .... "~C'. . .... ........ :.'-'.' . . . . . . . . . . . . !2 ),. "- ~ In- >. "- Figure 2. TI and explicit fracture models used in the finite difference simulations. The elastic moduli of the TI model were selected to match the static equivalent elastic properties of a medium containing 90 parallel fractures with 1<T = 1~ N - - 7.25x109 Pa/m. ISSN 0148-9062 To cite this paper: Int. J. R o c k M e c h . & Min. Sci. 34:3-4, p a p e r No. 349. C o p y r i g h t © 1997 Elsevier Science Ltd Paper 349, Figure 3. X-Companent -]-I Explicit Fracture Mo~el Nodel :., . •.. ... . ~--=,..~." ;-. ~::_~.~.- :.~. ...: it, :.'~ ~ ~ f > • • . ,1, ,, ~1.... , - ...: • :i [ ,==,~ . 1." l'. Z-Component ~'. .:%. " .~.. [...~ "." . • :. ..1 :.. ::.,,R>.'. x ' " : ~-,..3~:, i" 1..'" 1" .~.'.., i" ."'"," ~-..1. ,":• ., ." Figure 3. Snapshots of the horizontal (x) and vertical (z) component of particle velocity for the TI model and explicit fracture model (90-fracture) at 16 msec. ISSN 0 1 4 8 - 9 0 6 2 To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 349. Copyright © 1997 Elsevier Science Ltd Paper 349, Figure 4. P- Wave Group Velocity P-Wave- 2~cx,, •~e'~' Ampl itude 1 50 ~ [ .-,,I~TI ~,,,,::,~. ~ [ "-I!,'~ i fc'.3-:e 1,2- L: x p h "-- "~.'.¢':3:'.:fr" •. Cr;. 326"0 .= 3120,0 5 75 .3.. 9.5C ,'2 '12 i11 i- b_ ;.4~-3 --...._ % 22_'-5 OCt. -9.25 ,~02,5 ,~.%. I [ I I I I 1.5 31; -'5 r~ 75 g~-i ,C,n~lJe ~ I r..T.irlr:!l!:l ! (a) L~.2[= -'5.5:': .... r I I 1 i 1 "5 ~0 4',.. ¢2~; .~l~ :)~: a.rl{::l e o f I-i# der.~.e "..-"eg "~e', (b). Figure 4. (a) P-wave group velocity and (b) normalized vertical component of the P-wave amplitude for the TI model and the explicit model. ISSN 0148-9062 To cite this paper: Int. J. Rock Mech. & Mm. Sci. 34:3-4, paper No. 349. Copyright © 1997 Elsevier Science Ltd P a p e r 349, Figure 5. P-Wave Spectral Density g.5 1C s I 3 1.~" i 2.5 IC ~ L ~,loJ~,I E × c.li.cil 11 "r::):t<;. 2 1.~" 1.5 I,,','~ 1 1C,'~ 5 1(-" ~ -5 10"~J 20~ ,103 F.Ic:.~'IC? 600 8-2~ 1000 (' Iz', Figure 5. Spectral density for P-wave traveling perpendicular to fractures in the TI model and the explicit fracture model. TABLES P a p e r 349, Table 1. Table 1. ANISOTROPY RATIO FOR TI AND EXPLICIT FRACTURE MODELS Model Peak amplitude in Peak amplitude in Anisotropy vertical direction horizontal direction ratio (rn/s) Explicit TI 5.9x 10.5 2.0 xl0 "2 , (m/s) ,. 1.4 x 10.2 1.9 xl0 -2 0.041 1.05 t,,, ISSN 0148-9062 To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 349. Copyright © 1997 Elsevier Science Ltd References References Coates R. T., Schoenberg M. Finite-difference modeling of faults and fractures, Geophys., 60:5, 1514-1526. Frazer L.N. 1995. Dynamic elasticity of microbedded and fractured rocks, d. Geophys. Res., 95:B4, 4821-4831. Levander A. R 1988, Fourth-order, finite-difference P-SV seismograms, Geophys., 53, 1425-1436. Helbig K. 1966. A graphical method for the construction of rays and traveltimes in spherical layered media: Part 2: Anisotropic case, theoretical considerations, Bull. Seis. Soc. Am., 57, 527-559. Helbig K. 1983. Elliptical anisotropy- Its significance and meaning, Geophys., 48, 825-832. Mavko G., Mukerji T., Dvorkin J. 1993 Rock Physics Formulas, Stanford University Pyrak-Nolte L.J., Myer L.R., Cook N.G.W. 1990a. Transmission of seismic waves across single natural fractures, d. Geophys. Res., 95:B6, 8617-8638. Pyrak-Nolte L.J., Myer L.R., Cook N.G.W. 1990b. Anisotropy in seismic velocities and amplitudes from multiple parallel fractures, d. Geophys. Res., 95:B7, 11345-11358. Pyrak-Nolte L.J., Xu J., Haley G.M 1992. Elastic interface waves propagating in a fracture, Phys. Rev. Lett., 68, 3650-3653. Roy S., Pyrak-Nolte L.J. 1995. Interface waves propagating along tensile fractures in dolomite, Geophys. Res. Lett., 22:20, 2773-2776. Schoenberg M. 1980 Elastic wave behavior across linear slip interfaces, d. Acoust. Soc. Am., 68:5, 1516-1521. Schoenberg M., Muir F. 1989 A calculus for finely layered anisotropic media, Geophys., 54:5, 581-589. Schoenberg M., Sayers C.M. 1995. Seismic anisotropy of fractured rock, Geophys., 60:1,204-211. Virieux J. 1986. P-SV-Wave propagation in heterogeneous media: Velocity-stress, finite-difference method, Geophys., 51:4, 889-901. White J. E. 1983 Underground Sound, Elsevier, New York, USA ISSN 0148-9062