Microseismic source inversionin anisotropic mediabyW. Scott LeaneyB.Sc., The University of Manitoba, 1983M.Sc., The University of British Columbia, 1987A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Geophysics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2014c© W. Scott Leaney 2014AbstractSedimentary rocks and shales in particular are known to be anisotropic, sometimesstrongly so, and hydraulic fracturing is now common practice in shale plays to en-hance the extraction of hydrocarbons. One of the ways to understand the hydraulicfracturing process is through the micro–earthquakes that it generates; it is thereforeof interest to study the impact that anisotropy may have on hydraulically inducedseismicity. This thesis is concerned with the inclusion of anisotropy into the geo-physical forward and inverse problems of microseismic sources.Ray theory is used for the forward problem — dynamic ray tracing in a mediumcomposed of homogeneous layers with vertical transverse isotropy (VTI) is usedand the possibility of qSv triplication is considered. Novel approaches to the in-verse problem are introduced, including waveform fitting in the frequency domainto recover the source function as well as moment tensor. Uncertainties in estimatedmoment tensor components are quantified with multi–variate normal sampling uti-lizing the full covariance matrix from the linear inversion. A new decomposition ofthe moment tensor is described that removes the distortions due to anisotropy localto the source and a new way to visualize an earthquake source is also introduced.The impact of anisotropy on moment tensor inversion and decomposition is shownto be significant.Field data recorded by two downhole arrays of multicomponent receivers are an-alyzed using the new techniques. The event collection suggests a source mechanismdominated by cylindrical dilatation. This is unexpected but is supported by resultsfrom another code. Future directions include extensions to lower symmetry formsof anisotropy and application to surface arrays. Preliminary analysis of downholerecordings of aftershocks of the 2008 Mw=7.9 Wenchuan earthquake is shown.iiPrefaceSince beginning my PhD in September 2007 my research program has sought theanswers to questions surrounding the impact of anisotropy on earthquake sources,how to model micro–earthquake sources caused by hydraulic stimulation, the use ofray–theory to describe propagation from source to receiver, the inverse problem ofrecovering the moment tensor from vector waveform data recorded downhole andthe decomposition of the moment tensor in anisotropic media to help understand thesource mechanism. Several articles and conference proceedings have been publishedrelated to the work described in this thesis.The fundamental ideas of the frequency domain inversion described in Chapter4 first appeared in the SEG 2008 conference proceedings (Leaney [57]). The analysisof Wenchuan aftershocks using the procedures described in [57] appears at the endof Chapter 7, the abstract for which appears in the 2008 AGU proceedings (Leaneyet al. [69]). The figures presented at that meeting are published here for the firsttime. My co–authors pointed me towards the data set and assisted with the dataformatting and preprocessing.Chapters 2 and 5 contain elements of work published in the EAGE 2010 con-ference proceedings (Leaney and Chapman [59]), while the most important partsof chapters 4 and 6 concerning the inversion and its application to field data haveappeared in a sequence of conference proceedings (Leaney et al. [62], Leaney andChapman [60], Leaney and Chapman [61], Leaney et al. [70], Leaney et al. [64] andLeaney et al. [63]).Chapter 6 is based largely on a publication of which I was not the primary author(Chapman and Leaney [22]), but contains new results which appear here for the firsttime. Parts of other published work of which I was a second author (Michaud andLeaney [79], Mizuno et al. [81]) are described as necessary in chapter 7 in the contextof field data application.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Industry overview: hydraulic stimulation . . . . . . . . . . . . . . . 11.2 Microseismic monitoring basics . . . . . . . . . . . . . . . . . . . . . 51.3 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 The Forward Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 Stress glut source . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2 Particle displacements . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3 Source function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.4 Discussion on other source models . . . . . . . . . . . . . . . . . . . 152.4.1 Angle–dependent corner frequency . . . . . . . . . . . . . . . 152.4.2 Tensile crack P and S corner frequency . . . . . . . . . . . . 162.4.3 Extended sources . . . . . . . . . . . . . . . . . . . . . . . . 162.5 Particle displacements in the frequency domain . . . . . . . . . . . . 172.5.1 A practical approach to Q . . . . . . . . . . . . . . . . . . . 172.5.2 Receiver response . . . . . . . . . . . . . . . . . . . . . . . . 18ivTable of Contents2.5.3 Final forward equation . . . . . . . . . . . . . . . . . . . . . 182.6 Moment tensor construction . . . . . . . . . . . . . . . . . . . . . . 192.6.1 The displacement discontinuity source . . . . . . . . . . . . . 202.6.2 The pressure part . . . . . . . . . . . . . . . . . . . . . . . . 212.6.3 Anisotropic moment tensor . . . . . . . . . . . . . . . . . . . 212.6.4 Angle convention . . . . . . . . . . . . . . . . . . . . . . . . 212.6.5 Constructing D . . . . . . . . . . . . . . . . . . . . . . . . . 222.6.6 Ascribing potency . . . . . . . . . . . . . . . . . . . . . . . . 232.6.7 Scalar moment . . . . . . . . . . . . . . . . . . . . . . . . . . 232.6.8 Source parameter set . . . . . . . . . . . . . . . . . . . . . . 242.6.9 Full potency tensor . . . . . . . . . . . . . . . . . . . . . . . 242.7 Radiation patterns . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.7.1 Radiation pattern computation . . . . . . . . . . . . . . . . 252.7.2 Radiation pattern examples . . . . . . . . . . . . . . . . . . 262.7.3 Basic ISO and TC sources . . . . . . . . . . . . . . . . . . . 272.8 Ray trace synthetics . . . . . . . . . . . . . . . . . . . . . . . . . . . 313 Layered VTI Dynamic Ray–tracing . . . . . . . . . . . . . . . . . . 393.1 Homogeneous, general anisotropy . . . . . . . . . . . . . . . . . . . 393.2 VTI parameterizations . . . . . . . . . . . . . . . . . . . . . . . . . 413.3 Example: Algerian shale, strong VTI anisotropy . . . . . . . . . . . 433.4 Example: Earth inner core anisotropy . . . . . . . . . . . . . . . . . 443.5 Layered VTI models . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.6 Two point VTI ray tracing . . . . . . . . . . . . . . . . . . . . . . . 473.6.1 Phase and group velocity . . . . . . . . . . . . . . . . . . . . 483.7 Ray shooting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.8 Dynamic ray attributes . . . . . . . . . . . . . . . . . . . . . . . . . 493.8.1 Geometrical spreading . . . . . . . . . . . . . . . . . . . . . 493.8.2 Purely vertical and horizontal rays . . . . . . . . . . . . . . . 503.8.3 Effective medium spreading . . . . . . . . . . . . . . . . . . . 503.8.4 Polarizations . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.8.5 Transmission and reflection coefficients . . . . . . . . . . . . 523.8.6 Effective or average Q along a ray . . . . . . . . . . . . . . . 523.9 1D model preparation and unfolding — preraytrc . . . . . . . . . . 533.9.1 Near–field ray–straightening at source and receiver . . . . . 53vTable of Contents3.10 Head wave ray tracing . . . . . . . . . . . . . . . . . . . . . . . . . . 543.10.1 Dynamic properties of head waves . . . . . . . . . . . . . . . 543.11 Triplicating qSv . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 553.11.1 Triplication condition . . . . . . . . . . . . . . . . . . . . . . 553.11.2 Solution selection . . . . . . . . . . . . . . . . . . . . . . . . 563.11.3 qP+qSv wavefronts . . . . . . . . . . . . . . . . . . . . . . . 574 Moment Tensor Inversion . . . . . . . . . . . . . . . . . . . . . . . . 614.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.2 Weighted, damped least–squares inversion . . . . . . . . . . . . . . . 634.2.1 Interpretation of equation (4.4) for this application . . . . . 634.2.2 Estimation and covariance . . . . . . . . . . . . . . . . . . . 644.3 Posterior uncertainties . . . . . . . . . . . . . . . . . . . . . . . . . 644.3.1 Multivariate normal sampling . . . . . . . . . . . . . . . . . 654.4 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.4.1 Estimating the regularization parameter, . . . . . . . . . . 674.4.2 Constrained solutions . . . . . . . . . . . . . . . . . . . . . . 684.5 Estimating s′(t) from m(t) . . . . . . . . . . . . . . . . . . . . . . . 684.6 Estimating mˆ given s′(ω) and m(ω) . . . . . . . . . . . . . . . . . . 694.7 Scalar moment from s(ω) spectral fit . . . . . . . . . . . . . . . . . 704.7.1 Compute times . . . . . . . . . . . . . . . . . . . . . . . . . 714.8 Synthetic examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.8.1 Moment signals and source function recovery . . . . . . . . . 714.8.2 More noise and damping . . . . . . . . . . . . . . . . . . . . 724.8.3 Quality of fit metric . . . . . . . . . . . . . . . . . . . . . . . 724.8.4 Multiple noise realizations . . . . . . . . . . . . . . . . . . . 734.8.5 Multivariate posterior sampling . . . . . . . . . . . . . . . . 764.8.6 A pathological example . . . . . . . . . . . . . . . . . . . . . 764.8.7 Comments on moment tensor uncertainty . . . . . . . . . . . 795 Moment Tensor Decomposition . . . . . . . . . . . . . . . . . . . . . 825.1 Principal decomposition . . . . . . . . . . . . . . . . . . . . . . . . . 825.1.1 Angle extraction . . . . . . . . . . . . . . . . . . . . . . . . . 835.1.2 Solution ambiguity . . . . . . . . . . . . . . . . . . . . . . . 845.1.3 Hudson–Pearce–Rogers parameters . . . . . . . . . . . . . . 84viTable of Contents5.1.4 Vavryc˘uk parameters . . . . . . . . . . . . . . . . . . . . . . 855.1.5 Limitations of conventional decompositions in anisotropy . . 865.2 Anisotropic moment tensor decomposition . . . . . . . . . . . . . . 875.2.1 Biaxial decomposition . . . . . . . . . . . . . . . . . . . . . . 885.2.2 Biaxial potency tensor decomposition . . . . . . . . . . . . . 885.2.3 Source geometry parameters . . . . . . . . . . . . . . . . . . 895.2.4 Biaxial potency parameterization . . . . . . . . . . . . . . . 905.2.5 Normalized EOS . . . . . . . . . . . . . . . . . . . . . . . . 905.2.6 Impact of anisotropy . . . . . . . . . . . . . . . . . . . . . . 915.2.7 Special sources . . . . . . . . . . . . . . . . . . . . . . . . . . 925.3 Decompositions of noisy M . . . . . . . . . . . . . . . . . . . . . . . 945.4 Source visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . 965.4.1 Hockey pucks . . . . . . . . . . . . . . . . . . . . . . . . . . 966 Field Data Application . . . . . . . . . . . . . . . . . . . . . . . . . . 986.1 Stimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 986.2 Monitoring geometry and hardware . . . . . . . . . . . . . . . . . . 1026.3 Model building . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1096.3.1 Sonic logs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1096.3.2 Log blocking . . . . . . . . . . . . . . . . . . . . . . . . . . . 1096.4 VTI model calibration . . . . . . . . . . . . . . . . . . . . . . . . . . 1126.4.1 VTI calibration algorithm . . . . . . . . . . . . . . . . . . . 1136.4.2 Event locations . . . . . . . . . . . . . . . . . . . . . . . . . 1156.5 Residual statics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1166.5.1 Iterative, complex correlation statics . . . . . . . . . . . . . 1166.6 Receiver statics example . . . . . . . . . . . . . . . . . . . . . . . . 1176.7 Q model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1176.8 Moment tensor inversion example . . . . . . . . . . . . . . . . . . . 1176.9 Multi–event results . . . . . . . . . . . . . . . . . . . . . . . . . . . 1276.9.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1347 Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1387.1 Inversion improvements . . . . . . . . . . . . . . . . . . . . . . . . . 1387.1.1 Augmented arrivals model . . . . . . . . . . . . . . . . . . . 1387.1.2 Green function errors . . . . . . . . . . . . . . . . . . . . . . 138viiTable of Contents7.1.3 Other inversions . . . . . . . . . . . . . . . . . . . . . . . . . 1397.2 Extensions to lower symmetry . . . . . . . . . . . . . . . . . . . . . 1407.2.1 Salt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1437.3 Extended source model . . . . . . . . . . . . . . . . . . . . . . . . . 1447.4 Waveform fit model calibration . . . . . . . . . . . . . . . . . . . . . 1447.5 Surface array and shallow grid geometry . . . . . . . . . . . . . . . 1457.6 Wenchuan aftershocks . . . . . . . . . . . . . . . . . . . . . . . . . . 1478 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155viiiList of Figures1.1 Well site of first hydraulic fracturing job, Kansas, 1947. Source:Schlumberger corporate marketing, used with permission. . . . . . . 21.2 Distribution of resource plays in North America. Figure copyrightSchlumberger, used with permission. . . . . . . . . . . . . . . . . . . 41.3 Distribution of shale gas plays throughout the world, circa 2011. Fig-ure copyright Schlumberger, used with permission. . . . . . . . . . . 51.4 3D microseismic event distribution for a multi-stage frac in a horizon-tal well. A subset of this field data set will be analyzed in Chapter6. Data used with permission from Lightstream Resources. . . . . . 62.1 Time (left) and frequency domain representations of far–field particledisplacement source time functions. The n = 1 or Brune pulse isshown in red, the n = 2 pulse in green. . . . . . . . . . . . . . . . . . 142.2 Far–field particle velocity pulses corresponding to the far–field dis-placement pulses shown in 2.1 . . . . . . . . . . . . . . . . . . . . . . 152.3 Vector and angle conventions. Left: Polar coordinate angles show-ing azimuth increasing clockwise from north and dip increasing fromvertical–up. Right: strike (φ), dip (δ) and rake (λ) angles with slipvector [s] and fracture plane normal vector nˆ. The ENU directionsare represented by the unit vectors ıˆ, ˆ, kˆ. aˆ is the strike directionvector and bˆ is the dip direction vector. . . . . . . . . . . . . . . . . 222.4 PHV radiation at 60◦ phase angle for a pure slip source with N–Sstrike slip on a vertical plane in an isotropic medium (solid) and VTImedium (dashed). Color coding follows PHV =BRG=blue,red,green).Left: un–normalized moment tensors, right: normalized moment ten-sors. Normalization serves to highlight the effect of anisotropy on theray strain tensor part of the radiation pattern. The absolute valueof radiated amplitude is shown in these figures (i.e. polarity is notshown). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26ixList of Figures2.5 Equal area complete focal sphere radiation patterns for qP , Sh andqSv rays (PHV are plotted top to bottom). A pure slip source withSDR = (60◦, 40◦, 20◦) and equal potency A[d] has been placed inan isotropic medium (left) and a VTI medium (right). The VTImedium has Schoenberg parameters (0.25, 0.20, 0.30) (see eq. 3.9).The upper–hemisphere is contained within the black circle, whichrepresents the horizontal plane. In these and similar figures red islarge, blue is small and green is mid–way between. These plots areshown in true relative amplitude to highlight the larger radiated am-plitudes in the generally stiffer anisotropic medium. . . . . . . . . . . 282.6 Radiation patterns for pure slip (left) and pure tensile source withpositive opening angle (right). In both cases the medium is VTI.Each plot is normalized individually. . . . . . . . . . . . . . . . . . . 292.7 VTI P radiation as viewed from the south for a pure slip source (red)and a pure tensile source with positive opening angle (blue). Thezero line is shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302.8 Equal area complete focal sphere radiation pattern (left) for P rays.A composite source with slip, tensile opening and negative expansionis used. The source parameter set is (−2, 0.20, 60◦, 40◦, 20◦, 45◦). Onthe right is shown the radiated amplitude from the south with a zeroline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.9 Three array geometry in map(top) and section view used for simula-tion. Two vertical arrays and one horizontal array surround a sourceat a distance of about 565m. All arrays have 12 3C receivers spacedat 30m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322.10 Ray theory particle velocity synthetics in a homogeneous VTI modelwith DD source parameters SDRO = (60◦, 40◦, 20◦, 45◦), n=2 fc =100Hz source wavelet, geometrical spreading and additive Gaussiannoise with standard deviation equal to maxamp/100. The receiverarray consists of two vertical arrays and one horizontal array, eachwith 12 3C receivers. From left to right are groups of 36 traces forthe East, North and Up (ENU ) components. The three groups of 12traces can be seen as the different arrays. . . . . . . . . . . . . . . . 332.11 Same as figure 2.10 except including absorption and dispersion due toQ. Q = 100 was used, chosen as equal for PHV rays. fref = 200Hzwas used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34xList of Figures2.12 The simulated event recording from figure 2.10 rotated to scalar PHV(qP, Sh, qSv) arrivals based on the ray–trace receiver polarizations(top). Bottom: PHV radiation patterns with receiver locations show-ing focal sphere sampling. An approximate homogeneous VTI Greenfunction is included in the radiation patterns for direct comparisonwith waveform amplitudes. Radiation patterns are displayed in truerelative amplitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.13 The maximum PHV envelope function applied to the waveforms offigure 2.12 with arrival times shown according to the color codingPHV =BRG=blue,red,green). . . . . . . . . . . . . . . . . . . . . . . 373.1 A ray path and wavefront with slowness (p), group velocity (V) andpolarization (gˆ) vectors where for clarity the polarization vector isshown for a quasi–shear ray. The distance between wavefronts dTapart in time is |V|dT in the ray direction V, and dT/|p| in theslowness direction p. . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.2 Thomsen’s (red) and Schoenberg’s ellipticity parameter Ep (blue)versus a velocity ratio measure√A11/A33 − 1. . . . . . . . . . . . . 433.3 Phase slowness (left) and group velocity (right) for a strongly anisotropicshale. qP=blue, qSv=green and Sh=red. . . . . . . . . . . . . . . . 443.4 Phase slowness (left) and group velocity (right) polar plots for the TIanisotropy of the inner core. qP=blue, qSv=green and Sh=red. . . . 453.5 qSv horizontal phase slowness versus group angle for the triplicatingshale of figure 3.3 (left) and X vs p. . . . . . . . . . . . . . . . . . . 563.6 A 150 Hz Brune pulse and its Hilbert transform (dashed). . . . . . . 583.7 qP and qSv wavefronts at 100 ms for the triplicating shale of figure3.3 including geometrical spreading using a causal Brune pulse withcorner frequency fc = 100 Hz. Red is high amplitude. Left: with-out Hilbert transform on the reverted branch; right: with Hilberttransform on the reverted branch. . . . . . . . . . . . . . . . . . . . . 593.8 Left: qP and qSv wavefronts at 100ms for the triplicating shale offigure 3.3 using a causal Brune pulse with corner frequency fc = 100hzand using the maximum time selection strategy. Right: Sh wavefrontsshown using the same color amplitude scale but different horizontaland vertical scales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.9 Layered model, rays, times and qP + qSv wavefronts for a source inan isotropic layer embedded in the triplicating medium. . . . . . . . 603.10 Layered model, PHV rays, PHV times and qP + qSv wavefronts at100ms for a layered VTI model used in commercial processing. . . . 60xiList of Figures4.1 Three array synthetics with noise as in figure 2.10 with from leftto right groups of 36 input traces (3 arrays of 12) for ENU compo-nents followed by the data reconstructed after inversion to determinemoment tensor and source function. Colored arrival times correspon-dence is PHV=BRG. . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.2 Left: Estimated moment signals mij and source function s′(t) as(upper left) scaled particle velocity, (lower left) scaled particle dis-placements; right: source function amplitude spectrum |s(ω)| andspectrum fit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.3 Same as figure 4.1 but with five times as much noise. Top: no damp-ing, bottom: damping parameter from GCV, estimated at 0.27. . . . 744.4 VR (red) and arrival–weighted VR (green) fit metrics for 100 increas-ing noise realizations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 754.5 Moment tensor components inverted for 100 increasing noise realiza-tions (left): red=m11, green=m22, blue=m33, cyan=m23, magenta=m13,black=m12. Right: A relative error metric. . . . . . . . . . . . . . . 754.6 5000 multi–variate normal realizations of moment tensor componentshonouring the posterior covariance matrix. red=m11, green=m22,blue=m33, cyan=m23, magenta=m13, black=m12. . . . . . . . . . . . 774.7 5000 multi–variate normal realizations of moment tensor componentsm22 vs. m11 honouring the full posterior covariance matrix (red) andhonouring only the diagonal terms (green). . . . . . . . . . . . . . . 784.8 Model covariance matrices for the 3–array synthetic (top) and the 2–array synthetic with pathological event location (bottom), normalizedby the average of the two intermediate diagonal elements. . . . . . . 794.9 Multivariate normal histograms from 5000 samples of moment ten-sor components honouring the posterior covariance matrix for the3–array (top) and the pathological 2–array case (bottom). red=m11,green=m22, blue=m33, cyan=m23, magenta=m13, black=m12. . . . . 805.1 The basic structure of a Hudson (k vs. τ) or ISO vs. CLVD di-amond crossplot. The positions of various special sources are indi-cated. TC=tensile crack, LVD=linear vector dipole, CD=cylindricaldilation, IP=isotropic plane. . . . . . . . . . . . . . . . . . . . . . . . 855.2 Normalized moment tensors and their eigen decomposition (eigenval-ues are shown above the TNP eigenvectors) for a slip source withSDR angles (5◦, 30◦, 30◦). Results for the same isotropic and VTImedia as for figure 2.4 are shown. . . . . . . . . . . . . . . . . . . . . 86xiiList of Figures5.3 Left: radiation pattern for a composite source with equal parts in–plane and normal displacement. The fault plane strikes N5E with 30degrees dip to the west. In–plane slip has rake= 30◦. Right: Theeigen analysis of the moment tensor, M, and of the source tensor, D. 875.4 Normalized source parameter crossplots for DD sources with openingangle χ ranging from −90◦ to +90◦. Left: isotropic medium; right:VTI medium. SDR angles are (60◦, 40◦, 20◦). Green: Hudson k vs.τ ; blue: Vavryc˘uk ISO vs. CLV D; red: Ê vs. Ô. . . . . . . . . . . . 925.5 The opening angle χ was increased from−90◦ to +90◦ in the construc-tion of M in the anisotropic medium. On the left are shown SDROangles recovered using the biaxes decomposition of DDD while onthe right they are shown obtained from M conventionally assumingisotropy. S=red, D=green, R=blue, O=cyan. . . . . . . . . . . . . . 925.6 Histograms for multivariate sampled moment tensor elements (upperleft) and histograms of normalized source parameters. Hudson (upperright), Vavryc˘uk (lower right) and Ê, Ô, Ŝ (lower left). . . . . . . . . 955.7 Decomposition parameters on a diamond crossplot for the multi–variate normal posterior samples of M using full covariance and diagonal–only (on top). Hudson (green and magenta); Vavryc˘uk (blue andyellow), Ê, Ô, Ŝ (red and cyan). The true model plots at (Ô, Ê) =(0.5, 0). On the right are SDRO angles from EOS (solid) and fromM assuming isotropy. Vertical lines show true angles. . . . . . . . . 965.8 New 3D graphical representation of a composite source showing ori-ented EOS effective volumes or potencies. Left: map view, right:side view from south. This synthetic example has significant −E =−Expansion (blue wire frame) and an opening angle χ = 45◦ tosplit A[d] into +O = Opening (thickness of red pucks) and S = Slip(pucks displacement). Strike, dip and rake parameters are 60◦, 70◦and 10◦, respectively. Of the two solutions, the one closest to a priororientation (in this case the true solution) has been chosen for display. 976.1 Graphical depictions of plug–and–perf (top), sliding sleeve completion(middle) and a schematic of the ball–seat operation to open frac portsfor sliding sleeves (bottom). . . . . . . . . . . . . . . . . . . . . . . . 100xiiiList of Figures6.2 Pumping data and microseismic event counts for stage 8. The slurry(fluid) rate is in blue with scale on the right as Flowrate in m3/min.;top rig pressure is in red with scale on the left in MPa; proppantconcentration is in green with scale on the right in kg/m3; event rateis in double green bars with scale on the left. Bars indicate eventcount in 2 minute windows. The early build–up of pressure (red)followed by a sudden decline indicates the sliding of the sleeve andopening of the port to the formation. . . . . . . . . . . . . . . . . . . 1016.3 Map view of the wells and event locations color-coded by stage. Thegrids are 100m on a side. . . . . . . . . . . . . . . . . . . . . . . . . . 1036.4 Side view from the southwest the wells and event locations color-coded by stage. The grids are 100m on a side. . . . . . . . . . . . . . 1046.5 Upper left: schematics of Schlumberger’s VSI downhole 3C shuttleshowing the detached sensor package design. Lower left: a photo ofa VSI shuttle with anchor arm extended and sensor package. Upperright: VSI-4 dimensions including telemetry cartridge and individualshuttle; Lower right: a VSI-40 array in the shop with wireline inter-connects, getting ready for deep offshore Gulf of Mexico deployment.Thanks to Bill Borland for the graphics. . . . . . . . . . . . . . . . . 1056.6 Photos of rigging up a VSI array for hydraulic fracture monitoringoperations (top) and a screen grab of the acquisition software for aVSI-12 (see tool lay-out on right). Thanks again to Bill Borland forthe images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1066.7 Oriented ENU waveforms after merging the data from the two mon-itoring arrays for a large detected and located event (left) and theamplitude spectra (right) shown in dB (trace-normalized) with themedian spectrum in black. . . . . . . . . . . . . . . . . . . . . . . . . 1076.8 Ray-rotated PHV waveforms for the large microseism event of figure6.7 and for a sliding sleeve event. This microseism has roughly tentimes the amplitude of the sleeve event. The thin coloured linesrepresent the PHV arrival times (PHV =Blue,Red,Green) using thecalibrated VTI model and event locations. . . . . . . . . . . . . . . . 1086.9 Elastic logs used for initial model building. The empirically deriveddensity is in units of gm/cc, Vp and Vs are in m/s (note the factor102). The reservoir is located near the minimum of the Vp/Vs ratiolog at about 840m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110xivList of Figures6.10 Left: Vp log and extended log using conjugate extrapolation to handleend effects prior to smoothing. This approach preserves the end-zonegradient. Right: Vp log and three layered models corresponding todifferent averaging criteria. The green coloured model correspondingto the transmission smoothing criterion is carried forward for VTImodel calibration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1126.11 Schlumberger’s microseismic VTI model calibration software interfaceshowing the data selection window (upper left), the picking window(lower right), the inversion parameter window (Upper right) and themodel window (lower left). . . . . . . . . . . . . . . . . . . . . . . . . 1146.12 Calibrated VTI model (left) and geometry with P rays from 9 cali-bration shots (right) together with color–coded Vp. Depth–dependentanisotropy follows a soft rock physics constraint. This layered, VTImodel reproduces PHV times picked on 9 sleeve “shot” waveformswith an rms residual of 3.3ms. . . . . . . . . . . . . . . . . . . . . . . 1146.13 Vp, Vp/Vs and Thomsen anisotropy parameters for the calibrated VTImodel (left), PHV rays and arrival times to 12 receivers located 300mfrom the event, which is located at the base of the reservoir. . . . . . 1156.14 Calibrated VTI model time-flattened PHV arrivals (left) and afterresidual statics corrections (right) using the complex correlation al-gorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1176.15 Top: Input ENU data (left three panels) and reconstructed ENUdata after moment tensor inversion and forward modelling for eventreconstruction. The arrival-weighted normalized energy fit is .664.Bottom: The same as top except the waveforms have been alignedbased on the P-arrival time (including statics). . . . . . . . . . . . . 1196.16 Left: Inverted moment signals and estimated source function; right:source function spectrum and spectral fit to extract scalar momentand corner frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1206.17 Radiation patterns including homogeneous VTI Green function forthe inverted moment tensor with receiver array focal sphere sampling.Also shown is the puck–glyph viewed from above to visualize thedecomposed source. The green arrow points North. . . . . . . . . . . 1216.18 Event #36 from the 182 event collection of stage 7, shown after ray–rotation to scalar PHV amplitudes and after time–alignment basedon calibrated VTI model times and complex correlation statics. . . . 1226.19 Event #36 from the 182 event collection of stage 7, shown in inputENU coordinates (left) in the pass–band 23-290Hz and reconstructedafter inversion (right three panels). Bottom: the same as top buttime–aligned on the direct P times including statics. . . . . . . . . . 123xvList of Figures6.20 Event #36, PHV with inverted radiation patterns shown upper-hemisphere equal–area with relative receiver locations indicated show-ing the sampling of the focal sphere. . . . . . . . . . . . . . . . . . . 1246.21 Event #36 estimated moment signals and extracted source func-tion (left) where the ordering of moment signals is now, from topto bottom: m11(t),m22(t),m33(t),m23(t),m13(t),m12(t), s(t). Right:source function amplitude spectrum and spectral fit for scalar mo-ment and corner frequency. . . . . . . . . . . . . . . . . . . . . . . . 1256.22 Event #36 3D glyph using the EOS decomposition. The slip termindicates right–lateral slip on a sub–vertical plane with a signifi-cant closing (−O component) and positive pressure change (red wireframe). The poles of the red wire frame are oriented along the T–axis. 1256.23 Event #36 normalized covariance matrix. Ordering from upper leftis m11,m22,m33,m23,m13,m12. . . . . . . . . . . . . . . . . . . . . . 1256.24 Event #36 multivariate normal sampling of M (upper left; as his-tograms (lower left); and in diamond crossplot space (right). . . . . . 1266.25 Stage 7 event locations in depth versus Northing with logs of Vp,Vp/Vs, model and receivers. The well separation is 884m. . . . . . . 1286.26 Stage 7 event locations in map view (1km x 1km) along with receiverlocations (red dots). The magenta event is event # 36 studied in theprevious section. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1296.27 Condition number versus event number for the stage 7 event collection(left); estimated regularization parameter (right). . . . . . . . . . . . 1296.28 Inverted source function amplitude spectra for 182 events from stage7 (left) and the spectral fits using a ω−3 decay model. . . . . . . . . 1306.29 Source function spectral fit parameters. Left: corner frequency fromthe ω−3 decay model; Right: moment magnitude Mw derived fromscalar moment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1306.30 Fit quality and parameter uncertainty metrics versus event number.Left: conventional energy fit (red) and arrival-weighted energy fit(green); maximum relative moment tensor error (right). . . . . . . . 1316.31 Decomposition parameters versus event number. Left: Hudson k, τ, δ(red, green, blue); right: normalized EOS parameters (EOS=Red,Green,Blue).1316.32 Normalized parameter crossplot for the 182 events of stage 7 showingHudson (green), Vavryc˘uk (blue) and normalized E vs. O (red). Left:optimized damping, right: intermediate damping. . . . . . . . . . . . 1326.33 Normalized parameter crossplots as in figure 6.32 for a reduced (left)and an increased damping parameter. . . . . . . . . . . . . . . . . . 133xviList of Figures6.34 SDRO (Strike, co-dip, rake, opening) angles from the DD part of thepotency tensor coming from the EOS decomposition. Color corre-spondence is SDRO=Red,Green,Blue,Cyan. . . . . . . . . . . . . . . 1336.35 3D display from Petrel of the inversion results shown using the EOSdecomposition and “pucks” glyph. The view is from the south–westand above. Also shown are the two vertical monitor wells with somereceiver locations (red and green) and the lateral treatment well witha yellow icon at the location of the stage 7 port. . . . . . . . . . . . 1346.36 Hudson parameters shown for an optimized regularization parametermulti–event run. Crossplot locations for special sources are shown(for a Poisson medium). . . . . . . . . . . . . . . . . . . . . . . . . . 1356.37 Hudson parameters shown color-coded by time and scaled by VR (aquality fit metric) for two different moment tensor inversion codes.The display was made with Schlumberger’s Petrel package. Left:the frequency domain anisotropic code described in this thesis; right:Schlumberger’s commercial time–domain code. . . . . . . . . . . . . 1356.38 Sleeve event ENU data (first three panels) and inversion (right threepanels) for different source models. The simpler force source model(three parameters) inversion (left) compared to moment tensor inver-sion (six parameters, right) fits the data significantly better. . . . . . 1376.39 VR fit quality versus event number for a collection of the largestevents from stages 4-13. The FSI model (green) with three fewerdegrees of freedom is able to fit the waveforms better than MTI (red)for four events (circled) known with confidence to be sleeve events. . 1377.1 qP , qSh and qSv group velocities versus group inclination (0 − 90◦)and azimuth shown in an equal area projection with vertical–up inthe center. Top: qP , lower left: qSh, lower right: qSv. qP and qSare normalized separately; shear velocities are plotted with the samecolor scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1427.2 Halite qP group velocities versus group inclination (0− 90◦) and az-imuth shown in an equal area projection with vertical–up in the center(left). Right: qP group velocities versus group inclination angle forall azimuths and those for the closest VTI model (black line). . . . . 1447.3 Top: P and S velocities (left) and rays (right) for a source at depthand surface receivers out to 10,000ft offset. Bottom: Simulated Pamplitudes for a double–couple source as recorded by a dense gridof receivers. Left: scalar amplitude for ray take–off angles; right:amplitudes including the effects of spreading, transmission loss, Q(100, for dominant frequency 150Hz), and projection onto the verticalcomponent. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146xviiList of Figures7.4 The location of the magnitude 7.9 May, 12th, 2008 earthquake andaftershocks slip map with the location of the Schlumberger downholearray from 5-7 June, 2008. . . . . . . . . . . . . . . . . . . . . . . . . 1487.5 Oriented ENU components for the largest event during stage 4 record-ing. The horizontal axis is seconds. The boxes simply highlight firstP and S phases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1487.6 Amplitude spectra for 66 minutes of recording during stage 4 of thefailed treatment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1497.7 1D elastic model and P and S head wave rays used for processing. . 1507.8 Input ENU waveforms, reconstructed ENU waveforms, residual ENUwaveforms and separated scalar PHV waveforms over the P arrivalwindow (top) and the S arrival window (bottom). Vertical axis is inseconds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151xviiiList of AcronymsSome commonly used acronyms are included here for quick reference.ENU : (East, North, Up) coordinate systemPHV : qP , Sh, qSv direct arrivals, or waveform data rotated from ENU basedon ray–polarizations at receivers.SDRO:S = Strike, measured clockwise from North;D = co-dip = 90◦ − dip where dip is measured from horizontal, thus D ismeasured from vertical and is 0 for a vertical fracture plane;R = Rake, defined counter–clockwise from horizontal as the motion of the hang-ing wall relative to the footwall;O = Opening angle, meaning the angle the displacement vector makes with thefault or fracture plane — the amount of tensile component.EOS: Expansion, Opening and Slip, components of the new moment tensordecomposition for anisotropic media.VSI: Schlumberger’s Versatile Seismic Imager tool.VTI: Vertical Transverse Isotropy, also referred to as TIV, meaning polaranisotropy such as that produced by a stack of thin horizontal layers with vary-ing elastic constants.MTI: Moment Tensor InversionFSI: Force Source InversionxixAcknowledgementsIn addition to my supervisor, Tad Ulrych, and my committee, I need to acknowledgeChris Chapman, Schlumberger Scientific Advisor and recipient of the Royal Astro-nomical Society Gold Medal (Geophysics). Personal communications with Chrishave been almost daily and are too numerous to cite, only a select few have beenincluded. I have developed computer code based largely on the theory developedthrough our collaboration. I would also like to acknowledge several colleagues inSchlumberger: Doug Miller (former Schlumberger), the late Mike Schoenberg, Cen-giz Esmersoy, Takashi Mizuno, Xin Yu, Les Bennett, Jim Rutledge, Shawn Maxwell(former Schlumberger), Bill Borland, Gwe`nola Michaud, Nobuyasu Hirabayashi, LesNutt, Henry Menkiti and Andrew Fryer.I also thank the exam committee: Michael Bostock, Andrew Calvert (SimonFraser University), Ronald Clowes, Matthew Yedlin and external examiner Mirkovan der Baan (U. of Alberta) who greatly improved the readability of this thesis. Aspecial thanks to Michael who took over so kindly as “surrogate supervisor” whenTad took ill.I gratefully acknowledge John Duhault and Lightstream Resources for permissionto use and show the field data.xxDedicationFor my mother, Violet, and my late father, Donald.xxiChapter 1Introduction1.1 Industry overview: hydraulic stimulationOld rocks tend to be hard, fast and tight in that their porosity is low and what isin the pore space doesn’t flow easily. Continental rocks tend to be old, while thoseoffshore, for example under continental shelves, are younger, less compacted andmore porous. A consequence of continental rocks being older and tighter is that it ismore difficult to extract hydrocarbons from the pore space – wells drilled into stratacontaining natural gas in the pore space, for example, must be stimulated in someway to encourage the hydrocarbons to flow to the well. One of the oldest techniquesto encourage the flow of hydrocarbons to wells is that of hydraulic stimulation.Hydraulic stimulation of wells is a technique that dates back to 1947 (Figure 1.1).The idea is simple: pump water (slurry) into the well under high enough pressure tofracture the formation at depth, thereby creating new pathways for the hydrocarbonsto flow to the well. Since the content of the pore space at depth is under pressure,when the pumps are turned off and the well is allowed to flow back, the increasedpathways to the well lead to an increase in hydrocarbon flow rates. Increases canbe significant, exceeding 100% and lasting for years. The process may even berepeated. For an extensive review of hydraulic fracturing, including environmentalimpact considerations, see King [53].The size of the hydraulic stimulation industry is surprising. As of 2009 an esti-mated 20 billion dollars per year is spent in North America on hydraulic stimulation,distributed over roughly 60,000 wells (numbers are based on Schlumberger’s num-bers and estimated market share), largely to tap into what’s called “tight gas” or“unconventionals”. Since work on this thesis was started the number of hydrauli-cally stimulated wells has grown to 89,000 and the focus has moved from gas to oil,with the Bakken formation centered in North Dakota leading the way.The advent of horizontal drilling over the last twenty years or so has furtherexpanded the tight gas market, since with a horizontal well and hydraulic fracturingcarried out at stages along its reach, a single well can drain a much larger volumeof a single formation.For schematics showing what a hydrofrac well–site entails, the interested readermay look at web sites such as safewatermovement.org, and www.neb-one.gc.ca.11.1. Industry overview: hydraulic stimulationFigure 1.1: Well site of first hydraulic fracturing job, Kansas, 1947. Source: Schlum-berger corporate marketing, used with permission.The areal requirements for the equipment can be large. Often a temporaryreservoir is constructed for water, sometimes with convoys of trucks to fill it. Thereare also containers for fluid additives, pumps, fuel supply, worker trailers etc. — abig investment. The environmental impact has led to controversy, not only for thesurface environmental impact but for concerns about groundwater contaminationand induced earthquake risk. EPA investigations are ongoing. But the industrygrows so evidently the large investments are paying off.Historically the types of reservoirs thought amenable to hydraulic stimulationwere tight sandstones, but in the last few years it has been discovered that forma-tions associated with shales, previously thought to be source rocks or seals but notreservoirs, are amenable to hydraulic stimulation or “hydraulic fracturing”. Theformations referred at as shale reservoirs are actually not shales at all and neednot even contain any clay minerals (Prasad, 2013, CSEG Recorder [84]). Some arecarbonate muds while others are siliciclastic. The industry thus suffers from aninexact terminology and alternate terms have been used such as “unconventionals”,“self-resourcing rocks”, “organic-rich rocks”, “mudstones”, etc., but the term shaleappears to have stuck. These “shales” tend to be low porosity (3–8%) and can beexceedingly tight, having permeabilities on the order of micro–darcies, sometimesmany times less permeable than old, tight sandstones. Importantly, the unconven-tional reservoirs may be anisotropic, and even if they have low clay content they are21.1. Industry overview: hydraulic stimulationcommonly encased in true shales which are usually strongly anisotropic. Anisotropyis a complicating factor and is an important aspect of this thesis.The discovery that some shales may be considered hydrocarbon reservoirs, com-bined with the technologies associated with hydraulic fracturing and horizontaldrilling, has lead to a boom in the North American oil and gas industry. Lead-ing the boom in the North American gas-shale market was the Barnett shale ofnorth–east Texas, essentially centered under the Dallas–Ft. Worth (DFW) airport.Subsequently several gas-shale plays have been identified and have been exploredand exploited. Figure 1.2 shows the main resource plays in North America andfigure 1.3 shows shale plays throughout the world, circa 2011. For up–to–datemaps one may visit the U.S. Energy Information Administration: www.eia.gov, orhttp://shalewatchblog.com/. Shale-oil or “unconventional liquids” plays, exempli-fied by the Bakken shale centred under western North Dakota are now hot properties.In spite of the success of hydraulic fracturing it is generally regarded as beingpoorly understood. For years the only information available at the surface was thewell–head pressure and pump rates. Geomechanical models were used to predictfracture growth, but the evolution and extent of fracturing remained difficult toquantify.During a frac job, pumps increase fluid (slurry) pressure until breakthroughoccurs and a fracture begins to propagate away from the wellbore. Top hole pressuresare known but bottom hole pressures are usually calculated, not measured, and theyare preferred for analysis. Slurry rates can exceed 12m3 or about 100 barrels perminute. The evolution of pressure and slurry rate is monitored to decide when toadd proppants, essentially sand grains, to the slurry to hold open the fractures whenthe pressure is reduced.After the pumps are turned off and the well is allowed to flow back, hydrocarbonproduction is monitored to see how successful the hydraulic stimulation has been.From the pressure and production data alone it is often difficult to determine whyone stimulation might be superior to another. Sometimes production actually getsworse, with water being produced after stimulation. Why? If production logs in-dicate that certain intervals are producing more than others, is this due to spatialdifferences in fracture penetration? If a map or image of fracturing were availablecould problems be avoided and stimulation optimized?Many of the techniques being developed to answer these and other impor-tant questions are classified as techniques for “Hydraulic Fracture Monitoring” or“HFM”, and one of the most popular techniques involves monitoring hydraulic frac-ture stimulations by the seismic method. This involves recording seismic signalswith a network of seismometers or geophones located downhole, in shallow well-bores or on the surface. Downhole arrays may be deployed temporarily in monitorwells or they may be permanently cemented in place. The analysis of such wave-form recordings to recover the information about the nature of the source is the31.1. Industry overview: hydraulic stimulationFigure 1.2: Distribution of resource plays in North America. Figure copyrightSchlumberger, used with permission.general focus of this thesis. The impact of anisotropy on waveform recordings andanalysis techniques, and recovering the source mechanism from the directional vari-ation in radiated amplitude are specifically the goal. The link of the inferred sourcemechanism to the evolving fracture network is the problem of interpretation, anda very difficult one tied up in fracture geomechanics (see, for example, Fisher andGuest [42]). The goal here is to improve the geophysics side of the equation so thatinterpretations may be better informed and more reliable.41.2. Microseismic monitoring basicsFigure 1.3: Distribution of shale gas plays throughout the world, circa 2011. Figurecopyright Schlumberger, used with permission.1.2 Microseismic monitoring basicsWhen fractures are created by hydraulic stimulation microearthquakes can occur,potentially revealing information about the fracturing process. Naturally occurringearthquakes that can be felt by humans are generally considered to have a magnitudeof at leastM = 3. Microearthquakes or microseisms related directly to newly createdfractures usually fall in the magnitude range from −3 to 0 while events in the rangefrom 0 to 3 are usually related to re-activation of movement on pre-existing faults.Maxwell (2014) [74] has used the analogy for the size of microseisms as from acan of soda to a carton of milk hitting the floor from counter height, so relativelyspeaking many of these events truly are microearthquake events. Typical depthswhere reservoirs are hydraulically stimulated are from 800m to 3500m below groundlevel. At these depths such small events can be, not surprisingly, difficult to recordon the surface, although some regions produce much larger events than others andcommercial success has been reported (e.g. Snelling et al. [102], microseismic.com).The most common type of microseismic monitoring is done with downhole seis-mometers, a variant on the borehole seismic (BHS) technique. In this approach amonitor well, usually a well previously used for production and hence cased withsteel, is dedicated to the borehole seismic tool. The borehole seismic tool comprisesan array of “shuttles” containing three component (3C) geophones clamped to theborehole wall either mechanically or electromagnetically. Current commercial arraysvary from eight shuttles to forty or more, connected by standard heptacable wire-51.2. Microseismic monitoring basicsFigure 1.4: 3D microseismic event distribution for a multi-stage frac in a horizontalwell. A subset of this field data set will be analyzed in Chapter 6. Data used withpermission from Lightstream Resources.line or optical fiber cable. Recorded signals are digitized downhole by a cartridgeand transmitted uphole by telemetry. Telemetry rates allow continuous recordingat temporal sampling rates of 0.25 or 0.50 milliseconds, depending primarily on thenumber of shuttles and the cable. Shuttle spacing is usually 15m (50ft) or 30m(100ft), a compromise between total array aperture and cost. Very long arrays pro-vide diminishing returns because of loss of signal for the shallowest receivers; moredense arrays are preferred but costly. Monitoring from horizontal wells has beendone but is relatively rare as it is more costly and poses operational and processingchallenges; monitoring from two or more vertical or sub-vertical wells simultaneouslyis more common.During the course of a borehole microseismic monitoring survey there is a stan-dard sequence of events. First the BHS tool is deployed in the monitor well(s). Thisis done first in order to record the perforation shots that are fired by specially de-signed perforation guns at the desired depth(s) to put holes through the steel casing,the annular cement and into the formation. These holes will provide the conduitfor the fluid. The “perf shots” are needed for seismic monitoring as they providea relatively large source at a known location for calibrations. Perf shots are usedboth to orient the 3C geophones which are free to rotate during deployment and tocalibrate a velocity model built from sonic logs. The geophone orientations are de-termined using polarization vectors determined from picked P wave arrivals. Thesepicked P and S times from all available perf shots may be used in a constrainedinversion for a transversely isotropic velocity model (e.g. Mizuno et al. [81]). This61.3. Thesis outlinemodel will produce perf shot locations as close as possible to the known locationsgiven uncertainties in the time picks and starting model. Given the required vectorgeophone rotations and a calibrated velocity model, microseismic monitoring maybegin. Pumps are turned on and microseismic events are detected and located us-ing an automated algorithm (e.g. Michaud and Leaney [79]) within 30 seconds oftheir time of occurrence. These event locations may be viewed at the well site ortransmitted to a remote location (i.e. a company office) and viewed in a 3D inte-grated interpretation environment in virtual near real time. Figure 1.4 shows finalevent locations for a multi-stage lateral well treatment, color-coded by stage. (Asubset of this field data set will be analyzed using the new techniques described inthis thesis.) Based on the evolving cloud of microseisms, decisions can be taken tomodify the stimulation program. If events are seen to occur out of zone this may in-dicate unwanted fracturing of an adjacent formation and possible future productionof water from that formation. Similarly, events occurring distant from the activestage perforations may indicate the need for remedial action. There is an increasingvariety of stimulation technologies available in the industry to localize and controlhydraulic fracturing, and microseismic monitoring, always used in conjunction withpressure data, is seen as an important technology to control stimulation processesand decisions.As impressive as some of the commercial successes of hydraulic stimulation havebeen, there remain many unanswered questions and there is a wealth of informationin microseismic data that is not being tapped — improvements are needed. Centralto everything to do with hydraulic stimulation is the geological setting and theproperties of the reservoir rocks and those surrounding the reservoirs. Since it iswell known that shales are anisotropic the question of the impact of anisotropy onthe microseismic technique is an important one. The influence of anisotropy on thecorrect recovery of microseismic source mechanisms is the main topic of this thesis.1.3 Thesis outlineIn chapter 2 the forward problem is described, starting with the expression for amoment tensor in terms of the time and volume integral of the stress–glut ratedue to Backus and Mulcahy (1976), [7], [8]. The anisotropic ray theory formulafor time–domain particle displacements in the far–field is given, the particulars ofthe source function are discussed and the formula for the source function used isgiven. Since the inversion is carried out in the frequency–domain, the final forwardequation is in that domain, including a practical approach to incorporate anelasticabsorption due to Q and the receiver response. The impact of anisotropy on thevarious terms in the forward equation is summarized, and the forward problemof constructing the moment tensor in an anisotropic medium is then discussed.Moment tensor construction requires certain angles, so the coordinate system and71.3. Thesis outlineangle conventions used throughout are given. Radiations patterns are computedand compared for isotropic and anisotropic media, for some basic and compositesources. Finally, three component (3C) synthetic seismograms are computed for ageometry with three downhole arrays of three–component (3C) receivers oriented inthe East,North,Up (ENU) coordinate frame. The medium is homogeneous VTI —transversely isotropic with a vertical axis of symmetry. The P, Sh and Sv (PHV )arrivals are shown rotated to a scalar amplitude based on their receiver polarizationvectors for a convenient inspection of synthetics together with the radiation patterns,which are displayed in an upper hemisphere equal–area projection including receiverlocations.In chapter 3 the theory and implementation details are given for the VTI raytracer which is the forward modelling engine and tool for the computation of theray theory Green function. The starting point is the Christoffel equation and thecase of homogeneous general anisotropy. VTI parameterizations are given prior todiscussing the layered VTI case. The ray tracer is exact and dynamic. The specialcases for geometrical spreading for purely vertical and horizontal rays are givenand effective ray length spreading is proposed as an unconditionally stable formfor spreading. The two–point ray tracing algorithm and critically refracted (headwave) ray tracing are described. The issue of qSv triplication is discussed, as is thestrategy adopted to handle it. Wavefront snapshots are computed for a triplicatingmedium and the Hilbert transform of the source function is included on the revertedbranch for the triplicating qSv arrival.Chapter 4 covers the frequency domain inversion for frequency–dependent mo-ment tensor, which includes data covariance weighting and regularization parame-ter estimation. Uncertainties are quantified following linear inverse theory, and forsingle events the full posterior covariance matrix is used for multi-variate normalsampling. The case of a pathological geometry is considered. For multiple eventsa parameter uncertainty metric is used and the synthetic with multiple noise real-izations is used to study the response due to noise. The source function is obtainedfrom the frequency or time–dependent moment tensor using complex principal com-ponent analysis, and a robust estimate for the real moment tensor is obtained usinga weighted deconvolution procedure. Scalar moment and corner frequency are ob-tained by a parametric fit to the source function amplitude spectrum.Chapter 5 deals with moment tensor decomposition and visualization. Con-ventional decomposition is reviewed and a new decomposition valid for generalanisotropic media is described. The distorting effect of anisotropy is shown byexample and the success of the new decomposition is demonstrated in recoveringfrom the moment tensor the correct source geometry, free form the distortions dueto anisotropy local to the source. Histograms are given for the noisy moment tensorssampled and shown in chapter 4. New 3D visualizations for the new decompositionare shown, which have been implemented in PetrelTM , Schlumberger’s commercialinterpretation platform.81.3. Thesis outlineReal data are analyzed in chapter 6. The real data come from the West Pembinafield in central Alberta and come courtesy of Lightstream Resources. The formationis the tight Cardium sandstone. The stimulation was done in a horizontal welland monitored by two vertical wells where 12 3C receivers were deployed. Resultsare shown for selected events from a central stage, and a force source inversion isshown to fit the sliding sleeve events better than a moment tensor inversion. Theresults show an unexpected collection of source mechanisms, but are consistent withSchlumberger’s commercial code run on events from the same stage.In chapter 7 several possible future directions are covered. Suggestions of howthe inversion can be improved are made, and the extension to lower symmetriesof anisotropy is considered. Group velocities are computed for an orthorhombicmedium and for a cubic (salt) medium. The extended source model is revisited,waveform fit model calibration is discussed and simulation results for a surface arrayrecording geometry are shown. Finally, data analysis of aftershocks of the great 7.9Wenchuan mega–thrust event of May 12, 2008 is shown. It is conjectured that thefrequency domain approach presented in this thesis could be ideal to analyze suchlong–coda data if multiple events are considered simultaneously.In chapter 8 the conclusions that may be drawn from this research are putforward.9Chapter 2The Forward ProblemIn this chapter the anisotropic ray–theoretical results for a point source as repre-sented by a moment tensor are reviewed, the chosen source–time function is dis-cussed and the method for moment tensor construction is given. The influence ofanisotropy on radiation patterns is examined and synthetics are computed for a3-well downhole geometry.2.1 Stress glut sourceThe seminal papers by Backus and Mulcahy [7], [8] introduced the general conceptof the stress glut, σglut, being the difference between the continuum stress σC ascalculated in the elastic model (according to Hooke’s Law) and the true stress inthe real Earth. That is,σglut = σC − σtrue. (2.1)They defined rigorously the (zeroth–order, static) seismic moment as the volume andtime integrals of the stress–glut rate or the volume integral of the stress–glut change,describing any small seismic source (small compared with seismic wavelengths andpropagation distances). That is, the second–order moment tensor M is defined asM =∫ ∫∂σglut∂tdtdV =∫[σglut]dV, (2.2)where the integrals are over the duration and spatial support of the source wherethe stress–glut rate is non–zero. If the source is of short duration compared to theperiods of interest, the time-dependent moment tensor can be written:M(t) = MH(t), (2.3)where H(t) is the Heaviside function.The stress–glut concept of a source is adequate to describe internal sources, i.e.where net forces and torques must be zero by Newton’s 3rd law. In general it maybe necessary to extend this to include unbalanced forces, i.e. external or extendedsources. This will include force and torque sources which are used to represent non–indigenous sources where the reactive mass is considered outside the Earth or ata distance from the experimental system. Examples of such external force sources102.2. Particle displacementsare vibroseis, hammer, weight drop, etc., represented as simple, unbalanced forcesacting on a plane normal to their direction. The topic of extended source will bereturned to later.The label “moment” in “moment tensor” comes from the unit of the momenttensor — force times distance = moment. The moment tensor has units [ML2T−2].As seen later in the construction of M, these units are also made up of the productof volume times stiffness, or [(L3)(ML−3L2T−2)].A moment tensor source can be represented as a set of equivalent forces (Akiand Richards [1], section 3.3) with three components of force and three possible armdirections producing nine generalized force couples (Aki and Richards [1] fig 3.7),the nine elements of the 3 x 3 moment tensor. Simple couples have force and armin the same direction and are sometimes called vector dipoles. Because the sourceis internal there is no net force or torque and the moment tensor is symmetric withsix independent elements.2.2 Particle displacementsThe advantage of the moment tensor description of a source is that the response interms of particle displacement, u, is a simple linear function of the moment tensor:u(t,x) = u(t,x; xS) : ∗M˙(t), (2.4)where u is the third–order strain Green function, the source is located at xS, M˙(t)is the time rate of change of M(t) (2.3), being a time domain δ–function scaledby M as defined in (2.2), : is the contraction over two dimensions between tensorsand ∗ represents convolution. Given M, which will be constructed later, particledisplacements may be calculated using a choice for u, for example finite difference(FD) modeling, spectral element (SEM) modeling, reflectivity modeling or ray the-ory. For the remainder of this thesis the ray theory approximation (e.g. Chapman[19]) will be used for u.In a homogeneous, isotropic medium the third order ray–strain Green functionisu(t,x; x S) ≈∑raysδ(t− T )4piρc3RgˆE, (2.5)where R = |x − xS| is the ray length, T = R/c is traveltime, c is velocity and ρ isdensity. The unit polarization vector of the ray is gˆ and E is the ray strain tensor:E =12(pˆgˆT + gˆpˆT ) (2.6)where pˆ is a unit vector in the slowness direction. The second–order scalar productin (2.4) with substitution (2.5) isE : M = EijMij (2.7)112.2. Particle displacementswith the Einstein summation convention. It is noted that (2.6) is a general vector–dyadic form that leads, in isotropic media, to the trigonometric terms in the formulasfor radiation pattern (e.g. Chapman [19] 4.6.18 – 4.6.20).The anisotropic ray theory form of (2.5) isu(t,x; x S) ≈∑raysgˆRRe(ΠT ∆(t− T)e−iσpi/2)4pi(ρRρ SVRV SS)1/2cSE, (2.8)where ΠT indicates the product of all transmission / reflection (T/R) coefficientsalong the ray, ∆(t) = δ(t)− i/pit is the analytic delta function, T is travel time, σ isthe KMAH index counting caustics along the ray and V and c are group and phasevelocities, respectively, for source S and receiver R. Considering the layered raytracing code that will be used later there will be no caustics, but as the ray tracerdoes return complex transmission/reflection coefficients and they may be turned on,the analytic delta function is retained for the time domain version of the ray–strainGreen function.To obtain (2.5) from (2.8), replace the analytic delta function with δ(t− T ) (nocaustics or T/R coefficients) and note that c = cS = VR = V S for the homogeneousisotropic case and that spreading S1/2 = cR. The details of the terms in (2.8) aregiven in chapter 3.At this point it is useful to consider a restricted version of (2.4), simplified forthe case of a pure shear dislocation in a homogeneous isotropic medium recorded inthe far–field (Aki and Richards [1] eq. 4.32):u(t,x) =Rθφ4piρβ3RµA ˙¯u(t−R/β). (2.9)Here Rθφ represents the shear wave radiation pattern for polar angle (inclination)θ and azimuth angle φ, µ is the shear modulus at the source, A is the area of thefault across which dislocation has occurred, R/β is the event arrival time and ˙¯u(t)is the time–derivative of the displacement at the source, which is u¯(t). The reasonfor the time derivative is that the pulse shape of displacement in the far field isproportional to particle velocities at the source (Aki and Richards [1] eq. 4.29).The total displacement is depicted as u¯(∞), also called the saltus of displacement.The saltus of displacement will appear as [d] later. The scalar moment of the sheardislocation is µA[d] = µAu¯(∞). The form of ˙¯u(t) is known as the source–timefunction or simply source function.The particle displacements given by (2.4) and (2.8) are for an impulsive, infinites-imal source, i.e. small in both temporal and spatial domains. In reality, the sourcehas some finite (but small) duration and extent, so M contains a source function122.3. Source functionand should really be written as M(t), making particle displacementsu(t,x) = u(t,x; x0) : ∗M˙(t) (2.10)≈∑raysgˆRRe(ΠT ∆(t− T))4pi(ρRρ SVRV SS)1/2c SE : ∗M˙(t). (2.11)This admits a different source–time function for every element of M, but it is com-mon to assume that the geometry of the source can be treated separately from itstime variation, so a single source function, s(t) is usually ascribed to M˙(t).2.3 Source functionExtensive theoretical work has been done on the study of earthquake rupture, andthe important results are summarized in Aki and Richards [1] chapters 10 and 11.The source–time function used from now on is defined in terms of the moment tensorof (2.2) asM˙(t) = Ms(t), (2.12)where s(t) is the time–dependent far–field particle displacement pulse. A simple andconvenient form for s(t) that is widely used as it agrees with the dominant featuresof many recorded earthquake source spectra is the so–called Brune pulse model [14].In its simplest form (e.g. Madariaga [73]) the Brune pulse isb(t) = M0ω2c te−ωctH(t), (2.13)where M0 is the zero frequency asymptote (the scalar moment) and ωc = 2pifc withfc the corner frequency. This wavelet is causal due to the Heaviside function H(t)and rises with time t until the negative exponential takes over. Its peak occurs at1/ωc and it has an inflection point at 2/ωc. An example of a Brune pulse is shownin figure 2.1. The amplitude spectrum of (2.13) isB(ω) =M01 + (ω/ωc)2 , (2.14)which shows that it decays as ω−2, in accordance with many teleseismic body waveobservations. A problem with (2.13) is that, in terms of far–field particle velocitiesit is discontinuous at t = 0. Beresnev and Atkinson [10] considered a general formof far–field displacement parameterized by an integer power n in addition to the risetime or corner frequency parameter τ = 1/ωc. Following the notation in Beresnevand Atkinson [10], with displacement at the source u¯(t), the far–field displacementsareu¯′(t) =u¯(∞)n!τ( tτ)ne−t/τ , (2.15)132.3. Source functionFigure 2.1: Time (left) and frequency domain representations of far–field particledisplacement source time functions. The n = 1 or Brune pulse is shown in red, then = 2 pulse in green.which, with the zero frequency asymptote u¯(∞) = M0 leads to the formula for thesource function used for simulations:s(t) =1n!τ( tτ)ne−t/τH(t). (2.16)Notice in (2.16) the parabolic increase with time for n = 2 and consequentiallya linear increase in far–field particle velocity. Figure 2.1 shows n = 1, 2 pulses using(2.16) and their amplitude spectra. The corresponding far–field particle velocitiesare shown in figure 2.2. The wavelet for n = 2 particle velocity is preferred forsimulation and processing as it rises linearly from 0 without the need for filtering aswould the n = 1 Brune pulse. However, the spectral decay of the n = 2 displacementpulse follows ω−3 rather than the preferred teleseismic result of ω−2. But since themeasurable spectral decay is subject to much uncertainty as it is influenced by manyfactors (e.g. Q, scattering, interference), violating the ω−2 decay does not seem asimportant as it is to have a well–behaved far–field particle velocity wavelet, so then = 2 version of (2.15) is used for simulations. The formula for the n = 2 particlevelocity wavelet of figure 2.2 isv(t) = ω3c t(1− ωct/2)e−ωctH(t). (2.17)A popular extension of the Brune spectral model is due to Boatwright [12]. Inthat model the spectrum includes a term for absorption due to Q and an extra142.4. Discussion on other source modelsFigure 2.2: Far–field particle velocity pulses corresponding to the far–field displace-ment pulses shown in 2.1.fall–off parameter, γ. It is given byBb(ω) =M0e−pift/Q[1 + (f/fc)nγ ]1/γ, (2.18)which is equivalent to the Brune model with infinite Q, n = 2 and γ = 1. Someauthors (e.g. Tomic et al. [111]) claim a superior fit to data with γ = 2. In thefrequency domain inversion that follows both attenuation and dispersion due to Qabsorption are explicitly included in the forward model, so the Q term in (2.18) isnot needed. Given the many other factors affecting the source pulse (see section2.4), the n = 2 pulse of Beresnev and Atkinson [10] will be used for the simulationsand inversions that follow. Nevertheless a few issues about source spectra deservefurther discussion.2.4 Discussion on other source models2.4.1 Angle–dependent corner frequencyMost models of fault dislocation are based on slip in the plane of the fault wheredisplacement occurs at some rupture velocity, vr. If displacement proceeds untilsome time T then the corner frequency is roughly 1/T . The model of Sato andHirasawa (1973) [94] predicts a dependence of corner frequency on angle from the152.4. Discussion on other source modelsfault, different for radiated P and S waves. Madariaga (1976) [72] computed finitedifference simulations and also found an angle dependence of P and S corner fre-quencies. An inversion of azimuth–dependent pulse width for rupture speed wasdeveloped and applied to large California earthquakes by Ammon et al. [4] and in-terestingly that technique was applied to a set of azimuth–dependent shear pulsesin reservoir–induced earthquakes by Tomic et al. [111], but few examples exist andmany other phenomena (e.g. anisotropic Q) may explain such observations. Whilean angle–dependent source pulse differing for P and S is easy enough to include inray theory synthetics, the general lack of azimuthal receiver sampling makes an in-version for such dependence unlikely to succeed and so it a complication consideredbut not pursued.2.4.2 Tensile crack P and S corner frequencyWalter and Brune [125] considered the case of a tensile crack rather than a slipsource and found that the ratio of P and S corner frequency was different for theirtensile crack model than for the pure slip model. Eaton [38] studied their model andadvocated using the S/P spectral ratio as a diagnostic for detecting tensile crackevents. As with the angle–dependent source frequency, differing P and S sourcefrequencies tied to the proportion of tensile vs. slip would be easy to include insynthetics, but including such parameters as free in an inversion is another compli-cation not pursued. However, some observations of much lower frequency S than Parrivals have been made that do not seem to be explainable with Q, so there may besome merit in considering the above complications in the future. Since the inversionprocedure described later recovers a different time dependence and hence frequencydependence for each element of M, this could conceivably aid in the interpretationof microseismic event spectral properties.2.4.3 Extended sourcesAs mentioned earlier in section 2.1, it is possible to consider a more general, extendedsource model where unbalanced forces are added to the balanced force model ofa moment tensor source. This would include a simple force source and an anti–symmetric part of the moment tensor to allow for sources with non–zero net torque.An example of a source with non–zero net torque might be a shear vibroseis onthe surface with imperfect lateral shaking or a downhole vibrator with a rotationalcomponent. The modified version of (2.10), augmented for an extended source model162.5. Particle displacements in the frequency domainisu(t,x) ≈∑raysgˆRRe(ΠT ∆(t− T))4pi(ρRρSVRV SS)1/2(2.19)∗[gˆS · f(t) +1cSE : M˙(t)−1c SΞ : M˙A(t)], (2.20)where f(t) is a time–dependent force function, Ξ is a source rotation tensor andMA(t) is the asymmetric moment tensor. Note that this source model has six extraunknowns, three for the force vector and three for the anti–symmetric part of themoment tensor. The model of a moment tensor source together with a force sourcehas been considered for earthquakes. For example Chouet et al. [24] considered sucha source as a model for volcanic events. The model of a simple force source will beconsidered later when analyzing sliding sleeve events in real data. The inversion ofthe complete extended source model given by (2.19) is an area for future research.In the frequency domain inversion method that follows, no prior assumption aboutthe source–time function will be made, but a fit to the estimated source spectrumwill be made assuming a ω−3 source model to recover the source attributes of scalarmoment and corner frequency. The ω−2 or Brune source spectrum fit is an option.2.5 Particle displacements in the frequency domainThe frequency domain equivalent of (2.10) isu(ω,x S,xj) ≈[rays∑kgˆkjGk(ω,x S,xj)Ek]: −iωM(ω), (2.21)where now the subscript k has been used to denote different ray types, receivers havebeen indexed by j and the scalar propagation terms in (2.10) have been taken upin Gk. Note the term −iω for the time derivative of M. −iωM(ω) = M when theHeaviside form of M(t) 2.3 is true. In the frequency domain the scalar propagationterms becomeGk(ω) =ΠTkeiωTk4pi(ρRρSVRV SS)1/2cS. (2.22)2.5.1 A practical approach to QAs a practical addition to the forward formulation an absorption term due to Q isincluded, aQ(ω), so that Gk(ω) is replaced byGQk (ω) = aQ(ω)Gk(ω), (2.23)172.5. Particle displacements in the frequency domainwhere (Varela et al. [119], eq. 6)aQ(ω) = e−pitkQk(f+i 2fpi ln(ffref)), (2.24)which includes both amplitude and phase terms. fref is the reference frequency,the frequency at which there is no dispersion. Here Qk represents the harmonicaverage of tik/Qik for ray k in layer i, the effective Q, and tk =∑i tik. In the raytrace treatment Q is independent of ray angle (i.e. isotropic) but may be differentfor the three different ray types (P, Sv, Sh) and different within each layer. Theharmonically averaged or effective Q, Qeff , is returned by the ray tracer.2.5.2 Receiver responsePrevious equations have described far–field recordings of particle displacements, butseismometers or geophones generally record particle velocity or acceleration. Thusto simulate actual recordings the receiver response function must be included. Thisis represented by a term rj(ω) where j is the receiver index. Receiver responsesare generally available from the manufacturer. Schlumberger’s GAC-D (GeophoneACcelerometer, D-model) response is proprietary and is used later to correct realdata examples to particle velocity. For simulations particle velocity is used witha roll–off below 3Hz and an anti–alias roll–off at 0.9fnyq. Sensitivity numbers areincluded in the simulation of real recordings so that correct scalar moment andmoment magnitudes may be recovered.While it is admitted that in practice there may be different receiver responsefunctions – coupling characteristics – for different components of a three–component(3C) receiver, in what follows they are assumed to be equal. Leaney et al. [67]discussed a station–consistent deconvolution method to address the variability ofcomponent–dependent receiver coupling which may be applied as a processing stepprior to further data analysis.2.5.3 Final forward equationIncluding the propagation term due to Q and the receiver response function, thefinal, frequency domain forward model to be used for modeling and source inversionof vector recordings v isv(ω,xS,xj) ≈ r(ω)[rays∑kgˆkjGQk (ω,x S,xj)E˜kj]· −iωM˜(ω), (2.25)where now the vector form for E and M has been used, as all operations on secondorder, symmetric tensors can be replaced by vector operations on 6–vectors — an182.6. Moment tensor constructionisomorphism (Silver and Jordan [101]). The modified form is used, requiring somefactors of√2:E˜ = (E11, E22, E33,√2E23,√2E13,√2E12)T(2.26)andM˜ = (M11,M22,M33,√2M23,√2M13,√2M12)T. (2.27)Using vector notation, the radiation pattern (2.7) is now written as a dot productE˜ · M˜. Again, all terms inside the [] brackets come from the ray tracer, and while(2.25) is valid for general anisotropy, only the layered VTI case will be considered.The impact of anisotropy can be found in every term of (2.25). In an anisotropicmedium the polarization at receiver j for ray k, gˆkj , is not simply related to the rayvector as in an isotropic medium (i.e. parallel for P and perpendicular for S). Infact, phase, group and polarization vectors are all different in a VTI medium exceptfor purely vertical and horizontal propagation. The exception is the Sh wave in aVTI medium, which remains polarized in the horizontal plane. It is for this reasonthat the quasi= q qualifier on P and Sv — qP and qSv — is sometimes droppedfor Sh waves in VTI media.Continuing with the discussion of the presence of anisotropy in the terms of(2.25), transmission and reflection coefficients are affected, and times are obviouslyaffected, and through them the effective Q along each ray. The velocities that appearin G — the denominator of (2.8) — group at source and receiver and phase at source,are angle dependent; geometrical spreading depends on ray angle; and the ray straintensor is impacted as it is made up of the dyadic involving the take–off phase andpolarization vectors from the source. Last but not least the moment tensor itself isaffected. With so many factors it is difficult to draw any general conclusions aboutthe impact of anisotropy on recordings, but the problem can be broken up intoits major constituents, for example traveltime, geometrical spreading and radiationpattern. One of the great advantages of ray theory is that such components canbe studied separately and switched on or off. Later in this chapter the impact ofanisotropy on radiation patterns will be shown, but first the construction of M needsto be described.2.6 Moment tensor constructionIn this section the procedure for construction of a moment tensor in an anisotropicmedium will be described. The well known formula for a displacement discontinuitywill be used, augmented by a new form for the isotropic part due to a pressurechange in a cavity surrounded by the anisotropic medium. The coordinate systemand angle conventions used throughout will also be introduced.192.6. Moment tensor construction2.6.1 The displacement discontinuity sourceThe components of a moment tensor described by a displacement discontinuity (DD)across a fault in an anisotropic medium are given by (from Chapman [19] eq. 4.6.11)Mij = A[d]cijklnkdl (2.28)where nˆ = (n1, n2, n3)T is the unit normal to the fault surface and dˆ(d1, d2, d3)T isthe direction of the displacement, i.e. the displacement is[d] = [d]dˆ, (2.29)where [d] is the total displacement across the fault in the direction of dˆ. (2.28) canbe written in tensorial form asMDD = c : DDD, (2.30)where : is the contraction over two dimensions between the 4th order stiffness tensorc and a 2nd order tensor DDD that has been called the source tensor by Vavryc˘uk[121]. At this point it is noted that a matrix–vector operation can replace the 2–dimensional contraction between the fourth–order stiffness tensor and second ordertensor in (2.30). Thus for computations the stiffness tensor c is replaced by themodified Voigt matrix (Chapman [19], eq. 4.4.14) C˜, and the second–order tensorD is replaced by the 6-vector D˜.Expressed in terms of the normal and displacement unit vector dyadic, DDD isDDD =12A[d](nˆdˆT + dˆnˆT). (2.31)If the displacement is in the plane of the fault, normal to nˆ, then (2.31) describes apure slip source. When dˆ has a component normal to the fault then (2.31) describesa source with opening or closing in addition to slip. The angle between dˆ andthe plane normal to nˆ is an essential parameter to define and construct a generalmoment tensor. It is defined bysin(χ) = dˆ · nˆ, (2.32)where −90◦ ≤ χ ≤ +90◦.But (2.28) or equivalently (2.30) with (2.31) is an incomplete description ofthe moment tensor as only five parameters are needed to define it — A[d] andfour angles for the vectors nˆ and dˆ. While MDD contains an isotropic part ifthere is a component of opening or closing normal to the fault, it is missing a pureisotropic part due to a pressure change. For now it is assumed that the missing sixthparameter is the pure isotropic part; later in chapter 5 it will be seen that such acomposition is (almost) always complete and unique.202.6. Moment tensor construction2.6.2 The pressure partTo complete the construction of a general moment tensor in an anisotropic mediuman isotropic M must be added to MDD. This is defined as M˜iso = [P ]V where [P ]is the pressure change in a cavity of volume V . In Chapman and Leaney [22] it isshown thatM˜iso = [V ]κembI (2.33)where [V ] = V ′[P ]/κemb is the volume change (+/– expansion) due to a pressurechange [P ] where V ′ is the effective volume of the cavity. In Aki and Richards[1], examples 3.7 and 3.8, the effective volume is given as V ′ = piα2R3/β2 fora spherical cavity in an isotropic medium. κemb is the embedded bulk modulus,known in isotropic media to be κemb = λ + 2µ but in anisotropic media it is morecomplicated. In Chapman and Leaney [22], appendix C4, κemb is derived for generalanisotropic media (eq. C26) and also for weakly anisotropic media which leads to apractical result involving the spherically averaged qP squared velocity:κemb '[3(C11 + C22 + C33) + 4(C44 + C55 + C66) (2.34)+ 2(C23 + C13 + C12)]/15. (2.35)This approximate form for κemb is used in moment tensor construction and decom-position.2.6.3 Anisotropic moment tensorThe final construction formula for a moment tensor in a general anisotropic mediais then (Chapman and Leaney [22] eq. 86)M = M˜iso + MDD = [V ]κembI +12A[d]c : [nˆdˆT + dˆnˆT ]. (2.36)Notice that the terms [V ] and A[d] have units [L]3 or volume. The ∼ is writtenover Miso to signify that it may not be the total isotropic part of M — there maybe an isotropic part of M due to an opening or closing part in MDD. Before using(2.36) to construct M, angle conventions need to be specified for the calculation ofthe direction vectors nˆ and dˆ.2.6.4 Angle conventionThroughout this thesis the (East,North,Up) or ENU coordinate system is used. Itstems primarily from the fact that Schlumberger downhole tools are manufacturedaccording to this convention and the 3C software is designed to orient data to satisfythat geographical coordinate system. This is contrasted with the NED system (e.g.Aki and Richards [1]) commonly used in earthquake seismology. Using the ENU212.6. Moment tensor constructionFigure 2.3: Vector and angle conventions. Left: Polar coordinate angles showingazimuth increasing clockwise from north and dip increasing from vertical–up. Right:strike (φ), dip (δ) and rake (λ) angles with slip vector [s] and fracture plane normalvector nˆ. The ENU directions are represented by the unit vectors ıˆ, ˆ, kˆ. aˆ is thestrike direction vector and bˆ is the dip direction vector.system, angles are specified so that a vector pointing vertically upwards has θ = 0◦,horizontal is θ = 90◦ and down is θ = 180◦. Azimuth increases clockwise fromnorth at φ = 0◦. Note that this is not the polar azimuth used by mathematicians inspherical polar coordinates which is measured from x, (East) in a counter–clockwisedirection, whereas θ is the same as in spherical polar coordinates.In terms of the needed strike, dip and rake angles the symbols φ, δ and λ willbe used. Using the convention adopted, the strike is defined as a vector aˆ makingan angle with the northerly direction. The dip angle is defined such that lookingin the strike direction a dipping fault plane dips down to the right and the planenormal, nˆ, is taken to point upwards. δ = 0◦ is a horizontal plane, δ = 90◦ is avertical plane. Thus the dip vector is bˆ = aˆ× nˆ. The rake is defined in the plane ofthe fault so that λ = 0◦ corresponds to pure left–lateral strike slip motion; λ = 180◦means the motion is pure right lateral. If λ < 0 the hanging wall has moved downrelative to the footwall and the fault is normal; if λ > 0 it is a reverse fault. Interms of the slip vector sˆ, dip vector bˆ and strike vector aˆ, the rake angle is definedusing the atan2 function (atan2d for degrees) as λ = atan2d(−sˆ · bˆ, sˆ · aˆ). Figure2.3 illustrates the vectors and angles using the adopted conventions.2.6.5 Constructing DTo construct DDD using (2.31), the fault normal nˆ, and the displacement vectordˆ are defined using four angles. nˆ requires the angles of strike (φ) and dip (δ)222.6. Moment tensor constructionwhile dˆ also requires the slip direction in the fault plane (the rake angle, λ) andthe previously defined opening angle χ. Given these conventions the fault normal isgiven bynˆ =cosφ sin δ− sinφ sin δcos δ (2.37)and the slip vector is given bysˆ =cosλ sinφ− cos δ sinλ cosφcosλ cosφ+ cos δ sinλ sinφsinλ sin δ . (2.38)Given these two directions and the opening angle, χ, defined in (2.32), the finaldisplacement unit vector is thendˆ = sˆ cosχ+ nˆ sinχ. (2.39)To complete the terms needed in (2.31) the term A[d] is needed. This is in units ofvolume but as seen later this term has historically been called “potency”.2.6.6 Ascribing potencyAlternatively to defining A[d] and [V ] separately, the total potency, |[V ]|+A[d] canbe specified along with a ratio, ν:ν = [V ]/(|[V ]|+A[d]), (2.40)where −1 ≤ ν ≤ 1. The total potency can then be defined, Vtot = |[V ]|+A[d], or itcan be set to unity, defining it implicitly later through the scalar moment, which isthe approach taken.2.6.7 Scalar momentIn an isotropic medium a pure slip source is a double couple with scalar momentmagnitudeM = 2−1/2‖M‖, (2.41)where ‖M‖ = (M : M)1/2 is the Frobenius norm of M. The factor of 2−1/2 isincluded for historical compatibility with the original expression for the momenttensor of a fault–slip source (Burridge and Knopoff [17]; Silver and Jordan [101])M = M(sˆnˆT + nˆsˆT ). (2.42)In an isotropic medium this produces a scalar moment magnitude M = µA[d] where[d] is the amount of slip. It is noted that other definitions of scalar moment haveappeared, e.g. Bowers and Hudson [13].232.6. Moment tensor constructionPractically, the scalar moment is usually a large number and can be unwieldy. Itis common to use the moment magnitude, Mw, to describe the size of an earthquake,usually taken as (Wikipedia)Mw =23(logM − 9.1) (2.43)where M is in units of N-m. Thus after the vectors nˆ and dˆ have been determinedthrough specification of the four angles, and given the relative proportions of A[d]and [V ], given the medium stiffness tensor and hence κemb, a moment tensor can becalculated. It is then scaled to have the desired scalar moment given input momentmagnitude (2.43).2.6.8 Source parameter setGiven the above conventions and definitions a final parameter set can be given interms of the six parameters stemming from the original M: (Mw, ν, φ, δ, λ, χ). As anod to hydraulic fracturing where fractures tend to be vertical, dip will sometimesbe specified and results presented using co–dip rather than dip, that is D = 90− δ.The symbol O will be used for opening angle, O = χ, providing an angle set SDROso that a left–lateral strike–slip event on a vertical plane striking N–S has onlyone non–zero parameter, Mw, and all angles equal to zero. As will be seen laterin chapter 5, these angles can be recovered from M, and using a new approach,without the distortions due to anisotropy.2.6.9 Full potency tensorThe term “potency” has been used in conjunction with the terms A[d] and [V ] in(2.31) and (2.36). The history behind this terminology is described in section 1 ofChapman and Leaney [22]. Chapman and Leaney [22] write the potency tensor as(eq. 87)D = s : M = s : MISO + s : MDD = [V ]κembK +12A[d](dˆnˆT + nˆdˆT), (2.44)where K is the hydrostatic pressure modulus tensor, K = s : I (Chapman andLeaney [22], eq. A7) and s is the compliance tensor. Note the additional termadded to what Vavryc˘uk calls the source tensor (2.31), describing the expansion byan amount [V ]. In an isotropic medium K is diagonal with equal elements 1/(3λ+2µ)(Chapman and Leaney [22] eqs. 88 and C1). In a VTI medium the K33 value willgenerally be greater than the other two diagonal elements. For general anisotropyK may have non–zero off–diagonal elements. Thus an isotropic pressure source withonly MISO 6= 0 can have a complicated potency tensor. Equally, a simple potencytensor can give rise to a moment tensor with non–zero off–diagonal elements. Some242.7. Radiation patternssources are simpler in potency and some are simpler in moment (see for example thespecial sources in Appendix B of Chapman and Leaney [22]). The potency tensorcan also be written in matrix–vector form asd˜ = C˜−1d˜ = S˜m˜ (2.45)where S˜ = C˜−1 is the inverse of the modified Voigt matrix.2.7 Radiation patternsIn this section several examples will be shown of radiation patterns for qP , Sh andqSv rays, in isotropic and VTI media, focusing on pure slip, pure opening or tensileand a composite slip+opening source.2.7.1 Radiation pattern computationAs mentioned in section 2.2, equation (2.7), the direction–dependent radiated am-plitude for ray type k for a moment tensor M is given by Ek : M (2.7), where E isthe second order ray–strain tensor at the source and : represents the scalar productbetween second order tensors. In an isotropic medium the ray direction is definedsimply by the take–off angles (θ, φ). But in an anisotropic medium, one has a choice.It is most convenient to index computation of the radiation pattern by phase anglessince the phase direction vector is input to the Christoffel equation (seen later, eq.3.2) to get polarization vectors gˆ which are both needed for E. While one couldplot radiation patterns according to phase direction, it is more representative toplot them according to the ray or group direction (3.1). (Having said that, figure2.4 uses phase angle to illustrate a point.) In the later section 3.5 the simplifiedexpressions for group velocity in VTI media are given.Psenc˘ik and Teles (1996) [87] and Leaney and Chapman (2010) [59] investigatedthe influence of anisotropy on radiation patterns. In Leaney and Chapman [59] itwas found that anisotropy can have a significant impact. Firstly, for a given area–displacement or potency, A[d], the scalar moment is generally larger in an anisotropicmedium simply because the medium is stiffer (the moduli in c are larger). It shouldbe noted that the isotropic medium assumed here corresponds to the vertical Pand S velocities, a faster isotropic medium can be found that fits the sphericallyaveraged velocities of the medium (Chapman [19], page 132, exercise 4.6). Usingsuch isotropic velocities would reduce the difference in magnitude but retain thedifference in form. Secondly, the radiation pattern is distorted due to influences ofthe stiffness tensor on M and the misalignment of phase and polarization vectorsthat go into E. For the canonical case of a pure slip source where dˆ is in theplane of the fault, anisotropy can cause non–DC (double couple) components in M(.e.g. Julian et al. [50]). This means that if not accounted for in the decomposition252.7. Radiation patternsFigure 2.4: PHV radiation at 60◦ phase angle for a pure slip source with N–S strikeslip on a vertical plane in an isotropic medium (solid) and VTI medium (dashed).Color coding follows PHV =BRG=blue,red,green). Left: un–normalized momenttensors, right: normalized moment tensors. Normalization serves to highlight theeffect of anisotropy on the ray strain tensor part of the radiation pattern. Theabsolute value of radiated amplitude is shown in these figures (i.e. polarity is notshown).of the moment tensor, false tensile components would be interpreted. Since thisinformation is of great interest in hydraulic stimulations, clearly this distortion dueto anisotropy is not wanted. A decomposition that removes the distortions due toanisotropy is described in chapter 5.2.7.2 Radiation pattern examplesLeaney and Chapman [59] investigated PHV radiation patterns for a pure slip sourcein ISO and VTI media. The impact of the ray strain tensor E on the radiationpatterns can be highlighted by normalizing the ISO and VTI moment tensors tohave equal norms ‖M‖, a result also shown in figure 2.4. These plots only weremade for a constant phase inclination angle (instead of group angle) of θ = 60◦,hence there is no distortion of the Sh radiation pattern because VTI anisotropydoes not alter the Sh polarization vector, so E is the same for both ISO and VTIcases.In figure 2.5 PHV (recall PHV stands for qP, Sh and qSv) radiation patterns areshown for a pure slip source on the entire focal sphere using an equal area projection.In this projection the group inclination angle θ is re–indexed to 2 sin θ/2. Since inmicroseismic studies the sensors are generally above the sources, preference is given262.7. Radiation patternsto the upper hemisphere so vertical–up is plotted in the center of the circle. Thesedisplays are plotted in true relative amplitude, showing that the same potency orA[d] produces larger radiated amplitudes (brighter colors) for the case of anisotropy.There is also a subtle but significant rotation or distortion. As will be seen later,the characteristic SDRO angles describing the orientation of the source are in errorwhen conventional moment tensor decompositions are used on moment tensors fromanisotropic media.In figure 2.6 the radiation patterns for a pure opening or tensile crack event withthe same fracture plane geometry is compared to that of the pure slip event in theVTI medium, but it is noted that these are normalized individually.To appreciate the differences in radiated qP amplitude the two cases are plottedtogether as viewed from the south in figure 2.7. Note the positive qP amplitudesfor the +TC source.A final radiation pattern example has equal parts slip and opening with χ = 45◦and an additional negative isotropic expansion through non–zero [V ] with ν = −0.20.The full parameter set (2.6.8) is (−2,−0.20, 60◦, 40◦, 20◦, 45◦). The qP radiation isshown in figure 2.8.2.7.3 Basic ISO and TC sourcesAlthough more will be said later on special cases of sources when moment tensordecomposition is discussed, a short digression seems appropriate here. Given the ap-proach taken of moment tensor construction given by (2.36) and the approximationdescribed in 2.6.2, a source with ν = ±1 has a diagonal M with equal values andradiates only P , equal in all directions with polarity determined by the sign of ν.A source with ν = 0 and χ = ±90 is a pure tensile crack (TC) source with openingor closing normal to an infinite plane. A TC source has a simple potency tensor(Chapman and Leaney [22], eq. B10b) and in an isotropic medium such a sourcehas a moment tensor with a dipole added to an isotropic explosion (Chapman andLeaney [22], eq. B10b). Equivalently, Aki and Richards [1], page 51, give the mo-ment tensor for a TC source in an isotropic medium as diag(M) = A[d](λ, λ, λ+2µ)where displacement is normal to the x3 = 0 plane. Although a TC source has aparticularly simple form of potency tensor some sources are simpler in the momentdomain. For example, a dipole has only one non–zero moment element but threenon–zero elements of its potency tensor (Chapman and Leaney [22], B4¯).272.7. Radiation patternsFigure 2.5: Equal area complete focal sphere radiation patterns for qP , Sh andqSv rays (PHV are plotted top to bottom). A pure slip source with SDR =(60◦, 40◦, 20◦) and equal potency A[d] has been placed in an isotropic medium(left) and a VTI medium (right). The VTI medium has Schoenberg parameters(0.25, 0.20, 0.30) (see eq. 3.9). The upper–hemisphere is contained within the blackcircle, which represents the horizontal plane. In these and similar figures red islarge, blue is small and green is mid–way between. These plots are shown in truerelative amplitude to highlight the larger radiated amplitudes in the generally stifferanisotropic medium.282.7. Radiation patternsFigure 2.6: Radiation patterns for pure slip (left) and pure tensile source withpositive opening angle (right). In both cases the medium is VTI. Each plot isnormalized individually.292.7. Radiation patternsFigure 2.7: VTI P radiation as viewed from the south for a pure slip source (red)and a pure tensile source with positive opening angle (blue). The zero line is shown.302.8. Ray trace syntheticsFigure 2.8: Equal area complete focal sphere radiation pattern (left) for P rays.A composite source with slip, tensile opening and negative expansion is used. Thesource parameter set is (−2, 0.20, 60◦, 40◦, 20◦, 45◦). On the right is shown the radi-ated amplitude from the south with a zero line.2.8 Ray trace syntheticsThe previously developed theory for the forward problem is now illustrated with anexample of ray theory synthetics. Leaney et al. 2011 [62] constructed an illustrativeVTI model to show the impact of VTI anisotropy on the recovery of potency tensorattributes. Here a similar approach is taken to bring together aspects of the forwardproblem — ray–strain Green function, source function and receiver response function— into synthetic three component seismograms for a selected downhole geometry.The source and receiver geometry is shown in figure 2.9). Rays are not shown asthe medium is homogeneous so all rays are simply straight lines. The details of theVTI ray tracer are discussed in the following chapter 3. Three component (ENU)seismograms are shown in figure 2.10. The receiver array traces are plotted groupedby East, then North, then Up. The model is homogeneous VTI with parametersVp = 4000 m/s, Vs = 2300 m/s, ρ = 2500 kg/m3 and Schoenberg anisotropy Ep =0.20, Ea = 0.25, Es = 0.30 (anisotropy and parameterizations are discussed in thenext chapter). The n = 2 pulse (2.17) is used with fc = 100Hz. Geometricalspreading is included but no Q. Gaussian noise has been added with standarddeviation equal to 1% of the maximum amplitude of all the ENU traces.For comparison, 2.10 is shown with QP = QH = QV = 100 included in thecomputation in figure 2.11. Reference frequency, fref in equation (2.24) was set to200Hz.A useful way of viewing waveform event data is as scalar PHV arrivals. While itwould be easy to obtain scalar PHV amplitudes by simply leaving out the receiverpolarization vector so as not to project arrivals onto the ENU components, with312.8. Ray trace syntheticsFigure 2.9: Three array geometry in map(top) and section view used for simulation.Two vertical arrays and one horizontal array surround a source at a distance ofabout 565m. All arrays have 12 3C receivers spaced at 30m.322.8. Ray trace syntheticsFigure 2.10: Ray theory particle velocity synthetics in a homogeneous VTI modelwith DD source parameters SDRO = (60◦, 40◦, 20◦, 45◦), n=2 fc = 100Hz sourcewavelet, geometrical spreading and additive Gaussian noise with standard deviationequal to maxamp/100. The receiver array consists of two vertical arrays and onehorizontal array, each with 12 3C receivers. From left to right are groups of 36 tracesfor the East, North and Up (ENU ) components. The three groups of 12 traces canbe seen as the different arrays.332.8. Ray trace syntheticsFigure 2.11: Same as figure 2.10 except including absorption and dispersion due toQ. Q = 100 was used, chosen as equal for PHV rays. fref = 200Hz was used.342.8. Ray trace syntheticsPHV separation accomplished by simply not summing over rays, an approach forreal data is needed. Here the approach taken is to project the data according to thereceiver polarizations coming from the ray tracer — a model–based projection. Theresult of this operation is shown in figure 2.12. The Sh arrival is perfectly isolated asit is polarized in the horizontal plane and is orthogonal to both qP and qSv; there issome cross–talk evident between qP and qSv owing to their non–orthogonal polar-ization vectors in a VTI medium. Leaney (2008) [57] described a frequency–spaceleast squares inversion to separate arrivals that minimizes the cross–talk problem(and which will be revisited later in the discussion on nonlinear inversions in section7.4). Using PHV projection as in figure 2.12 the variation of arrival amplitude andpolarity can now be more easily appreciated both between arrays and within arraysdue to radiation pattern changes with azimuth and inclination.Figure 2.12 also shows the PHV radiation patterns with the receivers added tothe focal sphere projection in order to show the sampling afforded by their locationsrelative to the source. Included in these radiation patterns are the scalar terms ofthe homogeneous Green function using the approximation of effective ray lengthspreading (next chapter, eq. 3.40) for efficiency. The scalar propagation terms areincluded here as the previously computed radiation patterns alone (E : M) aredeceptive of the true radiated amplitude. This is because the velocities appearing inthe denominator of (2.8) serve to scale up the radiated shear amplitudes. The PHVamplitudes seen in the waveforms in figure 2.12 can now be directly located on theradiation patterns at the locations of the black dots for receivers: the vertical arrayto the northwest is plotted as the first set of waveforms, then the SW vertical array,then the horizontal array to the east. Notice that the P radiation pattern has noblue colour and is everywhere positive due to the positive opening angle.Another useful display is to view all components simultaneously by making use ofthe signal envelope, computed using the Hilbert transform (e.g. Ulrych et al. [115]).This can be done by simply summing the individual envelopes of un–oriented XYZdata, or the oriented ENU components. The advantage of envelope summationis that it can be done to facilitate time picking without needing to know receiverorientation (or well deviation). This technique is described for arrival onset timepicking (Leaney [55]). Here something slightly different is done. Here the sample–by–sample maximum of PHV waveform envelopes is used. This is shown applied tothe PHV waveforms in figure 2.13.Also shown in figure 2.13 are the PHV arrivals times in color. Notice the sig-nificant time difference between the two shear arrivals. This is primarily due to theEs parameter, which when positive means that the horizontally travelling Sh willarrive earlier than the qSv. For more vertical angles the Ea parameter makes theqSv faster while the Sh becomes slower. It is not uncommon to observe a crossoverbetween Sh and qSv on real data. It should be said that the source–receiver dis-tance of about 565m was chosen for this simulation to assure clear separation ofshear arrivals. A more typical source–receiver separation value for downhole data352.8. Ray trace syntheticsFigure 2.12: The simulated event recording from figure 2.10 rotated to scalar PHV(qP, Sh, qSv) arrivals based on the ray–trace receiver polarizations (top). Bottom:PHV radiation patterns with receiver locations showing focal sphere sampling. Anapproximate homogeneous VTI Green function is included in the radiation patternsfor direct comparison with waveform amplitudes. Radiation patterns are displayedin true relative amplitude.362.8. Ray trace syntheticsFigure 2.13: The maximum PHV envelope function applied to the wave-forms of figure 2.12 with arrival times shown according to the color codingPHV =BRG=blue,red,green).372.8. Ray trace syntheticswould be about 300m.Before describing the inversion of waveforms such as those shown in figure 2.10 or2.11 for source function and moment tensor, the kernel of the forward and inverseschemes — layered VTI ray tracing — needs description.38Chapter 3Layered VTI DynamicRay–tracingIn this chapter the theory and implementation are described for a dynamic ray tracerfor layered media where the layers are homogeneous and transversely isotropic witha vertical axis of symmetry (so–called VTI media). While ray theory has a longhistory and no discussion of anisotropic ray tracing should exclude the work ofC˘erveny´ [26], C˘erveny´ [27], the majority of what follows comes from Chapman [19].The ray tracer is the kernel of most of the forward computations used in this thesis.To begin, a homogeneous, general anisotropic medium is considered.3.1 Homogeneous, general anisotropyBefore considering layered VTI models, it is instructive to consider rays in a homo-geneous medium. In an isotropic medium the ray (=group) and phase directions arethe same and polarizations are simply related to the ray (e.g. P polarizations arealong the ray and shear polarizations are normal to the ray), but in an anisotropicmedium this is not the case. Figure 3.1, adapted from figure 5.7 in Chapman [19]shows the geometry of ray path and wavefront; slowness, group velocity and polar-ization vectors.In a homogeneous anisotropic medium the ray connecting source and receiveris still the obvious geometrical straight line but the travel time for that ray is nottrivial to calculate. This is because, while T = L/V where T is time and L is theray length, the velocity needed in the anisotropic case is the magnitude of the groupvelocity vector in the direction of the ray. In an anisotropic medium the componentsof the group velocity vector are given by (Chapman [19] 5.3.20)Vi = aijklpkgˆj gˆl, (3.1)where aijkl is the density-normalized 4th order stiffness tensor, gˆj and gˆl are com-ponents of the normalized polarization vector gˆ, pk is the k − th component of thephase slowness vector p and summation is assumed. To determine gˆ and p thefollowing eigen (Christoffel) equation (Chapman [19] 5.3.26) must be solvedΓ̂gˆI = c2I gˆI (3.2)393.1. Homogeneous, general anisotropyFigure 3.1: A ray path and wavefront with slowness (p), group velocity (V) andpolarization (gˆ) vectors where for clarity the polarization vector is shown for aquasi–shear ray. The distance between wavefronts dT apart in time is |V|dT in theray direction V, and dT/|p| in the slowness direction p.for the eigenvalues c2I and eigenvectors gˆI , I = 1, 2, 3 for qP , qS1 and qS2 codes. Γ̂is the Christoffel matrix given by (Chapman [19] 5.3.15)Γ̂ = pˆj pˆkajk (3.3)where ajk are a collection of 3 x 3 sub-matrices built up from the density-normalizedVoigt moduli (Chapman [19] eq. 4.4.39). Thus, having specified a phase directionpˆ = (pˆ1, pˆ2, pˆ3) for (3.3), having solved (3.2) for phase velocities cI and polarizationvectors gˆI and setting p = pˆ/c, the above three equations allow the group velocityvector or ray direction to be determined.It is noted here that the group velocity vector gives the direction of energypropagation, and equation (3.1) belies certain complexities. The Christoffel equationis for a given phase direction or plane wave, along which amplitude is uniform, soto get the direction of energy propagation one needs a spread of slownesses, i.e. acurved wavefront to define the group direction which is normal to the phase slownesssurface (see figure 3.3 for an example of a phase slowness surface). Thus a derivativeof the slowness surface is needed to determine the group velocity. Thomsen (1986)[109], eq. 14* cites Berryman (1979) [11] for the following equation relating groupand phase velocities given the phase angle, θ in a VTI medium:V 2(φ(θ)) = v2(θ) + (dv/dθ)2, (3.4)where φ(θ) is an expression for group angle given phase angle found in Thomsen[109]. Equation (3.4) shows how the derivative of phase velocity is needed to obtain403.2. VTI parameterizationsgroup velocity. Equation (3.1) is a compact form for group velocity where certainderivatives simplify (see the discussion after Chapman [19] 5.3.20 and Exercise 5.1).It is important to realize that the ray connecting two points through the ho-mogeneous anisotropic medium is a straight line but a nonlinear problem must besolved to determine the phase direction that will produce the correct group direc-tion to connect the ray end points. For a general anisotropic homogeneous mediumtwo take-off angles, which define the initial phase vector from the source, must besought to minimize the distance from the ray to the target (receiver). Commonlya distance measure in the wavefront is used in 3D ray tracers; in VTI media raysare constrained to a vertical plane so only a single angle or parameter is needed andthe distance measure involves horizontal range only. The process of determiningparameter(s) to connect two points in a medium with a ray is known as two-pointray tracing and many techniques have been developed to do it efficiently, even incomplicated 3D models. Complications can arise even for some relatively simpleanisotropic media - for example multiple solutions are possible for the qSv ray intriplicating VTI media such as that shown in figure 3.3.Having found the ray of interest certain dynamic (amplitude) attributes such asgeometrical spreading and polarization vectors can be calculated. In a model con-taining interfaces that represent boundaries between regions with contrasting elasticproperties, transmission and reflection coefficients are also dynamic ray attributesof interest to be computed and included in the ray signature. As will be seen in thefollowing sections, all is greatly simplified (and efficient) in layered VTI models.3.2 VTI parameterizationsSince VTI symmetry is ubiquitous and the moduli of the Voigt matrix do not providemuch intuition about the angle-dependent nature of velocities, anisotropy parame-ters have been defined to quantify the amount of anisotropy. The exploration seismicindustry has adopted Thomsen’s (1986) parameters ([109]) to define a VTI medium: =A11 −A332A33(3.5)δ =(A13 +A55)2 − (A33 −A55)22A33(A33 −A55)(3.6)γ =A66 −A552A55(3.7)σ =V 2pV 2s(− δ) (3.8)where Aij = Cij/ρ are the density–normalized stiffness moduli. is often (erro-neously) referred to as a “percentage” P anisotropy, δ is a parameter controlling413.2. VTI parameterizationssmall offset or sub-vertical wavefront curvature or moveout, γ is the Sh equivalentof and σ is a parameter for qSv. Another parameter, η = ( − δ)/(1 + 2δ) dueto Alkhalifah and Tsvankin [3], is also used to describe P -wave moveout. η is alsoknown as an anellipticity parameter as it quantifies the deviation of the qP slownesssurface from an ellipse. Note the similarity of η and σ, Thomsen’s qSv parameter.Another useful VTI parameterization with a more rigorous normalization andsome beneficial properties is due to Mike Schoenberg, first published by Carrion et.al. (1992) [18] and then renamed in Schoenberg and Helbig (1997) [98] and again inSchoenberg and de Hoop (2001) [97]. In fact Schoenberg had written “an answer” toThomsen’s 1986 parameters in 1987 (pers.comm. 1989) but it remained unpublishedand unused until the Carrion paper which came out of their time together in Brazil.Schoenberg’s VTI parameters, with a unified naming convention areEp =A11 −A33A11 +A33(3.9)Ea =(A11 −A55)(A33 −A55)− (A13 +A55)2(A11 −A55)(A33 −A55)(3.10)Es =A66 −A55A66 +A55. (3.11)In place of Thomsen’s parameters and γ are now Ep and Es, normalized to bebounded on [-1,1] for “mild anisotropy”. In place of σ or indeed η or δ is now Ea,the anellipticity parameter. It uniquely describes the behaviour of qSv anisotropy— the deviation from an isotropic circle — but at the same time describes thedeviation of qP phase slownesses from an ellipse, hence the name and the subscript.It is also normalized to be bounded on [−1, 1] under conditions of mild anisotropy,but is generally positive for sedimentary rocks. For weak anisotropy the differencebetween and Ep is small and both remain close to an intuitive ratio-type measure ofhorizontal/vertical velocity (√A11/A33− 1) but for stronger anisotropy, as possiblefor shales, Thomsen’s parameter quickly grows and even exceeds 1. Figure 3.2 shows and Ep versus a ratio, which shows just how misleading can be for strongeranisotropy.Not surprisingly, Schoenberg’s parameters have been found superior for inver-sions. In the literature δ is often used with for qP and qP − qSv data but δis generally smaller, may be negative and is strongly correlated with . σ or ηwould be a better inversion parameter than δ but this parameterization is rarelyused. Schoenberg’s parameters are preferred, particularly for inversion; both willbe used when describing models. The mapping between Schoenberg and Thomsen423.3. Example: Algerian shale, strong VTI anisotropyFigure 3.2: Thomsen’s (red) and Schoenberg’s ellipticity parameter Ep (blue)versus a velocity ratio measure√A11/A33 − 1.parameters isEp =1 + (3.12)Ea =2A33A11 −A55(− δ) (3.13)Es =γ1 + γ. (3.14)A note here that 55 is usually used for the modulus controlling vertical shearpropagation in a VTI medium, being as it controls propagation in the 1− 3 or x− zplane, but since in a VTI medium A44 = A55 one can choose either. Chapman [19]uses A44.3.3 Example: Algerian shale, strong VTI anisotropyFigure 3.3 shows phase slowness and group velocity plots for qP , qSv and Sh for thestrongly anisotropic medium shown in Leaney [58]. The parameters were determinedfrom an inversion of qP and qSv phase slownesses and polarizations described inLeaney and Hornby [66]. This dramatic example of strong anisotropy comes fromAlgeria and is a shale of Silurian age. The VTI parameters for this shale were433.4. Example: Earth inner core anisotropydetermined to be Vp(0) = 3500 m/s, Vs(0) = 1530 m/s, Ep = 0.37, Ea = 0.64 withequivalent Thomsen parameters = 0.58, δ = −0.05. To generate the Sh curves Eshas been set equal to Ep; equivalently γ has been set equal to .Figure 3.3: Phase slowness (left) and group velocity (right) for a strongly anisotropicshale. qP=blue, qSv=green and Sh=red.Note the triplication in the qSv group velocity curve. This is discussed later insection 3.11.3.4 Example: Earth inner core anisotropyWhile the sedimentary rocks of the upper crust are often strongly anisotropic, globalseismology has revealed that other regions of the earth are also anisotropic. Forexample, the inner core is thought to be anisotropic with TI symmetry with a fastdirection oriented towards the poles (Song and Richards [104]). Buffett discussescore anisotropy in Physics Today, November 2013 [16]. Curiously, core anisotropy isthought to be stronger in the western hemisphere than the eastern. Using a densityfor iron of 13 g/cm3 the TI moduli of the core in km2/s2 g/cm3 are C11 = 1801,C12 = 865, C13 = 810, C33 = 1919, C55 = 445 and C66 = (C11−C12)/2 = 468. Thisleads to standard parameters (velocities in km/s): V p(0) = 12.15, V s(0) = 5.85,Ep = −0.032, Ea = 0.212, Es = 0.025. The phase and group velocity polar plotsare shown in figure 3.4.443.5. Layered VTI modelsFigure 3.4: Phase slowness (left) and group velocity (right) polar plots for the TIanisotropy of the inner core. qP=blue, qSv=green and Sh=red.3.5 Layered VTI modelsAs discussed above, to trace a ray in an anisotropic medium the group velocityvector is needed. In a 1D layered medium composed of VTI layers with symmetryaxis in the x3 direction the horizontal component of phase slowness is constant.Following Chapman [19] the sans serif font is used for the horizontal parts of vectorsso p = (p1, p2)T is the horizontal part of p. And without loss of generality due torotational invariance p will be referred to equivalently as p1. The horizontal distanceor range X( p) and travel time T ( p) are given by (Chapman [19] 5.7.3 and 5.7.4)X( p) =∑lV(l)V3(l)dx3(l) (3.15)andT ( p) =∑l1V3(l)dx3(l) (3.16)where l is the layer index and dx3 is the vertical thickness of the layer l and sum-mation continues until some specified maximum vertical distance Z (i.e. depth) isattained.Having specified the horizontal phase slowness p1, the horizontal and verticalcomponents of the group velocity are needed to solve the above ray equations. Twoof the solutions to (3.2) have polarizations in the plane of propagation, they are theqP and qSv rays, so-called because as opposed to the isotropic case where the Ppolarization is in the ray direction it is now not necessarily, hence “quasi − P” or453.5. Layered VTI modelssimply “qP”, similarly for qSv. A quadratic arises from the VTI simplification ofthe Christoffel equation, for qP − qSv so that (Chapman [19] 5.7.32)p23 =B ∓ [B2 − 4A33A44(A11p21 − 1)(A44p21 − 1)]1/22A33A44(3.17)where (Chapman [19] 5.7.33)B = A33 +A44 −Ap21 (3.18)with (Chapman [19] 5.7.31)A = A11A33 +A244 − (A13 +A44)2. (3.19)By definition the upper (–ve) sign in (3.17) corresponds to the qP and the othersolution is referred to as qSv.Coefficients of the form Aij above are the density-normalized moduli of theelastic stiffness tensor following the Voigt indexing convention.Having determined p3 from p1 and the medium moduli, the horizontal and ver-tical components of group velocity are (Chapman [19] 5.7.36)V1 = p1(2A11A44p21 +Ap23 −A11 −A44)/D (3.20)and (Chapman [19] 5.7.37)V3 = p3(2A33A44p23 +Ap21 −A33 −A44)/D (3.21)with (Chapman [19] 5.7.38)D = (A11 +A44p21 + (A33 +A44)p23 − 2. (3.22)A qP or qSv ray, or a ray containing a sequence of qP and qSv segments cannow be traced, where the ray advances for an increment of vertical thickness ∆x3by a horizontal range increment ∆x1 = V1/V3∆x3 and the increments are summedas in (3.15). In a homogeneous layer the incremental range may be written as∆x1 =p1p3(2A11A44p21 +Ap23 −A11 −A44)(2A33A44p23 +Ap21 −A33 −A44)∆x3. (3.23)Similarly incremental travel time for use in (3.16) may be written as∆T =D(2A33A44p23 +Ap21 −A33 −A44)∆x3. (3.24)The third solution of the VTI Christoffel equation is simpler and has polariza-tion in the transverse, horizontal direction and is referred to as the Sh ray. (In a463.6. Two point VTI ray tracinghorizontally layered VTI medium the Sh ray always has a horizontal polarizationnormal to the ray direction so the ”q” qualifier is not used.) Given the horizontalcomponent of phase slowness p1 the vertical component for the Sh ray is (Chapman[19] 5.7.24)p3 = (1A44−A66A44p21)1/2(3.25)and the horizontal and vertical components of group velocity are (Chapman [19]5.7.25)V1 = A66p1 (3.26)and (Chapman [19] 5.7.26)V3 = A44p3 (3.27)which again are used in (3.15) to trace the ray.3.6 Two point VTI ray tracingSince iteration is required for two-point ray tracing and dynamic attributes areusually wanted only after the ray has been found, basic VTI ray tracing ideas andthe two-point algorithm are discussed here, leaving dynamic attributes for followingsections. The descriptions will at times refer to obvious parameters and routinenames in the code so that this chapter may serve as a sort of doc file for the raytracing subroutines.Given a layered VTI model, the goal is to trace a ray from the top to thebottom of the stack of N layers. Since the ray path depends only on the horizontalslowness p1 it is reasonable to ask what bounds may be placed on that parameter.For a vertical ray there is no horizontal component to the ray so p = 0. For apurely horizontal ray p3 = 0 and the above equations can be used to show thatthe maximum horizontal slowness in any layer is p(max) = 1/V1(max) = 1/Vh orp(max) = 1/√A11 for qP , p(max) = 1/√A44 for qSv and p(max) = 1/√A66 forSh. Since the horizontal slowness anywhere along the ray (i.e. for all layers) cannotexceed the maximum allowed in any layer, the maximum slowness for a stack oflayers isp(max) = 1/max(Vh(l)) (3.28)which can be determined prior to any ray tracing. Note that in the case of a raywith a mixture of qP and qSv rays (i.e. containing mode conversions) Vh for theray type in layer index l must be used in (3.28).A Newton-Raphson line search algorithm may be used (e.g. Press et al. [85],section 9.4) to update an initial guess at p. An initial guess at p for iteration k = 1could be half of the maximum or p(k) = p(max)/2. The basic equation for the kth473.6. Two point VTI ray tracingp update ispk+1 = pk − (Xk − X)/(∂ Xk∂ pk). (3.29)The correct p is found when |Xk − X| < tol where tol is machine precision.The derivative ∂X/∂p will be given later for geometrical spreading computa-tions, but it has been recognized (Costa et al. [25]) and experience has taught thatnumerical problems can occur for sub-horizontal rays, so the above Newton-Raphsonalgorithm with semi-analytical ∂X/∂p is not used. Rather, a simple secant method(e.g. Numerical recipes, section 9.2 Press et al. [85]) is used with a fall-back to a bi-section if an update produces a pk > p(max). It is interesting to note that the secantmethod can be thought of as a finite difference version of Newton’s method but his-torical evidence indicates that it predates Newton by almost 3000 years (Wikipedia”Secant method”, Papakonstantinou [83]). The secant method is initialized with twotrial values chosen from within the pre-established bounds: p1 = sin(θ0)/Vh(max)where θ0 = atan(X3/Xr) and p2 = p(max)(1 − tol), yielding ranges X1 and X2.The trial value for p1 was chosen as above since in a homogeneous isotropic mediumthat is the solution, iteration is not needed and the convergence test will be satisfied.The update for k + 1 using the secant method is thenpk+1 = pk + (X− Xk)/(Xk − Xk−1)( pk − pk−1). (3.30)Certain fail-safe conditions have been found necessary in practice. First of all, thevalue p(opt) = p that produced the minimum |Xk− X| is kept. If a secant iterationproduces pk+1 > p(max) then the next trial value is set to pk+1 = ( p(opt) +p(max))/2. If |Xk − X| > |Xk−1 − X| then divergence is occurring and the nexttrial value is set half-way to the best value pk+1 = ( pk + p(opt))/2. If there hasbeen no change in the range then iterations have stagnated and a new value isassigned: pk+1 = ( pk + p(max))/2. These last two fail-safe conditions are rarelyneeded but have been found necessary for exceptionally thin, fast layers where thesource or receiver are close to an interface and the range increment loses numericalprecision. The final fail-safe is to limit the number of iterations itmax, returning aflag of no convergence, and in any case p(opt) is carried forward to compute dynamicattributes.3.6.1 Phase and group velocityAt this point it is noted that the phase velocity, c, can be evaluated at any pointalong the ray, for example at the source, cS . This is needed for the strain Greenfunction given in chapter 2, equation (2.8). Phase velocity is simplyc =1√p21 + p23. (3.31)483.7. Ray shootingSimilarly group velocity is simplyV =√V 21 + V23 , (3.32)also needed at source and receiver in the strain Green function.3.7 Ray shootingTwo-point ray tracing may not be needed, so a ray shooting option is included.This is signified by setting the maximum number of two-point iterations itmax tounity, in which case the horizontal slowness rayp normally output by two-point raytracing is now used as input and a ray is traced according to that ray parameter.Ray attributes are returned as before for two-point ray tracing with the exceptionthat the target range is now the range returned.3.8 Dynamic ray attributes3.8.1 Geometrical spreadingThe general form of the spreading function is derived from the transport equation(Chapman [19] 5.4.9) and is given in Chapman [19] as 5.4.19. Ursin and Hokstad(2003) [117] studied the special case of layered VTI media. The formula for spreadingisSR =∣∣∣kˆ · VˆR∣∣∣∣∣∣kˆ · Vˆ S∣∣∣∣∣∣∣∂x1∂p1∂x2∂p2∣∣∣∣ . (3.33)where Vˆ is the group velocity vector, meaning that the terms |kˆ · Vˆ| are the cosineof the ray angle with respect to vertical (kˆ). In a 1D VTI medium the out-of-planeterm, ∂x2/∂p2 = X/ p (Chapman [19] 5.7.15) so (3.33) becomesSR =∣∣∣kˆ · VˆR∣∣∣∣∣∣kˆ · Vˆ S∣∣∣∣∣∣∣Xp∂x1∂p1∣∣∣∣ . (3.34)The qP-qSv in-plane term ∂x1∂p1 in any layer l is∂(∆xl)∂p1=∆x1p1(1 +p1p3∆x1∆x3)+4p3∆x3A11A44p21∆x23 −Ap1p3∆x1∆x3 +A33A44p23∆x21Ap21 + 2A33A44p23 −A33 −A44. (3.35)The equivalent expression for the Sh ray is∂(∆xl)∂p1=A66∆x3A244p23=∆x1p1(1 +p1p3∆x1∆x3). (3.36)493.8. Dynamic ray attributesThese “dxdp” terms are accumulated in a loop over layers, then the product(3.34) is formed.3.8.2 Purely vertical and horizontal raysPurely vertical and purely horizontal rays require special handling. As x→ 0 p1 → 0and hence X/ p is undefined. In this case l’Hopital’s rule is used to obtain∂x2∂p2=x1p1=Xp→∂x1∂p1(3.37)and so ∂(∆x1)/∂p1 replaces X/ p in (3.33) and ∂(∆x1)/∂p1 is effectively squared.Spreading for purely horizontal qP and qSv is complicated, the final result comesfrom Chapman (pers.comm., 2011, “Horizontal ray”), while the needed terms listedbelow come from Chapman [20]. The spreading S for qP or qSv rays is given bySqP−qSv(horz) =−X22V 2h p1∂(V 23 )∂p1(3.38)∂(V 23 )∂p1=(d3D)2 [∂(p23)∂p1+ 2p23(1d3∂d3∂p1−1D∂D∂p1)]D = (A11 +A44) p2 + (A33 +A44)p23 − 2d3 = A p2 + 2A33A44p23 −A33 −A44∂d3∂p1= 2Ap1 + aA33A44∂(p23)∂p1∂(p23)∂p1=(∂B∂p1∓12C1/2∂C∂p1)/(2A33A44)∂C∂p1= 2B∂B∂p1+ 8A33a44(A11 +A44 − 2A11A44 p2)p1∂B∂p1= − 2Ap1C = B2 − 4A33A44(A11 p2 − 1)(A44 p2 − 1)B = A33 +A44 −A p2Purely horizontal Sh is handled specially with spreadingSSH(horz) = X2A44. (3.39)3.8.3 Effective medium spreadingGeometrical spreading is not always a stable computation. For example a ray thatcontains a long horizontal segment will have a large change in X for a vanishing503.8. Dynamic ray attributeschange in p1 and ∂X∂p blows up. Consequently alternatives need to be considered.An effective ray length spreading (Chapman [19] 5.4.25) is given asS1/2 = RcS . (3.40)It is equivalent to (3.33) in the homogeneous isotropic case with R the geometricalray length; cS is the phase velocity at the source. (Notice the intuition about unitsthat effective ray length spreading provides — (3.40) has units L2T−1 or area/time— the amount the wavefront increases in area per unit time.) Given the intuitiveform of (3.40) one might consider another approximate form for spreading. For ex-ample, the actual ray length (the sum of ray segment lengths in layers) could be usedfor R together with the phase velocity at the source. This simple approach has beenimplemented as the default option as it has been found to produce stable inversionswhere sub–horizontal rays are involved. Another possibility might be to replace thephase velocity at the source with an average or rms velocity for the ray. Yet anotherapproach would use exact VTI spreading but in a new homogeneous medium derivedfrom the layered medium. Such an approximation can be obtained using Backusaveraging [6] for a stack of isotropic or VTI layers and would be generally differentfor each receiver. While the theory for calculating long-wavelength effective mediaapplies strictly for layer thicknesses of less than 1/10th of a wavelength, it serves thepurpose of replacing the layered VTI medium with an average VTI medium for thepurposes of a stable spreading computation. It is noted that Schoenberg and Muir(1989) [99] presented the effective medium results for arbitrary anisotropy. Otherpublications have used effective models where the true moveout or arrival time curveis approximated by a parametric curve, for example a hyperbola plus a quartic term,for efficient “no ray tracing” time computations. Aldridge et al. [2] used a shiftedhyperbola approach for downhole microseismic event location; Leaney (2000) [54]used a non-hyperbolic approach for layered VTI models in a walkaway VSP velocityinversion; Leaney and Hope [65] included effective medium non-hyperbolic spreadingin long-offset AVO processing of surface seismic data.3.8.4 PolarizationsWhile the eigenvalues of the Christoffel matrix correspond to phase velocity squared(3.2), the eigenvectors in (3.2) are the polarization vectors. The analytical resultfor qP − qSv in a VTI medium isgˆ = ±sgn2apˆ1pˆ30(A33 −A44)pˆ23 − (A11 −A44)pˆ21±[((A33 −A44)pˆ23 − (A11 −A44)pˆ21)2+ 4pˆ21pˆ23a2]1/2, (3.41)where sgn(x) = x/|x|. As discussed later, inside the ray tracing routines all raysare considered as downward, with p3 positive. Upward travelling rays and their513.8. Dynamic ray attributesassociated attributes are handled outside the ray tracing routines so that, for ex-ample, only the 3-component would become negative for qP polarization computedusing (3.41). ForqSv, the 3–component is continuous across the equator, requiredfor smoothly varying radiation patterns. The Sh polarization vector lies in the hor-izontal plane and since it depends only on the ray azimuth it is set outside the raytracing routines, and is discussed later.3.8.5 Transmission and reflection coefficientsReflection and transmission coefficients at an interface are commonly referred toas “Zoeppritz” coefficients due to the first solution of the problem for isotropicmedia (Zoeppritz, 1919 [127]). Considering the case of a welded interface as weexpect given the confining pressures deep within the earth, the boundary conditionsare that (a) particle velocity should be continuous and (b) the normal componentof traction should be continuous (Chapman [19] section 4.3). Zoeppritz (1919)[127] considered coefficients normalized with respect to particle displacement butother normalizations are possible. The coefficients used here are normalized withrespect to anisotropic eigenvectors (Chapman [19] 6.3.29), consistent with the raytheory Green function (2.8). The theory for anisotropic coefficients is describedin Chapman [19] 6.3.2 and simplifications for the case of transverse isotropy in6.3.4. The routines called by ”tizoep” in the VTI ray tracer follow the approach ofSchoenberg and Protazio (1992) [100] for media with up-down symmetry, developedoriginally for plane-wave TI reflectivity modelling (Folstad and Schoenberg, 1992[43]). The routines called by tizoep represent the only externally coded routines usedin this thesis. Those routines have been tested against another external isotropicZoeppritz code, the Crewes applet (crewes.org applet) and also recently developedcode for general anisotropic media (Chapman [21]).The ray tracer accumulates the product of complex transmission and reflectioncoefficients as needed and returns the total transmission coefficient for the selectedray code as well as the complex reflection coefficient for a given indexed interface inthe model (see section 3.9) that was passed to the ray tracing routine.3.8.6 Effective or average Q along a rayA fully anisotropic treatment of Q is complicated, requiring as many Q valuesas anisotropic moduli (five for the VTI case). Here a practical approach is takenwherein three different Q values are specifiable, one each for Qp, QSH and QSV ineach layer. The ray tracer then accumulates the effective Q along the ray accordingtoTQeff=∑ltlQl(3.42)523.9. 1D model preparation and unfolding — preraytrcwhere T =∑l tl is the total time along the ray and tl is the time in layer l for thespecified ray type in that layer. The ray tracer returns the value Qeff for the raywhich may be used in the modified scalar propagation terms of the Green function(equation 2.23).3.9 1D model preparation and unfolding — preraytrcThe low level ray tracing routines raytrc for qP − qSv and raytrc sh for Sh tracerays in a single direction – the direction of increasing depth, from the top of layer#1 to the bottom of layer #N . To handle the cases of rays with reflections, rayswith mode conversions or upgoing rays, a routine “preraytrc” is called to extract asecond 1D model given the source depth, receiver depth and reflector depth. If thereflector depth is equal to the receiver depth then there is no reflection and incidentphase angle is set equal to the phase angle at the receiver. The layer index of theinterface closest to the reflector depth is returned by preraytrc and the propertiesof the layer beneath the reflector are stored in an appended layer Nlay + 1 for lateruse by the ray tracer in a call to tizoep.Also set in preraytrc is an integer array to indicate the ray code in each layer - qP(1) or qSv (−1) – that is passed to raytrc. Certain flags are set to facilitate settingmore common ray code types, namely the case of surface (water) layer multipleswhere mode indicates the number mode of surface reflections; the case of a reflectedP − sv mode conversion (mode = 2); pure qSv propagation (mode = 3); P −SV − sv mode conversion with reflection (mode = 4); P − SV − sv− p transmittedmode conversion at the same depth on downward and upward legs with shear–shearreflection (mode = 5). If mode = 4 or 5 then the depth of mode conversion is savedin the location of the layer depth z(nlay+2) and the nearest model interface to thisdepth is found.3.9.1 Near–field ray–straightening at source and receiverIt is well known that very thin, high velocity layers pose problems for wide angle(sub-horizontal) rays. Washbourne et al. (2008) [126] propose a novel techniquethat they call “wave tracing” which straightens a ray to make it better approximatethe travel path of lower frequency waves. In their approach the minimum timeFermat ray returned by two point ray tracing is straightened until the travel time hasincreased by an amount related to the Fresnel zone for a given frequency. Adaptingthis idea, a Fresnel distance criterion may be used to modify the model whenever thelayer thickness of the first and/or last layers (containing the source and receiver)is less than a Fresnel distance. If this criterion is met the properties of the endlayers are modified to be a distance-weighted average of the end layer and the layernext to it. This has the effect of turning the last contrasts of the ray path into a533.10. Head wave ray tracinggradient, effectively smoothing the ray attributes with depth and range, particularlyend point ray attributes such as polarization angles. This has been found to have avery positive effect on focal sphere take-off angle sampling and concomitant inversioncondition number when the source is near a boundary.A Fresnel frequency parameter ffresnel is set to 100Hz and the vertical P velocityis used in a Fresnel distance zfresnel = vp(0)/ffresnel. Then parameter ξ is simplymodified if layer thickness ∆zl < zfresnel with α = ∆zl/zfresnel asξfresnel = (1− α)ξl−1 + αξl (3.43)where parameter ξ represents one of five VTI parameters.3.10 Head wave ray tracingSince high velocity carbonates are common above and below unconventional reser-voirs it is of interest to consider the existence of head wave arrivals. A head wavemay exist if there is a horizontal velocity in the model that exceeds that of all layersbetween the depths of source and receiver. There may be more than one, and theymay have their horizontal ray segment either above or below the depth zone definedby source and receiver. In this case the head wave with the earliest arrival time ischosen and returned.The algorithm for head wave ray tracing makes use of the results of reflectionray shooting, where a ray is shot and reflected from candidate head wave interfacesusing a horizontal slowness determined from the maximum horizontal velocity inthe reflected ray model. The reflected ray model is constructed only if there existsa layer with a horizontal velocity Vh = 1/pmax greater than that of the modelbetween source and receiver depths. If the range of the reflected ray is less thanthe target range (source–receiver horizontal distance), then a head wave exists andits horizontal ray segment, xhor, must be added to the reflected ray along with itshorizontal travel time, thor = pmaxxhor.3.10.1 Dynamic properties of head wavesAki and Richards [1] discuss the dynamic properties of head waves in Box 6.4 p.205which is based on a high frequency asymptotic approximate result given as equation6.26. A similar expression can be found in Chapman [19] as eq. 9.1.52. Head waveamplitude attenuates with distance as r−1/2L−3/2 where r is the total horizontalrange and L is the length of the horizontal segment along the interface. It is notedthat this behaves like r−2 for r >> |z + z0| — a mostly horizontal ray path — adecay that is much stronger than r as for a direct wave in a homogeneous isotropicmedium. It is also pointed out that this blows up as L → 0. Given these twopoints several approximate spreading formulas have been tried, with the goal being543.11. Triplicating qSvto always be defined and greater than the homo–iso case. By analogy with theeffective ray length (3.40) and given that the head wave ray length will alwaysbe longer than the homo–iso ray, the total head wave ray length ||xhead|| is used,together with the horizontal velocity of the horizontal segment so thatS1/2head = ||xhead||/p1(head). (3.44)Phase and group velocities, polarizations and effective Q for head wave rays aredetermined as described above for direct and reflected rays.A second feature of head waves is that they should exhibit a factor i/ω corre-sponding to a time integration, meaning that they should be lower frequency andcontain a 90-degree phase rotation. Part of this effect may be handled by rotatingthe complex transmission coefficient; the 1/ω effect may be handled external to theray tracer at the wavelet convolution stage.3.11 Triplicating qSvAs seen in Figure 3.3 strong anisotropy may give rise to multiple qSv arrivals in agiven ray (group) direction. For oblique directions this will pose obvious problemsfor 2-point ray tracing, which seeks a single ray. Mathematically what’s happeningis that roots of the range equation (3.23) are repeated so that more than one phasedirection gives rise to the same group direction with different group velocities.3.11.1 Triplication conditionSchoenberg and Helbig (1997) [98] provide the condition for triplication in a VTImedium as the positive root of a cubic, and Schoenberg and Daley [96] investigatefurther details of qSv triplications, showing that, for the ubiquitous situation whereA11 > A33 triplication occurs closer to the 3 − axis than the 1 − axis. Defining aparameter g1g1 ≡4A55(A11 +A33)/2(A11 −A55)(A33 −A55)> 0 (3.45)Schoenberg and Helbig [98] provide the following cubic equation in trip3trip + (2g1 + 3)2trip + g1(g1 + 2)trip − g21A11A33[(A11 +A33)/2]2 = 0, (3.46)the positive root of which provides a test for triplication: if the anellipticity of themedium exceeds this value; if Ea > trip, then there is triplication. Thomsen andDellinger [110] provide a simpler, approximate test criterion for off-axis or obliquetriplication.553.11. Triplicating qSv3.11.2 Solution selectionWhat to do about the possibility of qSv triplication ? Since 2-point ray tracing seeksto find the phase direction that will produce a group direction such that the ray willland on the receiver, and since a triplicating medium means that there are multiplevalues of p that will produce the same group direction, 2-point ray tracing may (anddoes) jump between solutions (not shown). Figure 3.5 illustrates the problem. Forthis example a source was placed in the medium of figure 3.3 and a line of receiverswas placed 150m above it. Shown is p versus group angle and X vs p.Figure 3.5: qSv horizontal phase slowness versus group angle for the triplicatingshale of figure 3.3 (left) and X vs p.Since there can be only one group angle that will connect source and receiverthere are three choices for p. For practical applications it is desirable to have asingle qSv arrival, so the question becomes, which should one choose? The earliestarrival would seem an obvious choice but this will lead to a discontinuous wave-front and possibly unstable amplitudes (singularity in spreading) near the cusps oftriplication. The innermost wavefront surface, the slowest arrival, is a choice thatprovides a continuous wavefront and stable spreading. Making this choice, the fol-lowing strategy is adopted to find the ray corresponding to this arrival. First, atest is made on the ray signature in every layer through which the ray passes in themodel — the ray must contain an Sv leg in a layer. Then for all layers containingSv the triplication test is made and if any layer is a triplicating medium then a newroutine is called to determine the ray parameter corresponding to the maximum timearrival. In this routine rays are shot through the model for a densely sampled rangeof ray parameters (< pmax). Figure 3.5 shows the X(p) computation for the rayexample described above. Triplication occurs where the derivative of this function,∂X/∂p < 0. The target range is subtracted from this range function and zeros are563.11. Triplicating qSvdetected. If there are multiple zeroes (i.e. three) the zero is selected correspondingto the maximum time. The p-value is interpolated. This ray parameter is then usedfor dynamic ray computations.3.11.3 qP+qSv wavefrontsWhat about amplitudes ? Ray theory breaks down at the cusps in figure 3.3,dX/dp = 0, spreading becomes singular and amplitudes become infinite. Finitefrequency (e.g. finite difference) simulation is required (with carefully selected pa-rameters) to get amplitudes right. Extensions to ray theory (Maslov theory, Chap-man [19] Ch. 10) provide an approximate solution. Between the cusps on thereverted branch theory calls for a Hilbert transform. Using a causal Brune pulseb(t) = ω2c te−tωcH(t) its Hilbert transform isb¯(t) =ωcpi−ω2c te−ωctpiEi(ωct) (3.47)which tends to zero as t → ±∞. Ei(x) is the exponential integral function whichmay be computed for positive arguments using a power series for small x and anasymptotic series for large x. For negative arguments a different approach is needed.Ei(−x) = −E1(x) − ipi/x. Numerical recipes provides the functions ei(x) and ex-pint(n,x) which with n=1 is used to evaluate E1 and Ei(−x), taking the real part.It is noted that while the Brune pulse is causal its Hilbert transform is not, andthat to preserve causality cancellations occur with the tail of the qP arrival. While(3.47) tends to zero, this happens slowly, so practically a taper is used. For thecomputations shown in Figure 3.7 a taper from |t| = 5/ωc to |t| = 10/ωc has beenapplied. Figure 3.6 shows a 150Hz Brune pulse and its Hilbert transform. Figure 3.7shows qP and qSv wavefronts for a Brune pulse with corner frequency fc = 100Hz,without and with the Hilbert transform on the reverted branch. There is no di-rectionality at the source due to radiation pattern, only geometrical spreading isincluded to see how this causes amplitudes to vary with direction. The variationalong the wavefronts is significant and the cusp singularities on the qSv arrival areevident.Considering the maximum time selection criterion discussed in section 3.11.2, fig-ure 3.8 shows this pruned qSv wavefront. While the maximum time arrival changesabruptly in amplitude, it remains of consistent phase and does not include theanomalous behaviour near the cusps. Furthermore, the maximum time arrival actu-ally is the maximum amplitude arrival over most of the group angles. These anglescorrespond to more horizontal rays and so, for the downhole geometry a discontin-uous jump in qSv amplitudes is not expected.Is all the above trouble worth the worry? Maybe not. Thick layers of triplicatingmaterial are probably rare, and such material must dominate the summation of(3.35) in order for dX/dp < 0 and multiple zeros of the range function to occur.573.11. Triplicating qSvFigure 3.6: A 150 Hz Brune pulse and its Hilbert transform (dashed).Figure 3.9 shows a case where the source is in an isotropic layer and embedded inthe triplicating medium. The qSv wavefront has just begun triplicating at the timeof this snapshot (t = 80ms).Figure 3.10 shows qP and qSv wavefronts for a layered VTI model used for com-mercial processing. While a certain amount of focusing occurs, none of the layersis triplicating and so no complications arise. A clear example of qSv triplicationremains to be seen in real microseismic waveform data, although given some mea-surements from walkaway VSPs (e.g. Miller et al. [80], Leaney [58]) and knowing ofsome models calibrated using microseismic perforation shot times (unpublished), itmust occur. The observation of cusps in qSv waveform arrivals has been reportedin offset VSP data (Hake et al. [46]). The problem of handling qSv triplicationhas been dealt with for the sake of completeness, intellectual curiosity and under-standing, more than because it has been observed as a problem in real data. Itmay be observed in the near future, because fracture monitoring is beginning in thefield where the measurement of figure 3.3 was made, so its handling needs to beunderstood.Figures 3.9 and 3.10 illustrate an important aspect of exact VTI layered ray-tracing - the shear rays, qSv and Sh, are very different. This is due to the verydifferent angle-dependent response that they have to anisotropy. A weak anisotropicapproximation would not show such a ray path difference, as both shear rays wouldhave the same isotropic path, and perturbations about the same ray path would beused to compute ray attributes.583.11. Triplicating qSvFigure 3.7: qP and qSv wavefronts at 100 ms for the triplicating shale of figure 3.3including geometrical spreading using a causal Brune pulse with corner frequencyfc = 100 Hz. Red is high amplitude. Left: without Hilbert transform on the revertedbranch; right: with Hilbert transform on the reverted branch.Figure 3.8: Left: qP and qSv wavefronts at 100ms for the triplicating shale offigure 3.3 using a causal Brune pulse with corner frequency fc = 100hz and usingthe maximum time selection strategy. Right: Sh wavefronts shown using the samecolor amplitude scale but different horizontal and vertical scales.593.11. Triplicating qSvFigure 3.9: Layered model, rays, times and qP + qSv wavefronts for a source in anisotropic layer embedded in the triplicating medium.Figure 3.10: Layered model, PHV rays, PHV times and qP + qSv wavefronts at100ms for a layered VTI model used in commercial processing.60Chapter 4Moment Tensor InversionWhile many approaches for moment tensor inversion have been published, to myknowledge only the work of Ro¨ssler and colleagues (Ro¨ssler et al. [93], Ro¨ssler et al.[92]) incorporated anisotropy in the Green function, and then only in a limited way.In Ro¨ssler et al. [93] the theory developed by Ps˘enc˘´ıkand Teles [87] is used togetherwith the inversion algorithm of Dahm [28], Dahm et al. [29], but approximationsare made and the polarization vectors are not for the anisotropic medium. Theinversion is linearized and a multitude of starting models is used to initiate theinversion. In Ro¨ssler et al. [92] the anisotropic ray tracer ANRAY (Ps˘enc˘´ık [86])is used in a time–domain inversion for the six components of the displacement–discontinuity or point dislocation potency tensor. As mentioned previously thisis an incomplete description of a moment tensor stress–glut source. The modelsused were either weakly anisotropic with vertical heterogeneity or homogeneousand transversely isotropic. The inversion was applied to a swarm of earthquakesfrom West Bohemia. Moment tensor inversion applied to hydraulic stimulation inthe oilfield has been in the time domain and using isotropic models (Urbancic et al.[116], Trifu et al. [113], Nolen-Hoeksema and Ruff [82]) even if travel times have beencomputed in an effective VTI medium (Eisner et al. [39]). Song and Toksoz [103]described an isotropic waveform inversion of downhole data, reporting the recoveryof non–double couple mechanisms. They also claimed the complete moment tensorcould be inverted from a single vertical array using near–field measurements. Awaveform inversion was also developed utilizing sparse representation theory forefficiency Rodriguez et al. [91], but also in isotropic media.Many moment tensor inversion schemes have appeared in the literature, imple-mented in the time domain or frequency domain, usually assuming a common formof source function for all moment tensor elements. Julian et al. [50] (and referencestherein) provides an overview in the context of recovering non–DC components ofthe moment tensor.In this chapter a linear, waveform inversion for the moment tensor in the fre-quency domain is described. It combines the efficiency of ray tracing for the Greenfunctions with the potential accuracy of waveform fitting. Returned by the inversionis the moment tensor, each element possessing independent variation with frequencyor time. To satisfy a condition that all elements of the moment tensor should havethe same time dependence, complex principal components analysis is used to extract614.1. Overviewa common signal as the source function from which corner frequency and scalar mo-ment are recovered. Posterior uncertainties are estimated with multivariate normalsampling using the Cholesky decomposition of the model covariance matrix. Theinversion is illustrated on two synthetic data sets using the forward computations de-scribed in the previous two chapters. One of the synthetics illustrates a pathologicalgeometry limitation common in downhole microseismic monitoring.4.1 OverviewRecall equation (2.25) from chapter 2 for vector recordings at frequency ω, restatedfor convenience here:v(ω,x S,xj) ≈ rj(ω)[rays∑kgˆkjGQk (ω,x S,xj)E˜kj]· −iωM˜(ω), (4.1)with all terms now explained in the previous two chapters. The inverse of thereceiver response term will be required to transform recordings to displacementunits, but this will be done later so that inversion and data reconstruction canbe applied in native recorded units. Likewise the factor −iω is absorbed into thefrequency-dependent moment tensor, leading to an inversion for a new parametervector m(ω) = −iωr(ω)M˜(ω). Equation (4.1) is therefore re–written asd(ω,x S,xj) =[rays∑kgˆkjGQk (ω,x S,xj)E˜kj]·m(ω). (4.2)Equation (4.2) says that the recorded three component data at frequency ω areassumed to be equal to a summation of arrivals approximated by ray theory as asum over rays. The recorded data are assumed to be oriented in the ENU coordinatesystem which in practice is usually accomplished via a prior processing step. The lo-cation of the source is also assumed known, also from a previous processing step andthe model too is known from previous work. These steps will be discussed when realdata are considered. As a matter of order, the components of the data are orderedE1, N1, U1;E2, N2, U2; ... as are the components of the receiver polarization vectorin (4.1) for ray k: gˆj = (gE , gN , gU )Tj . Thus the data vector contains 3M elementswhere M is the number of receivers in the network. With this ordering in mind,and remembering the sum over rays for each component of receiver polarization,equation (4.2) is written in matrix–vector form asd(ω) = G(ω)m(ω), (4.3)where the dependence of G on ω enters in the travel time in the exponential (2.22)and from the inclusion of absorption and dispersion due to Q (2.23–2.24). Impor-tantly, it is noted that at this point no assumption has been made about the source624.2. Weighted, damped least–squares inversionfunction and each element of m is free to have a different source function. Forcingm(ω) = Ms′(ω) is a nonlinear problem so the issue of estimating s′(ω) from m(ω)is saved for later. The reason for using s′(ω) rather than s(ω) at this point is be-cause the estimate is in recorded data units. s(ω) will be recovered from s′(ω) whenthe transformation to particle displacements is applied. For now the linear inverseproblem to solve for the six elements of m at each frequency is addressed.4.2 Weighted, damped least–squares inversionLinear inverse theory provides the damped least–squares solution to problems of theform d = Gm (e.g. Menke [77])m =[GTG + 2I]−1GTd (4.4)where is a regularization or damping parameter, sometimes called the Tikhonovregularization parameter (e.g. Aster et al. [5]).4.2.1 Interpretation of equation (4.4) for this applicationSince the above inversion is happening at each frequency, it is understood that theT transpose symbol means complex–conjugate transpose or Hermitian transpose,where the imaginary part is multiplied by −1. This effectively reverses the signmultiplying time in the argument of the complex exponential in (2.22), leading tothe time–reversal view of this inversion (Leaney [57]). Kawakatsu and Montagner[51] also viewed moment tensor inversion in the context of time–reversal imaging.The scalar amplitude terms in G and the subsequent normalization by GTG aftercorrelation of the data with the transposed strain Green function can also be inter-preted as elastic imaging with an approximate (ray theory) one–way wave equation,but instead of determining, for example, elastic contrasts as with two–way propaga-tion and elastic scattering, the components of the moment tensor are determined.With this insight the above can be viewed as a deconvolution imaging condition(Sacchi, pers.comm., 2008). Another aspect of the least–squares inversion of (4.4)is that it involves vector beam forming, being essentially a weighted stack, and soattenuates random noise proportional to√N where N is the number of channels ortraces. The difference from conventional beam–forming is that the time trajectoriesthrough the data are not simple linear or parameterized curves but rather arbitraryaccording to the modeled arrival times, with polarizations connecting the vectordata, essentially tripling the length of the array as shown for the anisotropic VSPwavefield separation problem (e.g. Leaney [56]).634.3. Posterior uncertainties4.2.2 Estimation and covarianceThe damped least–squares result (4.4) can also be derived as a special case of themore general Bayesian point of view as considered by Tarantola ([105]). If the dataerrors have a multivariate normal distribution and the prior distribution for themodel parameters is also multivariate normal then the posterior probability densityis also Gaussian (Tarantola [105], eq. 1.103) with mean (Tarantola [105], eq. 1.106)m = (GTC−1D G + C˜−1M )−1(GTC−1D dobs + C˜−1M mprior) (4.5)and covarianceCM = (GTC−1D G + C˜−1M )−1. (4.6)For the case of diagonal data and model covariance matrices, CD = σ2I and C˜M =α2I, respectively, the minimization problem becomes equivalent to zero–th orderTikhonov regularization or damped least–squares with = α/σ. In the limit whereα is extremely large (total ignorance about the parameters in m), all terms involvingC˜M disappear and the above expressions becomem = (GTC−1D G + I)−1GTC−1D dobs (4.7)andCM = (GTC−1D G + I)−1 (4.8)which are the expressions that will be used for inversion. Note that a very smallregularization term has been retained, admitting that this will bias slightly thesolution but may be needed for stability. This is the weighted damped least–squaressolution (Menke, [76],[77],[78]).As shown by Tarantola ([105],[106],[107]) the covariance matrix CD should reallybe viewed as representing observational uncertainties rather than just data uncer-tainties, so (Tarantola [105], eq. 1.101)CD = Cd + CT (4.9)where Cd is the data covariance matrix and CT represents errors due to an inexacttheory — “modelization” uncertainties. In what follows Cd is taken as diagonalbut not constant. In other words the data error (noise) is assumed uncorrelatedbut variable from receiver to receiver. CT is a term that may be added whereverray theory is thought to break down, for example for nearly horizontal rays wheretrapped modes may be present in the data.4.3 Posterior uncertaintiesApproximate parameter standard deviations can be obtained from the square root ofthe diagonal entries of CM , since those elements are the variances of the posterior644.3. Posterior uncertaintiesmultivariate normal distribution of model parameters (e.g. Press et al. [85] eq.15.4.15). The standard deviation for moment tensor parameter i is:σi =√CMii . (4.10)If CM were purely diagonal the parameters would be linearly independent andresolved by the experimental set–up. This is never the case however, and in practiceCM contains non–zero off–diagonal entries (parameters are coupled and co–vary)as a result of deficiencies in experimental design. In the case of the moment tensorinverse problem, it is the result of insufficient sampling of the focal sphere by raysleaving the source, connecting to receivers in the available network through themodel. Moment tensor parameter resolution from borehole data has been studiednotably by Vavryc˘uk [122], Eaton [37] and Rodriguez et al. [90].In practice data covariance–weighting is used for parameter estimation but forparameter uncertainties the data misfit is used and the data covariance is diagonal.This approach effectively includes modeling errors in the parameter uncertainties.Thus the posterior covariance matrix will be computed usingCM = σ2d(GTG + I)−1 (4.11)where now G contains no weighting and the prior data variances are replaced withthe data misfit variance (e.g. Aster et al. [5], ch.2, page 27):σ2d =1NN∑i(di −Gim)2 (4.12)where N is the number of data, in this case three times the number of geophonecomponents.It is of interest when viewing results from multiple events to compute a singlenumber that is representative of parameter uncertainty. For this purpose the max-imum standard deviation is used, normalized by the scalar moment. The relativemoment tensor uncertainty metric is defined asσm = max(σi)/M. (4.13)4.3.1 Multivariate normal samplingIn what follows a simple metric is used for multiple events as defined above (4.13),but a more complete approach is used for single events. As described by Aster et al.[5], having determined the expected parameter values (4.5), also called the MAPor Maximum A Posteriori values, mMAP , and the associated posterior covariancematrix, which together fully describe the posterior multivariate normal (MVN) dis-tribution, it is straightforward to generate model realizations that honour that MVN654.4. Implementationdistribution. The method (Aster et al. [5], section 11.3 and example B.10) uses theCholesky factorization of the covariance matrixCM = LLT (4.14)followed by generation of a random solution according tom = mMAP + Lx (4.15)where x is a vector of independent and normally distributed random numbers withzero mean and standard deviation one. More will be said on the topic of momenttensor (MT) uncertainties when it comes time to look at MT decomposition.4.4 ImplementationAt each frequency a complex SVD (singular value decomposition) routine is used tosolve d = Gm with inverse data covariance weighting, so the problem is rewrittenwith Gw = C−1/2G and dw = C−1/2d as dw = Gwm. The SVD of Gw is (e.g.Menke [78], eq. 7.32)Gw = UΛVT , (4.16)where the matrix U is a matrix of eigenvectors that span the data space and V is amatrix of eigenvectors that span model space. Λ is a diagonal matrix of non–negativesingular values, ordered by decreasing size. For the nominally over–determinedproblem dealt with here where there are more data than unknowns only the firstp = 6 entries of diag(Λ) are non–zero. The inversion formula for the weighteddamped least–squares problem using the SVD ismest = V{Λ[Λ2 + 2I]−1}UTdw (4.17)and the model covariance matrix isCM = V{Λ2[Λ2 + 2I]−2}VT (4.18)where it is recalled that the matrix of singular values, Λ, is diagonal. The conditionnumber, the ratio of maximum to minimum singular values, tells how much noisein the data is mapped to model parameters and is a useful “invertibility” metric. Itis used by some as the measure of the well–posedness of the inverse problem, withlower numbers being better, but one really needs to compute the complete modelcovariance matrix to understand how data noise and geometry map to parameterestimation. Written out in components the posterior covariance matrix for theM = 6 moment tensor elements readsCM = Cov(mi,mj) =M∑i=1(VjiVijλ2i(λ2i + 2)2). (4.19)664.4. Implementation4.4.1 Estimating the regularization parameter, In the damped least–squares solution of (4.4) the prediction error (data misfit) isminimized subject to an a priori constraint that the solution length should re-main small. The weight given to the solution length is determined by the size of .Increasing decreases the length or size of the model vector, stabilizing the underde-termined part but at the expense of fitting the data. If the data noise variance (thedata covariance matrix) is known a priori then the value of that reduces the misfitbelow some fraction of the noise can be determined by trying decreasing values for until that criterion is satisfied. This is known as the discrepancy principle (e.g.Farquharson and Oldenburg [40]). Two other approaches attempt to automate theprocess — the L–curve method and the generalized cross validation method.L–curve methodIn the L–curve method the data misfit, ‖d−Gm‖ and model size, ‖m‖ are plottedagainst each other for varying values of . The optimum trade–off is taken as thepoint of maximum curvature, computed numerically using finite differences. Efficientimplementations are not generally possible because of local extrema in the L–curvecurvature, so a rather exhaustive line search is needed. A disadvantage is thatseveral functional evaluations (inversions, at least 3) are required to compute thecurvature, so extra evaluations are required near the end points.Generalized cross–validation (GCV) methodThe GCV method is based on the leave–one–out principle (Wahba [124]), whereinthe prediction error is summed over each left out datum and minimized. A gooddescription is given in Aster et al. [5], section 5.7. The final, efficient, approximateformula for the GCV functional to be minimized for isGCV () =m‖Gm − d‖2Tr(I−GG−g)2, (4.20)where it is recognized that in the denominator is the (regularized) data resolutionmatrix which is evaluated using the SVD as GG−g = UUT and the superscript −grepresents the regularized generalized inverse, e.g. Menke [78], p. 69. GCV has theadvantage that no end–point evaluations are needed to compute a finite differencecurvature. Both of the above methods have been implemented and tend to yieldsimilar values for . In practice, both are evaluated at the same time and the lesserof the two values is used, but as estimating is expensive it is usually done onlyonce or a few times (for one or a few events). It can be a critical parameter whenthe experimental set–up makes parameter estimates sensitive to noise.674.5. Estimating s′(t) from m(t)4.4.2 Constrained solutionsA common linear constraint is to force the isotropic part of the moment tensor, itstrace, to be zero: M11 + M22 + M33 = 0. While this deviatoric constraint is easyto implement it is of little interest as it disallows a DD source and results in onlya source with DC+CLVD which seems to have little physical meaning. Dufumierand Rivera (1997) [35] use the deviatoric constraint as M33 = −M11−M22 becauseM33 is the most poorly resolved parameter using surface waves, but this and indeedall linear constraints are not useful for the microseismic problem. A more usefulconstraint would be to allow only a DC or slip source, or to permit only a DDsource, but these are nonlinear constraints so one might as well pose a fully nonlinearproblem. Something will be said later on nonlinear inversions, and some work hasbeen done in that regard, but they are not the focus of this thesis.4.5 Estimating s′(t) from m(t)Given the regularization parameter, the inversion proceeds, looping over frequencyto provide the estimated six–vector m(ω) and after Fourier transform m(t), butas mentioned previously this admits a possibly different source function for everyelement of m. There may be new and interesting applications for a time-dependentmoment tensor, but it is of immediate interest to extract a single source function,s′(t). In this section a novel approach to estimate a single source function s′(t) ors′(ω) and m (or ultimately M) is described.In section 4.1 it was stated that enforcing all six components of M to havethe same source function was a nonlinear constraint. Usually, moment tensor in-version either assumes a source function or a simple parameterized form for thesource function and solves a larger linear problem (e.g. Menke [78], section 13.2).The approach taken here is to allow each component of M to have a different andarbitrarily complicated source function, then the goal of obtaining a single sourcefunction is accomplished by solving a second, linear estimation problem. This isdone using principal component analysis or PCA.While the technique of PCA is rooted in statistical data analysis, the geophysicalslant on it will be followed here. The goal is to extract the signal from an ensemble ofsignals that has the most in common with all input signals. As derived, for example,in Freire and Ulrych [44], this is the signal constructed from the input signals thathas maximum variance or energy. It turns out to be derived from the eigen analysisof the covariance matrix of the input signals, and specifically it is the projection ofthe data onto the eigenvector corresponding to the largest eigenvalue. This is alsocalled the first principal component. The input data are reconstructed exactly if alleigenvectors are used and can be filtered by reconstructing with a subset.Many different variations on this theme have appeared in the geophysical litera-684.6. Estimating mˆ given s′(ω) and m(ω)ture (e.g. Trickett [112], Ulrych et al. [114])). The variation used here to extract s(ω)from m(ω) is that of complex PCA. Complex PCA can be applied frequency–by–frequency or in the time domain after inverse Fourier transform utilizing the complextrace. Both have been implemented but the time domain approach has been foundto yield superior results for noisy data, so the previously estimated moment signalsm(ω) are transformed to the time domain and s′(t) is estimated from complex tracesm(t). The reason for using complex traces is that both polarities need to be allowedfor each mj(t) and complex PCA handles this through the complex eigenvectors.In fact the phase rotation needed to produce the maximum variance stack is givenby the ratio of imaginary to real components of the eigenvectors. Complex, timedomain PCA will be used again later on real data to compute residual arrival staticswhere the arrival polarities may be opposite for different receivers. Complex PCAcomputation proceeds as follows.The covariance matrix C components are constructed from the complex signalsm(t):Cjk =∑imj(ti)m∗k(ti), (4.21)the eigen analysis of C is performed:C = UΛUT , (4.22)where U is a matrix of eigenvectors and Λ is a diagonal matrix of eigenvalues. Thesource function s(t) is now taken as the first principal component, constructed attime sample i ass′(ti) = UT1 m(ti). (4.23)Efficient routines are available that do not compute the entire eigen analysis ofa matrix, but can only compute the first eigenvector and eigenvalue.Having estimated the scaled source function s′(t) it is transformed back to thefrequency domain to estimate the six normalized real components of the momenttensor.4.6 Estimating mˆ given s′(ω) and m(ω)Having an estimate of s′(t) and hence s′(ω) the equation m(ω) = mˆs′(ω) is rec-ognized as a deterministic deconvolution problem. (The reason for the hat over mis that its elements are normalized, the magnitude of m now being absorbed intos′(ω).) It is easily solved for an estimate of mˆ at frequency i asmˆi = Real[sT (ωi)m(ωi)sT (ωi)s(ωi) + ](4.24)694.7. Scalar moment from s(ω) spectral fitwhere is a number added to avoid zeros in the source spectrum, being a smallnumber multiplied by the peak source power.This provides as many estimates of mˆ as there are frequencies used, and robustweighted estimation techniques can be considered to improve the reliability of afinal estimate for mˆ in the presence of noise. One technique that has been founduseful is that of the weighted median estimate, where the weight at each frequencyis the power |s′(ω)|2. The weighted median then proceeds by replicating estimatesmi according to an integer normalized weight index. The integer normalized weightindex is calculated using a number nw = |s′(ωi)|/max(|s′(ω)|) ∗ Nf/20 where i isthe frequency index and Nf is the number of frequencies. nw is the number of timesthe j − th moment tensor element at frequency i, mj(ωi), is repeated. The medianof this new sequence is the weighted median.4.7 Scalar moment from s(ω) spectral fitThe above procedure provides the normalized elements of the moment tensor, mˆor M̂, and as mentioned, s′(ω) contains the magnitude of m, ‖m‖. Recalling thesource function discussion in 2.3 wherein the equation for the spectrum of the Brunepulse was given (2.14), the amplitude spectrum |s(ω)| will be fit with a function ofthe form (Beresnev and Atkinson [10], eq. 11)|s(ω)| =S0[1 + (ω/ωc)2](n+1)/2(4.25)which with n = 2 has two free parameters, S0 and fc = ωc/2pi. Prior to doing this,however, the estimated source function, s′(ω), has applied to it the transform fromrecorded units to particle displacement units. This is done by applying the inverse ofthe receiver response function, r(ω), which maps recorded units to particle velocity,followed by a temporal integration to displacement:s(ω) =1iω[r(ω)]−1s′(ω). (4.26)Removing the receiver response function is a deconvolution problem taken care ofby a routine based on proprietary parameters for Schlumberger’s GAC-D geophone.For some other geophones less information is available and a constant scale factoris used.The parameters S0 and fc are determined using a simple algorithm that callsa simplex routine, amoeba (Press et al. [85]). Note that under this source modelS0 = ‖M‖ = M0√2 where M0 is the scalar moment and the factor√2 comes fromthe historical definition of scalar moment for a double–couple. The final momenttensor isM = M̂‖M‖/√2 = M0M̂. (4.27)704.8. Synthetic examplesThe procedure of removing the receiver response function and integration to dis-placement is also applied to the estimates m(t) to recover the band–limited momentsignals, M˙(t). The above procedures will be illustrated on synthetic data.4.7.1 Compute timesThe major part of the inversion is linear in a loop over frequency with an additionalsmall step of complex eigen analysis of a 6 × 6 matrix to estimate the commonsource function. It takes generally less than two seconds on a standard laptop, andthis includes the small nonlinear inversions to compute arrival statics (discussedlater) and source spectrum fitting. The layered VTI ray tracing is also fast andcompute–times to store the Green functions for all receivers are included in the abovetime estimate. Times increase linearly with bandwidth and super–linearly with thenumber of receivers. Thanks to solving the (mostly) linear problem hundreds ofevents can be processed in minutes and parameters can be tested on event collectionsrather than on single events only.4.8 Synthetic examplesThe 3–array synthetic from the previous chapter will be used to illustrate the inver-sion described in the previous sections. Figure 4.1 shows again the synthesized eventas recorded by three linear receiver arrays with additive noise, it’s the same as figure2.10. Recall that the model was homogeneous VTI, that the source wavelet was ann = 2 pulse with corner frequency fc = 100Hz, geometrical spreading was on, thesource mechanism was composite with equal parts slip and opening and simulatedrecordings were in particle velocity. Shown in figure 4.1 are the reconstructed ENUdata after inversion using the forward equation (2.25). The beam–forming aspectof the inversion is recognized, as the additive random noise is attenuated in thereconstructed input data.4.8.1 Moment signals and source function recoveryThe estimated moment signals corresponding to the inversion of figure 4.1 are showntransformed to the time domain in figure 4.2 along with the estimated source func-tion using complex PCA (the real part of the first principal component). Alsoshown is the estimated source amplitude spectrum and the spectral fit from whichthe scalar moment, M0, the zero frequency asymptote, is obtained.714.8. Synthetic examplesFigure 4.1: Three array synthetics with noise as in figure 2.10 with from left toright groups of 36 input traces (3 arrays of 12) for ENU components followed by thedata reconstructed after inversion to determine moment tensor and source function.Colored arrival times correspondence is PHV=BRG.4.8.2 More noise and dampingThe inversion is shown in figure 4.3 for more additive noise. Now the beam–formingaspect of the inversion is illustrated more dramatically. Also shown in this figureis the inversion result using the GCV regularization parameter estimate ( = .27).While the reconstruction noise is reduced with damping, especially visible on thereconstructed U component, the arrival amplitudes are also noticeably reduced.Damping has allowed less noise through from data space into model space but hasreduced the size of the model vector, biasing the result away from the true solution.In practice regularization has been found necessary when the focal sphere is poorlysampled by the experimental set–up and interpretation of results should be done inconjunction with a careful review of this critical parameter.4.8.3 Quality of fit metricA commonly used metric to quantify the quality of fit normalizes the sum of squaredresiduals by the input data variance and is sometimes called “VR” or “variancereduction”: V R = 1 − ‖d − Gm‖/‖d‖. For a waveform inversion this may beimplemented in the time or frequency domain, but since in the present problemthere are only three (PHV) direct arrivals expected to dominate the waveforms itseems reasonable to weight the fit by arrival time windows so as not to give undueweight to noise. This is accomplished automatically by weighting the residuals by724.8. Synthetic examplesFigure 4.2: Left: Estimated moment signals mij and source function s′(t) as (upperleft) scaled particle velocity, (lower left) scaled particle displacements; right: sourcefunction amplitude spectrum |s(ω)| and spectrum fit.the 3C envelope (see 2.13) of the reconstructed data so that for noise realization kthe arrival–weighted fit quality metric isq(mk) = 1−∑iw(ti)[(d(ti)−Gm(ti))]2∑iw(ti)[d(ti)]2 . (4.28)Such a weighted fit metric is used as the objective function in an optional three–parameter nonlinear inversion for an updated event location, determined using thesimplex algorithm.4.8.4 Multiple noise realizationsFigure 4.4 shows the VR metric and the arrival–weighted VR metric versus noiserealization for additive Gaussian noise with standard deviation ranging from 0 to .2times the maximum amplitude of the noise–free synthetic. Noise standard deviationis incremented linearly in steps of .002. The arrival–weighted VR is seen to decreasemore linearly with noise.734.8. Synthetic examplesFigure 4.3: Same as figure 4.1 but with five times as much noise. Top: no damping,bottom: damping parameter from GCV, estimated at 0.27.Figure 4.5 shows inverted moment tensor components for increasing noise. Thesensitivity to noise is seen to be worst for m11, m22 and m33; this will be seenagain in the multivariate normal posterior sampling results. Also shown is a relativemoment error, defined as the maximum σi (eq. 4.10) divided by the maximummoment component. This provides a single, conservative number of moment tensoruncertainty.744.8. Synthetic examplesFigure 4.4: VR (red) and arrival–weighted VR (green) fit metrics for 100 increasingnoise realizations.Figure 4.5: Moment tensor components inverted for 100 increasing noise realizations(left): red=m11, green=m22, blue=m33, cyan=m23, magenta=m13, black=m12.Right: A relative error metric.754.8. Synthetic examples4.8.5 Multivariate posterior samplingUsing the full posterior covariance matrix from eq.(4.8) for the inversion of the noisydata result (undamped solution) shown in figure 4.3 (which corresponds to noise re-alization number 25 in figures 4.4 and 4.5) and the multivariate normal samplingscheme described in section 4.3.1, the posterior distribution of moment tensor el-ements are sampled. This approach captures the increase in uncertainty due tothe increase in noise, the relative variations in parameter uncertainty (the variationalong the diagonal of CM ) and parameter coupling or covariance (the off–diagonalelements of CM ). Figure 4.8 shows the model covariance matrix, normalized by theaverage of intermediate diagonal elements. From this it can be seen that the trace(M11,M22,M33) or isotropic part of the moment tensor is the least well resolved,and the diagonal elements are linearly dependent or coupled for these parameters asthe off–diagonal terms are significant. Figure 4.6 shows 5000 multi–variate normal(MVN) model realizations of moment tensor elements. To see the effect of includ-ing the off–diagonal elements of CM , figure 4.7 shows m22 vs. m11 for the case offull covariance sampling and diagonal–only sampling. As expected, including theoff–diagonal terms stretches the distribution due to the strong covariance betweenthese parameters. Figure 4.9 shows histograms computed from the sampled MVNdistribution honoring the full posterior covariance matrix.4.8.6 A pathological exampleAs mentioned previously, two vertical arrays of receivers do not provide sufficient fo-cal sphere coverage if the event lies in the plane of the two arrays. In this setting theforce couple normal to this plane cannot be resolved because its radiation pattern isnull towards the receivers. As an example, using vertical receiver arrays in a verticalN–S plane and a source near this plane means the m11 component is not resolved.Some regularization is required in this case as the problem is rank–deficient. Figure4.8 compares the covariance matrix for this case versus the well–resolved 3–arraysituation. Both matrices are normalized by the average of intermediate diagonalvalues. Notice the much larger condition number for the 2–array case. The actualvalue for CM (11) is somewhat arbitrary as it depends on the amount of damping se-lected for this ill–posed case. The variance can be reduced at the expense of biasingthe solution. Figure 4.9 compares multivariate normal histograms for the momenttensor components for the well–posed 3–array case and this case. As expected thepdf is much broader for the m11 component, indicating that it is an unresolvedparameter given the geometry and event position.764.8. Synthetic examplesFigure 4.6: 5000 multi–variate normal realizations of moment tensor componentshonouring the posterior covariance matrix. red=m11, green=m22, blue=m33,cyan=m23, magenta=m13, black=m12.774.8. Synthetic examplesFigure 4.7: 5000 multi–variate normal realizations of moment tensor componentsm22 vs. m11 honouring the full posterior covariance matrix (red) and honouringonly the diagonal terms (green).784.8. Synthetic examplesFigure 4.8: Model covariance matrices for the 3–array synthetic (top) and the 2–array synthetic with pathological event location (bottom), normalized by the averageof the two intermediate diagonal elements.4.8.7 Comments on moment tensor uncertaintyThe approaches taken here provide a look at posterior uncertainties in the estimatedmoment tensor using linear inverse theory. While errors due to the model and limi-tations of ray theory (CT ) will be partially captured in the case of real data by usingthe normalized squared error or arrival–weighted normalized squared error for theconstant diagonal of Cd, (σ2d), better approaches can be envisaged. Non–Gaussiannoise and model errors would be better quantified using a nonlinear inversion anda Monte Carlo sampling of free parameters which may include an event locationupdate (∆E,∆N,∆U), model parameter perturbations, receiver location perturba-tions, etc. Bootstrap approaches to confidence intervals (e.g. Menke [78]) can alsobe considered.But the recovery of the elements of the moment tensor is not the ultimate goalas the moment tensor is non–intuitive and difficult to interpret on its own. The truegoal is to recover information about the source geometry, derived from the momenttensor and which may be interpreted and used in other applications. Techniques formoment tensor interpretation as discussed in the next section attempt to extractmore intuitive parameters from M, but the question then becomes, what are theuncertainties in those parameters ? Error propagation from the covariance of M tothe covariances of new parameters based on the eigenvectors of M is not a simpleproblem(e.g.Du and Warpinski [34], Han et al. [47]). The approach taken and shownin the next chapter will be to look at decomposition parameters for each of the MVN794.8. Synthetic examplesFigure 4.9: Multivariate normal histograms from 5000 samples of moment tensorcomponents honouring the posterior covariance matrix for the 3–array (top) and thepathological 2–array case (bottom). red=m11, green=m22, blue=m33, cyan=m23,magenta=m13, black=m12.804.8. Synthetic examplesrealizations of M as shown in figure 4.9.Having estimated the elements of M with associated uncertainties the task nowis to try to recover the parameters that went into its construction as described inchapter 2.81Chapter 5Moment Tensor DecompositionIn this chapter a new moment tensor decomposition will be described, valid for gen-eral anisotropic media. Prior to that some basics are reviewed, as are the decompo-sitions of Hudson et al. [48] Vavryc˘uk [120],Vavryc˘uk [123]. Then the multivariatenormal realizations of moment tensors from the previous chapter are decomposedto produce histograms in source parameter space. Finally, new source visualizationsare shown for composite sources that make use of the new decomposition. Much ofthis chapter comes from Chapman and Leaney [22], wherein an extensive review ofmoment tensor decompositions is given.5.1 Principal decompositionCentral to moment tensor interpretation is the principal decomposition which de-composes M into its principal (eigen) vectors and values,M = B̂ΛB̂T (5.1)where B̂ is the matrix of orthonormal eigen vectors and Λ is a diagonal matrix ofeigenvalues. It is customary to order the eigen vectors according to their associatedeigenvalues and to denote them with T,N, P with ΛT ≥ ΛN ≥ ΛP . The choice ofletters comes from T =Tension, N =Neutral and P =Pressure. The matrix B̂ isthenB̂ = (T̂ N̂ P̂). (5.2)The principal vectors (5.2) form a complete, orthonormal basis and it is straightfor-ward to rewrite expansion (5.1) in terms of elementary sources, for exampleM = ΛT T̂T̂T + ΛNN̂N̂T + ΛP P̂P̂T , (5.3)but there are many other ways to write (5.3) and the decomposition into elemen-tary sources is completely non–unique. It is common to define an isotropic sourcerepresenting volume change leaving the remainder as deviatoric parts,Λi = Λi − ΛISO, (5.4)whereΛISO = (ΛT + ΛN + ΛP )/3 = tr(M)/3. (5.5)825.1. Principal decompositionOne decomposition into elementary sources isM = ΛISOI+(ΛN −ΛP )(T̂T̂T −P̂P̂T )+13(ΛT +ΛP −2ΛN )(2T̂T̂T −P̂P̂T −N̂N̂T ),(5.6)where the first term is identified with an isotropic, pressure (explosive or implosive)source (ISO), the second and third terms have no volume change with the secondterm a double–couple (DC) and the third term a compensated linear vector dipole(CLVD). Decompositions in the literature make different choices for normalizationand orientation of the CLVD, but contain the same idea of a decomposition intothree basic sources of ISO, DC and CLVD. In isotropic media, the DC is the classicfault–slip source with normal and displacement or slip vector sˆ in the fault planeobtained from the eigen vectors asnˆ =1√(2(T̂ + P̂), sˆ =1√(2(T̂− P̂) (5.7)or vice versa.5.1.1 Angle extractionAssuming the moment tensor was the result of a pure slip (i.e. DC) source thenormal and slip vectors can be used to calculate strike, dip and rake angles, but theabove equations for nˆ and sˆ (5.7) bely a problem. Since the eigenvectors T̂ and P̂are known only to within a sign there are actually eight possibilities for (nˆ, sˆ) pairs.For each of four possible nˆ vectors there are two possible orthogonal sˆ vectors (e.g.+(T+P ) has ±(T−P )). To whittle the number down to two solutions a sequence ofrules are used. First, the candidate vector for normal (e.g. ±T+P ) is selected basedon a positive vertical component so the normal points upward to follow the adoptedconvention (see section 2.6.4). Then, for that normal the two choices for sˆ are usedto construct a tensor as (nˆsˆT + sˆnˆT )/2. The scalar product is computed with theinput moment tensor and the choice for sˆ that yields a positive scalar product isselected as the slip vector. The process is repeated for candidate normals ±T − P .Given the two solutions strike and dip angles are calculated according toφ = atan2(−n2, n1), θ = acos(n3), (5.8)from which the strike vector aˆ = (sinφ, cosφ, 0) and dip vector bˆ = aˆ× nˆ then yieldthe rake angle: δ = atan2(−sˆ · bˆ, sˆ · aˆ).Isotropic DD caseFor the case of a displacement discontinuity (DD) source in an isotropic medium,Vavryc˘uk [123] gives expressions for the normal and displacement vectors in terms835.1. Principal decompositionof the sorted eigenvalues and eigenvectors of M. Defining coefficients from theeigenvalues in TNP ordered notationa1 =√ΛT − ΛNΛT − ΛP, a3 =√ΛP − ΛNΛP − ΛT(5.9)the normal and displacement vectors arenˆ = a1T̂ + a3P̂, d = a1T̂− a3P̂ (5.10)or vice versa. The slip vector may be obtained by projecting the displacementvector onto the fault plane: sˆ = dˆ− (dˆ · nˆ)nˆ which can be used as described aboveto calculate strike, dip and rake angles. In an isotropic medium the opening angle ofthe displacement vector and the fault plane may be determined from the eigenvaluesas (Vavryc˘uk [123])sinα =(ΛT + ΛP − 2ΛN )(ΛT − ΛP ). (5.11)5.1.2 Solution ambiguityIt is obvious from the way in which a moment tensor is constructed that it isimpossible to distinguish the fault normal nˆ from the displacement vector, dˆ. Recallthe term (nˆdˆT+dˆnˆT )/2. To choose a single solution, the solution with normal closestto a prior direction is selected.5.1.3 Hudson–Pearce–Rogers parametersHudson et al. [48] introduced a commonly used parameter set that makes use of anL1 normalization: |k| + |τ | + δ = 1. These parameters may be written (Chapmanand Leaney [22] eq. 42b)k =1mΛISOτ =1m2Λmin (5.12)δ =1m(Λmax − 2|Λmin|)with Λmax and Λmin being the maximum and minimum deviatoric eigenvalues (5.4)and m = |ΛISO| + Λmax. With the L1 normalization a k (ISO) versus τ (CLVD)crossplot maps to a diamond with ±ISO at the extremes of the ordinate and∓CLV D at the ends of the abscissa so pure DC plots at the origin. Figure 5.1shows a Hudson–type diamond crossplot. The special sources with labels TC, LVD,CD and IP will be discussed later. In Hudson et al. [48] an additional parameter setwas proposed that skews the diamond. This complication is based on an assumed845.1. Principal decompositionFigure 5.1: The basic structure of a Hudson (k vs. τ) or ISO vs. CLVD diamondcrossplot. The positions of various special sources are indicated. TC=tensile crack,LVD=linear vector dipole, CD=cylindrical dilation, IP=isotropic plane.uniform distribution of eigenvalues, three independent variables, which seems anuntenable assumption. While in commercial use (Baig and Urbancic [9]), the jus-tification for the skewed diamond is questionable so those parameters will not begiven or used here.5.1.4 Vavryc˘uk parametersVavryc˘uk (2001,2011) [120], [123] proposed similar parameters to Hudson et al. [48]with a different normalization:MISO = (ΛT + ΛN + ΛP )/3MCLV D = 2(ΛT + ΛP − 2ΛN )/3 (5.13)MDC = (ΛT − ΛP + |ΛT + ΛP − 2ΛN |)/2with normalization M = |MISO|+ |MCLV D|+MDC for parametersCISO = MISO/MCCLV D = MCLV D/M (5.14)CDC = MDC/M.A crossplot of Vavryc˘uk’s CISO versus CCLV D produces a similar diamond areaas the Hudson parameters except that a +TC source has positive ISO and CLV Dparameters as opposed to a negative CLV D parameter as for the Hudson parame-ters.855.1. Principal decompositionFigure 5.2: Normalized moment tensors and their eigen decomposition (eigenval-ues are shown above the TNP eigenvectors) for a slip source with SDR angles(5◦, 30◦, 30◦). Results for the same isotropic and VTI media as for figure 2.4 areshown.5.1.5 Limitations of conventional decompositions in anisotropyLeaney and Chapman [59] discussed the computation of radiation patterns from mo-ment tensors in anisotropic media but also demonstrated how anisotropy can distortthe underlying geometry of the source. By way of illustration they constructed amoment tensor in a VTI medium and showed how basic sources derived from theeigenvalues show non–DC components for a pure slip source. They also showed howthe eigenvectors are rotated and hence how the strike,dip,rake angles would be inerror. Figure 5.2 shows the eigen decomposition for slip–source moment tensors inisotropic and anisotropic media. Notice the zero intermediate eigenvalue for theisotropic case. Using a percentage measure of DC (e.g. Jost and Herrmann [49]) thetensor in the anisotropic medium is seen to possess only 73% DC. The difference inTNP vectors from which SDR angles are obtained is also notable.Similarly to Vavryc˘uk [121], in Leaney and Chapman [59] it was shown howthe correct source type and geometry is recoverable once anisotropy is taken intoaccount. Essentially the eigen analysis is carried out on the potency tensor or whatVavryc˘uk calls the source tensor, D = s : M, instead of on M. Figure 5.3 showsan example of a source with slip plus opening. Note now the zero intermediateeigenvalue of D, a characteristic of a DD source.However, as mentioned in chapter two, a displacement discontinuity (DD) orpoint dislocation source is an incomplete moment tensor as it is missing a pureisotropic part, requiring only five parameters for its construction (four angles anda size, A[d]). Thus, the decomposition of s : M as used by Vavryc˘uk [121] and asdiscussed in Leaney and Chapman [59] is an incomplete solution as it decomposes865.2. Anisotropic moment tensor decompositionFigure 5.3: Left: radiation pattern for a composite source with equal parts in–planeand normal displacement. The fault plane strikes N5E with 30 degrees dip to thewest. In–plane slip has rake= 30◦. Right: The eigen analysis of the moment tensor,M, and of the source tensor, D.an incomplete representation of a moment tensor.A complete moment tensor was constructed earlier (eq. 2.36); its decomposi-tion is now discussed, showing how the essential source geometry can be recovered,free from anisotropic distortion. This will complete the process forward from ge-ometric source parameters, through the anisotropic medium at the source to themoment tensor with source function, through the layered VTI medium to receivers(chapters 2 and 3), back to the moment tensor and source function at the sourcethrough waveform inversion (chapter 4), and now finally back to geometric sourceparameters.5.2 Anisotropic moment tensor decompositionTo begin, the equation from chapter 2 for moment tensor construction is restated:M = [V ]κembI +12A[d]c : [nˆdˆT + dˆnˆT ]. (5.15)The task is to recover from M the parameters of interest free from distortion due toanisotropy. These parameters are the fracture normal and displacement vectors nˆand dˆ, the size or potency of the DD part of the source, A[d], and [V ], the volumechange due to a pressure change [P ] in a cavity of volume V . It is recognized thatthe vectors nˆ and dˆ require four angles for their representation, for example theSDRO angles, making six parameters to recover from M in all.875.2. Anisotropic moment tensor decomposition5.2.1 Biaxial decompositionThe decomposition of the DD dyadic part of an isotropic moment tensor (the secondterm including A[d] in (2.36) above), seems to first appear in the appendix of Du-fumier and Rivera [35] without citation. Vavryc˘uk [121] develops the same essentialtheory without citing any prior work. While the impact of anisotropy on the DDpart of M is removed in Vavryc˘uk [121] by contraction of M with the compliancetensor, s, there is no isotropic term included. The same basic equations as thosegiven in Dufumier and Rivera [35] and Vavryc˘uk [121] can be found in a book byFederov (1965) translated from Russian (Federov, 1968, [41]) wherein can be foundan alternative decomposition for second order tensors. In Chapman and Leaney[22] it is termed the “biaxial” decomposition and its mathematical properties aredescribed in section 3.1 of that paper.The biaxial decomposition for a second order tensor X is written (Chapman andLeaney [22]) eq. 60) asX = x2I +x1 − x32(Φ̂−Φ̂T+ + Φ̂+Φ̂T−), (5.16)where x1 ≥ x2 ≥ x3 are the ordered eigenvalues of X and the unit vectors Φ̂±, thebiaxes, are given byΦ̂± = cosφ xˆ1 ± sinφ xˆ3, (5.17)with eigenvectors xˆ corresponding to the ordered eigenvalues. The angle φ is givenin terms of the eigenvalues ascosφ =∣∣∣∣√x1 − x3x1 − x2∣∣∣∣ , sinφ =∣∣∣∣√x2 − x3x1 − x3∣∣∣∣ (5.18)from which the opening angle is χ = pi/2− 2φ.5.2.2 Biaxial potency tensor decompositionNote the similarity in form of (5.16) and (5.15) or the potency tensor (2.44), restatedhere:D = s : M = s : M˜ISO + s : MDD = [V ]κembK +12A[d](dˆnˆT + nˆdˆT). (5.19)If the biaxial decomposition is applied to the potency tensor the result is a dis-placement discontinuity (DD) source plus an isotropic potency tensor. However, inan anisotropic medium an isotropic dilatation would correspond to a complicatedsource mechanism. Rather, an isotropic stress source is preferred, a change of pres-sure such as an explosion or implosion. This suggests an algorithm wherein theamount of M˜ISO is determined such that the intermediate eigenvalue of the remain-ing potency tensor is zero, meaning that the remaining source is a DD source, free885.2. Anisotropic moment tensor decompositionfrom distortions due to anisotropy. In a general anisotropic medium, M˜ISO requiresthe solution of a cubic (Chapman and Leaney [22] eq. A14), with the three rootscorresponding to zeroing the three eigenvalues, so one root zeros the required inter-mediate eigenvalue. In practice a simple line search is done using the secant methodwhere the intermediate eigenvalue is the objective function in a search for M˜ISO.This maybe be written asΛ2[s :(M− M˜ISOI)]= 0, (5.20)where now Λ1 ≥ Λ2 ≥ Λ3 are the ordered eigenvalues.UniquenessTape and Tape (2013, pers. comm.) pointed out that the above approach doesnot always provide a unique decomposition of M. Chapman and Leaney [23] havedetermined the necessary condition for uniqueness in terms of the eigenvalues ofthe hydrostatic pressure modulus, K. Considering a VTI medium the necessarycondition for the decomposition to be unique is that Thomsen’s parameter δ mustmeet the following condition:δ <2C44C33. (5.21)Looking at many published measurements of VTI anisotropy Chapman and Leaney[23] found only one that violated this condition, and as that measurement hada large negative anellipticity it can be discounted as an unreliable measurement.(Schoenberg (1994) [95] demonstrates that a stack of isotropic layers must have along–wavelength equivalent medium with positive anellipticity.) Although criterion(5.21) for VTI media is necessary for the above decomposition to be unique, it isnot a sufficient condition and a unique decomposition may still exist even if (5.21)is violated; it depends on the potency tensor D. No problems have been found inpractice to recover the correct geometric parameters for any realistic medium.5.2.3 Source geometry parametersHaving determined the value for M˜ISO the parameters of interest are recovered asfollows:[V ] = M˜ISO/κembA[d] = Λ1 − Λ2 (5.22)dˆ = ±Φ̂nˆ = ∓Φ̂,having applied the biaxial decomposition to DDD = s : (M− M˜ISOI).895.2. Anisotropic moment tensor decompositionIt is noted that the angle χ seen previously (2.32) is expressed in terms of thebiaxes angle φ as χ = pi/2 − 2φ (Chapman and Leaney [22], eq. 89), which isequal to Vavryc˘uk’s (2011) angle α given previously in (5.11) with the intermediateeigenvalue, ΛN , now zero. Equivalent expressions for the normal and displacementvectors can also be found in terms of the eigenvectors of DDD (Vavryc˘uk [121], eq.47 and 48), rather than the biaxes (5.17) which are given in terms of the angle φ.5.2.4 Biaxial potency parameterizationThe potency tensor for the DD isDDD =12A[d](dˆnˆT + nˆdˆT)= A[d] sinχnˆnˆT +12A[d] cosχ(sˆnˆT + nˆsˆT ) (5.23)= DN + DSwhere the vector sˆ is the unit displacement vector projected onto the fracture plane— the slip vector. The first term of (5.23), with displacements in the normal di-rection nˆ, is recognized as the tensile crack source while the second dyad is theslip source. The moment tensor has thus been decomposed into basic sources withintuitive meanings: a tensile crack, TC or “Opening” source, DN and a slip source,DS . As will be seen shortly, the TC source DN does not require two parameters forits representation as with the Hudson k−τ or Vavryc˘uk CISO−CCLV D parameters.The total potency magnitude isDDD = ‖D‖ = A[d](1 + sin2 χ)1/2(5.24)with individual potency magnitudesDN = 21/2A[d] sinχ (5.25)DS = A[d] cosχ (5.26)which satisfy the norm relationship D2N + D2S = D2DD. The parameters [V ], DNand DS have been called the EOS parameters (Leaney et al. [62]), Leaney et al.[64]) standing for Expansion, Opening and Slip. Note again that these parametersare in units of volume. “EOS” may be used to refer both to the parameters andthe SDRO angles derived from the DD potency tensor. Normalized versions of theEOS parameters are of interest for cross–plotting purposes and for comparison withother parameterizations.5.2.5 Normalized EOSBy analogy with the normalization of the Hudson parameters, and with the goal ofcomparison via a similar diamond crossplot as figure 5.1, in Chapman and Leaney905.2. Anisotropic moment tensor decomposition[22] two new parameters were proposed that make use of the biaxial potency de-composition. They are (Chapman and Leaney [22] eq. 98 and 99)ν =[V ]|[V ]|+A[d](5.27)n = −DNDDDA[d]|[V ]|+A[d], (5.28)where the negative sign on n has been included for cross–plotting consistency withthe Hudson τ parameter. A problem with this normalization is that a source withno [V ] and equal parts opening and slip does not plot halfway from the diamondcenter to the n vertex. This suggests a set of normalized EOS parametersÊ =[V ]|[V ]|+A[d](5.29)Ô = sign(χ) sin2 χA[d]|[V ]|+A[d](5.30)Ŝ = cos2 χA[d]|[V ]|+A[d](5.31)which also satisfy the L1 normalization |Ê|+ |Ô|+ Ŝ = 1.5.2.6 Impact of anisotropyThe impact of anisotropy can be viewed in crossplot space by synthesizing a momenttensor in a VTI medium and comparing decomposition parameters on the diamondcrossplot. Staying with the anisotropy parameters of the homogeneous VTI mediumused thus far for synthetics, and the SDR parameters (60◦, 40◦, 20◦), moment ten-sors are constructed for increasing opening angle χ increasing from −90◦ to 90◦.In an isotropic medium the Hudson k vs τ and Vavryc˘uk ISO vs CLV D param-eters trace out a straight line (recall the opposite sign in CLV D between Hudsonand Vavryc˘uk), but in a VTI medium there is significant distortion. For both me-dia the Ê, Ô, Ŝ parameters trace a horizontal line from Ô = −1 to Ô = +1, sincethe distorting effect of the anisotropy has been removed. This is shown in figure5.4. Figure 5.5 shows the recovered SDRO (Strike,Dip,Rake,Opening) angles us-ing the biaxes potency tensor decomposition and a decomposition of M assumingisotropy (using equations (5.9)–(5.11) and procedures described in 5.1.1). Note thecorrect recovery of SDR angles (60◦, 40◦, 20◦) and the opening angle χ using EOS,but significant errors due to anisotropy when M is decomposed assuming isotropy,including polarity flips for the rake direction in this case when χ approaches 90◦(pure opening).915.2. Anisotropic moment tensor decompositionFigure 5.4: Normalized source parameter crossplots for DD sources with openingangle χ ranging from −90◦ to +90◦. Left: isotropic medium; right: VTI medium.SDR angles are (60◦, 40◦, 20◦). Green: Hudson k vs. τ ; blue: Vavryc˘uk ISO vs.CLV D; red: Ê vs. Ô.Figure 5.5: The opening angle χ was increased from−90◦ to +90◦ in the constructionof M in the anisotropic medium. On the left are shown SDRO angles recovered usingthe biaxes decomposition of DDD while on the right they are shown obtained fromM conventionally assuming isotropy. S=red, D=green, R=blue, O=cyan.5.2.7 Special sourcesIt is useful to consider several special sources and where they plot on the diamond.925.2. Anisotropic moment tensor decompositionTensile crack, TCA tensile crack (TC) source, the opening or closing along an infinite plane, producesan isotropic term with an additional dipole term in the direction of opening, re-quiring both k and τ or ISO and CLV D for its representation. Since δ = 0 (thereis no DC), it maps to the perimeter of the diamond. In a Poisson medium, wherePoisson’s ratio=1/4 and Lame´ moduli are such that λ = µ so that the ratio ofcompressional to shear wave speeds is Vp/Vs =√3, Hudson parameters are k = 5/9and τ = −4/9 for a TC source as plotted in figure 5.1. The special case of a verti-cally oriented TC source was also considered by Aki and Richards [1] (page 51 andfigure 3.8), who show that a TC source is the superposition of three vector dipoleswith magnitudes in the ratio λ : λ : λ + 2µ. In the potency domain a TC sourcehas a particularly simple form, with corresponding ratios 0 : 0 : 1 (Chapman andLeaney [22], appendix B5, eq. B10b). Some sources are simpler in potency, some inmoment.Linear vector dipole, LVDThe tensile crack is a source type thought to be representative of crack opening asmight be expected due to fluid injection, but as the TC model is for an infinite planeit is clearly an end member of crack opening, the other end member being a simpledipole or LVD (linear vector dipole), where the crack opens at a single point. AnLVD has a simple moment tensor with ratios 0 : 0 : 1 and in a Poisson mediumits potency tensor has ratios −σ : −σ : 1 where σ is Poisson’s ratio (Chapman andLeaney [22], appendix B4). An LVD source has a smaller ISO part and Hudsonparameters k = 1/3 and τ = −2/3. Note that while there is no transverse motion fora TC source, the LVD contracts laterally as can be seen from the ratios −σ : −σ : 1(Chapman and Leaney [22], appendix B4, eq. B8b).Isotropic plane, IPThe isotropic plane source is thought to be a model for a perforation shot, whichfires a sequence of bullets to perforate the casing, arranged at equal angular spacingaround the 360◦ circumference of the cased borehole. Angular spacings are com-monly 120◦ distributed along a spiral trajectory over a distance of a few feet. Bulletsare fired in rapid succession following the burn–rate of prima–cord; essentially in-stantaneously. There needn’t be an even number of or symmetric arrangement ofbullets but an isotropic plane or “2D explosion” is one simple model for a perfora-tion shot. An isotropic plane (IP) has a simple moment tensor with ratios 1 : 1 : 0and potency tensor ratios 1− σ : 1− σ : −2σ (Chapman and Leaney [22], appendixB6). Hudson parameters for an IP sources are k = 1/2 and τ = 1/2.Other special sources are considered in appendix B of Chapman and Leaney [22]935.3. Decompositions of noisy Mfor the case of an isotropic medium with the properties of a Poisson solid. The CD(cylindrical dilatation) source is to the IP source as the TC source is to the LVD.In other words, CD and TC are simpler in potency; IP and LVD are simpler inmoment. The CD and IP special sources will figure in the discussion on real dataresults. For the case of general anisotropy special sources should be defined in thesimplest domain (moment or potency), rotated using second order tensor rotation,and converted to moment for simulation as needed using M = c : D.5.3 Decompositions of noisy MThe results of the noisy simulations and multi–variate normal posterior computa-tions of the previous chapter may now be viewed in the space of parameters ofgreater interest and intuitive meaning. Figure 5.6 shows the histograms for theelements of M from the previous chapter along with normalized source parameterhistograms. The Ô and Ŝ distributions are slightly skewed although centered onthe correct values (0.5,0.5). The broader histograms for Ê, Ô, Ŝ parameters andparticularly for the Ô and Ŝ parameters is attributed to the nonlinearity of theEOS decomposition and the fact that the realized moment tensors are essentiallycontracted with the compliance tensor (s : M) resulting in distortion. The Hudsonand Vavryc˘uk parameters are more simply related to M. It is possible to invertdirectly for D as explored by Leaney and Chapman [60]; it is still a linear problem,but the decomposition requires an isotropic pressure (MISO) so an inversion for Mis preferred.It is also interesting to look at noisy moment tensors plotted as decomposedparameters in diamond crossplot space, and also the characteristic source angles.This is shown in figure 5.7, where the results are also shown for the MVN samplingof a diagonal covariance matrix. The histograms of the SDRO angles show that thecorrect angles are recovered as expected with EOS; angles recovered from M underan isotropic assumption are significantly in error. Again, the somewhat sharperdistributions for the angles from M̂ are conjectured to be due to the mapping fromM̂ to D through the compliance tensor, D = s : M.Source parameters can be a useful aid in understanding mechanisms of collectionsof events, and they can be color–coded by a third attribute, for example time (e.g.Baig and Urbancic [9]), but to intuitively visualize the geometry of a compositesource can be challenging. In the next section a new way of source visualization isdescribed that makes use of the EOS decomposition.945.3. Decompositions of noisy MFigure 5.6: Histograms for multivariate sampled moment tensor elements (up-per left) and histograms of normalized source parameters. Hudson (upper right),Vavryc˘uk (lower right) and Ê, Ô, Ŝ (lower left).955.4. Source visualizationFigure 5.7: Decomposition parameters on a diamond crossplot for the multi–variatenormal posterior samples of M using full covariance and diagonal–only (on top).Hudson (green and magenta); Vavryc˘uk (blue and yellow), Ê, Ô, Ŝ (red and cyan).The true model plots at (Ô, Ê) = (0.5, 0). On the right are SDRO angles fromEOS (solid) and from M assuming isotropy. Vertical lines show true angles.5.4 Source visualizationIt is common to display the nature of a source via its radiation pattern in an equal–area display (e.g. figure 2.5) and a lower hemisphere display of the P radiation hascome to be known as the “beach ball” representation. For simple DC sources usershave become adept at inferring the source mechanism because changes in sign inthe radiation pattern immediately lead to the fault plane, but for more complicatedsources it is virtually impossible to interpret such displays. Some users might alsolike the 3D P radiation pattern, for example as shown in figure 17 of Chapman andLeaney [22], but considering that the new EOS parameters are in units of volumethe idea of depicting a decomposed source as a 3D graphical object suggests itself.In Chapman and Leaney [22] this idea is developed so that [V ] is equated with asphere, the fracture plane is represented by a disk or plate with a black colorednormal vector nˆ, and another vector is used for dˆ, emanating from the center of theobject. Color is used to depict sign, red for positive (expansion or opening); bluefor negative (contraction or closing).5.4.1 Hockey pucksAnother visualization choice (Leaney et al. [64]), which has come to be known asthe “hockey puck” display, ascribes the DD part of the potency — A[d] — to thevolume of a cylinder. The cylinder is then cut in half and displaced according to the965.4. Source visualizationFigure 5.8: New 3D graphical representation of a composite source showing orientedEOS effective volumes or potencies. Left: map view, right: side view from south.This synthetic example has significant −E = −Expansion (blue wire frame) andan opening angle χ = 45◦ to split A[d] into +O = Opening (thickness of red pucks)and S = Slip (pucks displacement). Strike, dip and rake parameters are 60◦, 70◦and 10◦, respectively. Of the two solutions, the one closest to a prior orientation (inthis case the true solution) has been chosen for display.amount of slip. Specifically, the following associations are made to the componentsof the 3D object.The total potency is ascribed to the volume of a sphere of radiusr0 =3√3(|[V ]|+A[d])/(4pi).A sphere is used to depict [V ] with radius rv = r0[V ]/(|[V ]| + A[d]) and a cylinderwith radius set to rDD = r0A[d]/(|[V ]|+A[d]) is used to depict the opening and slipcomponents. The height of the cylinder is set to h = 2rDD| sinχ| and the cylinder issliced in half so as to portray slip as the displacement between two “hockey pucks”with slip equal to s = rDD cosχ. Practically, the limiting case of pure S and no Orequires two disks with some minimum thickness so as to be seen from all angles.Colour is used for the sign of E and O, red for expansion and opening, blue forcontraction (negative pressure change) and closing. Additionally, the sphere for [V ]is drawn with a wire frame with poles along the T axis and the pucks are made semi–transparent so the wireframe sphere can be seen inside. Figure 5.8 shows two viewsof the hockey puck glyph for a specially chosen source constructed with −[V ] = A[d]and an opening angle χ = 45◦ to produce equal parts opening and slip.All 3D source display options discussed here have been engineered into Schlum-berger’s PetrelTM interpretation system which will be used to display real datainversion results.97Chapter 6Field Data ApplicationIn industry, getting permission to publish with field or “real” data examples can beproblematic as the oil and gas company owns the data and special legal permission isrequired. So while those of us in the service side of the business see many fascinatingdata sets only a small fraction ever get seen by the public. Ultimately, it usuallytakes a benevolent customer to get permission. In this case acknowledgement goesto John Duhault and Lightstream Resources (Calgary) for granting permission touse and publish with this real data set. While it is a two well monitoring data set(three or four would have been preferred), the waveforms are of excellent quality.The location is the West Pembina field of west-central Alberta, Cardium forma-tion. This oil field currently consists of 20 pools that have produced over 1.1 billionbbls from an original oil in place of 7.46 billion bbls (Alberta Geological Societywebsite). Duhault [36] provides a review of the geology and utilization of microseis-mic monitoring of the hydraulic stimulation of the Cardium formation. Figure ??shows the location of the Cardium with respect to other unconventional basins inWestern Canada.The Cardium is a marine sandstone of late Cretaceous (Turorian) age that over-lies the mudstones of the Blackstone formation and includes repeated successionsof silty mudstones through siltstones to fine– grained sandstones. These units areunconformably overlain by chert-pebble conglomerates which in turn are overlainby the marine mudstones of the Colorado (Wapaibi). The thicker, highly permeableand porous conventional reservoirs have been exploited since the 50’s, having beendiscovered by Mobil Oil in 1953. The current focus is on the margins where thegross reservoir is up to twelve meters thick with net sandstone thickness varyingbetween three and seven meters. In the West Pembina field the porosity averagesless than 10%, hence hydraulic fracturing is used to stimulate the formation to easethe flow of oil and increase production.6.1 StimulationExploitation has evolved to where it is now done by drilling lateral wells. Thefundamental idea is to drill roughly in the direction of minimum horizontal stress,as this is the direction in which fractures will tend to open. Well stimulations for thisdata set have been done using slotted-sleeve hardware, in contrast to a “plug–and–986.1. Stimulationperf” operation. (See figure 6.1.) In plug–and–perf, the well is cased and cementedand plugs are installed to isolate a zone to be perforated, charges are fired to createperforations in the casing, and then fluid pressure is increased to initiate fracturing.Another plug is set and the process is repeated, moving from the toe to the heel ofthe well.Maxwell et al. [75] describe the process and seismic signals generated when aslotted sleeve is used. In a slotted–sleeve or “ball–seat” operation, packers are usedto isolate ports and a sleeve is designed into the tubing such that the ports areopened when the sleeve slides. The sleeve slides to open the ports for a stage bydropping a ceramic ball into the well to eventually plug a hole on the end of thesleeve. This causes pressure to build up, a shear pin fails, the sleeve slides and theports are opened to the formation. Balls of increasing diameter are dropped to allowports to be sequentially isolated and stimulated. Figure 6.1 shows schematics forThe slurry is often “energized” with nitrogen. This is said to help in initiating thefrac into the formation through the ports; it also provides for an energetic calibrationshot with special characteristics as will be seen. Top-hole pressure is monitored foran abrupt pressure increase which signifies the seating of the ball. The slurry rateis then ramped up, in this case to a maximum of about 5.5 m3/min. or about 45bbl/min. An abrupt drop in pressure indicates that the sleeve has slid. After someminutes of pumping slurry, proppant is added in doses while pressure, slurry andproppant are monitored at the surface.Figure 6.2 shows pumping data with microseismic event counts for a stage ofthe West Pembina treatment. Notice the early peak in the red curve indicating ballseat and abrupt decline indicating sleeve sliding and port opening. The timing ofthis event correlates with an energetic, sometimes complicated source as recorded inthe monitoring wells. Eighteen such stages of treatment were carried out along thelateral; ten stages (4–13) were located between the vertical monitor wells. Noticein figure 6.3 how the ports were designed with knowledge of the expected fractureazimuth trend so as not to frac into the monitor wells. Data analysis focus will beon the data best recorded in both monitor wells coming from central stages 7 and8.996.1. StimulationFigure 6.1: Graphical depictions of plug–and–perf (top), sliding sleeve completion(middle) and a schematic of the ball–seat operation to open frac ports for slidingsleeves (bottom).1006.1. StimulationFigure 6.2: Pumping data and microseismic event counts for stage 8. The slurry(fluid) rate is in blue with scale on the right as Flowrate in m3/min.; top rig pressureis in red with scale on the left in MPa; proppant concentration is in green with scaleon the right in kg/m3; event rate is in double green bars with scale on the left. Barsindicate event count in 2 minute windows. The early build–up of pressure (red)followed by a sudden decline indicates the sliding of the sleeve and opening of theport to the formation.1016.2. Monitoring geometry and hardware6.2 Monitoring geometry and hardwareThe lateral treatment well was drilled from south to north to take advantage oftwo vintage vertical wells from 1960 which were used for the monitoring arrays ofthree component (3C) receivers. Figure 6.3 shows the map view of the survey layoutincluding all qualified microseismic event locations with the 18 stages color-coded.Figure 6.4 shows a view from the southwest with formation tops and receivers. Thereceiver spacing for these arrays is 30m, the top three receivers are not shown. Thesefigures come from the commercial processing report.In each monitor well Schlumberger’s VSI-12 borehole seismic tool was deployed.This tool is unique in the industry in that it has sensor packages that are releasedon springs to couple to the formation or casing, effectively de-coupling the geo-phones from the body of the tool. Opposite to the sensor package is a clampingarm that pushes the shuttle against the side of the borehole. The sensor package isequipped with three geophone-accelerometers (Schlumberger’s GAC-D sensors) ar-ranged orthogonally with the x-axis in the clamping direction and the z-axis alongthe well axis. The y-axis is normal to these two, transverse to the clamping di-rection. Because the deployed sensor package is so light the springs provide aneffective clamping force of 15kg. This isolated sensor package design is radicallydifferent from the standard design of a sensor package built as part of the body ofthe shuttle whose clamping force comes from the mechanical arm used to anchor thetool. It results in greater immunity to borehole modes such as tube waves, providessuperior vector fidelity and a generally lower noise floor.The GAC-D geophones are omni-tilt, damped accelerometers so no space andelectronics in the tool are needed for a gimballing device — geophone orientationis done in processing using a relative bearing angle measurement from the tool anddirect P wave polarization vectors. The GAC-D has a characteristic frequency of20Hz with a damping coefficient of .7. It has a flat response in acceleration downto 5Hz at 125C and up to over 400Hz.The signals from each 3C shuttle arrive at a cartridge at the top of the arraywhere digitization at .5ms sampling rate occurs prior to telemetric transmission tothe surface recording system on standard heptacable wireline. Figure 6.5 shows aschematic of a VSI shuttle, a photo of a shuttle, a drawing showing the dimensionsof a 4 shuttle VSI array with telemetry cartridge and a photo of a VSI-40 array toolin the shop with wireline interconnects, getting ready for a deep offshore Gulf ofMexico deployment. In figure 6.6 are photos of a rig-up of a VSI array for hydraulicfracture monitoring operations and a screen grab of the acquisition software in actionduring a VSI-12 HFM job.Acquisition proceeds in independent recording units at the two monitor well sitesand the data are merged and synchronized based on GPS time in the processingpackage. After processes of geometry QC, geophone orientation, event detection1026.2. Monitoring geometry and hardwareFigure 6.3: Map view of the wells and event locations color-coded by stage. Thegrids are 100m on a side.1036.2. Monitoring geometry and hardwareFigure 6.4: Side view from the southwest the wells and event locations color-codedby stage. The grids are 100m on a side.and location using a calibrated VTI model (next section), a subset of dual arrayevent files, one file for each oriented ENU component, is saved. Figure 6.7 showsa large event from a central stage 8 with no filtering and the amplitude spectra intrace-normalized dB along with the median spectrum. Broad-band signal is presentfrom 5-600Hz; the acquisition anti-alias roll-off can be seen above 800Hz. The dataquality is generally very good. All receivers worked continuously throughout thejob, providing high fidelity waveforms.Before moving on to model building and VTI calibration, figure 6.8 comparesray-rotated PHV waveforms for the event shown in figure 6.7 and a sliding sleeveevent. The differences in recorded radiation pattern can be seen for the event bothbetween the two arrays and within the arrays. An unmistakable delay is observedbetween the qSv and Sh arrivals. The sliding sleeve event shows a markedly differentrecording. It is characterized by a more complicated coda (an apparent doublet)and relatively weak shear. Importantly, the first–motion polarity of the P arrivalis opposite on the two arrays for the sleeve event compared to the microseismicevent. It is approximately 1/10th the amplitude of the microseismic event. As willbe seen, the sleeve events are more directional than the micro-earthquakes and arewell modelled by a horizontal force source pointing towards the heel of the well.Although relatively weak, the sleeve events provided enough pickable shear arrivalsfor a reliable VTI model calibration.1046.2. Monitoring geometry and hardwareFigure 6.5: Upper left: schematics of Schlumberger’s VSI downhole 3C shuttleshowing the detached sensor package design. Lower left: a photo of a VSI shuttlewith anchor arm extended and sensor package. Upper right: VSI-4 dimensionsincluding telemetry cartridge and individual shuttle; Lower right: a VSI-40 array inthe shop with wireline interconnects, getting ready for deep offshore Gulf of Mexicodeployment. Thanks to Bill Borland for the graphics.1056.2. Monitoring geometry and hardwareFigure 6.6: Photos of rigging up a VSI array for hydraulic fracture monitoringoperations (top) and a screen grab of the acquisition software for a VSI-12 (see toollay-out on right). Thanks again to Bill Borland for the images.1066.2. Monitoring geometry and hardwareFigure 6.7: Oriented ENU waveforms after merging the data from the two monitor-ing arrays for a large detected and located event (left) and the amplitude spectra(right) shown in dB (trace-normalized) with the median spectrum in black.1076.2. Monitoring geometry and hardwareFigure 6.8: Ray-rotated PHV waveforms for the large microseism event of figure6.7 and for a sliding sleeve event. This microseism has roughly ten times the ampli-tude of the sleeve event. The thin coloured lines represent the PHV arrival times(PHV =Blue,Red,Green) using the calibrated VTI model and event locations.1086.3. Model building6.3 Model building6.3.1 Sonic logsSonic logs are an essential data type for downhole microseismic monitoring as theyprovide the detailed vertical P and S velocities from which an initial isotropic modelis built. Sonic logs came from MW 14-19-48-10W5 (P and S ). These logs didnot go to the depth of the reservoir due to insufficient “rat-hole” (extra depth toaccommodate the bottom of the sonic tool) in the vintage wells, but an offset well(Bellatrix 16-19-48-11W5) had P and S sonic through the reservoir which was mergedwith the sonic from the MW. Unfortunately, a density log was not available so onewas constructed using an empirical relation from the compressional velocity. It wasvisually calibrated against the density log shown in figure 7 of Duhault [36]. Densityappears in the impedance denominator term (see, e.g. equation 2.25) and impactsthe scalar moment estimate. Figure 6.9 shows the elastic logs used for initial modelbuilding. They have had a 9,5 trimmed mean filter applied as a mild de-spike andmoving average conditioning. The 5 sample window corresponds to two feet, aboutthe vertical resolution of the sonic tool. The Vp/Vs ratio log is often a good proxyfor the volume of clay and this is also the case here. It is a minimum in the sandstonereservoir at 840m depth.6.3.2 Log blockingSince a layered ray tracer is going to be used for Green function computations the.5ft sampling of the logs needs to be upscaled. This is done using slowness smoothingfollowed by minimum thickness automatic log blocking. The log blocking algorithmis based on a cascaded or “compound” median filter (Leaney and Ulrych [68]). Theminimum layer thickness is determined based on a wavelength criterion. Liner andFei (2007, [71]) investigated equivalent media using Backus averaging theory [6]and 2D FD simulations, and came up with “The Backus number”. Based on theirrecommended scattering and transmission averaging length criteria, models can begenerated. The recommended upper bound for an averaging window is given interms of the minimum shear velocity and dominant frequency asL ≤a min(VS0)fdom, (6.1)where min(VS0) means the minimum vertical shear velocity, fdom is the dominantfrequency of the signal and a = 1/3 for scattering and a = 2 for transmission. Forrobustness the minimum VS0 is replaced with the 10th percentile, and the numberL is taken as the half–width of a raised cosine averaging window. Smoothing isapplied to the sonic slownesses (to preserve vertical traveltime) prior to determininglayer boundaries automatically with compound median filtering (Leaney and Ulrych1096.3. Model buildingFigure 6.9: Elastic logs used for initial model building. The empirically deriveddensity is in units of gm/cc, Vp and Vs are in m/s (note the factor 102). Thereservoir is located near the minimum of the Vp/Vs ratio log at about 840m.1106.3. Model building[68]). Using a dominant frequency fdom = 100Hz and the determined value V s(10) =1845m/s the transmission Backus length is LT = 37m (the scattering or reflectionBackus length is LS = 9.3m). While four samples in this window length wouldsuffice for adequate representation of velocity gradients ten samples has been founda better number, leading to a minimum layer thickness of 3.7m, which was usedfor automatic interface detection. Interface detection may be done using any log,for example Vp/Vs, ρVp, etc. The idea of selecting a log for blocking is to chooseone that represents changes in lithology. Blocking on Vp/Vs is a useful strategy fora sand–shale sequence as Vp/Vs is a good proxy for the volume of clay and henceshaliness; blocking on ρVp may be preferred for P–wave reflection imaging.End effects are handled in the log smoothing operation not by padding a halfwindow with log end values but by a mirror–flip strategy coined “conjugate extrap-olation”. In this approach the end values are used as mirror pivots and the log isextended with previous values in reverse order such that the gradient is preserved onall length scales of the filter. Figure 6.10 shows an example of this process appliedto the Vp log.As a brief personal aside, the conception of this idea dates to a memorable night inParis in 1992. Tad and I refer to this fondly as “the Johnny Walker method”. It isused commercially to this day. The original idea was to use it to extend time seriesprior to Fourier-based spectral analysis. This aspect has yet to be investigated.Figure 6.10 also shows the Vp log and layered models for LS , LT and LT+:LT+ = LT + (LT − LS)/2. The reason for extra smoothing is that errors areexpected due to logging, lateral heterogeneity and geometry so sharp interfacesrepresent an unrealistic precision, hence it is safer to be smooth and to remainuncommitted to what may be poorly known. It is noted that Backus averagingto upscale isotropic sonic layers to a VTI model is not done as it tends to vastlyunderestimate VTI anisotropy. This lower bound on VTI anisotropy is of limited useas picked calibration shot times will be used to calibrate the model and much largeranisotropy is generally required than that predicted by Backus averaging. For thisstudy the model with 150Hz transmission smoothing is carried forward for modelcalibration.1116.4. VTI model calibrationFigure 6.10: Left: Vp log and extended log using conjugate extrapolation to handleend effects prior to smoothing. This approach preserves the end-zone gradient.Right: Vp log and three layered models corresponding to different averaging criteria.The green coloured model corresponding to the transmission smoothing criterion iscarried forward for VTI model calibration.6.4 VTI model calibrationGiven the initial layered model the next step is to calibrate it for anisotropy. Theraw material used for this is calibration shot data from the treatment well. Thismay take several forms depending on the job. Commonly, perforation shots areavailable at the beginning of every stage, but as mentioned earlier these are notavailable for the slotted sleeve completion used at West Pembina. Sometimes stringshots are taken if perf shots are not available, more rarely a downhole source (eitherpiezo-electric or vibrator) may be deployed prior to stimulation. For this data setthe model calibration data were the slotted sleeve events identified in the data. Itis noted that to guarantee the availability of data for the geophone orientation stepthree surface vibroseis shots were acquired into the downhole tools. These shotsended up giving equivalent geophone orientations as the sleeve events. Had thesleeve events not been identifiable in the waveform data other strategies for modelcalibration would have been employed, for example using early events and assumingthey came from the sleeve port location.So the “ball-drop” or sleeve events were used for VTI model calibration. Thesewere distinguished from microseisms based on observed P polarity being oppositeon the two arrays. Nine such sleeve events were identified with confidence and wereused for VTI model calibration.1126.4. VTI model calibration6.4.1 VTI calibration algorithmThe method behind the VTI inversion is described by Mizuno, Leaney and Michaud(2010, [81]) and is part of Schlumberger’s commercial microseismic processing pack-age. It makes use of the same VTI ray tracer described previously and takes asinput PHV arrival times picked interactively. The inversion is a linearized inversionfor five free parameters: (∆Vp,∆Vs,∆,∆δ,∆γ) that makes use of numerical deriva-tives in a damped steepest descent optimization. The total traveltime residual forphase type i (i = P,H, V ) and receiver j is written in terms of partial derivativematrix and model update vector as∆Tij =∂T∂Vp∆Vp +∂T∂Vs∆Vs +∂T∂∆+∂T∂δ∆δ +∂T∂γ∆γ. (6.2)Since the inversion for VTI parameters in every one of many layers is impossiblyill-posed given the limited ray coverage provided by calibration shots, a soft, non-linear, rock-physics constraint is used wherein the magnitude of anisotropy is drivenby an auxiliary log to reduce the parameter space. The auxiliary log a(z) is a lognormalized to be on the interval [0, 1] and is used to drive the depth-dependence ofthe magnitude of an anisotropy parameter. Thus for anisotropy parameter ξ andscaling log a(z) the model is updated with depth-dependent change in anisotropyaccording to ∆ξ(z) = ∆ξa(z). In sand-shale systems using a(z) = Vp/Vs is a goodchoice as Vp/Vs is a good proxy for the volume of clay. In systems with carbonate,a(z) = 1/Vp is a better choice to drive anisotropy. Using such an auxiliary log con-straint keeps the parameter space small and the inverse problem well-posed whileyielding physically plausible depth-dependent anisotropy. Ideally the auxiliary logcomes from a sonic measurement of VTI anisotropy (Valero et al. [118]), but suchsonic anisotropy was not available for this study.Numerical computation of derivatives makes it easy to incorporate nonlinearconstraints. Other nonlinear constraints have also been implemented, for examplethe option to hold any parameter fixed at initial model values, within bounds or freewithin depth intervals. An option to invert for Schoenberg parameters, particularlyanellipticity, Ea instead of Thomsen’s δ has also been implemented. Using thisoption together with the auxiliary log to drive the depth-dependent magnitude ofanisotropy allows δ to be both positive and negative versus depth if needed. (This isnot possible starting from an isotropic model if δ itself is driven by the auxiliary logas it lies on the positive interval [0, 1].) Figure 6.11 shows the sleeve event for stage7 in Schlumberger’s “VMI” software. Figure 6.12 shows the results of the inversionalong with the geometry and rays in 3D for the nine sleeve event locations. Initialisotropic rms time residuals were 9ms, reduced to 3.3ms after inversion.Figure 6.13 shows the VTI model, PHV rays to receivers for an event at thebase of the reservoir 300m from the receivers and PHV times. The Sh and qSv raysare different due to the differing response to anisotropy.1136.4. VTI model calibrationFigure 6.11: Schlumberger’s microseismic VTI model calibration software interfaceshowing the data selection window (upper left), the picking window (lower right),the inversion parameter window (Upper right) and the model window (lower left).Figure 6.12: Calibrated VTI model (left) and geometry with P rays from 9 cali-bration shots (right) together with color–coded Vp. Depth–dependent anisotropyfollows a soft rock physics constraint. This layered, VTI model reproduces PHVtimes picked on 9 sleeve “shot” waveforms with an rms residual of 3.3ms.1146.4. VTI model calibrationFigure 6.13: Vp, Vp/Vs and Thomsen anisotropy parameters for the calibrated VTImodel (left), PHV rays and arrival times to 12 receivers located 300m from theevent, which is located at the base of the reservoir.6.4.2 Event locationsThe event locations are taken from Schlumberger’s commercial microseismic pro-cessing package. The original technique is described by Drew et al. [32], and also inMichaud and Leaney [79] and Drew et al. [33]. The approach makes use of a scan-ning algorithm that finds the location and origin time (x, y, z, t0) that maximizes acharacteristic function (CF) based on an onset measure (STA/LTA = Short-Time-Average / Long-Time-Average). The STA/LTA CF is computed on a 3C envelopefunction centred on P and Sh times computed given (x, y, z, t0), the VTI model andthe receiver locations. An additional weight is applied to prefer P and S arrivalsthat have orthogonal polarizations. The algorithm is efficient enough, with a multi-grid search, to be able to keep up with real-time acquisition for tens of receivers.This is the approach that produced the event locations shown in figures 6.3 and 6.4.Event sets for each stage or multiple stages are exported in a proprietary format(LDF) and are read by my software. As mentioned briefly in chapter 4, a nonlinearinversion to refine event relocations during moment tensor inversion is optional.1156.5. Residual statics6.5 Residual staticsIn spite of the best efforts and use of the sophisticated VTI model inversion de-scribed earlier, errors remain. No model is perfect. The hope is that ray trajec-tories, amplitudes and polarizations match the data after the model building andcalibration process but model inadequacies and geometry uncertainties cause errors.The true model may contain dip, heterogeneities within layers or a lower symmetryof anisotropy. Monitor wells don’t always have deviation surveys carried out andthis was the case with these vintage wells from 1960 — they are assumed perfectlyvertical. All such errors and uncertainties cause kinematic (arrival time) errors, andthey need to be accounted for to get a waveform fit that will allow for confidence inthe inverted moment tensor. The procedure developed to determine residual traveltime errors is based on an iterative, complex correlation technique.6.5.1 Iterative, complex correlation staticsRecall the ray-rotations of the synthetic data from chapter 2, shown in figure 2.12and for a real microseism and sleeve events in figure 6.8. Since statics will betreated separately for P, Sh and Sv arrivals, they are first projected according totheir receiver polarization vectors to maximize their amplitude and isolate themfrom each other — the PHV rotations. Due to the radiation patterns the polarityof the arrivals is likely to change with receiver direction relative to the source. Thisposes a problem for a correlation–based waveform alignment and this problem ishandled by using the complex trace.The algorithm proceeds by subtracting the predicted arrival time of the ray-rotated PHV arrivals. This time shift is done efficiently and accurately in thefrequency domain. At this point the arrival should be fairly well aligned, withresidual time shifts and possible polarity changes remaining. The complex trace iscomputed for each trace (e.g. Ulrych et al. [115]). The complex stack is computedas the first (complex) principal component s(t) = UT1 (t)d(t) (See the discussionin section 4.5, the same subroutine is used.) Time shifts relative to the stack aredetermined by correlating the complex stack with the individual arrival traces andinterpolating the peak of the envelope. The determined shifts are applied to eachtrace, a new complex stack is computed and the process is repeated. The stack andtime shifts generally cease changing after fewer than 8 iterations. As the relativetiming of PHV arrivals may have changed during this process, a final time shiftis determined for PHV stacks, choosing the one with maximum amplitude as thereference.1166.6. Receiver statics example6.6 Receiver statics exampleFigure 6.14 shows the microseism event from figure 6.8 after time alignment bysubtracting VTI model times from PHV rotated data and after residual staticsdetermination using complex correlation. The determined residual travel times areadded to modelled times to improve the kinematic model prior to moment tensorinversion. The improvement in waveform fit can exceed a factor of 2. This allowsa broader bandwidth to be used in inversion and reduces the uncertainty in theestimated moment tensor.Figure 6.14: Calibrated VTI model time-flattened PHV arrivals (left) and afterresidual statics corrections (right) using the complex correlation algorithm.6.7 Q modelThe default value for Q is 100, set to be equal for PHV arrivals and constant inall layers. The maximum gain allowed for Q compensation in the inversion is auser defined parameter, set to 5 in this case given the signal-to-noise ratio for anaverage event. The dominant frequency for the shear (S ) arrivals is observed to beless than that of the P. It is conjectured that a physical process may be responsiblefor the difference in S/P arrival frequency, perhaps one of those discussed in section2.4. In an attempt to compensate somewhat for the lower S frequency a lower Qwas used for S than for P. Values used after some trial and error are Qp = 90,Qsh = Qsv = 50. While some work has been done on a nonlinear waveform–fitinversion for Q, this remains an area for future work.6.8 Moment tensor inversion exampleHaving determined receiver statics the event is now inverted for its moment tensorsignals. Figure 6.15 shows input ENU data and reconstructed data after moment1176.8. Moment tensor inversion exampletensor inversion and forward modeling for event reconstruction. For this strong eventthe arrival–weighted normalized energy fit was 0.664, regularization or dampingwas set to 0.01. Also shown are the data after flattening on the P arrival time.The dominant trends in amplitude variations between arrays are captured, withinreceiver arrays some shear arrival variations are not captured correctly, particularlythose at minimum time which correspond to the most horizontal ray paths. Thiscould be due to a breakdown in ray theory where trapped modes are present in thedata but not in the forward model.Figure 6.16 shows the recovered moment signals, the estimated source function,the source function spectrum and the amplitude fit to extract scalar moment andcorner frequency. In the process of going from recording units of pseudo-accelerationto displacement, low frequency noise is blown up, so a low-cut filter below 15Hz hasbeen used. Variability in the estimated moment signals in terms of pulse durationcan be seen in figure 6.16, which probably stems from the difference in S and Pdominant frequencies in the data being mapped to different elements of the momenttensor. Obtaining different signals for each element of m no doubt contains infor-mation about the source, but interpreting what it means is open to conjecture andremains an area for future research.Figure 6.17 shows the inverted PHV radiation patterns with receiver array sam-pling of the focal sphere and the puck glyph to visualize the mechanism. Thisevent shows a predominant slip or double–couple type of mechanism, right–lateralon a sub–vertical plane. Of the two possible solutions the one selected has fracturenormal closest to the normal of the event trend (N45E, vertical plane).This single event result is rather exceptional and not representative of the wholedata set. The event shown in figure 6.18 is more representative, plotted after pro-jection to scalar PHV amplitudes and after time–alignment using calibrated VTImodel times and complex correlation statics. This event is characterized by a strongP arrival (in the direction of recordings relative to the source). Figure 6.19 showsthe results of the inversion, including a display after time–alignment to the directP arrival. The amplitude variations for the P and Sh arrivals are recovered verywell, the Sv not so. As the time–aligned PHV display indicates, the source waveletis significantly different for S compared to P, nevertheless the arrival–weighted VRfor this event is 0.736 so the majority of the radiated amplitude is explained bythe inversion. Optimum regularization for this event was determined using GCVto be 0.15. The PHV waveforms are shown again in figure 6.20 along with theinverted radiation patterns and receiver sampling of the focal sphere. It turns outthat relative to this event the receiver arrays were close to the Sh nodal planes.The inverted moment signals and extracted source function are shown in figure6.21 along with the source spectrum and spectral fit. The medium at the depth ofthis event (z=851m) is very weakly anisotropic (as it is in the sandstone reservoir).The decomposition parameters for this event are EOS = (0.426,−0.177, .397) andSDRO = (−156◦, 13◦,−174◦,−33.7◦), producing the 3D glyph shown in figure 6.22,1186.8. Moment tensor inversion exampleFigure 6.15: Top: Input ENU data (left three panels) and reconstructed ENU dataafter moment tensor inversion and forward modelling for event reconstruction. Thearrival-weighted normalized energy fit is .664. Bottom: The same as top except thewaveforms have been aligned based on the P-arrival time (including statics).1196.8. Moment tensor inversion exampleFigure 6.16: Left: Inverted moment signals and estimated source function; right:source function spectrum and spectral fit to extract scalar moment and corner fre-quency.which shows a composite mechanism with right lateral slip on a sub–vertical planewith a closing component (negative O angle) and positive pressure change (positiveE). The corresponding Hudson parameters are (k, τ, δ) = (0.527, 0.294, 0.178). Theinterpretation of composite mechanisms such as this is a topic of much debate (seediscussion below).The uncertainty in the moment tensor and decomposed parameters can be as-sessed using the multi–variate normal (MVN) sampling technique discussed andshown previously in chapters 4 and 5. The normalized covariance matrix (no damp-ing) for this event is shown in figure 6.23. The components along the trace of themoment tensor are the most poorly resolved, but there are significant off–diagonalelements indicating parameter coupling. 5000 samples honouring the un–normalizedcovariance matrix are shown in figure 6.24. The histograms of the MVN–sampledmoment tensors for this inverted event are also shown in this figure, as well as thenormalized parameter space. Notice where this event plots in normalized parameterspace — it is not in the center as for a pure slip or DC source, nor is it predominantlylike a tensile crack. It is closest to a cylindrical dilatation (+CD) or isotropic plane(+IP). Given that, the EOS decomposition which assumes a model for the sourceof DD + pressure does not bring much intuition about the source. Nevertheless,the EOS glyph is shown in figure 6.22. A CD source has a simple form of diagonal-ized potency tensor (Chapman and Leaney [22], eq. B15b), so if such sources areprevalent in hydraulically stimulated microseisms the DD+pressure interpretationmodel upon which the EOS decomposition is based needs to be reconsidered. Asource that is most like a cylindrical dilatation should be displayed as a cylinder,colored by sign and oriented in the direction of the eigenvector corresponding to theminimum eigenvalue (the long axis).1206.8. Moment tensor inversion exampleFigure 6.17: Radiation patterns including homogeneous VTI Green function for theinverted moment tensor with receiver array focal sphere sampling. Also shown isthe puck–glyph viewed from above to visualize the decomposed source. The greenarrow points North.1216.8. Moment tensor inversion exampleFigure 6.18: Event #36 from the 182 event collection of stage 7, shown after ray–rotation to scalar PHV amplitudes and after time–alignment based on calibratedVTI model times and complex correlation statics.1226.8. Moment tensor inversion exampleFigure 6.19: Event #36 from the 182 event collection of stage 7, shown in inputENU coordinates (left) in the pass–band 23-290Hz and reconstructed after inversion(right three panels). Bottom: the same as top but time–aligned on the direct P timesincluding statics.1236.8. Moment tensor inversion exampleFigure 6.20: Event #36, PHV with inverted radiation patterns shown upper-hemisphere equal–area with relative receiver locations indicated showing the sam-pling of the focal sphere.1246.8. Moment tensor inversion exampleFigure 6.21: Event #36 estimated moment signals and extracted source func-tion (left) where the ordering of moment signals is now, from top to bottom:m11(t),m22(t),m33(t),m23(t),m13(t),m12(t), s(t). Right: source function amplitudespectrum and spectral fit for scalar moment and corner frequency.Figure 6.22: Event #36 3D glyph using the EOS decomposition. The slip termindicates right–lateral slip on a sub–vertical plane with a significant closing (−Ocomponent) and positive pressure change (red wire frame). The poles of the redwire frame are oriented along the T–axis.Figure 6.23: Event #36 normalized covariance matrix. Ordering from upper left ism11,m22,m33,m23,m13,m12.1256.8. Moment tensor inversion exampleFigure 6.24: Event #36 multivariate normal sampling of M (upper left; as his-tograms (lower left); and in diamond crossplot space (right).1266.9. Multi–event results6.9 Multi–event resultsIn this section results for a collection of 182 events from stage 7 are shown. Fig-ure 6.25 shows 182 event locations for stage 7 in depth versus Northing with logs,model and receiver locations. The same events are shown in map view with receiverlocations in figure 6.26.For multiple events, several passes through the data are required to assess crit-ical aspects. Condition number is the first invertibility metric to look at, and isexpected to be high for this geometry near the plane of the monitor array wells.Based on a preliminary spectral analysis such as shown in figure 6.7 an initial band-width estimate was 10-550Hz. Using this bandwidth the inversion was run on all182 events with minimal damping (0.01). Subsequent inspection of estimated sourcespectra indicated a particle displacement signal bandwidth of 23-290Hz. Anotherpass was made with regularization parameter estimation. Figure 6.27 shows con-dition number and damping parameter versus event number. Spikes in conditionnumber correspond to events near the monitor well plane.Figure 6.28 shows inverted source function amplitude spectra and their spectralfits using the ω−3 source function model discussed in chapter 2, section 2.3 andchapter 4, section 4.7.The parameters of the spectral fits shown in figure 6.28 are shown in figure 6.29.They are the corner frequency for the ω−3 decay model and the moment magnitudederived from the scalar moment (the 0–frequency asymptote of the spectral fit).The corner frequency is seen to anti–correlate with moment magnitude, expectedin some source models as a lower frequency corresponds to a longer time duration,greater energy release and hence larger moment magnitude.The quality of fit metric, 1 − ‖d −Gm‖/‖d‖ and the arrival-weighted version(section 4.8.3) are plotted versus event number in figure 6.30. Also shown is themaximum relative moment uncertainty metric, being twice the standard deviationcorresponding to the moment tensor element with the highest posterior variance,divided by the estimated scalar moment.Figure 6.31 shows Hudson and normalized EOS parameters versus event numberwhen the regularization parameter is estimated differently for every event. The EOSparameters indicate that there is almost no O component indicating that dampingmay be over–estimated (median = 0.40). The Hudson, Vavryc˘uk and normalizedE–O parameters are shown in the diamond crossplot in figure 6.32 for “optimum”damping and for an intermediate damping value (0.10).Diamond crossplot results are shown for two other constant values of dampingparameter in figure 6.33, showing the criticality of this parameter when the exper-imental set–up is such that noise maps to model parameter estimates. The singleevent covariance matrix (see figure 6.23) indicated that the trace of M was the mostpoorly resolved, so noise is expected to map into those parameters. This can be1276.9. Multi–event resultsFigure 6.25: Stage 7 event locations in depth versus Northing with logs of Vp,Vp/Vs, model and receivers. The well separation is 884m.1286.9. Multi–event resultsFigure 6.26: Stage 7 event locations in map view (1km x 1km) along with receiverlocations (red dots). The magenta event is event # 36 studied in the previoussection.Figure 6.27: Condition number versus event number for the stage 7 event collection(left); estimated regularization parameter (right).1296.9. Multi–event resultsFigure 6.28: Inverted source function amplitude spectra for 182 events from stage 7(left) and the spectral fits using a ω−3 decay model.Figure 6.29: Source function spectral fit parameters. Left: corner frequency fromthe ω−3 decay model; Right: moment magnitude Mw derived from scalar moment.1306.9. Multi–event resultsFigure 6.30: Fit quality and parameter uncertainty metrics versus event number.Left: conventional energy fit (red) and arrival-weighted energy fit (green); maximumrelative moment tensor error (right).Figure 6.31: Decomposition parameters versus event number. Left: Hudson k, τ, δ(red, green, blue); right: normalized EOS parameters (EOS=Red,Green,Blue).1316.9. Multi–event resultsFigure 6.32: Normalized parameter crossplot for the 182 events of stage 7 showingHudson (green), Vavryc˘uk (blue) and normalized E vs. O (red). Left: optimizeddamping, right: intermediate damping.seen in the reduced damping result, as points have migrated towards the ISO partof the crossplot. Conversely, increased damping moves points away from the ISOvertices, biasing the result to a pure slip or DC source type. The optimum dampingparameter values determined previously have been used for all other multi–eventreal data results shown in this section.Figure 6.34 displays the derived SDRO angles from the DD part of the potencytensor coming from the EOS decomposition. Of the two solutions the chosen onehas fracture plane normal closest to a prior direction, namely the normal to theevent trend plane: N45E and vertical. Recall the D stands for co-dip, so a valueof 0 means a vertical fracture plane. Fracture plane dip is near vertical in bothcases, about 20◦; the strike directions are clustered about but slightly off the eventtrend; rake angles cluster near 0◦ and ±180◦, meaning strike–slip, both left– andright– lateral. Opening angles cluster around 0 but can exceed ±30◦ indicating asignificant tensile component.These results were exported to the Schlumberger Petrel system for visualizationusing the 3D glyphs described and shown in chapter 5. They are shown in figure 6.35in an oblique view that also shows the vertical monitor wells and lateral treatmentwell with the stage 7 port. It is difficult to capture the information in a singledisplay, and one really needs to play a movie along with the pumping data as shownin figure 6.2 to appreciate the spatio–temporal development and the correlation ofevent type with stimulation changes.1326.9. Multi–event resultsFigure 6.33: Normalized parameter crossplots as in figure 6.32 for a reduced (left)and an increased damping parameter.Figure 6.34: SDRO (Strike, co-dip, rake, opening) angles from the DD part ofthe potency tensor coming from the EOS decomposition. Color correspondence isSDRO=Red,Green,Blue,Cyan.1336.9. Multi–event resultsFigure 6.35: 3D display from Petrel of the inversion results shown using the EOSdecomposition and “pucks” glyph. The view is from the south–west and above.Also shown are the two vertical monitor wells with some receiver locations (red andgreen) and the lateral treatment well with a yellow icon at the location of the stage7 port.6.9.1 DiscussionThe results shown in the crossplots of figure 6.32 display an unexpected distribution.Referring to the guide provided by figure 5.1, which shows the location of somespecial sources, the events from stage 7 show a predominant CD or cylindricaldilatation type of source rather than a TC or LVD source type. This is shownclearly in figure 6.36. As a “reality check”, these results were compared to resultsfor events from the same stage using Schlumberger’s commercial MTI (MomentTensor Inversion) code, which is a completely different, time–domain, pick–basedinversion described in Leaney et al. [64]. The Hudson diamond comparison is shownin figure 6.37. The results have a different fit quality distribution (icon size) butthe remarkable distribution in the ±CD quadrants is present for both results. Oneis led to believe that the inversion result may be informative of something in theunderlying process.While one can imagine a dilating cylinder as an analogy for fluid injection caus-ing expansion and connection of a pre-existing fracture network, it is certainly nota commonly proposed interpretation model. The presence of a significant volumecomponent is generally attributed to tensile crack –type (±TC) events. But theconfinement above and below the reservoir by stiffer rock, and the development ofa “fracture fairway” within the reservoir by the stimulation process does evoke the1346.9. Multi–event resultsFigure 6.36: Hudson parameters shown for an optimized regularization parametermulti–event run. Crossplot locations for special sources are shown (for a Poissonmedium).Figure 6.37: Hudson parameters shown color-coded by time and scaled by VR (aquality fit metric) for two different moment tensor inversion codes. The display wasmade with Schlumberger’s Petrel package. Left: the frequency domain anisotropiccode described in this thesis; right: Schlumberger’s commercial time–domain code.1356.9. Multi–event resultsimage of a pipe–like conduit. Such a conduit has a horizontal axis of symmetry likean oriented CD or IP source mechanism. New interpretation models and techniqueswill be needed to properly assess this possibility, for example studying the orien-tation of the eigenvector corresponding to the minimum eigenvalue in moment (forIP) or potency (for CD) domains. Other further work suggested is to look for anydifference in mechanism between the reservoir events and the few events that grewupward. These may have been due to excessive pump rates (John Duhault, pers.comm.). As they occur in the overlying formation that is more clay–rich and moreanisotropic their mechanism deserves careful examination. Finally, as discussed ear-lier in section 2.4.3, it may be that the moment tensor is not the correct descriptionof a microseismic source due to hydraulic stimulation — additional forces may beneeded for a complete description. The full extended source model with twelve pa-rameters has not been explored in this thesis work, but the case of a simple forcesource has.Earlier in figure 6.8 the waveforms for a microseismic event and a sliding sleeveevent were shown. By modifying the Green function for a force instead of a momenttensor source (by removing a phase velocity from the denominator, replacing theray–strain tensor with the polarization vector at the source and reducing the modelvector to length 3), the same MTI (Moment Tensor Inversion) code can be used toinvert the sleeve events for a force source. Figure 6.38 compares the results of aforce source inversion (FSI) to those of a moment tensor inversion (MTI) for thesleeve event. In spite of the fewer degrees of freedom for the FSI it is able to fit thewaveforms significantly better than MTI, producing a VR (Variance Reduction orquality of fit metric) of .343 for FSI versus .069 for MTI. The inverted direction ofthe force was horizontal and pointing toward the south — toward the heel or up thetreatment well.To compare FSI and MTI further, a subset of the largest events for all stagesmonitored by the two receiver arrays was inverted and the quality of fit metric wascompared. Figure 6.39 shows the comparison and indicates where sleeve eventswere confidently identified. Clearly a moment tensor is a poor model of a slidingsleeve source. The non–indigeneous, external energy provided by surface pumpsmay likewise require an extended source model to fully represent and ultimatelyunderstand microseisms due to hydraulic stimulation.1366.9. Multi–event resultsFigure 6.38: Sleeve event ENU data (first three panels) and inversion (right threepanels) for different source models. The simpler force source model (three param-eters) inversion (left) compared to moment tensor inversion (six parameters, right)fits the data significantly better.Figure 6.39: VR fit quality versus event number for a collection of the largest eventsfrom stages 4-13. The FSI model (green) with three fewer degrees of freedom isable to fit the waveforms better than MTI (red) for four events (circled) known withconfidence to be sleeve events.137Chapter 7DirectionsIn this chapter improvements to the inversion and some future directions will bediscussed. Group velocity computations will be shown for an orthorhombic mediumand a cubic medium (salt). Aspects of surface array monitoring will be consideredand finally some data analysis of downhole recordings of aftershocks from the 2007Wenchuan mega–thrust earthquake will be shown.7.1 Inversion improvements7.1.1 Augmented arrivals modelDirect PHV arrivals certainly dominate microseismic recordings but reflections andmode–conversions are present. One can see a downward travelling shear reflectionfrom the top of the receiver array in figures 6.7 and 6.8. Including reflections, headwaves and mode conversions in the forward model is straight–forward but addsconsiderable complexity to the code. There are likely to be multiple occurrencesof such arrivals but one of each type is likely to dominate. Choosing the depthat which such a ray event occurs will require a model with sharper discontinuitiesthan that used for direct arrival computations. The amplitude of the reflection andmode conversion will call for R/T coefficients, placing further demands on having adetailed and accurate model. For these reasons an augmented arrivals model has notbeen pursued, although it would add information to constrain the source inversion.A better sampling of the focal sphere would be provided by the different take–offangles of reflection, mode conversion and head wave ray paths.7.1.2 Green function errorsAs discussed in the section on compute times (4.7.1), results can be generated quicklyand quality of fit and parameter uncertainty metrics can be used to assess results.However, there is no substitute for visual inspection of waveform fit displays suchas 6.15 which are inevitably enlightening. Even the beautiful example of figure6.15 exposes inadequacies. The amplitude variation within arrays is not always fitsatisfactorily, and this may be due to inaccurate spreading computation, incorrect Qcompensation or transmission loss. For stability reasons it has been found necessary1387.1. Inversion improvementsto switch off transmission losses as they are much too large (and often complex) forsub–horizontal rays. Likewise the effective ray–length spreading is used instead ofthe exact form as spreading is unstable for sub–horizontal rays and leads to obviouslyerroneous inversions. The idea put forward in section 3.8.3 to replace layered VTIspreading with homogeneous VTI spreading in an effective medium has not beendeveloped but is a definite possible improvement. Another possibility would beto solve for a parameterized form for spreading in a nonlinear step. Likewise, Qand location perturbations can be considered as free parameters to be determinedin a nonlinear inversion, all aimed at providing the best possible fit to waveformamplitude variations to reduce uncertainty in the estimated moment tensor.It is noted that Dahm (1996) [28] introduced the idea of empirical Green func-tions for moment tensor inversion. In that approach reference event or events withsimilar source mechanism are used as representative and other events are invertedby correlating their waveforms with the reference event. An interesting approachrooted in interferometry, it would avoid issues with Green function computation butnot for the reference event(s). It has not been pursued.7.1.3 Other inversionsIt was mentioned that inverting for the potency tensor rather than the momenttensor is also a linear problem, and it has been investigated (Leaney and Chapman[60]). But since the source model of an isotropic pressure source (moment) is sim-pler than an isotropic volume change (potency), and the subsequent biaxial potencytensor decomposition is simpler starting from moment, the inversion for momenttensor has been preferred. Nevertheless the inversion for potency tensor or a non-linear inversion for the parameters that are used to describe the potency tensor isof interest.Early in this PhD work, a fully nonlinear moment tensor inversion was pursuedin which the parameters were the parameter set described in section 2.6.8, but withthe waveform fit functional computation (wanting to stay away from pick–basedapproaches). It was frustratingly slow and was abandoned.Hybrid linear–nonlinear inversions are also possible. For example, the parame-ters of fracture strike and dip needed to specify the unit normal nˆ may be solvedfor in a 2–parameter nonlinear search with the remaining parameters solved withlinear inversion. Since the fracture plane orientation is often thought to be knowna priori with the highest confidence this may be an attractive direction.A multi-event inversion is attractive as events from different locations provideadditional sampling of the focal sphere. This makes single–well data amenable toMTI given sufficient event distribution. However, the events would all be assumedto have the same mechanism so some clustering criterion would be needed to groupsimilar events, and the differences in scalar moment would need to be accounted for.1397.2. Extensions to lower symmetryFinally, it is worth noting that while the events processed with the techniquesdeveloped in this thesis are so far of an impulsive nature with short duration, sincethe code is a frequency domain inversion it may be well–suited to the mechanismof long–period tremors, an active area of research in microseismics (e.g. Das andZoback [31], Das and Zoback [30]).7.2 Extensions to lower symmetryWhile the ray theory formulation of the forward problem (2.25) is valid for generalanisotropy, thus far all computations have been done for a horizontally layered VTImedium. A dipping layered VTI medium may be handled using the same ray tracerby rotating receiver coordinates about a reference point according to the (constant)geological dip and indeed this is what’s done in Schlumberger’s commercial modelcalibration, event location and moment tensor inversion applications. This addedcomplexity has not been incorporated into the modelling and inversion presented inthis thesis.VTI is without doubt the most important type of anisotropy to handle, but itis not unreasonable to expect an overprint of azimuthal anisotropy. The process ofhydraulic stimulation is, after all, creating fractures which at depth are most likelyvertical, so the next level of complexity or lower symmetry would logically be toadd fractures to the VTI layers. Work has also been done to use the phenomenonof shear wave splitting to invert for fracture compliance parameters assuming anHTI symmetry; Kendall et al. [52] provide a review of that approach. Grechkaand Yaskevich [45] investigated the problem of simultaneously estimating layeredanisotropy and event locations using perforation shot constraints, and Grechka andYaskevich [45] provide parameters for a layered medium with three triclinic layersand two VTI layers.As shown by Schoenberg and Helbig [98], fractures can be added to a VTI back-ground medium by adding compliance to the compliance tensor. When compliantvertical fractures are added to a VTI medium the medium becomes orthorhombic,with an additional symmetry axis normal to the fractures. Schoenberg and Helbig[98] use convenient, dimensionless parameters δN , δV and δH , representing excessnormal compliance and two excess shear compliances, vertical and horizontal, re-spectively. A background VTI stiffness tensor in Voigt notation with subscript b canbe made orthorhombic using the above dimensionless parameters following (Schoen-1407.2. Extensions to lower symmetryberg and Helbig, 1997):C11 = C11b(1− δN ) C12 = C12b(1− δN )C13 = C13b(1− δN ) C22 = C22b(1− δNC212bC211b)C23 = C23b(1− δNC12bC11b) C33 = C33b(1− δNC213bC211bC233b)C44 = C44b C55 = C44b(1− δV )C66 = C66b(1− δH).Such a medium has the fast (i.e. fracture) direction in the 2-direction. This ten-sor can be rotated using the modified Bond transform (Chapman [19], p. 93). Ifit is rotated about the vertical 3-axis it is referred to as “VOR” for “Vertical Or-thoRhombic”; a general rotation is referred to as “TOR” for “Tilted OrthoRhom-bic”. A stiffness tensor thus computed can be used with equations (3.1–3.3) tocompute phase slownesses, polarization vectors and group velocities for qP , qS1and qS2 rays given the two angles specifying a phase direction. By conventionqS1 is taken as faster than qS2 but other labelling strategies may be adopted. Acommon strategy is to label the shears as qSh for the one with most horizontalpolarization and qSv as the other. This works perfectly to handle the fast/slowcross–over for VTI shears since they have unambiguous polarizations and is the ap-proach taken here in view of VTI being the assumed dominant symmetry, but thereis no recipe that works for all symmetries. For example a purely HTI medium isbetter labelled using a fast/slow strategy. Figure 7.1 shows group velocities versusgroup angle in an equal–area projection. The medium is that used by Schoenbergand Helbig [98] with parameters Vp = 2.45 km/s, Vs = 1.41 km/s, ρ = 1 g/cc,Ep = 0.25, Ea = 0.37, Es = 0.20, δN = 0.10, δV = 0.20, δH = 0.27 rotated 45◦ clock-wise about the vertical axis. Vp increases with inclination due to the VTI and isfaster parallel to the fracture plane at N45E. The shear labelling problem is ap-parent. Around the edges of the circle where propagation is horizontal Vsv is fasteralong fracture strike; it is also faster near 45◦ due to the VTI anellipticity parameter(Ea). Around the edges Vsh is seen to be slow both along and normal to fracturestrike and fastest at 45◦ to these directions. In this plane the orthorhombic mediumlooks TI, so this is like the qSv for a VTI medium rotated 90◦. Thankfully, theSchoenberg and Helbig (1997) model is probably extreme, and crack complianceparameters of 1/2 or even 1/4 of these values are more likely in practice.As mentioned in chapter 3, section 3.1, 2–point ray tracing in such a medium(or in a medium with general anisotropy) is a two parameter nonlinear optimizationproblem since in general an additional azimuthal component of phase direction isneeded to determine the group velocity vector that will connect source and receiver.Two–point ray tracing in a homogeneous medium with general anisotropy has beenimplemented but ray tracing in layered general media has not (Snell’s Law requires1417.2. Extensions to lower symmetryFigure 7.1: qP , qSh and qSv group velocities versus group inclination (0 − 90◦)and azimuth shown in an equal area projection with vertical–up in the center. Top:qP , lower left: qSh, lower right: qSv. qP and qS are normalized separately; shearvelocities are plotted with the same color scale.1427.2. Extensions to lower symmetrythe solution of a sextic for general media and shear selection and labelling acrossinterfaces will be needed). The problem of estimating or calibrating such a layeredanisotropic model remains challenging based on the available traveltime samplingprovided by most survey geometries. Inevitably, prior information in the form of softrock physics constraints such as used for depth-dependent VTI model calibration willbe needed. Microseismic processing, including moment tensor inversion in generalanisotropic layers is an important future direction.7.2.1 SaltWhile not thus far prevalent in the stratigraphy of unconventional resources, saltbodies (e.g. Gulf of Mexico, GoM) and salt layering (South Atlantic margins) areprevalent in many economic geological sequences. There are also extensive evapor-ite layers throughout the Middle East that may come into play as unconventionalresources get developed there. The large allochthonous salt bodies in the GoM areknown to be predominantly halite (Raymer et al. [88]) which has cubic symmetry.Cubic symmetry is simpler than that of VTI, which is also called hexagonal, and yetcannot be described by the Schoenberg or Thomsen VTI parameterizations. In acubic medium A11 = A22 = A33, A44 = A55 = A66 and A12 = A13 = A23. A simplecubic anisotropy parameter can be defined based on the deviation from isotropy (M.Slawinski, pers. comm., 2011) as c = (A12−A11 + 2A44)/(A12 +A11− 2A44). Thisallows a medium with differing axial P and S velocities to be imbued with cubicanisotropy.Given the volume of salt in the GoM it is of interest to compute group veloc-ities for halite. This is shown in figure 7.2. qP group velocities are shown versusinclination and azimuth along with those for the closest TI medium (e.g. Bucataruand Slawinski [15]). Notice that, in contrast to a typical clastic sedimentary rock,velocity is less than or equal to the vertical velocity at all angles and azimuths.Even though the maximum difference is comparatively small (about 5% at 60◦), thelarge thickness of salt bodies means that this would have a significant impact onseismic imaging. Little attention has been paid to salt anisotropy and few measure-ments have been made, but those that have (Raymer et al. [89]) are consistent withthis model and it is expected that the emerging tilted orthorhombic imaging algo-rithms (e.g. Thomas et al. [108]) will yield improved seismic images when the trueanisotropy of salt is included. Given the presence of salts in “unconventionals” geo-logical sequences, general anisotropy may be needed in models used for microseismicprocessing as well.1437.3. Extended source modelFigure 7.2: Halite qP group velocities versus group inclination (0−90◦) and azimuthshown in an equal area projection with vertical–up in the center (left). Right: qPgroup velocities versus group inclination angle for all azimuths and those for theclosest VTI model (black line).7.3 Extended source modelThe ray theory formula for the forward problem including force, moment tensorand asymmetric moment tensor (torque) sources was given in section 2.4.3, andin Chapter 6 an example of a force–source inversion was shown to fit a slidingsleeve event better than a moment tensor source. It is possible that microseismicsources caused by hydraulic stimulation require an extended source model due tothe reaction mass being at such a distance (pumps on the surface). Some downholesources may also require a torque term to be included. Including a force source in amore general source model and studying resolution and covariance of such a modelis another future research direction.7.4 Waveform fit model calibrationEarly PhD work focused on this problem. Some of it is described in Leaney [57].A waveform–fit objective function, essentially a norm of the residuals of a linearinversion was minimized for either model unknowns or event location. In the case ofmodel unknowns, VTI parameters, vertical velocity scaling and model smoothnesswere inverted. As with the VTI model calibration approach described in chapter 6,the VTI parameters were constrained by lithology, for example using model Vp/Vsas a proxy for the volume of clay, and hence anisotropy. Another idea is to invert for1447.5. Surface array and shallow grid geometryoptimum model smoothness. This is possible if model dip is assumed known as thepolarization vector at the receivers carries information about the index of refractionacross the last interface. Not surprisingly, optimum model smoothness was foundto depend on bandwidth, with lower frequency data requiring smoother models.Due to the strong nonlinearities and trade–offs in model inversion several variantsof global search techniques were used, including Monte Carlo Markov chain and ahybrid annealing+neighborhood algorithm. These waveform–fit nonlinear inversiondirections were abandoned for being too costly but have merit, especially whencombined with an augmented arrivals model and more general anisotropy, as withmore free parameters to determine polarization information will help to constrainthe solution(s).7.5 Surface array and shallow grid geometryWhile downhole data have been considered throughout this thesis the surface record-ing geometry has been investigated and can be handled. Outstanding issues withsurface monitoring include the proper handling of the free surface and model cal-ibration for near–surface shear properties. The frequency domain, beam–formingapproach described in this thesis may be well–suited to low signal–to–noise sur-face environments provided surface statics and near–surface model calibration canbe achieved. A sparse network of 3C receivers cemented in shallow boreholes is apromising low–noise surface acquisition geometry.Figure 7.3 illustrates the impact of propagation losses and single (vertical) com-ponent recording on P amplitudes recorded on a surface grid due to a pure double–couple source. Clearly including such effects is critical for correct source mechanismrecovery.1457.5. Surface array and shallow grid geometryFigure 7.3: Top: P and S velocities (left) and rays (right) for a source at depthand surface receivers out to 10,000ft offset. Bottom: Simulated P amplitudes fora double–couple source as recorded by a dense grid of receivers. Left: scalar am-plitude for ray take–off angles; right: amplitudes including the effects of spreading,transmission loss, Q (100, for dominant frequency 150Hz), and projection onto thevertical component.1467.6. Wenchuan aftershocks7.6 Wenchuan aftershocksThe method described in Leaney [57] was applied to a unique data set from China(Leaney et al. [69]). Coincidentally, Schlumberger’s first–ever downhole monitoringjob in China was to take place in Sichuan province just as the great Wenchuan7.9 earthquake of May 12th, 2008 devastated the region. The job was delayed amonth, until June 5th. The VSI-7 tool was deployed in a vertical well at a depth ofabout 1500m below ground level. The completion was open–hole and no calibrationshots were acquired for geophone orientation, the thinking being that early eventscould be used to orient the tools and calibrate the model. But there were no eventsrecorded because pump rates were too low (20bpm). Consequently there were doubtsabout whether the equipment was functioning properly, but the processing engineernoticed periodic low frequency noise throughout the job. Subsequent data analysisindicated that the noise was aftershocks, our equipment was functioning properlyand Schlumberger got paid!Figure 7.4 shows the geographical setting — the location of the Wenchuan fault,a slip–history result due to A. Sladen along the fault and the location of the VSItool about 320km from the fault. To process the data, waveforms were decimatedfrom 1ms to 4ms sampling rate and concatenated into 120s records. Geophoneorientation was determined from two large events that were identified as being widelyseparated in azimuth (30◦) and assumed to be coming from either end of the faultregion. Figure 7.5 shows oriented ENU waveforms for the biggest event in stage 4recordings. The approximate time difference Ts − Tp is consistent with a range ofabout 300km. Amplitude spectra are shown for the 66 minutes of stage 4 recordingin figure 7.6, which contained 8 earthquake events. Bandwidth recorded by thedownhole accelerometers was 1-65Hz with a signal–to–noise of 40dB in the 3-10Hzband.A 1D elastic model was built from sonic logs over the receiver array and extendedto 90km depth based on a regional model. To approximate turning ray arrivals,the head wave ray tracing option described in chapter 3 was used (see figure 7.7).Locations were updated for the largest event recorded during stage 4/5 on June6th using the algorithm described in Leaney [57], which assumes an average S/Pradiation pattern ratio. This algorithm provides separated scalar PHV wavefieldsfrom the array of vector ENU data. Figure 7.8 shows results over the P arrival timewindow and over the S arrival time window. Input ENU data are reconstructedusing the inverted phase slowness and polarization of the dominant upward travellingPHV arrivals; residuals are seen to contain downward propagating reflections.A lengthy coda is observed in the separated scalar PHV waveforms. A receiverfunction was not computed, it is not known if crustal reflections are present in thedata. These data should be amenable to a multi–event moment tensor inversion,and because of the emergent arrivals with complicated coda, may be well–suited tothe waveform inversion presented here.1477.6. Wenchuan aftershocksFigure 7.4: The location of the magnitude 7.9 May, 12th, 2008 earthquake andaftershocks slip map with the location of the Schlumberger downhole array from 5-7June, 2008.Figure 7.5: Oriented ENU components for the largest event during stage 4 recording.The horizontal axis is seconds. The boxes simply highlight first P and S phases.1487.6. Wenchuan aftershocksFigure 7.6: Amplitude spectra for 66 minutes of recording during stage 4 of thefailed treatment.1497.6. Wenchuan aftershocksFigure 7.7: 1D elastic model and P and S head wave rays used for processing.1507.6. Wenchuan aftershocksFigure 7.8: Input ENU waveforms, reconstructed ENU waveforms, residual ENUwaveforms and separated scalar PHV waveforms over the P arrival window (top)and the S arrival window (bottom). Vertical axis is in seconds.151Chapter 8ConclusionsAnisotropy is important for seismic source inversion. In this thesis the theory wasdetailed for the forward problem of a moment tensor source and a frequency domain,waveform inversion was developed for moment tensor inversion using anisotropic raytheory. The approach combines the efficiency of ray tracing with the accuracy ofwaveform fitting. A new moment tensor decomposition was presented along with anew way to visualize a decomposed moment tensor. Results were shown for syntheticand field data from central Alberta.Synthetic microseismic data computations made use of dynamic ray tracing ina model composed of homogeneous vertical transverse isotropy (VTI) layers andincluded propagation losses due to geometrical spreading, transmission loss andanelastic absorption due to Q. The anisotropic Green function included impedancecoupling terms at source and receiver requiring group velocity in the ray direction,the phase velocity at the source, polarization vectors at source and receiver and thephase vector at source. Moment tensor construction was done using source geometryparameters to construct the potency tensor, contracted with the stiffness tensor atthe source to produce the moment tensor. Radiation patterns for a moment tensorconstructed in an anisotropic medium were computed and compared to those from anisotropic medium. Source–time function models were reviewed and one with an ω−3decay was used as the source wavelet for synthetics. A network of arbitrary receiverlocations is specifiable and three component vector recordings were simulated, fortwo vertical and one horizontal receiver array, including receiver response. Thesesynthetics were used to test and understand the inverse problem.The theory and method of the ray tracing computations were described for ho-mogeneous general anisotropy and a layered VTI model, and while the model is1D with homogeneous layer properties, the computations are exact, with no weakanisotropy assumption, and dynamic. Rays are qP , Sh or qSv and may containmode–conversions and reflections although only direct transmitted rays were con-sidered for the forward problem of this thesis. The problem of qSv triplication wasinvestigated, and to select a single arrival in the case of triple qSv arrivals the strat-egy of selecting the inner or slowest wavefront was adopted. Two–point ray tracingwas accomplished using a secant search method and a novel Fresnel–based criterionwas used for ray–straightening near source and receiver. This approach makes therays behave more like waves and improves the conditioning of the inverse problem by152Chapter 8. Conclusionsbetter sampling the focal sphere. The simplification of layered VTI models permitsfast, accurate ray tracing for efficient forward and inverse computations.Inversion was done in the frequency domain considering the linear problem forthe six elements of the moment tensor, frequency–by–frequency, using weighteddamped least–squares. Each moment tensor component was free to contain a dif-ferent source function, which may be informative for future source interpretationmodels. Regularization parameter estimation was done using both L–curve andgeneralized cross validation methods. A common source function was extractedfrom the moment tensor signals using the novel approach of complex principal com-ponent analysis. The real, normalized moment tensor was extracted by deconvolvingthe moment tensor signals by the estimated source function, also novel. The scalarmoment and corner frequency were determined from a fit to the source function spec-trum. Uncertainty analysis followed linear inverse theory but utilized multi–variatenormal sampling using the complete covariance matrix. Considering a single event,histograms were made for sampled moment tensors and for decomposed parametersfrom the sampled moment tensors. For multiple events an uncertainty metric wasused and the response of the inversion to increasing noise was shown. A pathologicalsurvey geometry with two vertical arrays was compared to a geometry with threelinear receiver arrays with one being horizontal.Moment tensor decomposition was reviewed and a new decomposition in termsof an isotropic pressure plus a displacement discontinuity (slip plus opening) inpotency — EOS — was described that is valid for general anisotropic media. Itenables the recovery of geometric source parameters free from the distortion dueto anisotropy at the source. Synthetic examples were shown to demonstrate theefficacy of the new decomposition compared to traditional decompositions such asthat due to Hudson et al. [48]. A new 3D graphical representation or “glyph” wasintroduced that makes use of the new decomposition.Field data were processed with the new techniques. The data set came fromwest central Alberta, courtesy of Lightstream Resources and comprised a dual wellmonitoring of a lateral stimulation. Initial model building from sonic logs madeuse of a dominant shear wavelength smoothing and blocking criterion. Automaticlog blocking made use of a cascaded or compound median filtering algorithm andconjugate extrapolation to handle end effects. A linearized inversion was used fordepth–dependent anisotropy estimation where qP, Sh and qSv arrival times werepicked on nine identified sliding sleeve events. Depth–dependent anisotropy was im-posed using a soft rock physics constraint via an auxiliary log (in this case the Vp/Vsratio). The calibrated VTI model was used for event location and subsequent mo-ment tensor inversion. Geophone orientation and event location processes were donein Schlumberger’s commercial microseismic processing package and waveform fileswere exported for moment tensor inversion. 182 events from stage 7 were inverted.Detailed results were shown for two high quality events and attributes were shownfor the event collection. The results do not conform to commonly accepted source153Chapter 8. Conclusionsmodels, dominated by a cylindrical dilatation or isotropic plane type of mechanism,but were checked against a commercial code and found to be consistent. The resultswere exported to Schlumberger’s Petrel interpretation software package where thenew 3D glyph was used to visualize results. The moment tensor inverse problem wasalso modified to consider a force source, and the force source inversion was foundto fit slotted sleeve events better than the moment tensor inversion, with the forcedirection pointing up the well towards the heel.Future directions to come out of this research include ideas to improve the for-ward model, anisotropy of lower symmetry, a new view on source type interpretationgiven the field data results, an extended source model including unbalanced forcesand the potential application of this frequency domain inversion to long duration,emergent events.Although the moment tensor inverse problem is linear, it is not necessarily well–posed, nor is it simple. The added complexity of including anisotropy was shown tobe important. It is hoped that this research will lead others to consider anisotropyin microseismic source inversion, so that the many mysteries that remain in micro-seismology may be explained.154Bibliography[1] Aki, K., and P. Richards, 2002, Quantitative Seismology, 2nd ed.: UniversityScience Books.[2] Aldridge, D. F., L. C. Bartel, N. P. Symons, and N. R. Warpinski, 2003, Gridsearch algorithm for 3D seismic source location: SEG Technical Program Ex-panded Abstracts, 730–733.[3] Alkhalifah, T., and I. Tsvankin, 1995, Velocity analysis in transversely isotropicmedia: Geophysics, 60, 1550–1566.[4] Ammon, C. J., A. A. Velasco, and T. Lay, 1993, Rapid estimation of rupturedirectivity: application to the Landers (Ms=7.4) and Cape Mendocino (Ms=7.2)California earthquakes: Geophys. Res. Lett., 20, 97–100.[5] Aster, R. C., B. Borchers, and C. H. Thurber, 2005, Parameter Estimation andInverse Problems: Elsevier Academic Press.[6] Backus, G. E., 1962, Long–wave elastic anisotropy produced by horizontal lay-ering: Journal of Geophysical Research, 67, 4427–4440.[7] Backus, G. E., and M. Mulcahy, 1976a, Moment tensors and other phenomeno-logical descriptions of seismic sources – I. continuous displacements: Geophys.J.R. Astr. Soc., 46, 341–361.[8] ——–, 1976b, Moment tensors and other phenomenological descriptions of seis-mic sources – II. discontinuous displacements: Geophys. J.R. Astr. Soc., 47,301–329.[9] Baig, A., and T. Urbancic, 2010, Microseismic moment tensors: a path to un-derstanding frac growth: The Leading Edge, 29, 320–324.[10] Beresnev, I. A., and G. M. Atkinson, 1997, Modeling finite-fault radiation fromthe ωn spectrum: Bulletin of the Seismological Society of America, 87, 67–84.[11] Berryman, J. G., 1979, Long–wave elastic elastic anisotropy in transverselyisotropic media: Geophysics, 44, 896–917.[12] Boatwright, J., 1980, A spectral theory for circular seismic sources; simpleestimates of source dimension, dynamic stress drop and radiated seismic energy:Bulletin of the Seismological Society of America, 70, 1–27.[13] Bowers, D., and J. Hudson, 1999, Defining the scalar moment of a seismicsource with a general moment tensor: Bulletin of the Seismological Society ofAmerica, 99, 2801–2814.[14] Brune, J., 1970, Tectonic stress and the spectra of seismic shear waves fromearthquakes: Journal of Geophysical Research, 75, 4997–5009.155Bibliography[15] Bucataru, I., and M. A. Slawinski, 2009, Invariant properties for finding dis-tance in space of elasticity tensors: Journal of Elasticity, 94, 97–114.[16] Buffett, B., 2013, Earth’s enigmatic inner core: Physics Today, 66, 37–41.[17] Burridge, R., and L. Knopoff, 1964, Body force equivalents for seismic disloca-tions: Bulletin of the Seismological Society of America, 54, 1875–1888.[18] Carrion, P., J. Costa, J. Pinheirho, and M. Schoenberg, 1992, Cross–boreholetomography in anisotropic media: Geophysics, 57, 1194–1198.[19] Chapman, C. H., 2004, Fundamentals of Seismic Wave Propagation: Cam-bridge University Press.[20] ——–, 2009, Ray tracing in anisotropic, 1D layered media: Technical report,Schlumberger Cambridge Research.[21] ——–, 2013, Reflection/transmission coefficients and their differentials: Tech-nical report, Schlumberger Gould Research.[22] Chapman, C. H., and W. S. Leaney, 2012, A new moment–tensor decompositionfor seismic events in anisotropic media: Geophysical Journal International, 188,343–370.[23] ——–, 2014, A correction to: A new moment–tensor decomposition for seismicevents in anisotropic media: Geophysical Journal International, in preparation.[24] Chouet, B., P. Dawson, T. Ohminato, M. Martini, G. Saccorotti, F. Giudi-cepietro, G. de Luca, G. Milana, and R. Scarpa, 2003, Source mechanisms ofexplosions at stromboli volcano, italy, determined from moment tensor inversionsof very-long-period data: J. Geophys. Res., 108, 2019–2030.[25] Costa, J., M. Schoenberg, and D. Miller, 1991, Ray tracing in layeredanisotropic media: Presented at the 2nd International Congress of the Brazil-ian Geophysical Society, Salvador.[26] C˘erveny´, V., 1972, Seismic rays and ray intensities in inhomogeneousanisotropic media: Geophys. J.R. Astron. Soc., 29, 1–13.[27] ——–, 2001, Seismic Ray Theory: Cambridge University Press.[28] Dahm, T., 1996, Relative moment tensor inversion based on ray–theory: theoryand synthetic tests: Geophys. J. Int., 124, 245–257.[29] Dahm, T., G. Manthei, and J. Eisenbla¨tter, 1999, Automated moment tensorinversion to estimate source mechanisms of hydraulically induced micro–seismicityin salt rock: Tectonophysics, 306, 1–17.[30] Das, I., and M. D. Zoback, 2013a, Long–period long–duration seismic eventsduring hydraulic stimulation of shale and tight–gas reservoirs — Part 2: Locationand mechanisms: Geophysics, 78, 118–117.[31] ——–, 2013b, Long-period long–duration seismic events during hydraulic stim-ulation of shale and tight–gas reservoirs — Part 1: Waveform characteristics:Geophysics, 78, 97–108.[32] Drew, J., D. Leslie, P. Armstrong, and G. Michaud, 2005, Automated micro-seismic event detection and location by continuous spatial mapping: SPE Pro-ceedings, 95513.156Bibliography[33] Drew, J., R. S. White, F. Tilmann, and J. Tarasewicz, 2013, Coalescence mi-croseismic mapping: Geophys. J. Int., 195, 1773–1785.[34] Du, J., and N. R. Warpinski, 2011, Uncertainty in FPSs from moment tensorinversion: Geophysics, 76, WC65–WC75.[35] Dufumier, H., and L. Rivera, 1997, On the resolution of the isotropic componentin moment tensor inversion: Geophys. J. Int., 131, 595–606.[36] Duhault, J., 2012, Cardium Microseismic West Central Alberta: A CaseHistory: CSEG Recorder, 37, 48–57.[37] Eaton, D. W., 2010, Resolution of microseismic moment tensors: SEG Ex-panded Abstracts, 2789–2793.[38] Eaton, D. W., and C. W. Hogan, 2012, Crack-tip stress field, Coulomb fail-ure, and the spectral characteristics of tensile rupture: Presented at the CSEGGeoconvention.[39] Eisner, L., S. Williams-Stroud, A. Hill, P. Duncan, and M. Thornton, 2010,Beyond the dots in the box: Microseismicity–constrained fracture models forreservoir simulation: The Leading Edge, 29, 326–333.[40] Farquharson, C. G., and D. W. Oldenburg, 2004, A comparison of automatictechniques for estimating the regularization parameter in non–linear inverse prob-lems: Geophys. J. Int., 156, 411–425.[41] Federov, F. D., 1968, Theory of elastic waves in crystals, Trans: Bradley,J.E.S.,(original in Russian, Nauka Press, Moscow, 1965): Plenum Press, NewYork.[42] Fisher, T., and A. Guest, 2011, Shear and tensile earthquakes caused by fluidinjection: Geophysical Research Letters, 38, 05307.[43] Folstad, P. G., and M. Schoenberg, 1992, Low frequency propagation throughfine layering: SEG Expanded Abstracts, 1279–1281.[44] Freire, S. L. M., and T. J. Ulrych, 1988, Application of singular value decom-position to vertical seismic profiling: Geophysics, 53, 778–785.[45] Grechka, V., and S. Yaskevich, 2013, Inversion of microseismic data for triclinicveolocity models: Geophysical Prospecting, 61, 1159–1170.[46] Hake, J. H., E. C. A. Gevers, C. M. van der Kolk, and B. W. Tichelar, 1998, Ashear experiment over the Natih field in Oman: pilot seismic and borehole data:Geophysical Prospecting, 46, 617–646.[47] Han, J. Y., B. H. W. van Gelder, and T. Soler, 2007, On covariance propagationof eigenparameters of symmetric n–D tensors: Geophysical Journal International,170, 503510.[48] Hudson, J. A., R. G. Pearce, and R. M. Rogers, 1989, Source type plot forinversion of the moment tensor: J. Geophys. Res., 94, 765–774.[49] Jost, M. L., and R. B. Herrmann, 1989, A student’s guide to and review ofmoment tensors: Seism. Res. Lett., 60, 37–57.[50] Julian, B. R., A. D. Miller, and G. R. Foulger, 1998, Non–double–couple earth-quakes: Theory: Reviews of Geophysics, 36, 525–549.157Bibliography[51] Kawakatsu, H., and J. Montagner, 2008, Time–reversal seismic source imagingand moment tensor inversion: Geophysical Journal International, 175, 686–688.[52] Kendall, J.-M., J. P. Verdon, and A. F. Baird, 2014, Evaluating fracture–induced anisotropy using borehole microseismic data: CSEG Recorder, 39, 56–63.[53] King, G., 2012, Hydraulic fracturing 101: What every representative, envi-ronmentalist, regulator, reporter, investor, university researcher, neighbor andengineer should know about estimating frac risk and improving frac performancein unconventional gas and oil wells: Presented at the SPE Proceedings.[54] Leaney, W. S., 2000, Look–ahead walkaways using effective VTI models: SEGTechnical Program Expanded Abstracts, 1743–1747.[55] ——–, 2001, Seismic method of performing the time picking step, WO patent01/88570 A1: Technical report, Schlumberger Limited.[56] ——–, 2002, Anisotropic vector plane wave decomposition for 3D VSP data:SEG Expanded Abstracts, 2369–2372.[57] ——–, 2008a, Microseismic inversion by least–squares time reversal and wave-form fitting: SEG Expanded Abstracts, 1347–1351.[58] ——–, 2008b, Polar anisotropy from walkaways: The Leading Edge, 27, 1242–1250.[59] Leaney, W. S., and C. H. Chapman, 2010, Microseismic sources in anisotropicmedia: EAGE Expanded Abstracts, F014.[60] ——–, 2013, A ray+waveform inversion for the potency tensor: Presented atthe CSEG geoConvention: Integration, Calgary.[61] ——–, 2014, Anisotropic ray–waveform moment tensor inversion: EAGE Ex-panded Abstracts, E108.[62] Leaney, W. S., C. H. Chapman, and T. J. Ulrych, 2011, Microseismic sourceinversion in anisotropic media: Presented at the CSEG Geoconvention: Recovery,Calgary.[63] Leaney, W. S., C. H. Chapman, and X. Yu, 2014a, Anisotropic moment tensorinversion and visualization: Presented at the SEG Annual Meeting, Denver.[64] Leaney, W. S., C. H. Chapman, X. Yu, L. Bennett, S. Maxwell, J. Rutledge,and J. Duhault, 2014b, Anisotropic moment tensor inversion and visualizationapplied to a dual well monitoring survey: CSEG Recorder, 39, 48–54.[65] Leaney, W. S., and R. R. Hope, 1998, Borehole–guided long offset AVO pro-cessing for improved lithology classification: SEG Technical Program ExpandedAbstracts, 230–233.[66] Leaney, W. S., and B. E. Hornby, 2007, Depth–dependent anisotropy fromsub–salt walkaway VSP data: EAGE Expanded Abstracts, H011.[67] Leaney, W. S., M. D. Sacchi, and T. J. Ulrych, 2009, Station–consistent decon-volution for multi-source borehole seismic data: Journal of Seismic Exploration,18, 229–237.[68] Leaney, W. S., and T. J. Ulrych, 1992, Compound median filtering and the me-dian decomposition: log processing applications: Journal of Seismic Exploration,158Bibliography1, 173–190.[69] Leaney, W. S., A. Voskamp, L. Bennett, M. Craven, and Y. J. Li, 2008, Highfrequency downhole recordings of sichuan aftershocks: Presented at the AmericanGeophysical Union, Fall Meeting 2008, abstract S11D-0.[70] Leaney, W. S., X. Yu, C. H. Chapman, L. Bennett, S. Maxwell, J. Rutledge,and J. Duhault, 2014c, Anisotropic moment tensor inversion and visualizationapplied to a dual well monitoring survey: Presented at the CSEG Geoconvention:Focus, Calgary.[71] Liner, C., and T. Fei, 2007, The Backus number: The Leading Edge, 26, 420–426.[72] Madariaga, R., 1976, Dynamics of an expanding circular fault: Bulletin of theSeismological Society of America, 66, 639–667.[73] ——–, 2007, 4.02 seismic source theory, in: Treatise on Geophysics, H.Kanamori Ed.: Elsevier.[74] Maxwell, S., 2014, Microseismic imaging of hydraulic fracturing: Improved en-gineering of unconventional shale reservoirs: SEG Distinguished Instructor SeriesNo. 17: Society of Exploration Geophysicists.[75] Maxwell, S., Z. Chen, I. Nizkous, R. Parker, Y. Rodionov, and M. Jones, 2011,Microseismic evaluation of stage isolation with a multiple–fracport, openhole com-pletion: Presented at the CSUG/SPE.[76] Menke, W., 1984, Geophysical Data Analysis: Discrete Inverse Theory:Academinc Press, Elsevier.[77] ——–, 1989, Geophysical Data Analysis: Discrete Inverse Theory, SecondEdition: Academinc Press, Elsevier.[78] ——–, 2012, Geophysical Data Analysis: Discrete Inverse Theory Matlab edi-tion, Third Edition: Academic Press, Elsevier.[79] Michaud, G., and W. S. Leaney, 2008, Continuous microseismic mapping forreal–time event detection and location: SEGExpandedAbstracts, 1357–1361.[80] Miller, D. E., W. S. Leaney, and W. H. Borland, 1994, An in situ estimation ofanisotropic elastic moduli for a submarine shale: Journal of Geophysical Research,99, 21659–21665.[81] Mizuno, T., W. S. Leaney, and G. Michaud, 2010, Anisotropic velocity modelinversion for imaging the microseismic cloud: EAGE Expanded Abstracts, F014.[82] Nolen-Hoeksema, R. C., and L. J. Ruff, 2001, Moment tensor inversion of micro-seisms from the b–sand propped hydrofracture, m-site, colorado: Tectonophysics,336, 163–181.[83] Papakonstantinou, J., 2007, The Historical development of the secant methodin 1–D: Presented at the Annual meeting of the Mathematical Association ofAmerica, The Fairmont Hotel, San Jose, CA.[84] Prasad, M., 2013, Introduction to unconventional gas: CSEG Recorder, 38,20–25.[85] Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, 1992,159BibliographyNumerical Recipes in Fortran 77, 2nd ed.: Cambridge University Press.[86] Ps˘enc˘´ık, I., 1998, ANRAY Package, version 4.10.: Dept. of Geophysics, CharlesUniversity, Prague, Czech Republic, 148, 403–404.[87] Ps˘enc˘´ık, I., and T. N. Teles, 1996, Point source radiation in inhomogeneousanisotropic structures: Pure Appl. Geophys., 148, 591–623.[88] Raymer, D., A. Tommasi, and J. M. Kendall, 2000a, Predicting the seismicimplications of salt anisotropy using numerical simulations of halite deformation:Geophysics, 65, 1272–1280.[89] Raymer, D. G., J. M. Kendall, D. Pedlar, R. R. Kendall, M. C. Mueller, andG. J. Beaudoin, 2000b, The significance of salt anisotropy in seismic imaging:SEG Technical Program Expanded Abstracts, 562–565.[90] Rodriguez, I. V., Y. Gu, and M. D. Sacchi, 2011a, Resolution of seismic–moment tensor inversions from a single array of receivers: Bull. Seism. Soc. Am.,101, 2634–2642.[91] Rodriguez, I. V., Y. J. Gu, and M. D. Sacchi, 2011b, Simultaneous recovery oforigin time, hypocenter location and seismic moment tensor using sparse repre-sentation theory: Geophys. J. Int., 188, 1188–1202.[92] Ro¨ssler, D., F. Kru¨ger, I. Ps˘enc˘´ık, and G. Ru¨mpker, 2007a, Retrieval of sourceparameters of an event of the 2000 West Bohemia earthquake swarm assumingan anisotropic crust: Stud. Geophys. Geod., 51, 231–254.[93] Ro¨ssler, D., F. Kru¨ger, and G. Ru¨mpker, 2007b, Retrieval of moment tensorsdue to dislocation point sources in anisotorpic media using standard techniques:Geophysical Journal International, 169, 136–148.[94] Sato, H., and T. Hirasawa, 1973, Body wave spectra from propagating shearcracks: Journal of the Physics of the Earth, 84, 829–841.[95] Schoenberg, M., 1994, Transversely isotropic media equivalent to thin isotropiclayers: Geophysical Prospecting, 42, 885–915.[96] Schoenberg, M., and T. Daley, 2003, qSv wavefront triplication in a transverselyisotropic material: SEG Expanded Abstracts, 137–140.[97] Schoenberg, M., and M. V. de Hoop, 2000, Approximate dispersion relationsfor qP-qSv–waves in transversely isotropic media: Geophysics, 65, 919–933.[98] Schoenberg, M., and K. Helbig, 1997, Orthorhombic media: Modeling elasticwave behavior in a vertically fractured earth: Geophysics, 62, 1954–1974.[99] Schoenberg, M., and F. Muir, 1989, A calculus for finely layered anisotropicmedia: Geophysics, 54, 581–589.[100] Schoenberg, M., and J. Protazio, 1992, ’Zoeppritz’ rationalized and generalizedto anisotropy: Journal of Seismic Exploration, 1, 125–144.[101] Silver, P. G., and T. H. Jordan, 1982, Optimal estimation of scalar seismicmoment: Geophys. J.R. Astr. Soc., 70, 755–787.[102] Snelling, P. E., M. de Groot, and K. Hwang, 2013, Characterizing hydraulicfracture behaviour in the HornRiverBasin with microseismic data: SEG TechnicalProgram Expanded Abstracts, 4502–4507.160Bibliography[103] Song, F., and M. N. Toksoz, 2011, Full waveform based complete momenttensor inversion and cource parameter estimation from downhole microseismicdata for hydrofracture monitoring: Geophysics, 76, 103–116.[104] Song, X., and P. G. Richards, 1996, Evidence of inner core anisotropy fromtravel time observations: Nature, 382, 221.[105] Tarantola, A., 2005, Inverse Problem Theory and methods for model param-eter estimation: SIAM.[106] Tarantola, A., and B. Valette, 1982a, Generalized non–linear inverse problemssolved using the least–squares criterion: Rev. Geophys. Space. Phys., 20, 219–232.[107] ——–, 1982b, Inverse problems = Quest for information: J. Geophys., 50,159–170.[108] Thomas, M., S. Mothi, and P. McGill, 2012, Improving Subsalt images us-ing tilted–orthorhombic RTM in Green Canyon, Gulf of Mexico: SEG TechnicalProgram Expanded Abstracts, 1–5.[109] Thomsen, L., 1986, Weak elastic anisotropy: Geophysics, 51, 1954–1966.[110] Thomsen, L., and J. Dellinger, 2003, On shear wave triplication in transverselyisotropic media: Journal of Applied Geophysics, 54, 289–296.[111] Tomic, J., R. E. Abercrombie, and A. F. do Nascimento, 2009, Source pa-rameters and rupture velocity of small M ≤ 2.1 reservoir induced earthquakes:Geophysical Journal International, 179, 1013–1023.[112] Trickett, S. R., 2003, F-xy eigenimage noise suppression: Geophysics, 68,751–759.[113] Trifu, C.-I., D. Angus, and V. Shumila, 2000, A fast evaluation of the seismicmoment tensor for induced seismicity: Bull. Seism. Soc. Am., 90, 1521–1527.[114] Ulrych, T. J., M. D. Sacchi, and S. L. M. Freire, 1999, Eigen image processingof seismic sections: Covariance analysis for seismic signal processing:: Soc. Expl.Geophys., Kirlin, R. L., and Done, W. J., Eds., 53, 241–274.[115] Ulrych, T. J., M. D. Sacchi, M. Graul, and T. Taner, 2007, Instantaneousattributes: the what and the how: Exploration Geophysics, 38, 213–219.[116] Urbancic, T. I., C.-I. Trifu, R. A. Mercer, A. J. Feustel, and J. A. G. Alexander,2000, Automatic time–domain calculation of source parameters for the analysisof induced seismicity: Bull. Seism. Soc. Am., 86, 16271633.[117] Ursin, B., and K. Hokstad, 2003, Geometrical spreading in a layered trans-versely isotropic medium with vertical axis of symmetry: Geophysics, 68, 2082–2091.[118] Valero, H. P., T. Ikegami, B. Sinha, S. Bose, and T. Plona, 2009, Sonic disper-sion curves identify TIV anisotropy in vertical wells: SEG Expanded Abstracts,381–385.[119] Varela, C. L., A. L. R. Rosa, and T. J. Ulrych, 1993, Modeling of attenuationand dispersion: Geophysics, 58, 1167–1173.[120] Vavryc˘uk, V., 2001, Inversion for parameters of tensile earthquakes: J. Geo-phys. Res., 106, 16339–16355.161Bibliography[121] ——–, 2005, Focal mechanisms in anisotropic media: Geophys. J. Int., 161,334–346.[122] ——–, 2007, On the retrieval of moment tensors from borehole data: Geophys.Prosp., 55, 381–391.[123] ——–, 2011, Tensile earthquakes: Theory, modeling and inversion: J. Geo-phys. Res., 116.[124] Wahba, G., 1990, Spline models for observational data: SIAM, Philadelphia.[125] Walter, W. R., and J. N. Brune, 1993, Spectra of seismic radiation from atensile crack: Journal of Geophysical Research, 98, 4449–4459.[126] Washbourne, J. K., K. P. Bube, P. Carillo, and C. Addington, 2008, Wavetracing: Ray tracing for the propagation of band-limited signals: Part 2 — appli-cations: Geophysics, 73, VE385–VE393.[127] Zoeppritz, K., 1919, U¨ber erdbebbenwellen vii b. U¨ber reflexion und durch-gang seismischer wellen durch unstetigkeitsfla¨chen: No¨chr. der Ko¨niglichen Gesell.Wiss. Go¨ttingen Math.-Phys. KI., 66–84.162
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Microseismic source inversion in anisotropic media
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Microseismic source inversion in anisotropic media Leaney, W. Scott 2014
pdf
Page Metadata
Item Metadata
Title | Microseismic source inversion in anisotropic media |
Creator |
Leaney, W. Scott |
Publisher | University of British Columbia |
Date Issued | 2014 |
Description | Sedimentary rocks and shales in particular are known to be anisotropic, sometimes strongly so, and hydraulic fracturing is now common practice in shale plays to enhance the extraction of hydrocarbons. One of the ways to understand the hydraulic fracturing process is through the micro--earthquakes that it generates; it is therefore of interest to study the impact that anisotropy may have on hydraulically induced seismicity. This thesis is concerned with the inclusion of anisotropy into the geophysical forward and inverse problems of microseismic sources. Ray theory is used for the forward problem --- dynamic ray tracing in a medium composed of homogeneous layers with vertical transverse isotropy (VTI) is used and the possibility of qSv triplication is considered. Novel approaches to the inverse problem are introduced, including waveform fitting in the frequency domain to recover the source function as well as moment tensor. Uncertainties in estimated moment tensor components are quantified with multi--variate normal sampling utilizing the full covariance matrix from the linear inversion. A new decomposition of the moment tensor is described that removes the distortions due to anisotropy local to the source and a new way to visualize an earthquake source is also introduced. The impact of anisotropy on moment tensor inversion and decomposition is shown to be significant. Field data recorded by two downhole arrays of multicomponent receivers are analyzed using the new techniques. The event collection suggests a source mechanism dominated by cylindrical dilatation. This is unexpected but is supported by results from another code. Future directions include extensions to lower symmetry forms of anisotropy and application to surface arrays. Preliminary analysis of downhole recordings of aftershocks of the 2008 Mw=7.9 Wenchuan earthquake is shown. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2014-08-19 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0165970 |
URI | http://hdl.handle.net/2429/50045 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2014-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2014_september_leaney_walter.pdf [ 23.31MB ]
- Metadata
- JSON: 24-1.0165970.json
- JSON-LD: 24-1.0165970-ld.json
- RDF/XML (Pretty): 24-1.0165970-rdf.xml
- RDF/JSON: 24-1.0165970-rdf.json
- Turtle: 24-1.0165970-turtle.txt
- N-Triples: 24-1.0165970-rdf-ntriples.txt
- Original Record: 24-1.0165970-source.json
- Full Text
- 24-1.0165970-fulltext.txt
- Citation
- 24-1.0165970.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0165970/manifest