Numerical investigation of fracture‐induced seismic anisotropy

Numerical Investigation of Fracture-Induced Seismic Anisotropy 1 1 2 1 Weidong Yi , Seiji Nakagawa , Kurt T. Nihei , James W. Rector , Larry R. Myer 2 & Neville G.W.Cook 1 1 Department of Materials Science & Mineral Engineering, 579 Evans Hall, University of California, Berkeley, CA 94720, USA 2 Earth Sciences Division, E.O. Lawrence Berkeley National Laboratory, 1 Cyclotron Road, MS 90-1116, Berkeley, CA 94720, USA SUMMARY A two-dimensional elastic finite difference code was developed to examine the characteristic of seismic waves propagation in a fractured rock containing multiple, aligned fractures. The displacement-discontinuity boundary conditions were used to model the fractures explicitly. The effect of fracture spacing on the wavefiled generated by an explosion source was examined by changing the spacing from 0.6 to 0.15 of a wavelength. A wavefiled was computed for a transversely anisotropic (TI) medium with elastic moduli equivalent to the effective static moduli of the most densely fractured system. The results showed significant differences between the amplitudes, velocities, and frequency content of the waves in the explicit and equivalent medium fracture models. 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 lead to an unusually strong velocity and amplitude anisotropy which cannot explained by the TI medium approximation. The characterization may prove useful in characterizing fractures in reservoir rock. 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 much 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 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 a zero-frequency effective medium approximation. The additional compliance of the fractures are lumped into normal and shear stiffnesses, and the overall stiffness of the fracture plus rock system is captured in effective elastic moduli (e.g., White, 1983; Schoenberg and Muir, 1989; Schoenberg and 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 and 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). Fractures in rock are thin, localized regions of low compliance. This compliance results in sharp jumps in the normal and tangential components of displacement across the fracture, the magnitude of the jump 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; PyrakNolte et al., 1990a) boundary conditions for the fracture are given by τ u x+ − u x− = zx κx (1) + = τ− τ zz zz + = τ− τ zx zx (2) Here, the stresses are represented by τij, κ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 displacementdiscontinuity boundary condition is a generalized boundary condition in the sense that it degenerates to the boundary conditions for a welded interface as κ→ ∞ and to those for a traction-free interface as κ→0. 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 and Pyrak-Nolte, 1995) have established the validity of the displacement-discontinuity model for describing the dynamic properties of fractures. the external boundaries of the models, boundary conditions were imposed. transmission coefficient as 0.99 for normal incident Pwave transmission. The fracture spacing was approximately as 1/8, 1/4/ and 1/2 of a wavelength. For the TI fracture model, the five elastic constants were selected to exactly represent the static effective properties of the 80-fracture system. Explicit Fracture Model 900 Receiver 2 0 Source 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), 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). 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. COMPARISON OF TI MODEL WITH EXPLICIT FRACTURE MODEL 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 . Figure 1 shows the domain of computation which consists of a square region of 12 wavelengths in size. The simulations were performed using a broadband explosion source (first-derivative of a Gaussian function, central frequency of 374 Hz) located at the center of the model. At absorbing The explicit fracture model consisted of 80, 40 and 20 horizontal fractures modeled as displacement-discontinuity boundary conditions with normal and shear stiffnesses, κT = κN = 6.8x1010 Pa/m, which yields the 12 λ t u z+ − u z− = zz kz 0 Receiver 1 12 λ Figure 1: The geometry of the TI and explict 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 80 10 6.8x10 Pa/m. parallel fractures with κT = κN = Snapshots of the horizontal and vertical components of particle velocity at 16 ms are shown in Figure 2 for the explicit fractures models with different number of fractures per wavelength. As the number of fractures per wavelength increases, the velocity of the waves propagation in the vertical direction decreases due to the accumulation of the group time delays across the fractures. The amplitude of the scattered incoherent waves becomes small compared with distinct coherent waves as the fractured system becomes as effectively homogeneous medium. The scattering and conversion of waves introduces slower propagating secondary wavefronts. For models with closely spaced fractures, these wavefronts become parabolic. Superposition of the two parabolic arcs forms a triplication in the wavefront, a characteristic of the wavefield observed for very strongly anisotropic TI medium. A snapshot at 16 ms for the model is shown in Figure 3. Compared with the results for the explicit model, the shape of the outer first wavefrontis more elliptic and the velocity anisotropy is weaker than to the explicit model. It is also noted that triplication is not observed in the inner second wavefront and scattered waves are not pertbertion the wavefield. In addition to triplication, it should also be noted that seismic waves energy appear focused along the vertical direction. A possible explanation for this phenomena is that seismic waves are multiply reflected, resulting in localized resonance in the vertical direction. Explicit Fracture Model X-Component N =20 fractures. Wave follow following the first part of waves also show ‘ringing’ which seems to decay for the 80fracture model, but still show significant amplitude. TI Model X-Component Z-Component Z-Component N =40 Fig 3 Snapshots of the horizontal (x) and vertical (z) component of particle velocity for the TI model with the static effective properties of the 80-fracture at 16 ms. Horizontal Component Traces from the System Containing 20, 40 and 80 Fractures at Receiver 1 (0 degree) 0.35 N = 80 N =80 Particle Velocity (m/s) 0.3 0.25 N = 40 0.2 0.15 N = 20 0.1 0.05 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Time (sec.) Fig 4 The seismic waveforms for the horizontal component received at Receiver 1 for explicit fracture model. Fig 2 Snapshots of the horizontal (x) and vertical (z) component of particle velocity for the explicit fracture model (80, 40 and 20fracture) at 16 ms. The horizontal component of observed waves at Receiver 1 and (Fig.1 ) and the vertical component at Receiver 2 are shown in Fig 4 and 5 respectively. Due to the symmetry of the problem, other components do not have significant amplitudes. Similar results for TI medium are shown in Figure 6. The first arrival of the waves observed in the horizontal direction (Figure 4) shows a very similar waveform followed by reflected or guided waves which vary for different fracture spacing. For 80-fracture case, the late arriving waves decay as the number of the fractures per wavelength increases , and the medium approach those of the effectively homogeneous TI medium. On the other hand, waves observed in the vertical direction show quite different behavior for differing numbers of fractures. Unlike waves in the horizontal direction, velocity, frequency, and the amplitude of the first arrival decrease the number of the fractures increases. This due to the frequency dependent filtering effect of the The most remarkable feature of the waves observed at Receiver 2 for explicit model is that the second wavefront amplitude increases significantly as the number of the fractures increases. One dimension al simulations have revealed that this second wave cannot be seen for plane waves propagating perpendicular to the fractures. From the snapshots shown in Figure 2, it can be seen that this second wave is the ‘pseudo S-wave’ with particle motion approximately parallel to the wavefront. This wave cannot be seen in Figure 6 because for TI medium used for this simulation, the wavefront becomes horizontal at Receiver 2, losing both vertical and horizontal components of the particle velocity. For the explicit models, triplication leads to the coexistence of two wavefronts which intersect at an angle on the vertical axis. When the particle motion of the two wavefront superposed, a large vertical motion results on the vertical axis. The amplitude of second increases as the triplication becomes more distinct for models with smaller fracture spacing. Vertical Component Traces from the System Containing 20, 40 and 80 Fractures at Receiver 2 (90 degrees) 0.35 0.35 N = 80 Particle Velocity (m/s) Particle Velocity (m/s) 0.25 N = 40 0.2 0.15 Receiver 2 X-Component Receiver 1 0.25 0.2 0.15 0.1 N = 20 0.1 Z-Component 0.3 0.3 0.05 0.05 0 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.005 0.04 Time (m/s) Fig 5 The seismic waveforms for the vertical component of particle velocity received at Receiver 2 . CONCLUSION This study has used finite difference simulations to investigate the effects of multiple, parallel fractures on seismic waves generated by a point explosion source. These simulations reveal that even when the stiffness of the fracture is high (i.e., P-wave normal incidence transmission coefficient of 0.99), the accumulation of frequency-dependent time delays and reflection losses arise because the effective medium fracture model is an inherently a static approximation that neglects frequencydependent time delays, reflections, and transmission filtering which occur as the wave crosses each fracture. A comparison of the explicit fracture model for fracture-towavelength ratios of 1.7, 3.3, and 6.7 (i.e., 20, 40, 80 and fractures, respectively) shows that as this ratio increases: (1) the amplitudes of multiply reflected waves propagating normal to the fractures and fracture channel waves propagating along the fractures decrease, and (2) the wavefronts develop triplication. This latter phenomena is a characteristic of strong anisotropy. No such triplication are present in the wavefront of the effective medium fracture model. In summary, this study has demonstrated that the dynamic properties of fractures can differ significantly from that predicted by static effective medium approximations, and that the frequency-dependent nature of fractures should not be ignored in the seismic characterization of fractured rock. ACKNOWLEDGMENTS 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-AC0376SF00098. 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Time (sec.) Fig. 6 The seismic waveforms received at Receiver 1 and 2 for the horizontal and vertical component of particle velocity for TI model. 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, J. 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. Pyrak-Nolte, L.J., Myer, L.R. & Cook, N.G.W. (1990a). Transmission of seismic waves across single natural fractures, J. 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, J. 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, J. 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