Developments in waveform tomography ofland seismic datawith applications in south-central British ColumbiabyBrendan Robert SmithymanB.Sc., Queen?s University, 2007A 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)October 2013? Brendan Robert Smithyman 2013AbstractWaveform tomography (WT) is an advanced class of seismic imaging methods that may be ap-plied to generate detailed models of subsurface velocity and attenuation by analysis of recordedfield seismograms. WT comprises two steps: 1) initial velocity model building by traveltimeinversion (raytracing) and 2) full-waveform inversion to update the model of velocity based onanalysis of refracted data waveforms. The nature of this method requires that the initial modeland simulated survey parameters reproduce closely those used for the field data acquisition. Theapplication of WT to on-land seismic data is challenged by the geometry irregularities inherentin crooked-line two-dimensional (2D) seismic data acquisition. This thesis develops two method-ologies that permit the inversion of on-land crooked-line vibroseis multi-channel seismic (MCS)reflection data by traveltime tomography (TT) and subsequent viscoacoustic full-waveform inver-sion (FWI) to generate 2D cross sections of P-wave velocity. The first method uses 2.5D TT andapplies a TT-derived static correction to enable the subsequent application of 2D FWI. The sec-ond method uses 2.5D or three-dimensional (3D) TT followed by 2.5D full-waveform inversionthat models 3D survey geometry (a new development). These technical developments are mo-tivated and tested by a multi-part case study, which analyses data from the 2008 GeoscienceBC Nechako Basin seismic survey. Both the static-correction method and 2.5D FWI method areapplied to a single acquisition line. The results are analysed and compared with the conclusionthat 2.5D WT provides a superior result, but at additional computational cost. This cost is stillsignificantly lower than that associated with 3D FWI. The 2.5D WT method is applied to threedatasets from the 2008 survey, and TT alone is applied to a fourth. The resulting models constrainthe positions and lithology of Eocene volcanic rock units, Cretaceous strata and Jurassic or earlierbasement rocks in the Nechako-Chilcotin Plateau region of south-central British Columbia. Theseresults are synthesized and interpreted in conjunction with colocated surficial-geology maps, ex-iiploration well logs and geophysical results by other workers. This interpretation enhances thespecific knowledge of lithostratigraphy along these acquisition profiles and informs general geo-logical and geophysical analyses of the region.iiiPrefaceAll of the work presented herein was carried out at The University of British Columbia byme, with advice, supervision and editorial authorship by my supervisor, R. M. Clowes. My thesiscommittee (M. Bostock, A. J. Calvert and L. Kennedy) provided editorial feedback on initial draftsof the published papers and contributed comments and suggestions during committee meetingsand in personal discussions. The technical thesis objectives and design of the research programwere determined by me in consultation with R. M. Clowes. These design choices were made inresponse to challenges posed in the processing and interpretation of the seismic data from thecase study (2008 Geoscience BC Nechako Basin dataset). I was not involved in the survey designof the case study.Chapters 1 and 2 were prepared by me and represent my original work with the exception ofFigure 1.3, which is reproduced unmodified from Calvert et al. (2011).Chapter 3 was published inGeophysics as Smithyman and Clowes (2012), entitled ?Waveformtomography of field vibroseis data using an approximate 2D geometry leads to improved velocitymodels?. It is reproduced unmodified, except for formatting changes and minor corrections forclarity. I was the lead author and was responsible for design of the research program, analysisof the data and composition of the paper. R. M. Clowes provided feedback, advice and editorialsupervision throughout the process. Figure 3.7 compares my results with those of others andreproduces modified images from their work (WesternGeCo MDIC, 2010; Calvert et al., 2011).Chapter 4 was published in the Journal of Geophysical Research: Solid Earth as Smithymanand Clowes (2013), entitled ?Waveform tomography in 2.5D: Parameterization for crooked-lineacquisition geometry?. It is reproduced unmodified, except for formatting changes and minorcorrections for clarity. I was the lead author and was responsible for design of the researchprogram, analysis of the data and composition of the paper. R. M. Clowes provided feedback,ivadvice and editorial supervision throughout the process. Some figures were modified from theprevious paper/chapter.Chapter 5 is my original work, with the exception of Figure 5.1 (modified from Calvert et al.,2011). I am the lead author and am responsible for design of the research program, analysis ofthe data and composition of the paper. R. M. Clowes provided feedback, advice and editorial su-pervision throughout the development of this chapter. Figures 5.5?5.9, 5.12?5.13 and 5.16?5.17reproduce PostSTM images from CGG Veritas processing of the case study data, first summarizedby Calvert et al. (2011), which I compare with my results. Figure 5.8 reproduces well-log inter-pretations fromMwenifumbo and Mwenifumbo (2010) and Riddell (2011), which I compare withmy results. Structural interpretations were prepared by me in consultation with E. Bordet [PhDstudent (Geology) in EOAS at UBC].Chapter 6 is my original work. R.M. Clowes provided advice and comments.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Research objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Traveltime inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 FAST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.3 Full-waveform inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.3.1 Wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.3.2 Forward modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.3.3 Velocity model inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.3.4 Logarithmic L2 misfit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121.3.5 Iterative Tikhonov regularization . . . . . . . . . . . . . . . . . . . . . . . . 141.3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151.4 Vibroseis acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15vi1.5 Difficulties posed by on-land seismic acquisition . . . . . . . . . . . . . . . . . . . 161.6 The Nechako Basin dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181.6.1 Geometry and acquisition parameters . . . . . . . . . . . . . . . . . . . . . 211.7 Geological context and background . . . . . . . . . . . . . . . . . . . . . . . . . . 212 Technical Developments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.1 Handling of coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.2 Traveltime tomography procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.3 First-arrival picking strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.4 Source signature estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313 Waveform Tomography of Field Vibroseis Data using an Approximate 2D Geometryleads to Improved Velocity Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.2 Technical background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.3 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.4 Case study: Nechako Basin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.4.1 Geological background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.4.2 Traveltime inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.4.3 Full-waveform inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.4.4 Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564 Waveform Tomography in 2.5D:Parameterization for Crooked-Line Acquisition Geometry . . . . . . . . . . . . . . . . 574.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 614.3.1 Traveltime inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61vii4.3.2 Full-waveform forward modeling . . . . . . . . . . . . . . . . . . . . . . . 624.3.3 Full-waveform inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.3.4 2.5D, AVO, source and attenuation considerations . . . . . . . . . . . . . 674.4 Case study: Nechako Basin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.4.1 Geology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.4.2 Survey parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.4.3 Model-building procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.4.5 Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 914.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 944.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 965 New Geophysical Models for Subsurface Velocity Structure in the Nechako?ChilcotinPlateau from 2.5D Waveform Tomography . . . . . . . . . . . . . . . . . . . . . . . . . 975.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 975.1.1 Technical background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.1.2 Seismic lines and related works . . . . . . . . . . . . . . . . . . . . . . . . 1015.2 Waveform tomography results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.2.1 Line 05: southern Nazko . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1055.2.2 Line 06: northern Baezaeko . . . . . . . . . . . . . . . . . . . . . . . . . . 1155.2.3 Line 10: southern Baezaeko . . . . . . . . . . . . . . . . . . . . . . . . . . 1165.2.4 Line 15: Tibbles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1245.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1315.3.1 Attenuation inversion and data fit . . . . . . . . . . . . . . . . . . . . . . . 1315.3.2 Fence diagram interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . 1325.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1396 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1416.1 Fulfillment of thesis objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141viii6.2 Assessment of research impact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1436.2.1 Methodological contributions . . . . . . . . . . . . . . . . . . . . . . . . . 1436.2.2 Geological contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1456.3 Future research directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1466.3.1 Geological constraints and regularization . . . . . . . . . . . . . . . . . . . 1466.3.2 Quality control and assessment . . . . . . . . . . . . . . . . . . . . . . . . 1476.3.3 Automation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149AppendicesAppendix A: Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156A.1 Analysis and publication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156A.2 2.5D full-waveform inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156A.3 Open source software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157A.3.1 pygeo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157A.3.2 waveserv . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157A.3.3 multiplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158ixList of Tables4.1 Relative reduction of objective function . . . . . . . . . . . . . . . . . . . . . . . . 89xList of Figures1.1 The idealized seismic cross-section . . . . . . . . . . . . . . . . . . . . . . . . . . 171.2 Two possible solutions are shown . . . . . . . . . . . . . . . . . . . . . . . . . . . 191.3 An overview of the Nechako Basin . . . . . . . . . . . . . . . . . . . . . . . . . . 202.1 Four simplified cases are presented . . . . . . . . . . . . . . . . . . . . . . . . . . 252.2 An example shot gather is presented . . . . . . . . . . . . . . . . . . . . . . . . . . 282.3 Panels of a) original first-arrival picks . . . . . . . . . . . . . . . . . . . . . . . . . 302.4 Several source wavelets are shown . . . . . . . . . . . . . . . . . . . . . . . . . . 313.1 Relevant geological terranes in central British Columbia . . . . . . . . . . . . . . . 403.2 Geometry of seismic line 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.3 A schematic stratigraphic column . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.4 Traveltime and residual plots are shown . . . . . . . . . . . . . . . . . . . . . . . . 453.5 Velocity models from traveltime inversion . . . . . . . . . . . . . . . . . . . . . . 463.6 Real (top) and synthetic (bottom) data . . . . . . . . . . . . . . . . . . . . . . . . . 473.7 (a) The velocity model from Figure 3.5b is shown . . . . . . . . . . . . . . . . . . 514.1 Overview of geological terranes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.2 Surface lithology and the presence of Quaternary cover . . . . . . . . . . . . . . . 674.3 The natural logarithm of data amplitude . . . . . . . . . . . . . . . . . . . . . . . . 734.4 (a) The starting velocity model from traveltime inversion . . . . . . . . . . . . . . . 764.5 (a) Velocity difference from 2D FWI . . . . . . . . . . . . . . . . . . . . . . . . . . 774.6 Sensitivity kernels are displayed for pairs . . . . . . . . . . . . . . . . . . . . . . . 824.7 A collection of source gathers from source 200 . . . . . . . . . . . . . . . . . . . . 844.8 A collection of source gathers from source 350 . . . . . . . . . . . . . . . . . . . . 85xi4.9 A collection of source gathers from source 600 . . . . . . . . . . . . . . . . . . . . 864.10 Frequency-domain panels showing amplitude . . . . . . . . . . . . . . . . . . . . 884.11 The data residual is shown for all interations . . . . . . . . . . . . . . . . . . . . . 904.12 Velocity and attenuation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 925.1 The surficial geology of the study area . . . . . . . . . . . . . . . . . . . . . . . . . 995.2 The 2008 Geoscience BC Nechako Basin seismic survey comprised . . . . . . . . 1035.3 Line 05; starting (a) and final (b) velocity models . . . . . . . . . . . . . . . . . . . 1065.4 Line 05; field data (black) and synthetic data (green) . . . . . . . . . . . . . . . . . 1075.5 Line 05; the PostSTM image from Calvert et al. (2009) is overlaid . . . . . . . . . . 1095.6 Line 05; the PostSTM image from Calvert et al. (2009) is overlaid . . . . . . . . . . 1105.7 Well sonic logs from exploration wells . . . . . . . . . . . . . . . . . . . . . . . . 1115.8 Well-log interpretations are reproduced . . . . . . . . . . . . . . . . . . . . . . . . 1125.9 Line 06; (a) The model of P-wave velocity . . . . . . . . . . . . . . . . . . . . . . . 1185.10 Line10; starting (a) and final (b) velocity models . . . . . . . . . . . . . . . . . . . 1195.11 Line 10; field data (black) and synthetic data (green) . . . . . . . . . . . . . . . . . 1205.12 Line 10; the PostSTM image from Calvert et al. (2009) is overlaid . . . . . . . . . . 1215.13 Line 10; the PostSTM image from Calvert et al. (2009) is overlaid . . . . . . . . . . 1225.14 Line 15; starting (a) and final (b) velocity models . . . . . . . . . . . . . . . . . . . 1255.15 Line 15; field data (black) and synthetic data (green) . . . . . . . . . . . . . . . . . 1275.16 Line 15; the PostSTM image from Calvert et al. (2009) is overlaid . . . . . . . . . . 1285.17 Line 15; the PostSTM image from Calvert et al. (2009) is overlaid . . . . . . . . . . 1295.18 The phase, amplitude and total data residuals . . . . . . . . . . . . . . . . . . . . . 1335.19 The phase, amplitude and total data residuals . . . . . . . . . . . . . . . . . . . . . 1345.20 Fence diagram depiction of waveform tomography velocity models . . . . . . . . . 1375.21 Fence diagram depiction of interpreted lithostratigraphy . . . . . . . . . . . . . . . 138xiiAcknowledgementsThanks first and foremost to Ron Clowes, my supervisor. It has been an absolute privilege.Thank you to my outstanding PhD committee, Michael Bostock, Andy Calvert and Lori Kennedyfor your support and guidance. Thank you to Esther Bordet for sharing your maps and expertise,which greatly improved the quality of my geophysical interpretation in Chapter 5. Thank youto Gerhard Pratt, Rie Kamei, Drew Brenders, and Michael Afanasiev at the University of West-ern Ontario for endless discussions and support; the technical work in this thesis is built on thefoundations provided by your work and the work of our colleagues in the field. I appreciate thehelp and feedback from many others in the research community with whom I have spoken atconferences and meetings. Thank you to the editors and reviewers at Geophysics and JGR: SolidEarth for constructive comments that greatly improved the quality of my manuscripts.Funding for the early part of this project was provided by Geoscience BC. This work was alsosupported by the Natural Sciences and Engineering Research Council of Canada (NSERC) througha Canada Graduate Scholarship (Doctoral) to me and a Discovery grant to Ron. UBC providedFour-Year Fellowship funding at the end of my NSERC term. Computational facilities for thisresearch were provided by WestGrid, part of Compute Canada, which was absolutely critical tothe timely completion of this research. I appreciate the license provided by GLOBE Claritas fortheir processing software, which was vital during the final stages.I have worked with an incredible (and incredibly diverse) group of people in EOS/EOAS atUBC, and I am sure to forget to thank someone. Some of the people with whom I?ve workedclosely are: Andrew Schaeffer, Ben Postlethwaite, Brett Gilley, Jounada Oueity, Mika McKinnon,Dave Semeniuk, Rebecca Taylor, Anna Hippmann, SashaWilson, Phil Hammer, Evan Smith, RalfHansen and Mark Jellinek. All of the people in the SLIM and GIF groups, including but not limitedto: Felix Herrmann, Ian Hanlon and Doug Oldenburg; Rajiv Kumar, Haneet Wason, Li Xiang andxiiiTristan van Leeuwen; Elliot Holtham, Dikun Yang, Dave Marchant and Rowan Cockett. I?d alsolike to thank our outstanding support staff; everyone here is startlingly competent, but I?ve hadthe pleasure to bug, pester or chat with some more than others, namely: John Amor, Cornel Pop,Michael MacIsaac, Nick Irvine and Sukhi Hundal; Teresa Woodley, Cary Thomson and CeciliaLi.Finally, on a personal note, I would like to thank my family for supporting me through thiswork and generally for ?putting up with my nonsense?. My wife Megan is incredible, and shedeserves at least half the credit for this thing, because I wouldn?t have managed any of it withouther support. My parents Debbie and Robert Smithyman and my grandparents Colin, Sylvia, Earland Helen; they deserve the other half of the credit, for getting me where I am today. My sisterKelsey: thanks for keeping me on my toes! My parents in-law, Marilyn and Nigel Blamire, andmy sister-in-law Emily (who also keeps me on my toes); I think they?re all awesome, and for someinexplicable reason the feeling is mutual. Eric: thanks for all the chats, because seriously, it?s awonder that anyone else puts up with either of us.xivTo Megan, without whom none of this would be possible.xvChapter 1IntroductionRefraction processing of surface vibroseis multichannel seismic (MCS) reflection data to de-termine near-surface velocity variations is typically limited to near-offset refraction statics usingtraveltime-based procedures. Waveform tomography combines inversion of first-arrival travel-time data with full-waveform inversion of densely sampled refracted waveforms (Pratt and Wor-thington, 1990; Pratt, 1999). Since inversion of the waveform amplitude and phase is not lim-ited by the ray-theory approximation (as is traveltime tomography), identification of low-velocityzones and small scattering targets is possible. In order to proceed with the full-waveform inversionof refraction data, it is important to have a suitably accurate starting velocity model. This is typi-cally a synthesis of velocity models from refraction statics processing and/or traveltime inversion(i.e., ray tracing) along with a priori geological information (Pratt and Goulty, 1991).This work involves the development and application of 2.5-dimensional (2.5D) waveform to-mography to refraction phases from vibroseis MCS reflection data collected in the Nechako Basinof south-central British Columbia. This represents a challenging application of 2.5D waveform to-mography because of the sub-basalt setting of the basin, the non-ideal geometry associated withland acquisition, and the absence of critical low-frequency information in the data. Detailedcharacterization of the top- and bottom-basalt interfaces can potentially improve significantly theresolution of deep reflectors. In addition, a proof of concept application of 2.5D waveform to-mography for conventional vibroseis field data may encourage adoption of this methodology infuture projects.11.1 Research objectivesThe research objectives of this thesis may be grouped into two main categories:1. The technical development of new processingmethods to enable waveform tomography inversionof on-land crustal exploration seismic reflection data.2. The development and interpretation of geophysical imaging results in the south-central intermon-tane belt of the western Canadian Cordillera using waveform tomography.In part, the latter objective drives the former, as research developments in the technical aspectsare motivated and tested by the challenges posed in producing the geophysical results. This thesisis structured around a series of case studies drawing from recent exploration seismic work in theNechako?Chilcotin plateau. Data were collected in 2008 (Calvert et al., 2009) that are processedin this thesis to produce geophysical models. These geophysical models and their interpreta-tions represent a contribution to the understanding of geological structures along several profiles,with benefits to both scientific understanding and potentially industrial development. Waveformtomography (including full-waveform inversion) is a complex geophysical method, and many ap-proximations are typically made in order to improve the tractability, feasibility and affordabilityof the approach. Technical developments in this thesis permit the application of 2.5D waveformtomography to datasets acquired in challenging on-land situations, and substantially improve theaffordability of the approach for academic researchers and practitioners in independent industry(viz., small-scale seismic contractors).1.2 Traveltime inversionTraveltime inversion is a class of methods for seismic velocity model inversion in which the first (orpossibly later) refracted and/or reflected wave arrivals from a known source are timed. The mea-sured arrival times are compared with synthetic arrival times from a presupposed starting model,and that model is updated to reduce their misfit. The misfit is typically measured using a weightedleast-squares criterion, although the use of alternative misfit criteria is gaining in popularity. The2synthetic arrival times are found using raytracing (either direct shooting or a wavefront marchingalgorithm), and the kernel functions are based on the backprojection or backpropagation of thetiming errors to the model. In first-arrival traveltime inversion, the least-time path between eachsource and receiver is considered.Layer velocity methods are well suited to seismic methods using near-vertical-incidence wave-forms because they represent the earth as a series of horizontal to sub-horizontal layers. Lateralvariations are accounted for by varying the depth and thickness of each layer along its length.The resulting model provides averaged interval velocities for each layer. Traveltime tomographymethods, in comparision, typically work on a gridded model and allow variation in vertical andlateral directions. In addition to being a viable standalone analysis technique, traveltime tomog-raphy methods are a useful means of generating a starting model for full waveform inversion (Prattand Goulty, 1991). These methods are based on ray theory, in which a particular wavepath isapproximated by a raypath. The ray integral,T =?ray1c(x)ds (1.1)expresses the traveltime T in terms of an integral along the raypath, through a time-independentvelocity field c(x). This results from the eikonal equation,|?? |2 = 1c2(x)drds(1.2)wherein ? is the wavefront time delay, and drds is a unit vector pointing along the ray. Onlyone arrival phase is synthesized per source-receiver pair, and in this simple implementation am-plitude is not considered. Inversion to generate an updated model may then be applied throughbackprojection of data residuals (e.g., Hole, 1992; Hole and Zelt, 1995) calculated by numericalforward modelling (Vidale, 1990).31.2.1 FASTFirst Arrival Seismic Tomography (FAST) by Colin Zelt is an iterative solver that computes anupdated velocity model based on a gridded input model (Zelt and Barton, 1998). This worksunder the approximation of the eikonal equation, and consequently is not sensitive to low-velocityanomalies under the size of the Fresnel zone.FAST provides means for controlling the model regularization in both x- and z-directions,when used in two dimensions. FAST additionally supports 3D operation and model updates withsmoothness and smallness constraints. As a result, this also permits operating on the crookedaspects of a 2D land survey. An optimally simple model can be found through the use of Tikhonovregularization, in which the weighting of the model objective function is gradually reduced untila target data misfit is achieved. The exit criteria of FAST is based on a presupposed noise levelfor each traveltime pick; however, it would be relatively straightforward to incorporate morecomplex logic due to the manner in which the results are stored at each iteration.The resolution of the traveltime tomography method is on the order of the diameter of the firstFresnel zone,??Lwherein ? is the relevant temporal frequency and L is the source-receiver pathlength (Williamson, 1991;Williamson and Worthington, 1993). Because of the effects of multiplescattering, anomalies smaller than this scale length will experience wavefront healing, i.e., thefirst-arrival time will not contain information about the anomaly. The effect of such anomalies onthe wavefield may be distinguishable to methods that consider waveform amplitude and phase,such as waveform inversion.1.3 Full-waveform inversionThe family of techniques known as full waveform inversion has grown out of the fields of seis-mology and non-linear inversion since the mid 1980s. In contrast with traveltime tomography orrefraction statics, waveform tomography is more computationally expensive, but provides bene-fits in terms of its ability to see beyond the first arrival and consider data phases that are invisibleto a traveltime inversion. In this document, I will use the term waveform tomography to refer4to a collection of similar methods for full waveform inversion of transmitted wave energy. Inparticular, the procedures and implementation I use refer to the Pratt Waveform Tomographymethod and its implementation as the code ?FULLWV? (Pratt, 1990, 1999). Where possible, Iwill make distinctions between the general field of full waveform inversion and this particularimplementation.Initial work by Lailly (1983) and Tarantola (1984a) showed that the inverse problem for amodel could be solved as iterative back propagation of data errors. Rather than an expensiveglobal search, the inverse solution could be formulated in terms of the correlation of forward-propagated synthetic waveforms and backward-propagated data residuals. This provides modelresiduals in a form compatible with a local-descent scheme.Because of this correlation mechanism, the resolution of the waveform tomography methodis controlled by convergence to a particular data phase. The starting model must be sufficientlyaccurate to reproduce the measured (field) data within one cycle, or a ?cycle-skip? may occur,which results in convergence towards a local minimum of the objective function. Consequently,the resolution of the waveform tomography output is on the order of ?2 , wherein ? is the wave-length at a particular frequency of operation.Waveform inversion differs from traveltime inversion because the former models and invertssome approximation to the wave equation. At minimum, this means that a waveform forward-modelling method will synthesize multiple data phases and their corresponding amplitudes. Thecorresponding inversion procedure will incorporate this information and in theory produce a su-perior model to traveltime inversion. The term full-waveform inversion usually refers to waveforminversion in which all relevant data phases are considered, though the usage is not completelyconsistent among researchers.The computational domain and level of simplification varies among full-waveform inver-sion implementations. Tarantola (1986) and Mora (1987) discuss full-waveform inversion ofacoustic data by time domain finite differences, whereas Pratt and Worthington (1990) and Liaoand McMechan (1996) use frequency-domain finite-difference implementations. Vigh and Starr(2008) discuss some of the benefits of each.5The chief benefits of a frequency-domain approach are: 1) improvedmodelling (and inversion)efficiency in 2D through direct factorization of the forward modelling terms; and 2) access tomulti-scale inversion benefits by inverting data from low temporal frequencies prior to data fromhigh temporal frequencies. A time-domain approach benefits from simpler implementation andpartitioning of the model for computational parallelism, and also offers the ability to time-windowdata without risk of aliasing. The memory requirements of direct matrix factorization methods in3D are in most cases prohibitively large, and so 3D implementations are more straightforward toimplement in the time domain (though notable counter-examples exist, e.g., Ben-Hadj-Ali et al.,2008).The work by Tarantola (1984a) laid the groundwork for applied waveform tomography. Heshowed that the iterative back propagation of data residuals could be used in a descent scheme,instead of an expensive global search through the space of the objective function. The operationis similar to reverse time migration, in which an image is formed by forward propagation of syn-thetic waveforms and back propagation of receiver waveforms onto a model grid. In reverse timemigration an image is formed directly from the data by either excitation-time or cross-correlationimaging conditions. Xu et al. (2010) provide an excellent recent review of the technique.In waveform tomography, instead of back propagating recorded data, the data residuals arebackpropagated and cross-correlated with the synthetic wavefield. This results in reverse-timemigration of the data errors into model space. Rather than forming an image directly, an image ofthe error in the model is produced. The error image acts as the gradient in a local descent scheme(preconditioned non-linear conjugate gradient), which iteratively updates the model in pursuit ofreduced data (and possibly model) misfit. Some methods precondition (scale) the gradient withadditional terms (e.g., an estimate of the Hessian, which results in a quasi-Newton method).1.3.1 Wave equationThe propagation of waves in the real earth is determined by the properties of the materials throughwhich the waves pass. In seismology, it is typical to assume that the physics of seismic waves arewell-modelled by the elastic wave equation, which is a partial differential equation in space and6time. This assumption is valid for waves that have an amplitude sufficiently small to make theresponse of the earth time-invariant, which is valid in practice for the seismic case outside of near-field earthquake-generated events. The displacement vector u is written in terms of displacementpotentials:u = ?? +??? (1.3)wherein ? is the scalar displacement potential and ? is the rotational (vector) displacement po-tential. For dilatational P-waves, the elastic wave equation may be written as in Equation 1.4:1?2?2?t2? = ?2? (1.4)For rotational S-waves, the wave equation is expressed in Equation 1.5:1?2?2?t2? = ?2? (1.5)The terms ? and ? are the P-wave and S-wave propagation velocities and:? =?K+ 43?? (1.6)? =??? (1.7)The P-wave velocity is dependent on the bulk modulus (K; incompressibility), shear modulus (?;rigidity) and the density (?) of the earth materials; the S-wave velocity is dependent on the shearmodulus (?) and the density (?). Hence, P-waves travel faster than S-waves by definition. TheP-wave equation may be written in terms of pressure under the assumption that the medium doesnot possess a shear strength (stiffness), which results in the acoustic wave equation:1c2?2?t2p = ?2p (1.8)7This is analogous to Equation 1.4, except that ? has been replaced by the acoustic velocity (c)and the equation is explicitly in terms of pressure (p). Equation 1.8 is rearranged and reposed asan inhomogeneous version with the addition of a source term, s:[?2 ? 1c2?2?t2]p = ?s (1.9)The addition of viscous damping leads to a modified viscoacoustic wave equation,[?2 ? ? ??t? 1c2?2?t2]p = ?s (1.10)in which the term ? relates to loss of energy due to intrinsic attenuation (viz., by Ohmic heating).1.3.2 Forward modellingPratt et al. (1998) describe the constraining partial differential equation for numerical forwardmodelling in matrix parlence, and I repeat part of their discussion here. Their Equation 1 is:M?2u?(t)?t2+ Ku?(t) = f?(t) (1.11)The termM is a so-calledmass matrix and K is the stiffness matrix. The time-domain vectors u? andf? are the solution wavefield and the source term (respectively) for the inhomogenous discretizedwave equation. Their Equation 2:M?2u?(t)?t2+ C?u?(t)?t+ Ku?(t) = f?(t) (1.12)describes the case in which viscous damping is added to the problem, and C is called the dampingmatrix. However, Pratt et al. (1998) take the temporal Fourier transform to lead to a modifiedversion of the problem,Ku(?) + i?Cu(?)? ?2Mu(?) = f(?) (1.13)8which is the discretized inhomogeneous Helmholtz equation with viscous damping andu(?) =? ???u?(t)e?i?tdt and f(?) =? ???f?(t)e?i?tdt (1.14)This problem is posed in the frequency domain, and the boldface vectors u and f defined inEquation 1.14 are the Fourier-transformed equivalents of the time-domain vectors u? and f?. Fornotational convenience, Pratt et al. (1998) rewrite Equation 1.13 as:Su = f or u = S?1f (1.15)wherein S is termed the impedance matrix. Pratt et al. (1998) point out that the second formin Equation 1.15 is only representational, since the method they employ (which I also use) isan implicit finite-difference numerical approach. In FULLWV, this solution is implemented withan implicit 9-point finite-difference star (Pratt and Worthington, 1990; ?tekl, 1997), and absorb-ing boundary layers are implemented using the perfectly matched layer approach Roecker et al.(2010).1.3.3 Velocity model inversionSection 1.3.2 shows how forward modelling is carried out for a source signature s, possibly com-posed of a distribution of sources, over both space and time. The forward-modelling code deter-mines the response of the medium at a particular location to a source perturbation in the wavefieldat another location, dependent on the pre-existing model. Of course, if the pre-existing model iscorrect, the problem has already been solved. In lieu of this happy situation, the goal becomes toimprove the model based on recorded data. To measure the degree of correctness in our model,we define an objective function E(m) (Equation 1.16) that consists of the L2 norm of the datamisfit calculated in a model m.E(m) = 12?dT ?d? (1.16)9The quantity ?d represents the data misfit column vector, the superscript ? indicates complexconjugation, and the superscript T the matrix transpose (Pratt, 1999). If the model perfectly re-constructs the earth response from the true model, we would expect the misfit function to beminimized; in fact, if the data were noise free, the error would tend towards zero.We represent the current model as m, and the starting model as m0. Each of these is rep-resented by a vector of complex velocities, dimensioned into a 2D array. Since we presumethe starting model is incorrect, but broadly similar to the true model, we would ideally like tofind a vector ?m such that m = m0 + ?m matches the true earth model mt. In a linear problemthis would allow the inversion procedure to converge successfully in only one iteration from anystarting model. However, the nonlinearity of the wave equation with respect to velocity meansthat the objective function in Equation 1.16 does not decrease monotonically towardsmt. Conse-quently, a local gradient operation is not guaranteed to find mt, and so the starting model is veryimportant.Using an estimated starting model (velocity and possibly attenuation in terms of 1Q ) and aknown source distribution, the forward modelling code provides a pressure field dr(r, t) at modellocations of interest. For each model location that coincides with our data (i.e. the location of areceiver at the time we were recording), we can make a comparison between dr(r, t) and dobs(r, t),the recorded data at r. This comparison gives us ?d = dr ? dobs, the data residual vector.To generate a model update that reduces the data misfit, a gradient direction is found by takingthe derivative of E(m) with respect to each model parameter:?mE =?E?m = R{J?T ?d?}(1.17)In this case, J?T is the transpose of the Jacobian J?, which contains the Fr?chet derivatives. Thismatrix is expensive to form explicitly, but can be computed efficiently using the adjoint-statemethod (Pratt, 1999; Plessix, 2006):?mE = R{J?T ?d?}= R{F?T[G??1]T?d?}(1.18)10wherein G??1 is the inverse of the seismic Green?s function G?, and F? is composed of so-calledvirtual source terms:f(i) = ??G??mid (1.19)In this way, the gradient is formed from the forward propagated synthetic wavefield and thebackward propagated data residuals (Pratt, 1999). The backward propagated wavefield could bethought of as a stencil on the portion of the model that the data could be sensitive to, given theform of G??1 (which is implicitly dependent on m for the seismic problem). F? acts to illuminatethe portions of the model that are inside the kernel of the forward modeling operator, based onthe estimated source signatures. By convolving these two fields, the gradient ?mE is formed(Equation 1.16) over the extent of the model space.The gradient is representative of the inconsistencies between the starting model and the familyof models supported by dobs. For some sufficiently small step-length ?, the model mk can beupdated asmk+1 = mk ? ??E (1.20)with a guarantee that, by the measurement norm used to define the error, the misfit willdecrease. The step-length ? may be computed by a line search.The waveform forward modelling solution can be implemented in either the time domainor the frequency domain. Each method benefits from certain characteristics that the other lacks.Since the inverse method relies on forward modelling operations to carry out its updates of thecomplex velocity model, the inversion automatically operates in the same domain.The time-domain implementation of acoustic forward modeling benefits from:? The ability to window input data without risk of aliasing.? Simplified parallelization of the solution for multiple processors.The frequency-domain implementation benefits from:? The option of non-linear operation in terms of the chosen frequencies.11? The ability to easily account for dispersion and attenuation (Pratt, 1999).? The convolution operation discussed in Section 1.3 is implemented as a simple element-wisearray multiplication.The Pratt Waveform Tomography method described in this document is a frequency-domainfinite-difference acoustic solver, operating in two dimensions. The forward propagator G? is de-composed by LU -decomposition (Song and Williamson, 1995; Forsythe and Moler, 1967) intolower- and upper-triangular matrices, which makes finding its inverse G??1 a relatively trivial oper-ation. In addition, the LU matrix factors can be stored for multiple applications of the propagator.A minimum of three forward modelling operations must be carried out:? The forward modelling to find the virtual source terms in F?? The forward modelling (backwards in time) to find the backpropagated wavefields[G??1]T?d?? At least one forward modelling operation in calculation of the step length1.3.4 Logarithmic L2 misfitIn Chapter 4 I discuss the use of a variant misfit criterion, the logarithmic L2 norm. This is theresult of work by Shin and Min (2006). As is discussed in Chapter 4, the use of the logarithmic L2norm instead of the conventional L2 norm may lead to improved numerical stability, particularlyin the presence of data outliers.S(?)u?(?) = f(z)gs(?)ei?s(?) (1.21)Shin and Min (2006) begin by defining a one-dimensional frequency-domain problem (our Equa-tion 1.21), which is roughly analogous to our Equation 1.15 (originally from Pratt et al., 1998;their Equation 5); however, note the notational difference wherein u? represents the frequency-domain wavefield and f(z) is representative only of the positional component of the source term.Also, note that for notational brevity I have removed the explicit dependence on ? from the terms,12but that all of these terms are implicitly frequency-dependent. Shin and Min (2006) defineg = gsei?s (1.22)as the temporal source signature with amplitude gs and phase ?s. The terms in the wavefield u?are expressed asu?j = gsAmj ei?mj +i?s . (1.23)The scaling of these terms is controlled by the amplitude of the source signature gs and the am-plitude of the Green?s function (the earth response) Amj at the jth depth location. The phase ofthe wavefield is a sum of the phase delays from the source term ?s and the Green?s function ?mj .The modelled wavefield at the surface isu?1 = gsAm1 ei?m1 +i?s (1.24)and the measured seismogram isd?1 = Af1ei?f1 (1.25)Shin and Min (2006) note that the phases are wrapped, but make the assumption that no cycleskip exists in the modelling from the initial model estimate. The misfit criterion is then:ln u?1d?1= ln gsAm1Af1+ i[?m1 + ?s ? ?f1](1.26)i.e., the logarithms of the synthetic wavefields are compared with the measured wavefields, result-ing in a separable misfit criterion over amplitude and phase. The objective function for a singlefrequency as defined by Shin and Min (2006) is:E = 12???[ln gsAm1Af1]2+[?m1 + ?s ? ?f1]2???(1.27)13and the gradient ?k with respect to the kth model parameter pk is:?k = ?pkE =?E?pk= ln(gsAm1Af1)(1Am1)?Am1?pk+[?m1 + ?s ? ?f1] ?m1?pk(1.28)In implementation, the gradient for this method is computed similarly to the case for theL2 norm. The logarithmic scaling on the gradient amplitudes acts as a preconditioner, whichcompresses the dynamic range of the gradient and results in greater weight being given to low-amplitude components of the misfit. This may be expected to be an improvement on the L2 normin cases when the error distribution of the data are well-modelled by the heavy-tailed probabilitydistribution of the logarithmic L2 norm, in comparison with the Gaussian distribution of the L2norm.1.3.5 Iterative Tikhonov regularizationIn some of the processing for Chapter 3 I apply a method known as iterative Tikhonov regular-ization to regularize the model updates during full-waveform inversion. In a linear inversionframework, the goal is to find the optimal model x that satisfies the system of equations Ax = bwith the forward modelling operator A for a data vector b. The least-squares misfit is defined as|Ax? b|2 . (1.29)Tikhonov regularization refers to a method for constrained optimization, in which the penaltyfunction is augmented by a misfit term that acts on the model vector itself as in|Ax? b|2 + ? |?(x? x0)|2 (1.30)The matrix ? may simply be an identity matrix, or some more complex weighting function, andthe term ? acts as a scale factor to control the relative weighting of the model constraint and thedata constraint in the misfit. For the special case when x0 (the reference model) is not zero and? ? I, this criteria will tend to penalize models different from the reference model.14Iterative Tikhonov regularization involves the addition of a vector (such as a scaled versionof x0, sign-reversed) to the gradient prior to computation of the step length during a linearizedinversion. This will tend to have a similar effect to Tikhonov regularization applied to the objectivefunction (for correct choice of the scaling factor), but is more straightforward to include in anexisting local-descent inversion framework.1.3.6 SummaryIn this project, I make use of full-waveform inversion in the frequency-domain under the 2D, vis-coacoustic approximation (Pratt andWorthington, 1990). This technique draws on previous workcarried out by myself and others (Smithyman et al., 2009; Brenders and Pratt, 2007b,a). The geom-etry issues posed by land seismic acquisition make processing the data from the Nechako Basincase study intractable with an unmodified 2D approach: Chapter 3 discusses static-correction ap-proaches that can improve the tractability of the 2D full-waveform inversion problem. Chapter 4extends the work of Song and Williamson (1995); Song et al. (1995); Zhou and Greenhalgh (2006)and others to formulate an improved 2.5D approach to solving the problem for full-waveform in-version of land seismic data.1.4 Vibroseis acquisitionVibroseis is a method for generating a seismic source signature at the ground surface without theuse of explosives (Crawford et al., 1960). Large vibrating trucks are triggered simultaneously bya continuous control signal with controllable frequency content. After acquisition, the controlsignal is correlated with the recorded data from each receiver location to produce seismic tracessimilar in appearance to data shot with a delta pulse.Active source seismology commonly makes use of several methods of generating a controlledshot. Explosives are commonly used to produce an impulsive source similar to a delta functionin time. An explosive source is close to minimum-phase (because of the minimum delay releaseof energy implicit in the reaction).15Because the cross correlation is band limited by the range of frequencies in the (commonlylinear) vibroseis sweep, the wavelet does not perfectly mimic a delta function: Both causal andacausal trailing components exist in the source term. However, the minimum phase earth reflec-tivity series also plays a role. The recorded seismogram x(t) is:x(t) = s(t) ? w(t) ? e(t), (1.31)wherein s(t) is the vibroseis sweep, w(t) is the wavelet containing system response charac-teristics, and e(t) is the earth?s reflectivity series (Yilmaz and Doherty, 1987). In correlation, thevibroseis sweep s(t) is replaced by the Klauder wavelet, formed from the autocorrelation of thevibroseis sweep itself. However, correlation with the earth reflectivity series and system response(assumed to be minimum phase) produces a mixed phase source signature. The mixed-phasenature of the source makes first-arrival picking difficult, since the zero-time phase of the sourceis not coincident with the first peak or trough.1.5 Difficulties posed by on-land seismic acquisitionOn-land 2D seismic acquisition has some disadvantages in comparison to offshore seismic acqui-sition in several key ways:? the signal-to-noise ratio in the data is typically lower? topographic effects create shadow zones and complex near-surface reflections that are difficultto model? the existence of ground roll makes it difficult to use data from near-offset traces? land use concerns may dictate the survey geometry, in particular leading to crooked-line acqui-sitionThe latter issue (crooked-line geometry) is of particular interest in this study, and is the motivationfor technical developments that follow. On-land acquisition also has its advantages, e.g., acqui-sition of long-offset and wide-azimuth data may be more effective because the available seismicaperture is wider without the presence of a water layer.16yxz(3-D) True Path(2-D) Projected PathAcquisition Line(xs, ys, zs)(xs, zs)(xr, yr, zr)(xr, zr)Figure 1.1: The idealized seismic cross-section for this 2D survey is presented in the (x, z) plane. The2D seismic line is depicted in red, and diverges either side of the idealized profile (constant y). Source-receiver paths are indicated in orange (3D) and green (2D), which represents a gross simplification ofthe problem for purposes of discussion (see text). In the 2D version, the y-coordinate is identical forsource and receiver (simple projection), and therefore the path length differs.Figure 1.1 provides a simplified overview of crooked-line 2D seismic acquisition. The orangeray represents the original path taken by the seismic energy from source to receiver. When thesedata are processed using a 3D method, the source and receiver positions may be represented inthree dimensions (x, y, z), and the source-receiver offset will be modelled approximately correctly.However, in 2D processing only two coordinates exist (x, z), and some projection or mapping17must be made to convert between the acquisition geometry and the simplified 2D geometry. Ifthe y-coordinate is simply ignored, the 2D source-receiver path length may be shorter than the3D path length.Figure 1.2 (cf. Figure 1.1) shows two methods for dealing with the geometry reprojection.Figure 1.2a depicts the result of formulating the model geometry in terms of a ?slalom line?, whichis a smoothed version of the line geometry produced by projection to CDP (Common Depth Point)bin centres. This configuration is typical in seismic reflection processing, because source andreceiver positions are preserved (though typically static-corrected to some datum). The model isdefined by the CDP locations directly, and a seismic image (stack, migration image, etc.) is drapedfrom the top of each CDP bin centre. Figure 1.2b presents an alternative method for ensuring path-length consistency, in which the source and receiver locations in the 2D profile are shifted alongthe line either side of the midpoint in order to simulate the longer raypath encountered in 3D.The formulation of waveform tomography discussed in this thesis is carried out in a 2D modelgeometry, and therefore some sort of geometric correction is necessary. Chapter 3 demonstratesone solution that uses static corrections from 3D raytracing to time-shift the field seismic data.Chapter 4 discusses a superior method that makes use of 2.5D full-waveform inversion to handlethe 3D geometry directly, with some necessary assumptions about model characteristics.1.6 The Nechako Basin datasetThe 2008 Geoscience BC Nechako Basin seismic survey comprises a total of 13 acquisition lines,the longest of which is in excess of 80 km (Calvert et al., 2009). Figure 1.3 is reproduced fromCalvert et al. (2011) and depicts the locations of both the new Geoscience BC lines and theolder Canadian Hunter seismic lines (not discussed here). Waveform tomography processing isapplied to the seismic data from Lines 05, 10 and 15 of the Geoscience BC dataset. Line 06 isalso processed using traveltime tomography, but full-waveform inversion is not applied to Line06 due to time and budget constraints.Initial development work (i.e., Chapters 3 and 4) is motivated by the challenges presented by18Model geometry reprojected to slalom line (CDP bin centres)Path-length corrected with constant midpointPathlength APathlength A(xs, ys, zs)(xs ? ?x, zs)(xr, yr, zr)(xr + ?x, zr)(xr, zr)(xs, zs)a)b)Figure 1.2: Two possible solutions are shown for reprojecting 3D source/receiver geometries in 2Dmodel space. The slalom (CDP) line approach (a) is commonly used in reflection seismology. An-other possibility is the reprojection of sources and receivers to adjusted locations, maintaining source-receiver offset and midpoint (b).Line 10. An excerpt from this line was chosen based on the relative straightness of the profile (seeChapter 3), with a goal of applying 2D waveform tomography processing. However, crooked-line seismic acquisition makes the application of 2D traveltime tomography and full-waveforminversion intractable.19Figure 1.3: An overview of the Nechako Basin seismic acquisition profiles in relation to surficialgeology and points of interest. Reproduced from Calvert et al. (2011), their Figure 2.20Two methods are developed with the goal of improving the tractability of full-waveform in-version for Line 10 (Chapters 3 and 4). Chapter 5 details the application of 2.5D waveformtomography to Lines 05, 10 and 15 of the Geoscience BC Nechako Basin seismic dataset, andtraveltime tomography to Line 06.1.6.1 Geometry and acquisition parametersCalvert et al. (2009) details acquisition parameters for the 2008 Geoscience BC Nechako Basinseismic survey. The acquisition was carried out with the primary goal of maximizing signal tonoise ratio, since this was a main detractor from previous seismic data acquired by CanadianHunter in the 1980s. Vibroseis sources were used, with a nominal source spacing of 40 m anda receiver group spacing of 20 m. Part of Line 06 was acquired with 20 m source spacing. Fourvibroseis trucks were triggered simultaneously a total of four times per vibration point, with aresulting diversity-stacked composite source approximately 60 m wide along-line. The sweepduration was 28 s, with bandwidth from 8?64 Hz and 0.9 s on and off tapers. Nominal fold was240.1.7 Geological context and backgroundThe region known as the Nechako Basin is loosely defined as the south-central intermontane beltof the western Canadian Cordillera. This term has been used inconsistently by workers in the field,leading to ambiguous definitions of the region?s extent. Riddell (2011) discusses this in detail. Forpurposes of present and future publications as of 2013, it has been suggested that researchers usethe term Nechako?Chilcotin plateau and/or specifically identify the geological and geographicalfeatures (personal communications; A. Calvert). Chapters 3 and 4 use the term Nechako Basin inorder to maintain the text of Smithyman and Clowes (2012) and Smithyman and Clowes (2013)intact, whereas Chapter 5 makes use of the term Nechako?Chilcotin plateau.The Nechako?Chilcotin plateau has been an occasional hydrocarbon exploration target since1931, with several wells drilled throughout the region. Recent efforts (since about 2006) have21been motivated by government-industry partnerships involving Natural Resources Canada, Geo-science BC, the British Columbia Ministry of Energy, Mines and Petroleum Resources and severalCanadian universities (Riddell, 2011). Massey et al. (2005) compiles unified geological informa-tion for the province of British Columbia, and acts as a basemap for most geophysical resultsin this area since its publication. Detailed mapping work and rock physics studies (e.g., Bordetet al., 2011; Kushnir et al., 2012) have greatly improved the understanding of geological featuresin the region. Geophysical surveys (e.g., Spratt and Craven, 2010; Calvert et al., 2011) have pro-vided valuable information about the structures in the Nechako?Chilcotin plateau region, whichmay be interpreted alongside rock physics information to tie together information from geologicalmapping.The south-central intermontane belt of the Canadian Cordillera is mainly underlain by Jurassic-aged rocks of the Stikine terrane, though the eastern part of the region is underlain by Jurassicand Triassic rocks of the Cache Creek terrane. Jurassic basement rocks are mainly arc volcanicand volcaniclastic rocks in the Stikine terrane, with more sedimentary rocks found in the base-ment to the east. A series of Cretaceous volcanic, volcaniclastic and clastic sedimentary rocksoverlie the basement. The conglomerates and sandstones in the Cretaceous succession are mostprospective for resource development, and are therefore chief exploration targets for industry ef-forts. The surficial geology is dominated by volcanic rocks that were erupted contemporaneouslywith Eocene extension. These dominantly felsic volcanic rocks are overlain by a thin cover ofNeogene basaltic rocks (of the Chilcotin group). The volcanic cover contains large P-wave veloc-ity contrasts, and has historically been a challenge to geophysical exploration. Characterizationof these volcanic rocks by waveform tomography is a primary goal of this thesis.22Chapter 2Technical DevelopmentsDue to the presentation of following material (Chapters 3 through 5) as self-contained researchpapers, several key elements are discussed only briefly. The present chapter provides some back-ground material that is not included in the individual research chapters to follow.2.1 Handling of coordinatesThe 2.5D waveform tomography method (Chapters 4 and 5) generates 3D wavefields using aFourier synthesis of multiple 2D wavefields. This yields benefits in terms of the simplicity of theproblem and reduces the computational scale in comparison to 3D processing. However, theinversion method is only valid under circumstances in which the simplification of the 3D waveequation is valid (see Bleistein, 1986; Song and Williamson, 1995). In particular, this means thatthe earth model is invariant in the direction of the Fourier synthesis operation, i.e., perpendicularto the grid of the 2D sub-problems. There is no ability to model cross-line effects or scatteringfrom reflectors out of the plane of the model; therefore, wave arrivals in field data from out-of-plane scattering will be (incorrectly) mapped to in-plane features. In practice, the validity of thisassumption will vary for the real 3D earth, depending on the local-scale variability of the site ofinvestigation. We discuss three cases, in which the elevation (z-axis position) of each source orreceiver remains unchanged.For earth materials that:1. are approximately 1D, it is sensible to align the forward-modelling and inversion grid to bestapproximate the true acquisition geometry (e.g., a best-fit line based on source and receiver po-23sitions). This is in lieu of a preferred orientation arising from geological structure considerations.Choosing to fit the line to the array geometry minimizes the number of plane-wave componentsnecessary to adequately simulate the 3D wavefield.2. possess 2D structure, it is most desirable to align the forward-modelling and inversion grid tocross geological strike. This maintains the validity of the 2.5D projection, which makes 2.5Dfull-waveform inversion attractive as an alternative to 3D methods.3. possess complicated 3D structures, it is advisable to use 3D methods. However, dependingon the length-scales of structures in different directions, it may be possible to apply 2D or 2.5Dprocessing with some success. In this case it is still desirable to align the grid to cross the dominantgeological strike as in the 2D case.Four cases are presented in Figure 2.1 that are more or less well-suited to analysis using 2.5Dwaveform tomography. In the case of 1D earth structures (Figure 2.1 a), 2.5D FWI is tractable,because the model is homogeneous in the direction perpendicular to the acquisition line. How-ever, in the presence of 2D structures (e.g., dipping strata), the orientation of the line matters. Ifthe acquisition line is along strike (Figure 2.1 b), the 2.5D approximation is ill-posed and the 3Dwavefield may not be synthesized correctly. If the acquisition line crosses geological strike (as inFigure 2.1 c), the 2.5D approximation is well-posed because the model is homogeneous in thecross-line direction. In practice, complex 3D structures will exist at some length scale in the realearth (e.g., Figure 2.1 d). If the length scale of out-of-plane geometry is sufficiently large then2.5D full-waveform inversion will be tractable. If the length scale of out-of-plane geometry issmall relative to the seismic offest then 3D full-waveform inversion is expected to yield a superiorresult.In order to accommodate geometries that are non-ideal, a simple coordinate rotation andtranslation can be applied to convert between acquisition coordinates (x, y, z) and modellingcoordinates (x?, y?, z?). This enables the methodology to adapt to acquisition geometries that donot cross geological strike perfectly, as long as the alignment error is moderate. This may beapplied in the event that a priori knowledge of the dominant geological strike is available.24??????x?y?z???????=??????cos ? sin ? 0? sin ? cos ? 00 0 1????????????x? x0y ? y0z ? z0??????(2.1)The coordinates from the original acquisition (e.g., x = UTM Easting and y = UTM Northing,z = Elevation above Mean Sea Level) are transformed to new modelling coordinates by firstsubtracting the location in the original coordinates (x0, y0, z0) that will serve as the origin of themodelling grid. A coordinate rotation about the z-axis is then applied (the angle controlled by ?)to align the x? coordinate with the horizontal coordinate of the 2D modelling grid.1D Structure2D Structure ? Line Crosses Strike 3D Structure2D Structure ? Line Along StrikeFigure 2.1: Four simplified cases are presented, showing the line of seismic acquisition relative tolarge-scale geological features. In these cartoons the seismic array is represented by a dashed line,and the coloured layers represent different geological units. These cases are discussed further in thetext.252.2 Traveltime tomography procedureTraveltime tomography may be applied as 2D or 3D, depending on the characteristics of theproblem geometry. Unlike with full-waveform inversion, geometric spreading and amplitudeeffects are not relevant to the Eikonal equation (Chapter 1.2). In choosing to apply 2D or 3Dtraveltime inversion one must consider several factors:1. Processing traveltimes in 3D incurs an additional computational cost that may be prohibitive forlarge problems.2. Processing traveltimes in 2D may simplify the geometry of the problem too much, resulting intraveltime modelling errors.3. Processing traveltimes in 3D permits out-of-plane model structures.Factors 1 and 2 are trade-off parameters, but in practice modern computers make it reasonableto process traveltimes in 3D if any substantial errors would result from 2D processing. Factor 3must be considered as a benefit or hinderance depending on the aim of the inversion project. Forsurvey geometries that provide sufficient coverage to interrogate 3D features, this last factor isof paramount importance. However, when traveltime inversion is used to develop a model thatwill feed into 2.5D (or, in principle, 2D) FWI, the situation is more complex. Because 2.5D FWImakes use of a 2D model of earth physical properties (viz., velocity and possibly attenuation), 3Dheterogeneity cannot be parameterised. The 2.5D methodology is predicated on the assumptionof homogeneity along the cross-line direction. It is necessary to transform the 3D model to somerepresentative 2D equivalent.One approach to the dimensionality reduction from 3D to 2D is to impose smoothing in thethird dimension during the inversion process, as discussed in Chapters 3 and 4. This has the virtueof developing features in an iterative manner that fit the 3D traveltime data using a nominally 2Dmodel. We term this methodology 2.5D traveltime inversion. Alternatively, traveltime inversionmay be carried out in 3D, with computation of 3Dmodel updates. A representative 2Dmodel maythen be extracted using a spatially-weighted average of the 3D model parameters. The naturalapproach is to derive this weighting from the raytracing hit counts in the 3D model. This is then26roughly equivalent to extracting the values in CDP bins. However, because this modifies the 3Dmodel, it is necessary to model traveltimes in the resulting 2.5D model (using 3D methods) toensure that the data fit is still acceptable. This technique is used in the processing of Lines 06 and15 of our case study (see Chapter 5), whereas the smoothing approach (discussed above) is usedfor Lines 05 and 10.2.3 First-arrival picking strategyIn order to carry out traveltime tomography, the one-way travel time of the seismic waves mustbe known for relevant arrivals. In the case of first-arrival tomography, each traveltime ?pick? rep-resents the earliest arriving P-wave signal for the corresponding source/receiver pair. In principle,it is desirable to generate these traveltime data automatically by computerized analysis of seismictraces (waveforms). In practice, this is seldom a tractable problem with field seismic data becauseof several factors, viz.:? the field seismic wavefield may include signals that arrive late in time and possess higher ampli-tude than the first arrival.? ambient seismic noise (either random or systematic) may make it difficult to distinguish the firstarrival.? the seismic wavelet may be mixed-phase, making the selection of the time-zero phase moredifficult.? the data may contain systematic precursor arrivals due to processing or acquisition artefacts thatare similar in appearance to the first arrival.? the direct and refracted waves may be shadowed by topography, causing the first arrival in theseismic data to conflict with the first arrival generated in a simplified synthetic model.In lieu of a functional automatic picking methodology, it is necessary to carry out manualfirst-arrival picking on the datasets seismic data. Examples from the Nechako Basin case studyare discussed in Chapters 3, 4 and 5. The field seismic waveforms are preprocessed for first-arrivalpicking to improve the visual appearance of the data and the coherencey of events, without fil-27050010001500-500 1 755Zero o!setReceiver numberReduced time (ms)High-amplitude late arrivalsAcausal sourceLow SNRFirst-arrival pickFigure 2.2: An example shot gather is presented in reduced time format with a reduction velocityof 4500 m/s (from Nechako Basin Line 15; see Chapter 5). First-arrival picks are indicated. Severalfeatures are indicated that affect the first-arrival picking strategy (see text).tering that would affect the waveform shape on a local scale. The seismic traces are debiased (toremove zero-frequency signal) and the samples in each trace are normalized by the maximumtrace amplitude. The initial picking is carried out on every 32nd shot gather, followed by exami-nation and possibly picking on every 32nd receiver gather, which allows visual quality control onthe first-arrival picks. The picking strategy then proceeds on every every 16th shot gather, and soon until the desired picking density is reached. Since the number of shots and receivers in a mod-ern seismic acquisition line is very large, some form of automatic or semi-automatic procedure isnecessary to reduce the picking time.GLOBE Claritas is a commercial seismic processing suite produced by a subsidiary of GNSScience. The interactive sv interface of GLOBE Claritas is employed for each stage of first-arrivalpicking. The initial pickmethod is based on amanual pick, but a semi-automatic cross-correlationautopicker is used to predictively estimate the appropriate pick for the next trace, recursivelythrough the gather (search window ?10ms, cross-correlation length 30 ms). These picks are thenupdated to fit the nearest positive inflection point in the trace on a local level. In this way, a shotgather of data waveforms that possess high signal-to-noise level may be picked in as little as 1?228minutes. Data that are noisier or that contain complex wave interference near the first arrival arepicked by hand without applying the predictive picking algorithm. As a quality control measure,seismic traces that contain systematic noise or exhibit a low signal to noise ratio are not picked.Figure 2.2 exemplifies some of these cases. The presence or absence of the first arrival pick mayalso be used in further waveform processing to exclude noisy data waveforms.Pick interpolationIt is not necessary in all cases to generate a first-arrival pick for every seismic trace. Due tothe smoothing regularization that is applied in traveltime inversion, it is possible to generate auseful velocity model with a subset of the total dataset. However, since the approximate first-arrival picks may be used in other portions of the data workflow (e.g., in selecting live traces forfull-waveform inversion), it is desirable to interpolate the portions of the dataset that have beenexcluded by subsampling, while preserving the masking of data regions that possess low signal-to-noise ratios. The picks may be interpolated to a higher density by use of a 2D spline interpolationor a simple nearest-neighbour interpolation (the latter is sufficient in our strategy); in this mannerit is possible to take a dataset that has been picked on alternate gathers and generate approximatepicks for the intermediate traces. Figure 2.3a shows a set of unmodified first-arrival picks from thepicking strategy outlined. Portions of the dataset are absent because the presence of acquisitionnoise makes picking unreliable (and would cause instabilities in FWI). The picks are interpolated(Figure 2.3b) to fill in the missing traces, but it is desirable to preserve the exclusion mask that ispresent on the original data from initial quality-control efforts. A procedure may be applied toupdate the pick mask for the interpolated dataset:1. A 1-bit image is generated based on the presence or absence of first-arrival picks in the originalimage.2. The 1-bit occupancy image is converted to floating point format in the range [0, 1] and smoothedby a 2D Gaussian kernel of a size appropriate to the initial sampling interval.3. A mask is calculated based on the value of each element in the smoothed occupancy imagecompared to an appropriate threshold value in the range (0, 1).29Source stationReceiver stationSource stationReceiver stationTraveltime (ms)Traveltime (ms)Original picksInterpolated picksa)b)Figure 2.3: Panels of a) original first-arrival picks and b) interpolated picks are displayed for an exam-ple dataset (from Nechako Basin Line 15; see Chapter 5). With the exception of those areas excludedbecause of low signal-to-noise ratios, first-arrival picks are shown corresponding to every second shotgather. The interpolated picks (b) are generated using nearest-neighbour interpolation between datavalues that exist in the original dataset. In order to preserve the exclusion of traces from quality controlefforts, the picks are masked based on a threshold calculated from a smoothed version of the originalpick mask (Gaussian ?1 = ?2 = 2.0, threshold = 0.6; see text).4. The mask is applied to the interpolated traces.The interpolation by nearest-neighbour method is too simple to be used reliably for traveltimeinversion directly, and therefore only the original picked traveltimes are employed to generate a300?4004000 20 40 45 50 55 60Iteration NumberTime (ms)Figure 2.4: Several source wavelets are shown estimated at different iterations during full-waveforminversion on Line 05 of the Nechako Basin Seismic Survey. The source signature is nominally zero-phase, with a dominant trough indicating the centre of the wavelet equivalent to a ?first arrival? in thevibroseis case.tomography velocity model. The interpolated picks are sufficient, however, to be used in com-putation of mute windows and Laplace-domain time damping terms for waveform preprocessing;the necessary precision for these uses is on the order of hundreds of milliseconds.2.4 Source signature estimationIn the vibroseis configuration, the seismic wavelet is defined as in Equation 1.31, reproducedhere:x(t) = s(t) ? w(t) ? e(t) (2.2)The wavelet x is produced by a convolution of the vibroseis sweep (s), the response characteristicsfrom the sources and geophones (w) and the earth response (e). In the case studies discussed inthis thesis, the source signatures are estimated during full-waveform inversion after the methoddiscussed by Pratt (1999)s = utd?utu? , (2.3)with s representing the complex source term for a single frequency, u representing the estimatedwavefield and d representing the field data. Hence, the estimated wavelet combines these termsand should converge towards the true wavelet x as the numerical earth model converges towards31the true earth model. Figure 2.4 shows several estimated source source signatures for differentiterations during processing for Line 05 of the Geoscience BC Nechako Basin seismic survey(discussed in Chapter 5). The signature changes little between the starting model (iteration 0) andthemodel from iteration 50. However, introduction of amplitude information in the FWI objectivefunction at iterations 51 and above results in subsequent changes in the attenuation model (seeChapter 5). The changes to the attenuation model affect how velocities are dispersed, and resultin the change in the wavelet shape seen at iterations 55 and 60 of Figure 2.4. The estimationof the source wavelet makes the full-waveform inversion process somewhat insensitive to signalsconvolved with the seismic source, and hence deconvolution is not applied prior to full-waveforminversion in my methodology (or the methodology of Pratt (1999) and related works).32Chapter 3Waveform Tomography of Field VibroseisData using an Approximate 2D Geometryleads to Improved Velocity ModelsSmithyman, B. R., and R. M. Clowes, Waveform tomography of field vibroseis datausing an approximate 2D geometry leads to improved velocity models, Geophysics,77(1), R33?R43, doi:10.1190/GEO2011-0076.1, 2012.Waveform tomography, a combination of traveltime tomography (or inversion) and waveforminversion, is applied to vibroseis first-arrival data to generate an interpretable model of P-wavevelocity for a site in the Nechako Basin, south-central British Columbia, Canada. We use con-strained 3D traveltime inversion followed by 2D full-waveform inversion to process long-offset(14.4 km) first-arrival refraction waveforms, resulting in a velocity model of significantly higherdetail than a conventional refraction-statics model generated for a processing workflow. Thecrooked-line acquisition of the data set makes 2D full-waveform inversion difficult. Thus, a proce-dure that improves the tractability of waveform tomography processing of vibroseis data recordedon crooked roads is developed to generate a near-surface (< 2 km) velocity model for the studyarea.333.1 IntroductionVibroseis multichannel seismic acquisition is typically designed to produce high-quality near-offset reflection-data with the goal of imaging reflectors at depth. The interrogation of the shal-lowest layers of a study region is a secondary goal in many seismic surveys. The applicationof refraction-statics, which typically involve a form of traveltime inversion, is normally used toaccount for near-surface heterogeneities when processing later arrivals. To a lesser extent, it issometimes used to produce velocity models that can be useful for near-surface geological inter-pretation. Due to the difficulty of picking first-arrival traveltimes at longer offsets, the full lateralextent of the geophone spread is not always used. In contrast, our traveltime and full-waveforminversion efforts are designed primarily to produce high-quality velocity models from vibroseisfirst-arrival data to aid in near-surface interpretation. To date, many successful applications ofwaveform tomography have been shown in synthetic data studies (e.g., Brenders and Pratt, 2007a;Ben-Hadj-Ali et al., 2008; Krebs et al., 2009; Shah et al., 2010) and field projects in crosshole(Pratt and Goulty, 1991; Pratt et al., 2005) and offshore settings (Dessa et al., 2004; Kamei andPratt, 2010). Onshore applications have been shown using explosive sources (e.g., Operto et al.,2004; Malinowski and Operto, 2008; Jaiswal et al., 2009), but typically include low-frequencyinformation that is not available in conventional vibroseis acquisition. An excellent overview ispresented by Virieux and Operto (2009). To the best of our knowledge, results from waveformtomography processing of band-limited on-shore vibroseis data acquired on a crooked-line havenot previously been published.In this study, we use the extended offset data available from the 2008 Geoscience BCNechakoBasin vibroseis seismic survey (Calvert et al., 2009) to demonstrate the efficacy of our approach.We develop a method that accounts for several of the limiting factors that might otherwise makewaveform tomography of vibroseis data intractable. The use of data from offsets up to 14.4 km inthe traveltime inversion process provides velocity models with depths of investigation on the orderof 2?3 km (dependent on local geology). Full-waveform inversion improves on the resolution andfidelity of the traveltime inversion velocity models by fitting the waveform amplitude and phase.34The waveform tomography algorithm that we applied is based on a 2D implementation (Pratt,1999) that could be used effectively with limited computational resources. However, the vibroseisdata were acquired on existing (crooked) roads and were not necessarily conducive to 2D analy-ses. The combination of time shifts due to out-of-plane geometry and data that are band-limiteddue to the vibroseis source provided challenges in taking advantage of the full 14.4 km horizontaloffsets of the data. Because of the crooked-line acquisition geometry used in the Nechako Basinsurvey, we found it necessary to carry out geometry corrections to make the waveform inversionproblem tractable for the 2D implementation.Improved velocity models from our studies (and the methods that produced them) have rele-vance in seismic interpretation and processing workflows for two main reasons: (1) interpretationcan provide valuable information about the near-surface region that is not well parameterizedby reflection methods; and (2) additional near-surface velocity information may be used to im-prove stacking and migration results in seismic-reflection workflows. Information from (1) can behelpful in identifying differing rock types and their relevance for further exploration or geologicalunderstanding. This is the example that is used in our study. By carrying out waveform inversionon first-arrival data, improved resolution may be achieved in the near-surface, i.e., the shallowest2?3 km (for 14 km offsets) that are normally mainly interrogated by refracted waves rather thanreflected ones. Improved velocity models also may enable better images of the subsurface whenprestack depth migration or other specialized processing procedures are carried out.3.2 Technical backgroundWe apply a technique known as waveform tomography, in which high-quality traveltime data areprocessed by traveltime (tomographic) inversion (Zelt and Barton, 1998) followed by frequency-domain, 2D acoustic full-waveform inversion of preconditioned waveform data (Pratt and Wor-thington, 1990; Song et al., 1995; Pratt, 1999; Brenders and Pratt, 2007a). This technique takesadvantage of the reduced nonlinearity of the traveltime inversion objective function comparedto the objective function found in full-waveform inversion. Ideally the velocity model resulting35from traveltime inversion is similar to a low-wavenumber model using wave-equation methods.By careful use of traveltime-inversion techniques, a starting model can be designed that is closeto the global minimum of the eventual solution.Tarantola (1984b) and Mora (1987) discuss full-waveform inversion of acoustic data by time-domain finite differences, whereas Pratt and Worthington (1990) and Liao and McMechan (1996)use frequency-domain finite-difference implementations. Vigh and Starr (2008) discuss some ofthe benefits of each. The chief advantages to a frequency-domain approach are: (1) improvedmodeling and inversion efficiency in 2D through factorization of the impedance matrix; and (2)access to multiscale inversion benefits by inverting low frequencies first. Because of its low com-putational cost, we employ 2D frequency-domain full-waveform inversion. The density modelis derived from an empirical relation using the P-wave velocity model (after e.g., Brenders andPratt, 2007a). Other workers are exploring 3D full-waveform inversion and/or waveform tomog-raphy (Plessix, 2009;Mika et al., 2010; Virieux et al., 2010;Warner et al., 2010), as well as elasticwave-equation implementations (Brossier et al., 2009a; Kamei and Pratt, 2010; Singh et al., 2010).The application of full-waveform inversion requires that the characteristics of the survey bereproduced accurately when generating synthetic data (forward modeling). This is the simplestin cases where geometry and survey characteristics are regular or easily controlled; examplesinclude synthetic studies, marine acquisition and crosshole experiments. Land acquisition ingeneral, and vibroseis acquisition in particular, is more difficult to simulate in a 2D full-waveformmodeling code because topographic features and land-use concerns often control the placementof source points and receiver groups. The incorporation of offline (y-direction) station offsetsis often unavoidable, and the irregular geometries produce effects in the data that cannot bemodeled correctly with a 2D processing workflow. The exact constraints quantifying how muchdeviation from 2D is acceptable will vary greatly depending on the processing methods used.For 2D full-waveform inversion, the traveltime errors due to incorrect encoding of offsets shouldbe much less than one half cycle at the frequency of interest, or equivalently the path-lengthdifference should be much less than one half wavelength. We have designed a method that makesthe problem tractable when these errors are small. In the absence of computational limitations,36however, the application of a 3D method (Plessix, 2009; Mika et al., 2010; Virieux et al., 2010;Warner et al., 2010) should always provide superior results.The characteristics of the vibroseis acquisition provide advantages and disadvantages for full-waveform inversion. Vibroseis commonly allows for high signal-to-noise ratios through highfoldsurveys and stacking of shots. The spatial sampling of vibration points is also typically much finerthan that found with explosive shots. This is beneficial, especially with horizontally travelingwaves, because the station spacing controls the maximum frequency that can be inverted withoutspatial aliasing. However, vibroseis seismic data are band-limited by the vibroseis sweep. Thelow-frequency signals (i.e., below 8 Hz in our case study) that are very beneficial in full-waveforminversion are typically not included in the vibroseis source.3.3 MethodWe develop a method for application of 2D waveform tomography to vibroseis data acquiredalong crooked roads. This approach incorporates approximations that are not necessary or de-sirable in a 3D full-waveform inversion workflow. Our goal is not to improve upon existing 3Dworkflows; rather, we aim to overcome obstacles that limit the effectiveness of 2D waveformtomography in cases that may not justify 3D processing.We make use of the traveltime inversion package FAST (First-Arrival Seismic Tomography;Zelt and Barton, 1998) to develop 3D velocity models that are constrained to be smooth inone dimension. Although we constrain the model updates to be homogeneous in one direc-tion, the raytracing methods applied are still dependent on the eikonal equation in a 3D model;the smoothness in the y-direction comes from model regularization constraints. Hence, this is a?2.5D? traveltime inversion method. This allows us to develop the initial model in the contextof accurate source and receiver locations in three dimensions, but to prepare for later process-ing to be carried out in 2D. We apply smoothing regularization to favor horizontal features withan anisotropic Gaussian filter. We also implement a maximum filter on the 3D model updatesin the y-direction (i.e., perpendicular to the 2D plane of interest) to reflect the sensitivity of the37eikonal equation to the maximum local velocity. This reduces the tendency of the traveltime in-version to develop narrowly supported unrealistic high-velocity features in 3D and improves thestability and convergence rate. The step length of each model update is controlled by iterativelyreweighted Tikhonov regularization.To partially account for the 3D acquisition geometry, we apply a static correction to thewaveform data before beginning full-waveform inversion. This is derived from the differencebetween the synthetic traveltime data using the 3D acquisition geometry and the desired 2Dprojected geometry.Full-waveform inversion is carried out for velocity using the program FULLWV (Pratt, 1999).We begin inversion by using successive groups of frequencies from the low-frequency part of thedata domain. The source signature is inverted based on the initial model (i.e., from traveltimeinversion), and updated at several points through the procedure (a linear best-fit source signatureestimate from the current model). The amplitude versus offset characteristics of the data are nor-malized in a log-linear fashion to fit the synthetic data (heuristically; this partially accounts forirregularities in applying 2D wave-equation backpropagation to field seismic data). The veloc-ity model updates are regularized via iterative Tikhonov regularization to be close to the initialmodel from traveltime inversion (i.e., the initial model is also the reference model). During theprocess, the reference model is updated on the basis of the earlier stages of inversion, beforeincorporating additional data frequencies. The velocity model updates are computed by a con-jugate gradient algorithm and preconditioned by 2D low-pass filtering in wavenumber, whichlimits the perturbation sharpness to the theoretical resolution supported by the maximum inver-sion frequency at each stage. The smoothing filter is mildly anisotropic (2:1) to favor horizontalfeatures (appropriate given the acquisition geometry and basin environment of the test case).3.4 Case study: Nechako BasinThe motivation and test case for our method comes from the Nechako Basin seismic explorationwork undertaken by CGG Veritas on behalf of Geoscience BC in July 2008. We carried out38waveform tomography inversion of vibroseis data collected along line 10 of the Nechako Basinseismic survey; refer to Calvert et al. (2009) for details on data collection. We used a subsetof the full line comprising 699 source points, each with 960 active receiver groups arranged ina split-spread configuration (720 to the west, 240 to the east). The source point spacing was40 m and receiver group interval was 20 m. The data were correlated vibroseis records fromdiversity stacking of four sweeps, each using four vibrators. The nominal source-group widthwas 60 m along-line. The vibroseis signal was a linear sweep from 8 to 64 Hz over 28 s (i.e.,2 Hz/s), including 0.9 s tapers. This line is located along the northern boundary of the mostprospective (south-eastern) region of the Nechako Basin, approximately 100 km west of Quesnel,south-central British Columbia (Figures 3.1 and 3.2).3.4.1 Geological backgroundThe Nechako Basin is a sedimentary basin in the Intermontane Belt of the western CanadianCordillera (Figure 3.1). This area has been characterized as prospective for hydrocarbon develop-ment. Hayes et al. (2003) identified the south-eastern portion of the basin as having the highestprospectivity.Figure 3.2 shows the study area in detail, highlighting the geometry of the seismic acquisitionline and the surface geological units. Figure 3.3 presents a simplified depositional/emplacementhistory for the study area. The structural geology of the site is generally not known sufficiently wellto produce an authoritative stratigraphic section; the closest well control comes from two wellseach about 50 km from seismic line 10, which differ greatly. Basaltic andesite and volcaniclasticrocks of the Hazelton Group (part of the Stikine Terrane) were deposited in the Early to MiddleJurassic. These were followed by mixed sedimentary rocks of the Bowser Lake Group, but theseare not observed in our area of study.A variety of terrigenous clastic sequences were deposited in the Cretaceous, exemplified bythe Middle Cretaceous Skeena Group rocks and the overlying Late Cretaceous sediments of theSustut Group. Sedimentary units are seen to be in excess of 2500m thick in some parts of the basin(Hannigan et al., 1994), but appear to be significantly thinner in our study area. Ootsa Lake rocks39Figure 3.1: Relevant geological terranes in central British Columbia, showing the location of the studyarea. The Nechako Basin is subdivided into two main regions; this study encompasses the northernedge of the south-eastern region, roughly at the boundary between rocks of the Stikine Terrane andthose of the basin proper (Massey et al., 2005).including mixed sediments and rhyolite were deposited before and contemporaneously with En-dako Group mixed volcanic rocks; the volcanic sequences were likely erupted in conjunctionwith Eocene extensional tectonics. The Chilcotin Group rocks are flood basalts typically 50 mthick (Calvert et al., 2009) erupted in the Neogene. These are covered locally by mixed Quater-40Figure 3.2: Geometry of seismic line 10 in relation to lithology and several other seismic lines from the2008 Geoscience BC Nechako Basin vibroseis survey; map area as indicated in Figure 3.1. The mapis projected to UTM coordinates (NAD 27). The surface of the central portion of line 10 is dominatedby the Ootsa Lake group, including rhyolite and other volcanic rocks erupted in the Eocene. Bothflanks are overprinted by the Chilcotin basalt (Massey et al., 2005), characterized by variable observedseismic velocities most likely due to variable vesicularity and brecciation. To the north, the HazeltonGroup volcanic rocks of the Stikine Terrane appear to plunge beneath line 10 toward the south. The2D geometry for the model profile used in waveform inversion is shown by the red line. The receiverarray for this study spans the model line; the part covered by the active source array is shown in blue.nary sediments that are quite variable over the study area. The prospective units for oil and gasexploration in the Nechako Basin are Middle Cretaceous clastic sedimentary rocks, mainly theSkeena Group, that overlie the Stikine Terrane (Hannigan et al., 1994). They contain anticlinalstructural traps formed by compression through the mid-Jurassic to mid-Eocene and normal faulttraps from subsequent extension. The overlying Late Cretaceous sedimentary rocks are classifiedas nonprospective based on the absence of oil and gas shows, despite the existence of structuressimilar to those in the Skeena Group rocks (Hannigan et al., 1994). They do not outcrop in the re-gion of our study. The Hazelton Group, part of the Stikine Terrane, outcrops near the seismic lineas shown in Figure 3.2 (Massey et al., 2005). The near-surface is dominated by Eocene volcanicrocks of the Ootsa Lake and Endako groups, and by the younger Neogene Chilcotin basalt. Thesevolcanic units also dominate the near-surface seismic response. Based on rock-physics results,the Chilcotin basalt was originally expected to show a higher bulk P-wave velocity than the Ootsa4120014513097652352452.6GDGDGDGDHazeltonGroupSkeenaGroupSustutGroupOotsa LakeGroupEndakoGroupChilcotinGroupQD JurassicCretaceousEoceneNeogeneP?OG0D (Ma)166Figure 3.3: A schematic stratigraphic column for the region of study, outlining the deposi-tional/emplacement sequences relevant to line 10 of the Nechako Basin seismic survey. D?dioriteplutons; GD?granodiorite plutons; OG?Oligocene; P?Paleocene; Q?Quaternary.42Lake and Endako units. However, recent work by Calvert et al. (2011) indicates that brecciationmay cause the Chilcotin basalt to possess lower P-wave velocities throughout much of the region.Quaternary deposits of differing types and varying thicknesses overlie the older rocks (Figures3.2 and 3.3). Recent work by others has increasingly found that the Nechako Basin is more ap-propriately thought of as a collection of related subbasins (J. Riddell, personal communication,2011).3.4.2 Traveltime inversionVelocitymodels of the shallow subsurface are typically developed as part of the reflection-seismologyworkflow to facilitate common depth point (CDP) stacking and migration; however, these modelsare often coarse and of limited use for interpretation. When applying waveform tomography, itis important to produce detailed velocity models at the traveltime inversion stage; this, in turn,requires precise traveltime picks. Without a sufficiently accurate starting model, waveform inver-sion methods cannot succeed. We picked first-arrivals over the full (14.4 km) offset range of 699shot gathers; the bandwidth present in the early arriving waveforms was approximately 8?16 Hz.We picked a particular phase in the band-limited data, hence the traveltime picking error is likelyon the order of 20?30 ms. The picking error was controlled by the dominant data frequency inthe first-arrivals, which remained relatively consistent over the offset range of interest. The as-signed error for the picks was chosen to be uniform across the data set; however, low-confidencepicks were manually removed before inversion. Traces with large traveltime data residuals wereexcluded from subsequent processing of the data waveforms. We found that very careful consid-eration was necessary due to significant shingling in the first-arrival waveforms. This was due inpart to the attenuation of high-frequency parts of the data waveforms and the (nominally) zero-phase source imposed by the vibroseis correlation.We initially processed these picks using FAST (Zelt and Barton, 1998) and GLI3D (a well-known commercial package; Hampson and Russell, 1984). Both codes were able to fit the first-arrival picks in 3D with an rms traveltime misfit of 20?25 ms. Because of the need for manualcontrol over the traveltime-inversion process, we primarily used FAST (an academic code). To43find the smoothest model that fit the traveltime data within the expected misfit range, the velocitymodel was regularized by an anisotropic 3DGaussian filter at each iteration. The filter was chosento enforce smoothness in the crossline direction and to encourage a smooth model in the 2Dplane of interest. A 5:1 horizontal to vertical smoothing ratio was chosen to reflect the expecteddominance of horizontal geological features and the refracted arrivals used in our procedure. Thesize of the smoothing filter was chosen empirically to guide the model toward an rms data misfitin the target range. Operating in constrained 2.5D (i.e., 3D geometry and 2D velocity model),FAST was able to fit the first-arrival data within an rms misfit of 27 ms. Given that dominantfrequencies of the first-arrival data were 8?16 Hz, the resulting velocity model was appropriateas a low-wavenumber starting model for full-waveform inversion.Figure 3.4a shows real and synthetic first-arrival data and residuals resulting from forwardmodeling in the constrained 3D FAST model. The resulting velocity model is shown in Figure3.5a. The data residuals are approximately normally distributed, with an rms misfit of 27 ms,i.e., comparable to the picking error. The true 3D geometry of the survey cannot be representedcorrectly in the 2D full-waveform inversion process. Figure 3.4b shows synthetic traveltime datacalculated in the 3D geometry (as in Figure 3.4a) and a best-fit 2D projected geometry. The corre-sponding residuals represent the traveltime error resulting from the 2D geometry approximation(under the infinite-frequency approximation and in the velocity model of Figure 3.5a). The magni-tude of these errors corresponds closely with the locations along the acquisition line in which theoffset error is high (from the 2D projection; compare with the line geometry visible in Figure 3.2).We used the traveltime residuals presented in Figure 3.4b to static correct the data waveformsbefore carrying out full-waveform inversion.44Figure 3.4: Traveltime and residual plots are shown for (a) real and 3D synthetic data, and (b) 2D and3D synthetic data. The horizontal axis represents source-receiver midpoint along the line for eachresidual. The 3D synthetic data (generated in a 2D model) predict the true data within an rms misfitof approximately 27 ms, and the residuals are well distributed about zero mean. The 2D syntheticdata differ from the 3D synthetic data due to the projection of the geometry onto a plane (strikingapproximately 106?).45Checkerboard testing result(a)(b)(c)(d)Figure 3.5: Velocity models from traveltime inversion (a) and subsequent full-waveform inversion (b). The difference (c) is shown to helpidentify features that are changed in the full-waveform inversion. The result of a corrugation testing procedure (d) is provided to help judgethe expected resolution. The waveform data-fit from traveltime inversion was excellent in the eastern portion of the model, leading to verysmall perturbations from full-waveform inversion.46Shot 225 Shot 475SyntheticReal3210Traveltime (s)14.4 E 4.8 W0.0Source-receiver o!set (km)3210Traveltime (s)14.4 E 4.8 W0.0Source-receiver o!set (km)3210Traveltime (s)14.4 E 4.8 W0.0Source-receiver o!set (km)3210Traveltime (s)14.4 E 4.8 W0.0Source-receiver o!set (km)Line noiseTime shiftfrom geometryShallowre"ectionsTime shiftdi#cult to modelRe"ectionsHigh amplitudeand velocityModelledaccurately Secondary arrivalsreproducedHigh-amplitude$rst arrivalTime shifted$rst arrivalsFigure 3.6: Real (top) and synthetic (bottom) data are shown for two shot gathers; locations shown in Figure 3.2. Shot 225 (left) is represen-tative of the poorer data fit seen in the western portion of the model, due primarily to problems in approximating the geometry. Shot 475(right) is representative of the high-quality data fit found in the eastern portions of the model. Annotations highlight some of the relevantdata features.473.4.3 Full-waveform inversionThe velocity model in Figure 3.5a was the input (starting) velocity model for full-waveform in-version. The goal of the full-waveform inversion process is to improve on the velocity modelprovided by traveltime inversion in two main ways:? Finer spatial resolution can be expected because of the migration-like generation of the modelupdates.? Ability to resolve low-velocity zones is possible due to the incorporation of multiple phases andamplitude information.The aforementioned static corrections were effective at improving the quality of the initialwaveform fit. To test this, initial synthetic waveforms (also used in inversion) were calculatedusing a delta function source signature. The initial correspondence between the phase of thereal data and synthetic data was high in the eastern portion of line 10 but lower toward the west.This is likely due to the large offline displacements of sources and receivers (relative to the ideal2D line) at the western end. Smaller displacements were well accounted for by the time-shiftingmethod.We carried out full-waveform inversion between 8 and 11 Hz, using an approach in whichwe iteratively increased the frequency content of the inversion by fractions of a Hertz. The upperend of the inverted frequency range was limited by the ability to fit the real data waveforms(discussed below). The initial stages used 8.0, 8.5 and 9.0 Hz frequencies; subsequent stagesincorporated data up to 11 Hz at 0.5 Hz intervals. At early stages, a delta function source wasused to approximate the response of the ideal broadband vibroseis source. Once frequenciesabove 9 Hz were incorporated, the inversion proceeded with a synthetic source signature deriveddirectly from the data. The resulting updated velocity model is presented in Figure 3.5b.Corrugation testing was carried out to aid in judging the spatial resolution in different parts ofthe velocity model. This was performed by multiplicatively perturbing the waveform tomographyvelocity model (Figure 3.5b) by +/?5% in an alternating checkerboard pattern. The maximumperturbation was 344 m/s, and each cell of the checkerboard was 1000 m wide by 500 m high.48Synthetic data were modeled using the perturbed velocity model. Figure 3.5d shows the recov-ered perturbations after inversion with the model from Figure 3.5b. The spatial resolution of themethod is finer in the near-surface where the angular coverage is the greatest and the propaga-tion velocities are the lowest. The wavepaths are predominantly horizontal at depth, which limitsthe maximum achievable horizontal resolution and depth of penetration with inversion of datafrom transmitted wave energy. Vertical resolution remains reasonable at a depth of 1?2 km, buthorizontal resolution is most likely coarser than 1 km. In regions where the acquisition line isstraight, fitting of reflection arrivals is expected to improve the depth of penetration. However,fitting reflections is less valid in parts of the model where the array contains large crossline offsetsand consequently the focusing at depth is less effective.Figure 3.6 shows a comparison of real and synthetic data from two representative shot gathersin the data set. This acts as a quality control on the success of the method, and in particularprovides valuable information about which features in the model are robust. The data werewindowed in time for use in the full-waveform-inversion process, and therefore we are particularlyinterested in the data fit within specific regions. Analysis of the frequency content of the earlyarrivals suggests that these data may support additional higher frequencies (up to?16 Hz) in someparts of the model (e.g., the well-characterized region near shot 475). The maximum supportedfrequency from the source array spacing is approximately 24 Hz, based upon a dominant near-surface velocity no less than 1500m/s. However, the early arriving data have very low amplitudesabove 16 Hz, which limits significantly the benefits of additional high-frequency processing.The quality of the data fit in the eastern regions of the model was relatively high from theinitial model, resulting in small perturbations, on the order of +/?50 m/s, to the velocity model.However, much stronger perturbations, which dominate the difference image in Figure 3.5c, wereproduced on the western end of the model. The numerical resolution of the inversion algorithmin the western portion of the model was lower due to reduced data fold west of the 12 km pointon the projected line, which is reflected in Figure 3.5d. The high-amplitude model perturbationsin the western part of the line can most likely be traced to problems in the initial data fit from thetraveltime tomographymodel, which in turn are due to unresolved geometry errors. The geometry49correction incompletely accounts for approximation errors in the full-waveform inversion stage,which are significant toward the western end of the line (see geometry in Figure 3.2). The staticcorrection improves the full-waveform inversion fidelity in the very near-surface (less than about 1km depth) across most of the seismic line, but substantial model artifacts remain from the seismictraces with the largest crossline source-receiver offsets. In the shallowest regions, low-velocityzones are resolved that were not identifiable in the model from tomography alone. These arehighlighted well in Figure 3.5c ? for example, the small high-velocity anomaly visible at 17 kmis resolved sharply in the waveform tomography model and the region surrounding it exhibitslower velocity; the square-sided feature at 22 km exhibits sharper edges than in the model fromtomography, with corresponding low-velocity features to either side.To apply advanced processing techniques using the velocity model we have created, addi-tional steps would be necessary. Since the tomographic arrivals we invert interrogate only theshallowest 2?3 km of the site, additional velocity analysis would be necessary to establish aninterval velocity model throughout the depth-range of interest. An examination of the anisotropycharacteristics of the site would also be necessary if applying the model from waveform tomogra-phy to a problem such as prestack depth migration of reflected waves. Full-waveform inversionin 2.5D or 3D would preserve the ability to invert reflection-data as well as the tomographicarrivals (though without necessarily accounting for anisotropy), and increase the useful depth ofinvestigation of the method; however, these methods come with additional computational costand are not part of the current study.3.4.4 InterpretationThe velocity models derived through traveltime tomography and full-waveform inversion can bedirectly interpreted to assist in developing a geological model of the region. Because this repre-sents a location where the structural geology is poorly understood, we compare our results withthose from other workers, who used different data-processing methodologies and geophysicalmethods. The region of highest confidence in our work is the eastern portion of the line, due tothe relatively low line-curvature. Figure 3.2 presents the surface geological features and Figure500123Depth (km)Quaternary Chilcotin basalt Ootsa Lake rhyolite Quaternary Chilcotin basaltABC CDEFG KDepth (km)012345MNOP QR210Time (s)Folded Chilcotinvolcanic rocks Resistive volcanic rocksResistive Hazelton basementHJL(a)(b)(c)Resistivity (Ohm?m)1021011031041500600045003000Velocity (m/s)Depth (km)012345MNOP QR0750050002500Velocity (m/s)PreSDM Velocity (WesternGeco, 2010)PreSDM Re%ectivity (WesternGeco, 2010)PostSTM re%ectivity and MT resistivity (Calvert et al., 2011)Figure 3.7: (a) The velocity model from Figure 3.5b is shown annotated with features of interest(described in main text). Color legend is identified by the right-angle arrow. (b) A portion of thepoststack time-migration image overlain with resistivity values (colors), based on resistivity valuesfrom magnetotelluric studies; labels (H, J, L) have been inserted (modified from Calvert et al., 2011).The horizontal and vertical resolution are variable between the two methods. These models sharesome complementary features, especially with regard to the identification of subbasins. However,the time stretch present in the result from CGG Veritas processing (Calvert et al., 2011) makes directcomparison of features challenging. Color legend is identified by the horizontal arrow. (c) Upper:velocity model (upper 5 km only) from prestack depth-migration processing by WesternGeCo MDIC(2010) for comparison with our velocity models. This velocity model is shown for completeness, buta direct comparison with the result in (a) is inappropriate, given the different intended uses of the twomodels. Lower: Prestack depth-migrated reflectivity image from WesternGeCo MDIC (2010). Lettersare added to facilitate discussion (see text). Dashed lines identify several interfaces interpreted by usfrom the velocity and reflectivity images.3.3 provides a simplified stratigraphic sequence. Figure 3.5 shows the velocity model derived forinterpretive purposes and any subsequent reflection-data reprocessing (to be done by others). Forinterpretation, we refer to labels in Figure 3.7; the velocity model in Figure 3.5b is reproducedin Figure 3.7a and 3.7b. Several of the features in (3.7b) the time-migration image are deempha-sized by the processing and/or depth conversion; however, common features can be identified51on (3.7a) the velocity model from waveform tomography and the PreSDM reflectivity image.The presence of the Chilcotin basalt in the eastern portion of the survey region may be re-sponsible for a local high-velocity anomaly (A). This corresponds with rock-physics informationsuggesting that the Chilcotin basalt is distinguishable by high seismic velocities relative to thenearby Ootsa Lake and Endako groups. However, recent work (Hayward and Calvert, 2009;Calvert et al., 2011) has found that the velocity of the Chilcotin basalt in situ is typically some-what lower than that of the nearby Eocene volcanic rocks. Thus, the high-velocity feature maybe due to a different unit. Note that this is at the far eastern extent of our model, and may notbe as well resolved. Immediately west of this, there is an apparent increase in the recoveredheterogeneity in the near-surface (B). This is most likely due to volcanic rocks underlying theQuaternary cover; however, the variability in velocity likely occurs at a finer scale than we areable to distinguish with the chosen model discretization interval. This feature appears to correlatewell with the distribution of recent surface sediments, extending for about 10 km over the middle-eastern portion of the line. Alternatively, the apparent heterogeneity may be due to near-surfaceor coupling effects that are too shallow for our model to resolve. Seismic high-velocity features atC are interpreted to be due to the Hazelton Group volcanic and volcaniclastic rocks, which out-crop immediately north of line 10. This unit most likely plunges southward beneath the Eocenevolcanic rocks that dominate the near-surface. This interpretation corresponds to the knowledgethat the basin is deeper toward the south (Hannigan et al., 1994), but it does not preclude the pos-sibility that Hazelton Group outcrops toward the north could be responsible for the high-velocityresponse. Because of the presence of this feature, we have not attempted to interpret subbasinstructures in this part of the model.A deepening of the low-velocity region (?3000?3500 m/s) to approximately 700 m (D) isinterpreted as a subbasin. This appears, from interpretation of velocities, to be infilled by OotsaLake rhyolite, suggesting that the formation of the basin structure is at least Eocene in age. Theambiguity of the velocity information (because of overlapping velocity ranges in different rocktypes) makes it difficult to determine the deeper basin rock types from seismic results alone. Asharply defined high-velocity feature at E could be a volcanic plug. Another subbasin structure52is seen farther west at F. However, the confidence of the full-waveform inversion is lower in thisregion, due to the poorer performance of the 2D geometry approximation (Figure 3.2). An 8?10km wide synformal structure at G was initially interpreted to be an artifact of the ray tracing usedin the initial model but subsequent comparisons with other studies have modified this view (seebelow).Calvert et al. (2011) presented a preliminary interpretation of the vibroseis seismic-reflectionsections collected on behalf of Geoscience BC (see Calvert et al., 2009; modified in our Figure3.7b), combined with magnetotelluric information (Spratt and Craven, 2010). The time-migratedreflection-data primarily represent deeper structures than we might expect to see from the first-arrival refraction data alone, but there is a region of overlap where results of the two methodscan be directly compared. Additionally, interpretations of rock types from magnetotelluric resultscomplement seismic interpretations. Calvert et al. (2011) interpret Eocene and younger extensionin the seismic-reflection section through the middle-western portion of line 10 (?3?27 km on thescale of Figure 3.7). They identify highly resistive Hazelton Group basement rocks from mag-netotelluric studies in a location that corresponds well with the high-velocity features identifiedabove (C). Furthermore, a full-graben structure is interpreted between 16 and 25 km (H), and ahalf-graben east of 3 km (J; distances relative to Figure 3.7). They find evidence for post-Eoceneextension in the Chilcotin basalt near the surface in the western portion of line 10 (Figure 3.2),leading to folding and brecciation. Additionally, recent reinterpretation of older seismic results(Hayward and Calvert, 2011) indicates strike-slip faulting perpendicular to the seismic line inquestion. This may explain regions of local variability and localized low-velocity zones.Considering results from Calvert et al. (2011), several additional features of interest may beidentified. We see a reasonable correspondence between the graben structure they interpretedat H and a deepening of the bedrock interface in our results (D, F; Figure 3.7). The geologicalunit at depth here is not well constrained, but we presume the presence of Cretaceous sedimentsbetween the younger volcanic rocks and Stikinia. The implication is that the graben structure wasinfilled by Eocene volcanic rocks after or contemporaneous with Eocene extension. However,we identify a high-velocity feature between D and F that we assume to bound the two subbasins.53Additionally, the subbasin we identify from velocity information at D is significantly shallowerthan the graben structure identified by Calvert et al. (2011). The feature we interpret at E is in aregion possessing variable resistivity (Figure 3.7), but it is not identified in the seismic-reflectioninterpretation.Our results show high-velocity, near-surface features and heterogeneities (compare the near-surface in subbasins F and D) that are spatially correlated with the brecciated Chilcotin basaltinterpreted by Calvert et al. (2011). However, the full-waveform inversion responsible for pro-ducing this heterogeneity in the model was negatively affected by geometry errors in this region,and the exact placement and extent of the heterogeneity is not well constrained. Likewise, we seea high-velocity region at F that cannot be well constrained, although the heterogeneity presentcorresponds to an interpreted fault from the reflection-seismic and magnetotelluric work. Thehigh-velocity feature at K seems to have the same orientation and location as a highly reflectivefolded package (L) in the results of Calvert et al. (2011). Resolution and reliability of our modelwest of the 10 km point (near K, Figure 3.7) is very degraded, due to the omission of back-shotsbeyond this point; however, we have included it for comparison with Calvert et al. (2011) andfuture work.Recent prestack depth migration (PreSDM) processing of the Geoscience BC line 10 data setbyWesternGeCo MDIC (2010) provides information about reflectivity in a geometry compatiblewith our velocity model. Velocity and density models used for PreSDMwere determined byWest-ernGeCo MDIC (2010) from joint inversion of magnetotelluric, gravity and seismic information;the upper part of their velocity model is reproduced for reference purposes in Figure 3.7c. Theirnear-surface velocity structure is similar to the traveltime tomography and waveform tomographymodels in this paper (Figure 3.5), but is significantly smoother. Information in the near-surfacewas developed primarily from maximum-amplitude traveltime inversion and gravity data inver-sion, whereas the deeper model was constrained mostly by iterated migration velocity analysisand magnetotelluric data inversion. To date, these results have not been interpreted by otherresearchers, and we attempt to form a joint interpretation alongside our result. It is inappropriateto compare directly the velocity model from waveform tomography (Figure 3.7a) and the velocity54model used in migration (Figure 3.7c, upper), because the intended applications of these resultsare very different. A direct comparison between the resulting migrated image and our waveformtomography velocity model (Figure 3.7c, lower and Figure 3.7a, respectively) is more appropriate.This comparison shows strong reflection events at several features we previously identified.A region of strong reflectivity is present in the shallow near-surface atM (Figure 3.7c) that co-incides spatially with a heterogeneous region identified at B in the waveform tomography velocityresult. This is also coincident with the presence of Quaternary cover on the interpreted surficialgeology (overlaid on Figure 3.7a). High reflectivity is also visible at the same lateral location in thetime-domain reflectivity section from Calvert et al. (2011), and likely represents a reflection off ofthe top boundary of Stikine Terrane basement rocks of the Hazelton group (i.e., C). Strong reflec-tivity is also present at N, which appears to delineate a continuous basin reflector encompassingthe same extent as features D and F in Figure 3.7a. Little indication is present of a feature corre-sponding to the local high-velocity anomaly at E; however, the occurrence of moderately strongnear-surface reflectivity near O corresponds spatially with velocity heterogeneities observed inthe waveform tomography velocity model (Figure 3.7a). Additionally, a reflection feature belowO resembles the local shallowing of the velocity interface at F and may be a feature that is dis-tinguishable from the interface continuing from N. A second subbasin is interpreted at P, whichagrees with a feature seen in our velocity model (G) and the half-graben (J) interpreted by Calvertet al. (2011); however, the western extent of this feature is not resolved well. Local high reflectiv-ity at Q and the presence of a high-velocity anomaly at K are interpreted as a shallowing of theflanking subbasins. However, this feature could be due to heterogeneities in the Chilcotin basaltat surface or a plutonic feature. A deeper reflector (at a depth of ?3 km) is observed at R and thepresence of linear features suggests that this may be a continuation of the interface interpreted atN. This may indicate volcanic in-filling of a deeper basin, which would explain observations ofshallower high-velocity features (Figure 3.7a) and high conductivity (Figure 3.7b).553.5 ConclusionsWe present results of waveform tomography applied to multichannel seismic data from the Geo-science BC Nechako Basin seismic survey. To successfully process these data, we tested severaltechniques to overcome the limitations of the 2D full-waveform inversion techniques we applied.Careful initial model building and static corrections can make this problem tractable. However,we found that geometry effects that might usually be ignored in MCS reflection processing cancause significant problems for full-waveform inversion (especially in the near-surface). Modelingand inverting data acquired in a 3D geometry using 2D techniques requires careful and time-intensive correction of the waveform data, and the assumptions used to make these correctionsaffect the result significantly. We built the initial model using 3D first-arrival traveltime inversion,which enabled static corrections to be applied to the waveform data; one result is that the abil-ity to incorporate nontomographic waveform data is hampered. Consequently, this technique ismostly useful for incorporating amplitude information from waveforms that interrogate a similarsubsurface region to the tomographic arrivals modeled by ray-theory methods.In general, 3D full-waveform inversion techniques do not suffer from the geometry approxi-mation errors that limited this workflow, and are therefore expected to produce a better result. Inproblems similar to our test case, in which offline offsets are not large, a 2.5D full-waveform in-version may present the best compromise between 2D and full 3D methods. This has the benefitof properly accounting for 3D geometry and geometric spreading effects, but with much of thesimplicity of a 2D approach.The detailed velocity information available from waveform tomography processing of seismicdata is valuable in making interpretations about features in complex near-surface geology. This isespecially true when rock-physics information is able to constrain the range of velocities expectedfor each rock type. However, the relative difficulty of interpreting velocity information alonemeans that joint interpretation along with other geophysical information (e.g., high-resolutionPreSDM images and other physical property measurements) is necessary to present a completegeological model.56Chapter 4Waveform Tomography in 2.5D:Parameterization for Crooked-LineAcquisition GeometrySmithyman, B. R., and R. M. Clowes, Waveform tomography in 2.5D: Parameteriza-tion for crooked-line acquisition geometry, JGR Solid Earth, doi: 10.1002/jgrb.50100,2013.A method for 2.5D viscoacoustic waveform tomography that can be applied to generate 2Dmodels of velocity and attenuation from inversion of refraction waveforms on land seismic reflec-tion data acquired along crooked roads is developed. It is particularly useful for typical crustalreflection surveys. First-arrival traveltime tomography is applied using a 3D method, but withconstraints on the intermediate 3D velocity model; the result is the starting model for the nextstep. A 2.5D frequency-domain full-waveform inversion stage parameterizes 3D geometry in theseismic source and receiver arrays, with the assumption that the velocity and attenuation mod-els are homogeneous in a prescribed out-of-plane direction. This approach results in superiorresults compared to a strictly 2D approach when the acquisition line is crooked, with a moderateincrease in computational cost. A case study using data acquired in the Nechako Basin in south-central British Columbia, Canada exemplifies and validates the procedure. The velocity modelderived from 2.5D waveform tomography is compared with that from a previous study in which2D waveform tomography was applied to the same dataset and with results from 3D traveltime to-57mography. The resolution and accuracy of the velocity model from 2.5D waveform tomographyare demonstrated to be greater than those from traveltime or 2D waveform tomography. A modelof viscoacoustic attenuation, which was not possible in the 2D case, is also generated. Thesemodels are interpreted jointly to highlight features of geological interest, such as a sedimentarybasin, basement rocks and faults, from surface to about 3 km depth.4.1 IntroductionWaveform tomography is a method for building seismic velocity models through multiple appli-cations of nonlinear inversion of seismic waveform (phase and amplitude) data. An initial stageof inversion is carried out via a conventional raytracing tomography approach, typically usingan iterative non-linear inversion framework that seeks to improve a simple starting model by in-version of picked first-arrival traveltime data. The model resulting from traveltime inversion isthen iteratively updated by full-waveform inversion of the seismic dataset to better predict thedata waveforms. The successful application of the combination of traveltime inversion and full-waveform inversion is predicated on the assumption that the model from the initial tomographystage is representative of a low-wavenumber version of the ultimate result.We implement a 2.5D method for the full-waveform inversion stage. It reproduces resultsequivalent to those of a 3D method with a specially constrained model. This leads to significantcomputational benefits, in terms of reduced processing time compared to a 3D workflow. Theuse of a 2.5D method allows for improved accuracy when modeling 3D geometric spreading ef-fects that cannot be reproduced in a 2D workflow. By constructing the equivalent 3D wavefieldin the receiver plane, the method accounts for a crooked-line acquisition geometry that wouldmake 2D processing unfeasible. Waveform tomography in 2.5D is thus particularly applicable tomost crustal reflection data, which are generally acquired in 2D geometry along existing, crookedroads. Waveform tomography cannot be applied to such data without some allowance for thecrooked line acquisition geometry, such as an approximate geometric correction by traveltimestatics (Smithyman and Clowes, 2012), a modification of the modeled source and receiver loca-58tions in 2D to account partially for 3D offset (not considered here), the approach discussed in thispaper, or application of full 3D waveform tomography.The first part of this paper provides a discussion of the 2.5D full-waveform inversion method,including special focus on the extension of the method to account for 3D geometry effects, whichwe believe represents a new development. We then provide a real-data case study to contrastthis approach with a more typical 2D approach based on reprocessing of previous work. Ourexamples derive from a dataset acquired in south-central British Columbia, Canada.4.2 BackgroundFull-waveform inversion may be implemented using a variety of approximations and simplifica-tions, to suit the computational strategy to the problem at hand. One common choice is the useof the acoustic (or viscoacoustic) wave equation in two spatial dimensions. This limits the size ofthe problem to a manageable number of model parameters, and the forward-modeling problemcan be efficiently formulated in the frequency domain as a boundary-value problem separableover frequencies. This approach is attractive for use with seismic data acquired using 2D seismicacquisition methods. However, the physical world is three-dimensional, and field seismic dataalways include some information that cannot be readily synthesized using a two-dimensional nu-merical code. The use of the acoustic approximation also leads to inaccuracies in modeling (andconsequently errors in inversion), which we discuss in Section 4.3.4.Active-source waveform tomography and full-waveform inversion have become trusted toolsin research (Takougang and Calvert, 2011; Malinowski et al., 2011; Kamei et al., 2012), and arebeing increasingly used in industry settings as well. A number of research groups and companieshave developed practical massively parallel full-waveform inversion programs that operate in 3D(e.g., Ben-Hadj-Ali et al., 2008; Plessix, 2009). The researchers responsible for these efforts aimto tackle not only the challenging problem of 3D full-waveform inversion, but also the significantcomputer science problems (and computational costs) relating to the scale of the science. Thesecomputational issues are in most cases inseparable from the geophysical problem. There exists,59however, a middle ground between costly 3D processing and comparatively simplistic 2D pro-cessing: a ?two and one half dimensional? method solves a simplified 3D problem in such a waythat the overall solution is formed from a series of solutions to smaller problems (Bleistein, 1986).Each of these smaller problems is a 2D problem of much lower complexity than the full 3D prob-lem. For certain classes of problems, the use of a 2.5D method can convey some of the benefitsof 3D processing, while reducing significantly the computational demands in comparison to a3D approach. We extend the work of other researchers in the field (Song and Williamson, 1995;Song et al., 1995; Pratt et al., 1996; Zhou and Greenhalgh, 2006) to include the ability to handle3D geometry.The use of a 2.5D processing methodology for full-waveform inversion is not, in and of itself,a new idea. The 2.5D approach to generating 3D-equivalent wavefields is used commonly in thefield of seismic migration. Song and Williamson (1995) discuss the use of a 2.5D method for for-ward modeling of acoustic waveforms, and Song et al. (1995) extend the method to full-waveforminversion in a companion paper. More recently, Zhou and Greenhalgh (2006) present a method-ology for efficiently sampling the transverse wavenumber component when synthesizing seismicdata. In the past decade, much research in this field has been towards full-waveform inversion in3D (Ben-Hadj-Ali et al., 2008; Plessix, 2009) and improved accounting for elastic wave propaga-tion (Sears and Barton, 2010). Because of the increased cost when compared to 2D processing(significant at the time), the 2.5D method of Song et al. (1995) did not enter into common use.However, there exists a plethora of seismic data collected along 2D acquisition profiles that donot contain cross-line information necessary to solve for a 3Dmodel, but do contain 3D geometrythat is irreconcilable with a 2D approach. Furthermore, while the computational cost of 2.5Dforward modelling (and inversion) is somewhat higher than the equivalent 2D approach, the costis still on a similar computational order; 3D processing is two computational orders more expen-sive than 2D processing. Significant savings can be made by choosing 2.5D processing over afull 3D implementation. This is often a controlling factor in academia, or for small companiesworking in the geophysics or engineering sectors.604.3 MethodologyWaveform tomography comprises a series of processing and inversion steps designed to producea model of velocity, and possibly viscoacoustic attenuation (Song et al., 1995; Brenders and Pratt,2007a), that optimally predicts a set of field seismic data. An initial 1D model is constructedof velocity variation with depth, based on t-x interpretation or known average properties. Thismodel is then updated through traveltime inversion (Section 4.3.1) to generate a gridded modelthat optimally predicts the seismic first-arrival data. Full-waveform inversion is applied (Section4.3.3) to update the intermediate model from traveltime inversion, imaging structures that areilluminated by wave arrivals and that are not well-predicted in the tomography stage (which islimited by asymptotic ray theory). Because of the use of the acoustic wave equation, isotropyassumptions, simple representation of the source and receiver arrays, and a variety of other ap-proximations implicit in the method, some heuristic or empirical corrections must be appliedbefore and after full-waveform inversion (Section 4.3.4). The result is assessed based on a priorigeological knowledge and comparison of synthetic data waveforms with the true data waveforms,for purposes of quality control.4.3.1 Traveltime inversionBecause full-waveform inversion is a highly nonlinear process, it is necessary to approach theinversion with a reasonably good estimate of the true model. This model (called the starting modelor initial model) is iteratively updated by the full-waveform inversion process to produce a modelthat better predicts the data waveforms. In the method of waveform tomography, the startingmodel for full-waveform inversion is typically built by a more conventional velocity analysisstage. This is less computationally intensive than the later full-waveform inversion stage, but noless important to the final result.We use the program FAST (First Arrival Seismic Tomography; Zelt and Barton, 1998) to pro-duce a starting model for full-waveform inversion. A first-arrival traveltime datum is picked (typi-cally by hand) for each of the seismic traces to be considered in the later full-waveform inversion.61This traveltime represents the arrival of the earliest waveform to reach a given part of the receiverarray from the selected source. FAST is initialized with a simple (near-1D) model of the estimatedvelocity structure of the site, and models these traveltimes using an eikonal equation method (i.e.,fast marching; Vidale, 1990; Hole and Zelt, 1995). The predicted traveltime for each trace iscompared to the measured traveltime (i.e., the pick) and the residual is weighted based on theinverse of the expected picking error. The weighted residuals are spread across the (smoothed)numerical raypaths by iterative backprojection (Hole, 1992)) to form a descent direction at eachiteration (equivalent to a gradient method of inversion). The descent vector is regularized by theaddition of scaled vectors representing the roughness in the vertical and horizontal directions,which encourages a smooth model. The initial model is perturbed with a series of modified ver-sions of the regularized gradient, scaled by a step length. The new model that optimally reducesthe data misfit is found, and the process is repeated using the new model as the starting modelfor the subsequent iteration. A conjugate gradient scheme is followed until the data are fit withinthe expected picking error.In order to produce a model compatible with the application of 2.5D full-waveform inversion,it is necessary to extract a 2D slice from the 3D volume representing the traveltime tomographymodel. To simplify this process, we implement two additions to the method of Zelt and Barton(1998): 1) the local grid of the model is oriented such that the x- and z-directions match those ofthe consequent 2Dmodel, and 2) the 3Dmodel from FAST is strongly smoothed in the y-directionat each iteration, to ensure that the model heterogeneity is nominally 2D. The resulting volumeis then sliced along the x-z plane and possibly re-sampled for use in full-waveform inversion.4.3.2 Full-waveform forward modelingThe numerical forward modeling in this study is carried out by the method of Pratt (1999), withmodifications discussed in the following text. We proceed in a similar fashion to the discussionin Pratt et al. (1996), for the case when y is not necessarily zero. The 3D acoustic wave equation(Equation 4.1) describes the propagation of acoustic waves in space and time, through a region62of acoustic velocity c (x, y, z).(?23 ? 1c2(x,y,z)?2?t2)p (x, y, z; t)= ?s (t) ? (x|xs) ? (y|ys) ? (z|zs) (4.1)However, the method of frequency-domain full-waveform synthesis solves the Helmholtz equa-tion (Equation 4.2), which represents the time-independent form after applying the method ofseparation of variables.(?23 + ?2c2(x,y,z))P (x, y, z;?)= ?S (?) ? (x|xs) ? (y|ys) ? (z|zs) (4.2)A spatial Fourier transform is applied to Equation 4.2 that separates the result into a number of2D wave equations as presented in Equation 4.3, with the additional assumption that the modeldoes not vary in the direction described by the y-coordinate (i.e., c(x, y, z) = c(x, z)). This repre-sents a decomposition over plane waves in the y-coordinate, and therefore leads directly to therequired assumption of infinite support in that axis (i.e., model homogeneity).(?22 + ?2c2(x,z) ? k2y)P? (x, ky, z;?)= ?S (?) ? (x|xs) ? (z|zs) e?iky(y?ys) (4.3)In Equation 4.3, ky is set separately for each Fourier component; this means that each solutionfor a different value of ky is independent. In order to reconstruct the wavefield equivalent to thesolution of the 3D wave equation, an inverse spatial Fourier transform is applied, i.e.:P (x, y, z;?) =?P? (x, ky, z;?) dky (4.4)It is convenient to solve the 2D Helmholtz equation (in Equation 4.3) for the case in the y = 0plane, i.e., with the omission of the transform of the ? (y|ys) term from Equation 4.2, as shown in63Equation 4.5.(?22 + ?2c2(x,z) ? k2y)P?0 (x, ky, z;?)= ?S (?) ? (x|xs) ? (z|zs) (4.5)In this case, P?0 is the wavefield in the plane of the source, and the Fourier synthesis in Equation4.4 is modified to contain a phase shift applied to P?0, which results in the same synthesizedwavefield P :P (x, y, z;?) =?P?0 (x, ky, z;?) eiky(y?ys)dky (4.6)We note that there is a specific region of interest along the integration axis of Equation 4.4:Zhou and Greenhalgh (2006) define kc as the critical value of the propagation medium (this pointis also made more obliquely by Song et al., 1995). This marks the threshold between propagatingwaves for values of ky inside the range of ?kc < ky < kc and evanescent waves characterizedby wavenumbers outside this range. The practical application of the transform in Equation 4.4can therefore be concerned chiefly with constructing the equivalent 3D wavefield from thosewavenumbers within ?kc of the zero wavenumber. Since the solutions to the wave equation ofinterest are symmetric about ky = 0, the discretized version may be simplified to a discrete cosinetransform over the reduced space (square brackets indicate discretized quantities):P [?] [x, z] =N?i=0aP?i [?] [x, z] cos (ky,i?y)?ky (4.7)with ky,i = (0, . . . , kc)and a =?????1 if i = 0;2 otherwise.The transform in Equation 4.8 represents a discretized version of the inverse Fourier trans-64form; the sampling of the spectrum in the transform over ky was the topic of work by Zhou andGreenhalgh (2006). In that case, synthetic wavefields were calculated for the governing waveequation with values of ky controlled by the abscissae of the Gauss-Legendre integration formula.We have incorporated this method, which results in somewhat improved numerical stability incomparison to regular sampling, for a constant number of integration components. In this case,the sampling of ky in Equation 4.8 is not regular, and ?ky is controlled by the Gauss-Legendreweights.4.3.3 Full-waveform inversionThe full-waveform inversion stage of the waveform tomography processing is based on themethodin Pratt (1999) and related works, but extends the work of Song et al. (1995) to account for 3Dgeometry (this comes from the forward modeling method discussed above). The inversion is car-ried out using the adjoint state method; i.e., the gradient descent vector for an iterative descentscheme is computed by multiple forward modeling steps (Lailly, 1983; Tarantola, 1984a; Prattet al., 1998). The objective function is based on the logarithmic l2 norm, in which the gradient ispreconditioned by taking the natural logarithm of the data residuals (Shin and Min, 2006; Kameiet al., 2012). The gradient points in the direction of fastest (local) increase of the objective func-tion; for any position in model space that is not a local extremum of the objective function, thesign-reversed gradient represents a local descent operator. There exists a perturbation of somesize along this vector that will reduce the logarithmic l2 norm of the data residuals. The precon-ditioned gradient vector indicates the direction of descent, and a perturbation that is a multipleof this vector is applied. If the data residual is not decreased, the step length is reduced and theprocess is repeated. The conjugate gradient preconditioning method is used to improve the rateof convergence (Polak, 1969), which modifies the gradient for iterations after the first.In order for the non-linear inversion process to converge towards the true earth model, theinitial model information (from tomography or earlier stages of FWI) must predict the traveltimesof a given arrival mode within one half cycle. Since the duration of one cycle is longer at lowfrequencies, the non-linear inversion scheme is more tolerant of model errors when processing65the low-frequency components of the data waveforms. The initial stages of inversion operate ona set of data including only the lowest frequency components. Subsequent stages of the inversioninclude information from higher frequency components (i.e., with reduced band-limiting), whichgradually provide improved resolution and accuracy. This enables a multi-scale approach toinversion, in which coarse model information is resolved prior to inverting for fine structure.Figure 4.1: Overview of geological terranes in central British Columbia. The study area is indicated.A series of geological mapping efforts and geophysical surveys have been carried out throughout thebasin, especially in the central and southeastern portions. Map data from Massey et al. (2005).66Figure 4.2: Surface lithology and the presence of Quaternary cover are shown for the area indicated inFigure 4.1. Much of the study area is thinly covered by the Neogene Chilcotin basalt, which is seldomthicker than 10 m (Bordet et al., 2011). The shallow seismic response is dominated by the Eocenevolcanic rocks of the Ootsa Lake and Endako groups. Stikine Terrane rocks of the post-accretionaryHazelton group outcrop north of the seismic line, and outcrops of Cretaceous sediments are alsovisible east of the line. The geometry of seismic line 10 is shown, with several other lines from theseismic survey indicated for reference purposes. The endpoints of the 2D seismic profile are indicatedalong with source and receiver array extents. Shots 200, 350 and 600 are the locations of the vibrationpoints for which waveform data are shown in Figures 4.7, 4.8 and 4.9, respectively. UTM coordinatesare used, with NAD 27 map datum.4.3.4 2.5D, AVO, source and attenuation considerationsThe discretized version of the Fourier transform (Equation 4.8) is synthesized from a finite num-ber of component wavefields. The limited number of components results in a periodicity of thesource y-coordinate, which manifests as numerical ghost sources. The representation of eachsource along the y-coordinate is a numerical approximation to a spatial delta function, whichtakes the form of a Dirichlet kernel with a shape controlled by the sampling in the discreteFourier transform. The spacing of the ghost sources in the out-of-plane direction is 2pi?ky (Songand Williamson, 1995). Incorporating additional components in the data synthesis increases thisrepeat distance. A similar phenomenon is seen in the synthesis of time-domain data from a seriesof Fourier-component data (Song and Williamson, 1995). In both cases, the solution of choice isto include additional frequency or wavenumber components. However, because of the implicit67periodicity resulting from band-limited Fourier synthesis, some time or space wrap-around ef-fects will still contaminate the data. The application of a frequency domain exponential damping(Sirgue and Pratt, 2003; Shin et al., 2010) selectively amplifies early-arriving waveforms, whiledamping late-arriving waveforms. This substantially attenuates the effects of the ghost sources(and time-periodic effects, equivalently) relative to the main arrivals. The use of this method hasalso been shown to improve convergence by restricting the inversion to fitting the early-arrivingwaveforms (Brenders and Pratt, 2007a), which avoids the later wave coda that may contain wavearrivals with elastic propagation modes. We do not model a free surface, instead extendingthe near-surface low-velocity region with a homogeneous (non-reflective) half-space of the samevelocity and background attenuation above the sources and receivers. Consequently, the down-weighting of free-surface multiples by the exponential damping factor is also beneficial, since wedo not expect to reproduce these features in the synthetic data.While the dominant factors that affect data amplitude variation with offset (AVO) behaviour?viz., geometric spreading effects?are reproduced accurately by 2.5D forward modeling, a num-ber of residual effects remain that may result in discrepancies between the field data and theequivalent synthetic seismogram. These include elastic effects that are not modeled by the vis-coacoustic code (e.g., direct S-waves, mode conversions between P-wave and S-waves), discrep-ancies in the modeled source and receiver radiation patterns, near-surface coupling effects, etc.In order to accommodate the bulk of these AVO effects, we carry out heuristic normalization ofthe data amplitudes to those present in forward-modeled data from an early (smooth) version ofthe model. This follows the method described by Brenders and Pratt (2007a): 1) the real and syn-thetic data are scanned to determine the bulk AVO response, 2) a least-squares linear regressionis applied to the logarithms of the amplitude response vs. absolute offset, and 3) the field data arescaled based on this factor such that their bulk AVO response approximates that of the syntheticdata. In Equation 4.8, the scaled data are computed from the original data by an offset-dependentmultiplicative scale factor. The offset is r. A0 is the intercept and A1 the slope of the log-linearregression line, which describes the bulk amplitude difference between the synthetic and realdata.68dscaled = dorig ? eA0+A1r (4.8)This approach has the benefit of mitigating the large-scale AVO effects of approximation errorswhile retaining small-scale variations in data amplitude that are produced by attenuation or scat-tering. However, the background attenuation model used in forward modeling affects the com-puted amplitudes that go into the AVO correction and consequently may affect the attenuationinformation extracted during inversion.The source signature for the seismic experiment may in some cases be known, but the apparentencoded source is also dependent on the near-surface response in the source region, which maynot be modeled accurately by a discretized finite-difference framework (especially in absenceof an accurate free surface). After the method of Pratt (1999), a source signature is estimatedfrom the data in a linearized fashion by finding the complex coefficients at each frequency (i.e.,amplitude and phase) that lead to the optimum fit between real and synthetic data. This avoidsthe need for a known source signature as input to the inversion, and since a unique series of co-efficients is found for each source, unmodeled near-surface effects on amplitude and phase canbe accommodated to some extent. Other methods for accommodating source and receiver cali-bration can be considered, e.g., a determination of a scaling factor for each source and receiverthat normalizes the amplitude of the data on a per-gather basis. These affect the weighting of thedifferent data elements in the inversion, and have recently been advocated by others in the field.However, additional re-weighting of the data amplitudes on a source and receiver basis was notcarried out in our case study.Other researchers have discussed the limitations of inverting for velocity and attenuation insome detail (e.g., Mulder and Hak, 2009; Hak and Mulder, 2011). For linearized inversion,especially in situations when the aperture of illumination of subsurface features is incomplete, anambiguity exists between the resolution of numerical intrinsic Q and velocity. This ambiguity isreduced when the target is illuminated from many angles (e.g., with near-offset and wide-anglerefraction and reflection wavepaths) and when strong multiple-scattering is modeled correctly,69which produces an equivalent effect. Non-linear inversion enables recovery of accurate high-resolution velocity and attenuation structures, and therefore is less susceptible to this ambiguityat late iterations. In order to mitigate remaining non-uniqueness in this case, we begin with ahomogeneous (fixed) model of numerical Q, and only update this model after the velocity modelhas already been improved using phase-only inversion. The model of attenuation that is invertedfrom data waveforms during full-waveform inversion is also not able to distinguish between effectsof intrinsic (P-wave) Q and coda Q from scattering. This coda Q component is also influenced bymode conversions from P- to S-waves, which are not modeled by the acoustic forward-modelingmethod. Consequently, unmodeled scattering effects and mode conversions may result in loss ofamplitude in the forward-scattered P-wave, which are imaged as high-attenuation anomalies inthe numerical model of inverse Q.4.4 Case study: Nechako BasinWe initially developed this 2.5D waveform tomography methodology for use in processing aseries of seismic reflection datasets acquired along crooked roads in the Nechako Basin, south-central British Columbia, Canada. In a previous work, we applied waveform tomography pro-cessing using the same starting model as this work; however the full-waveform inversion stepwas carried out on static-corrected data waveforms, using a 2D processing method (Smithymanand Clowes, 2012). The methodology described in the previous work improves the tractability ofthe full-waveform inversion step when compared to a methodology that does not account for de-viations from 2D geometry. However, the static-correction approach introduces additional errorsinto the data waveforms, and is only applicable for certain wave arrivals (i.e., those that are well-modeled by the initial traveltime tomography processing). Here we present an improved versionof the result from (Smithyman and Clowes, 2012); the new 2D methodology uses the logarithmicl2 norm rather than the l2 norm. It results in fundamentally different model perturbations that aremore geologically consistent than those seen when the l2 norm was used. This revised 2D resultallows us to make direct comparison with the 2.5D result that uses the logarithmic l2 norm (the70topic of this paper), and we detail the advantages offered by the 2.5D full-waveform inversion(FWI) step.The Nechako Basin is an area of interest for both hydrocarbon and mineral resource explo-ration. Geoscience BC, in conjunction with the British Columbia Department of Energy, Minesand Petroleum Resources, funded and organized a series of exploratory projects involving bothindustry and academia from about 2006 to 2012. As a result, a variety of geophysical and geo-logical mapping efforts have been applied, on both the regional and site-specific scales. Theseinclude detailed geological mapping, potential field surveys and a series of seismic reflectionprofiles along existing roadways (Ferri et al., 2007; Calvert et al., 2009, 2011).4.4.1 GeologyThe Nechako Basin is located in the Intermontane belt of the western Canadian Cordillera. Theregion comprises a series of sedimentary sub-basins. A regional map (Figure 4.1) depicts thelocation of our site of interest, and a map of the local surficial geology (Figure 4.2) provides con-text for this discussion. The underlying basement is composed of Stikine Terrane sedimentaryand volcaniclastic rocks deposited up to the Middle Jurassic. Middle Cretaceous sedimentaryrocks of the Skeena group overlie the Stikine Terrane rocks, followed by Late Cretaceous sedi-mentary rocks of the Sustut group (among others). Both are terrigenous; the former are classifiedas prospective for hydrocarbon development (Hannigan et al., 1994), while the latter are not?both units contain structural traps from compressional tectonics in the mid-Jurassic to Paleoceneand extensional tectonics in the middle to late Eocene. Eocene volcanic rocks of the Ootsa Lakeand Endako groups dominate the near-surface; these were erupted contemporaneously with mid-to late-Eocene extension. Neogene Chilcotin basaltic rocks overlie these in some areas (seldomthicker than 100 m, and often considerably thinner). Quaternary deposits of varying thicknessesand types overlie other units.Smithyman and Clowes (2012) present a combined interpretation of the geology along Line 10,using information gleaned from 2D waveform tomography processing of the same data examinedin this study. The results of that work are additionally compared with an image from prestack71depth migration (PreSDM) along the same line (WesternGeCoMDIC, 2010), and with a result thatincludes interpretations of both a post-stack time migrated seismic section and magnetotelluricdata (Calvert et al., 2011).4.4.2 Survey parametersAn industry-style multichannel seismic reflection survey was carried out in the summer of 2008by CGG Veritas, on behalf of Geoscience BC. This example case deals with 2.5D waveformtomography processing of a subset of data from seismic Line 10. The data comprise records from699 source points, each with 960 live receiver channels arranged in a split-spread configuration(720 receivers along-line to the west of the source, and 240 receivers to the east). The vibroseissource points were spaced 40 m apart, and the receiver-group interval was 20 m. Hence, thenominal along-line receiver offsets extended to approximately 14.4 km. However, due to linecurvature, the straight-line maximum offset was 12 km on average. The vibroseis sweep was 28seconds long (including 0.9 s on- and off-tapers), with frequencies ranging from 8 Hz to 64 Hz.The source consisted of four diversity-stacked sweeps from four vibrators each, distributed over60 m along-line (Calvert et al., 2009). The geometry of the seismic line is shown overlaid onsurficial geology in Figure 4.2.4.4.3 Model-building procedureThe initial model building process for this study is identical to that discussed in Smithyman andClowes (2012), with the exception that none of the steps carried out for the static correction werenecessary for the 2.5D processing (due to the accommodations for 3D geometry at the FWI stage).First-arrival traveltimes were picked for each trace in the 699 source gathers we examined.Picking was carried out on the direct and refracted waves, which possess a bandwidth from 8-16Hz. The data are generally good quality, but some traces were excluded from further processingbecause of the presence of noise. This includes random background noise (e.g., due to wind andother ambient conditions on the line swamping signal at long offsets) and systematic noise dueto logging trucks and equipment that were active during the data acquisition. Since picking was72Amplitude variation with osetFieldSynth0 2 4 6 8 10 12 14101613?5+5+0?10?15Oset (km)ln (A)? ln (A)Figure 4.3: The natural logarithm of data amplitude is shown for the original correlated seismic data(Field in the upper panel) and an early synthetic result (Synth) using a unit-amplitude delta functionsource. Colours represent the density of points (magenta indicates the highest density and red thelowest); the variation in data amplitude for a given offset range is higher far away from the source.The mean data amplitude was calculated for each dataset at a number of bins over the offset range(indicated by x). From this, a best-fit log-linear trend (black line) for each dataset was determined. Thedifference between these trends represents a scaling relation that adjusts the bulk characteristics of thefield data amplitudes to correspond to the unit-amplitude synthetic data (lower panel; A0, intercept;A1, slope).carried out on the full offset range of the data (nominally 14.4 km), we observed some shift in thephase of the waveform associated with the highest-amplitude arrivals. Unlike with the picking ofexplosive-source data, the acausal wavelet from vibroseis correlation makes determining the first-arrival phase difficult at widely-varying offset ranges (due to differences in the correlated waveletbetween near and far offsets, deconvolving with a uniform wavelet is of limited effectiveness).The picking error was between 20 and 30 ms, and was controlled by the ability to distinguish aparticular arrival phase. Additionally, the seismic records did not contain information before timezero, and as such some of the short-offset traces did not possess the full record of the correlated(acausal) source.73FAST (Zelt and Barton, 1998) was applied to generate a 3D model of seismic velocities basedon nonlinear inversion of the traveltime data. The model size was 70 km ? 4 km ? 6 km, with a100 m cell size. The x- and y-coordinates were oriented such that the y-coordinate was transverseto the line of best fit through the plan-view projection of the seismic array. This preserved thegeometry, but ensured that the largest components of the source-receiver offsets were representedby the x-coordinate offset. The model was iteratively updated until a RMS traveltime misfit of 27ms was reached. An asymmetric smoothing operator was applied to the model at each iterationin order to encourage the smoothest model that fit the traveltime data; the smoothing in the y-direction was sufficiently strong to ensure that the model was homogeneous in that coordinate.This also represents a simplification, with an assumption similar to the 2.5D approximation usedin the waveform code: the smoothing along the y-direction is strictly valid only if the acquisitionline crosses geological strike. The traveltime misfit was compatible with the expected pickingerrors in the data; this represented an average misfit of less than one quarter cycle at 8 Hz, whichindicated that the full-waveform inversion stage should be tractable. A slice in x, z coordinates ofdimensions 701 ? 61 was extracted, which represented the average model velocities along they-axis. This slice was sub-sampled to a cell size of 50 m ? 50 m, and cropped for purposes ofFWI to a model size of 901 ? 81 (i.e., dimensions of 45 km ? 4 km).The waveform data were windowed around the first-arrivals, in order to limit the inversion tofitting mainly the refractions. This had the added benefit of eliminating most contributions fromelastic propagation modes (S-waves and mode-converted arrivals), which generally arrive laterin the seismic record. However, loss of amplitude in the P-wave arrivals due to elastic modeconversion is still present in the refractions. To avoid contamination of the result by ground-rolland the air-wave, the data for receivers with offsets smaller than 1 km were not included in theinversion.In order to compare the results of 2D and 2.5D full-waveform inversion, we undertook areprocessing of the results from Smithyman and Clowes (2012) using modified preconditioningparameters (viz., logarithmic l2 norm and adjusted smoothing parameters). These parameterswere chosen to generate a 2D result with a workflow as close as possible to the 2.5D workflow74that is the topic of this study. Full-waveform inversion was carried out to improve the resolutionand fidelity of the model of seismic velocity, in each of the 2D and 2.5D cases. The data werepreconditioned using an exponential damping factor (? = 0.8 s) to emphasize early arrivals. Inearly iterations we used a small portion of the seismic data, band-limited within a few Hz of the 8Hz minimum frequency. We mainly attempted to fit the waveform phase, and did not attempt tosolve for complex velocity (and therefore attenuation, parameterized as 1/Q) in the initial stagesof inversion. This is important, because both velocity and attenuation models affect the amplitudeof arrivals, but the phase of the forward-scattered waveforms is mainly controlled by the velocitymodel (Brenders and Pratt, 2007a; Smithyman et al., 2009). A fixed model of Q = 200 was setfor the initial stages of inversion.In the final stages of 2.5D FWI we fit both waveform phase and amplitude, and additionallyallowed the gradient updates to affect the full complex c, and thereby to update the attenuationmodel. For the later stages of inversion including amplitude information, the data were scaledusing a log-linear correction to fit the bulk amplitude variation with offset as determined by anearlier stage of inversion (Figure 4.3; see also Section 4.3.4). This normalized the data amplitudesto the unit-amplitude delta-function source employed at early stages of inversion. We appliedFWI in several stages. In all cases, frequencies were accumulated in the inversion for maximumstability, and the gradient was formed from the frequency-domain data at 0.25 Hz intervals. Theintegration over y-wavenumbers (only present in the 2.5D case) contained 40 components, from0 to kc (see Equation 4.8), selected independently for each temporal frequency component usingthe Gauss-Legendre method of Zhou and Greenhalgh (2006).The following steps were carried out:1. 10 iterations using a constant model for Q = 200, including 5 frequencies from 8 to 9 Hz2. 5 iterations with Q = 200, including 13 frequencies from 8 to 11 Hz3. 5 iterations with Q = 200, including 21 frequencies from 8 to 13 Hz4. 8 iterations with inversion for c and Q, including 33 frequencies from 8 to 16 Hz. Updates mutedoutside 10 km < x < 35 km to preserve stability.We chose to accumulate frequencies (i.e., invert additional frequencies in combination with the750123Depth (km)0123Depth (km)0123Depth (km)0123Depth (km)0123Depth (km)0123Depth (km)1500700051673333Velocity (m/s)00.0670.0440.022Inverse Qa)b)c)d)e)f )Figure 4.4: (a) The starting velocity model from traveltime inversion. (b) The velocity model from 2Dwaveform tomography using the l2 norm in Smithyman and Clowes (2012). (c) The velocity modelfrom 2D waveform tomography using the logarithmic l2 norm. (d) The velocity model from 2.5Dwaveform tomography using the logarithmic l2 norm after 20 iterations (for direct comparison withc). The detail recovered is higher in the 2.5D result. (e) The velocity model from 2.5D waveformtomography after 28 iterations. Low-velocity features are identified more reliably than in (a), and depthof investigation is substantially better than with traveltime inversion alone. (f) Attenuation model from2.5D waveform tomography (accompanies e). The model updates were zeroed outside of 7.5 km <x < 37.5 km; background Q = 200. All panels are muted to indicate the portion of the modelinterrogated by the raypaths in the initial model; recovery of features is most reliable in this region.frequencies of the previous steps) in the inversion process, for purposes of stability and to reducethe number of potential factors contributing to modeling errors while assessing the capabilities ofthe 2.5D method; however, in principle significantly fewer simultaneous frequencies should berequired as long as appropriate damping and regularization are applied (e.g., using the methodof efficient waveform inversion described by Sirgue and Pratt, 2004).760123Depth (km)0123Depth (km)a)b)c)0123Depth (km)0123Depth (km)0123Depth (km)d)e)-100010000Velocity Di!erence (m/s)Figure 4.5: (a) Velocity difference from 2D FWI in Smithyman and Clowes (2012) with l2 norm (i.e.,Figure 4.4 b minus a). (b) Velocity difference from new 2D FWI in this paper with logarithmic l2norm (i.e., Figure 4.4 c minus a). (c) Velocity difference from 2.5D FWI (i.e., Figure 4.4 e minus a).(d) Checkerboard testing result for 1% perturbation on velocity in 2.5D FWI result. (e) Checkerboardtesting result for 1% perturbation on Q in 2.5D FWI result.4.4.4 ResultsThe initial model from FAST processing of the first-arrival data is shown in Figure 4.4 (a). Thismodel image is darkened in the portions not sampled by the raytracing algorithm used in trav-eltime inversion; this represents part of the null space of the traveltime inversion. The velocitymodel resulting from 2D (l2-norm) waveform tomography processing by Smithyman and Clowes(2012) is reproduced in Figure 4.4 (b). In order to facilitate comparison under the same normas the new 2.5D FWI result, a result from 2D full-waveform inversion of static-corrected wave-form data by reprocessing of the results of Smithyman and Clowes (2012) using the logarithmicl2 norm is presented in Figure 4.4 (c). The new result uses processing parameters equivalent tothe 2.5D workflow, and results in FWI model perturbations that are quite different from thosediscussed by Smithyman and Clowes (2012). The earlier 2D waveform tomography study was77based on similar processing methods to this study, but used a preconditioned conjugate gradientalgorithm to minimize the l2 norm of the waveform data residuals. In both 2D results, large dataerrors existed due to unresolved geometry errors in the western part of the seismic line. In addi-tion, noise was present in the survey data from nearby industry and environmental effects, andwas similarly strongest in the western part of the survey. The l2 norm is sensitive to data outliers,and these errors dominated the updates to the model in the previous result due to the l2 normused. The pattern of model updates from the previous study is visible in Figure 4.5 (a), and theareas with the strongest updates also correspond to the parts of the line with the largest residualgeometry errors due to deviations from a straight line (cf. line geometry in Figure 4.2). Additionalwork using the logarithmic l2 norm has yielded FWI updates that are more stable in the new 2Dresult (Figures 4.4 c and 4.5 b) and this suggests that the FWI updates to the model presented bySmithyman and Clowes (2012) were not imaging geological structure. The data fit in the easternpart of the line in the previous study was comparatively good when compared to the western part.Hence, the assessment of the error under the l2 norm meant that the updates from traces with anominally acceptable (but not excellent) data fit were not considered (or improved upon) in thepresence of larger errors due to limitations in the geometry static correction. The sensitivity of thel2 norm to data outliers (e.g., noise bursts or modeling errors) has been documented in researchby several workers, who have argued for a variety of treatments to mitigate the problem. Someof these approaches make use of alternate norms, such as the l1 norm (Brossier et al., 2009b), thelogarithmic l2 norm (Shin and Min, 2006), and others (Brossier et al., 2010; Aravkin et al., 2011)that reduce sensitivity to data outliers. Another approach involves processing of trace-normalizedseismic data under the l2 norm and considers only the phase residuals (e.g., Operto et al., 2004;Bleibinhaus and Hilberg, 2012), which may yield significant benefits for stability because of re-duced sensitivy to noise. The use of the logarithmic l2 norm (Shin and Min, 2006) reduces theeffect of large data residuals on the gradient update, and tends to increase the strength of theperturbations away from the source and receiver arrays (acting as a preconditioning operation onthe gradient), which improved the stability of the 2D and 2.5D full-waveform inversions in thepresent study. However, the difficulties associated with modeling out-of-plane geometry remain78for the 2D approach and preclude useful model updates at frequencies above about 13 Hz. Figure4.4 (d) shows an intermediate velocity model after 20 iterations of 2.5D FWI, for direct compar-ison with the result from 2D FWI using the logarithmic l2 norm (Figure 4.4 c). The final modelsof velocity and attenuation resulting from the new 2.5D waveform tomography processing arepresented in Figure 4.4 (e) and (f). This new result contains information from inversion of the fullbandwidth of the early-arriving (forward-scattered) waveforms (i.e., 8?16 Hz) and does not sufferfrom the inaccuracies introduced by a static correction using refraction traveltimes in the previ-ous 2D methodology. Note, however, that much of the structure produced in the velocity modelfrom 2.5D waveform tomography is developed in the first 20 iterations (which are comparable tothe 20 iterations carried out in the 2D logarithmic l2 norm processing method); compare Figure4.4 (c) and (d). The contribution to the velocity model from the band of frequencies between 13and 16 Hz is minor in the 2.5D result (see Figure 4.11 and surrounding discussion). The mainfeatures of the velocity model from 2.5D FWI are substantially established in the first 20 iterations(compare Figure 4.4 d and e, which are visually indistinguishable). However, the final 8 iterationsof the 2.5D FWI process constrain the model of attenuation.Figure 4.5 (a) shows the difference between the final velocity model from 2D full-waveforminversion in the result of Smithyman and Clowes (2012) (reproduced in Figure 4.4 b) and thestarting velocity model from traveltime inversion (Figure 4.4 a). Equivalently, Figure 4.5 (b) showsthe updates from the re-processing of of the 2D workflow using the logarithmic l2 norm (discussedabove). The FWI model update from the 2.5D result is shown in Figure 4.5 (c), i.e., the differencebetween Figure 4.4 (e) and Figure 4.4 (a). The increased stability from fitting the 3D geometryand from the improved accuracy in amplitude modeling allows for improved confidence in theimaging throughout the depth range. However, many of the features in Figure 4.4 (a?e) aresubstantially similar between the 2D log l2 norm and 2.5D results (or, for that matter, betweenthe two full-waveform inversion results and the starting model), since each stage of the inversionyields an incremental improvement in resolution and accuracy. The most significant changesin the character of the model updates (compared to the result of Smithyman and Clowes, 2012)come from the change in data norms. In several locations, the 2.5D full-waveform inversion code79smoothed or restructured features that were present in the initial model and were left unmodifiedin the 2D model, e.g., a raypath artifact visible at 1?2 km depth and approximately 23?25 kmlateral position (Figure 4.4). There are also some indications of ringing features in the 2D result(e.g., Figure 4.4 c, 12 km lateral position) that are not present in the 2.5D result or the startingmodel, though the same is true to a lesser extent in the 2.5D result (e.g., Figure 4.4 e, 9 kmlateral position). In both cases, the effect is most noticeable on the flanks of the model where theangles of target illumination are incomplete (see e.g., Figure 4.5 b and c). The areas in which theupdates from full-waveform inversion are most complex can be readily determined by examiningFigure 4.5 (b and c). In the case of the 2.5D inversion, several locations seem to show recoveryof complex structures with relatively high wavenumbers (e.g., 12 km along-line, shallow; 20 kmalong-line, middle depths; 25 km along-line, shallow).To estimate the interpretable resolution of the model, we carried out checkerboard perturba-tion testing on both the velocity and attenuation for the 2.5D FWI result (shown in Figure 4.5e and f). From checkerboard testing results it is possible to determine the sensitivity of the seis-mic experiment to perturbations on the final velocity and Q models. The velocity model andQ model were multiplicatively perturbed by a pattern of rectangular checkers, of amplitude 0.01(i.e., 1%). Adjacent checkers possessed opposite polarity, and each checker was 1000 m in widthand 500 m in height. Each parameter (velocity and Q) was perturbed independently, while hold-ing the other parameter constant. New (noise-free) data waveforms from the perturbed modelwere synthesized using the same forward modeling code as that used in the FWI. These syntheticwaveforms were processed using the same parameters as the final iterations of the real-data FWIstudy; however, the use of noise-free viscoacoustic synthetic data means that the resolution es-timate is optimistic, compared to the result that would be expected from noisy data recorded inthe field. Other limitations include the lack of a free surface, potential out-of-plane scattering,and unmodeled scattering from elastic effects; all of these will also impact practical resolution,and are not modeled in this resolution test. Checkerboard tests estimate the modeling resolutionbased on the recovered earth model, but do not test the validity of that earth model (and hencecannot be used as ground truth in a real data case study such as this). Examination of Figure804.5 indicates that the main perturbations in the seismic velocity and attenuation models occurredbetween lateral positions 10 km and 35 km. Because of instabilities in the formation of the gra-dient at 16 Hz along the flanks of the study (due to incomplete illumination from the seismicarray), only a subset of the model was updated in the last stages of the 2.5D real-data inversion.The ambiguities between recovery of velocity and attenuation models are also strongest in theflanks, where the angles of target illumination are limited. The checkerboard pattern is recoveredin the region of good illumination, but the updates were muted outside this region for the samereasons as in the field-data inversion. Resolution is highest in the shallowest 1 km to 1.5 km ofthe model (relative to the source/receiver array), and between about lateral positions 23 km and33 km. Outside this region, striped features indicate an imperfect recovery of small perturbations,with the effect of anisotropic smoothing of the images of the checkers. The checkerboard patternis recovered best in regions of the model with a strong velocity gradient. This is most-likely dueto improved reflection amplitude and signal-to-noise ratio, combined with shorter propagationdistances for the refracted modes that are important for determining the velocity and attenuation.In discussing resolution, it is also informative to consider the portions of the model to whicheach waveform component is most sensitive. Figure 4.6 presents numerical sensitivity kernelscomputed for single pairs of sources and receivers with an offset approximately half that of themaximum offset of the seismic survey. These kernels are presented for 8 Hz and 16 Hz frequen-cies, which represent the lowest and highest usable frequencies present in the forward-scatteredfield data. In most cases, the anomalies resolved by full-waveform inversion will be on the orderof the size of the sensitivity kernel, or some fraction thereof. However, the stacking of multi-ple kernels (for the set of all sources and receivers) can lead to high model resolution, espe-cially with multiple iterations of nonlinear inversion. The lowest data frequencies (close to 8 Hz)provide a resolution similar to traveltime tomography, but with improved physics that leads tolarge-wavenumber model updates at the earlier inversion stages. The highest data frequenciesinterrogate mainly the near-surface region, due to the rapid decay of amplitude with regards toboth offset and depth; however, the model updates resulting from the high-frequency portion ofthe gradient contain significantly greater information about fine structure than the low-frequency810123Depth (km)0123Depth (km)0123Depth (km)0123Depth (km)a)b)c)d)7.2 kmFigure 4.6: Sensitivity kernels are displayed for pairs of sources and receivers 7.2 km apart (i.e., at50% of maximum offset), calculated in the final model from 2.5D waveform tomography (Figure4.4 e and f). Kernels are shown for the 8 Hz wavefield (a and b) and the 16 Hz wavefield (c andd) corresponding to unit-valued point sources and receivers. Blue represents the positive portion ofthe kernel; the central blue band is analogous to the first Fresnel zone. Red represents the negative(anti-correlative) portion of the kernel. The wavefields displayed are tapered near the model edgesusing the same parameters as the gradient preconditioning during inversion; however, no additionalpreconditioning is applied. Hence, high-wavenumber anomalies can be seen in these images that aretypically removed by low-pass wavenumber filtering.wavefields. Additionally, the use of the log l2 norm can be seen as a preconditioner on the gra-dient (and therefore on the sensitivity kernels), which increases the strength of updates in thedeepest parts of the model at all inversion frequencies.In order to assess the quality of the inversion result, we compare the original (field) seismicdata and synthetic data generated using the initial velocity (Figure 4.4 a), final 2D logarithmic l2norm velocity (Figure 4.4 c), and final 2.5D velocity and attenuation models (Figures 4.4 e and f)for sources 200, 350 and 600 (Figures 4.7, 4.8 and 4.9 respectively). In each figure, the field dataand synthetic data from the 2.5D result are presented for each trace (Figures 4.7, 4.8 and 4.9; a andb in each figure). The latter three panels of each of Figures 4.7, 4.8 and 4.9 show overlay (wiggle)plots for the three synthetic datasets compared with the field data; however, these plots display82only a subset of traces to enable readability. These results are presented with trace-normalizedamplitudes and in reduced time with a reduction velocity of 4500 m/s. Both real and syntheticdata were windowed to 1000 ms records from 200 ms before the traveltime pick to 800 ms afterthe pick. As a result of the vibroseis correlation process, the source signature is acausal and partof the seismic record occurs before the nominal first-arrival traveltime (see technical discussion,Smithyman and Clowes, 2012). Real data were low-pass filtered to remove information above16 Hz prior to inversion (to attenuate noise and avoid aliasing; the first-arrival bandwidth wasnaturally limited to less than 16 Hz outside of very short offsets), and these band-limited data areshown in Figures 4.7 through 4.9 for comparison with synthetic data. The 2D synthetic data (inpanel d of each respective figure) are presented with a reversal of the static correction, which wasimplemented on the field data for the FWI, to allow for direct comparison with the uncorrecteddata waveforms.Figure 4.7 shows a series of source gathers corresponding to source 200, the location of whichis indicated in Figure 4.2. In general, the phase of the real-data arrivals is predicted more accu-rately by the result from 2.5D FWI (Figure 4.7 e) than in either the starting model (Figure 4.7 c) orthe model from 2D logarithmic l2 norm FWI of static-corrected data (Figure 4.7 d). In several lo-cations (labelled A), there are features on the first-arriving waveforms that are not well replicatedby the initial or 2D data but are closely recovered in the result from 2.5D processing. However,the result from 2D FWI also shows substantial improvement over the result from the initial model,and the fit is quite similar to that from 2.5D FWI at a number of points on the first arrival (e.g., la-bel B). In some locations (labelled C), there are strong diffractions present in the field data that arepartially recovered by 2.5D FWI but not in the initial model or result from 2D FWI. These featuresare probably recovered in part during the final 8 iterations of 2.5D FWI, when the bandwidth ofthe inversion is maximal and attenuation is imaged.Figure 4.8 shows data corresponding to source 350 (see Figure 4.2 for location). The fielddata (Figure 4.8 a) are quite complex in this region, as the westward-travelling waves mainlyinterrogate the (interpreted) basin structures (viz., the synformal features from about 12 km to 26km lateral position). The data waveforms are not recovered accurately in either the initial model83Field Data2.5D Synthetic DataOverlay: Initial model dataOverlay: 2D synthetic dataOverlay: 2.5D synthetic data01250Traveltime (ms)01250Traveltime (ms)01250Traveltime (ms)01250Traveltime (ms)01250Traveltime (ms)0 4.7 km E14.4 km W Source?Receiver O!setC CCAA B BBa)b)c)d)e)Figure 4.7: A collection of source gathers from source 200 of the Nechako Basin seismic survey. Dataare presented time-reduced, with a reduction velocity of 4500 m/s. (a) Field data. Letters refer to fea-tures discussed in the text. (b) Synthetic data generated from the final result from 2.5D waveformtomography. (c?e) Synthetic data (green) overlaid on field data (black); field data are subsampled toevery 5th channel and synthetic data are subsampled to every 10th channel. (c) Synthetic data over-lay generated from the initial model, using 2.5D modeling. (d) Synthetic data overlay generated frommodel resulting from 2D logarithmic l2 norm waveform tomography of static-corrected data wave-forms; modeled in 2D, with appropriate reversal of static correction for comparison with unshiftedfield data. (e) Synthetic data overlay generated from the final result from 2.5D waveform tomography.or in the result from 2D logarithmic l2 norm FWI (Figure 4.8 c and d respectively). However, the2D FWI result improves on the fit from the starting model over most of the offset range towardsthe west, at least for waveforms closely following the first arrival (those labelled D). The resultfrom 2.5D FWI does somewhat better (Figure 4.8 e), and reproduces several secondary arrivals84Field Data2.5D Synthetic DataOverlay: Initial model dataOverlay: 2D synthetic dataOverlay: 2.5D synthetic data01250Traveltime (ms)01250Traveltime (ms)01250Traveltime (ms)01250Traveltime (ms)01250Traveltime (ms)0 4.7 km E14.4 km W Source?Receiver O!setDE EFGa)b)c)d)e)Figure 4.8: A collection of source gathers from source 350 of the Nechako Basin seismic survey. Dataare presented time-reduced, with a reduction velocity of 4500 m/s. (a) Field data. Letters refer to fea-tures discussed in the text. (b) Synthetic data generated from the final result from 2.5D waveformtomography. (c?e) Synthetic data (green) overlaid on field data (black); field data are subsampled toevery 5th channel and synthetic data are subsampled to every 10th channel. (c) Synthetic data over-lay generated from the initial model, using 2.5D modeling. (d) Synthetic data overlay generated frommodel resulting from 2D logarithmic l2 norm waveform tomography of static-corrected data wave-forms; modeled in 2D, with appropriate reversal of static correction for comparison with unshiftedfield data. (e) Synthetic data overlay generated from the final result from 2.5D waveform tomography.(E) as well as generally fitting the first arrival with improved accuracy. A complex feature on thefirst-arrival (F) is only fit in the 2.5D result. In all cases, some phase shift is present on the easternend of the gather, though less than a cycle in the case of the 2.5D result, which fits the waveinterference from multiple arrivals (G) to some extent.85Field Data2.5D Synthetic DataOverlay: Initial model dataOverlay: 2D synthetic dataOverlay: 2.5D synthetic data01250Traveltime (ms)01250Traveltime (ms)01250Traveltime (ms)01250Traveltime (ms)01250Traveltime (ms)0 4.7 km E14.4 km W Source?Receiver O!setHIa)b)c)d)e)Figure 4.9: A collection of source gathers from source 600 of the Nechako Basin seismic survey. Dataare presented time-reduced, with a reduction velocity of 4500 m/s. (a) Field data. Letters refer to fea-tures discussed in the text. (b) Synthetic data generated from the final result from 2.5D waveformtomography. (c?e) Synthetic data (green) overlaid on field data (black); field data are subsampled toevery 5th channel and synthetic data are subsampled to every 10th channel. (c) Synthetic data over-lay generated from the initial model, using 2.5D modeling. (d) Synthetic data overlay generated frommodel resulting from 2D logarithmic l2 norm waveform tomography of static-corrected data wave-forms; modeled in 2D, with appropriate reversal of static correction for comparison with unshiftedfield data. (e) Synthetic data overlay generated from the final result from 2.5D waveform tomography.Figure 4.9 shows data corresponding to source 600 (see Figure 4.2 for location). Overall, thedata fit is quite good for synthetic data generated in both the 2D (Figure 4.9 d) and 2.5D (Figure4.9 e) models. The eastern-travelling waves are predicted accurately from the initial model, andare improved upon in both FWI results. For western-travelling waves, a partial cycle-skip is visible86at about 4.5 km offset (H). This indicates that a region of the model is recovered with a velocitythat is slightly too high (corresponding to the deeper portions of the high-velocity region visible at20 km lateral position in Figure 4.4 a?e, where some numerical artifacts are visible at about 2 kmdepth). The first-arrival is modeled more accurately from 2.5D FWI than from the initial model or2D FWI throughout the longer-offset portion of the gather (about 5.5 km onward beginning at I),which is representative of different propagation paths that do not interrogate the region affectedby numerical artifacts.The data waveforms are predicted more accurately in part because of the incorporation of thefull 3D geometry information by 2.5D waveform tomography; in several cases there are notice-able perturbations in the shape of the traveltime curve which are not accurately modeled with a2D projected geometry (see Smithyman and Clowes, 2012 for discussion of traveltime data andexample gathers without reversal of the static correction). Waveform tomography also benefitsfrom fitting secondary arrivals that are not visible to traveltime-inversion methods. Low-moveoutemergent arrivals from deeper layers, which are modeled more accurately by the waveform to-mography process, are observed at several points in Figures 4.7 through 4.9. The emergent wave-forms are representative of wavepaths that pass through comparatively fast earth materials, andcontrol the first-arrival traveltimes at offsets greater than the crossover point. The expression ofthese waveforms at offsets shorter than the crossover point is exemplary of a phenomenon that isnot considered by traveltime inversion, but is fit by the full-waveform inversion process.Comparison of time-domain waveforms alone may not be sufficient to assess the success ofthe non-linear waveform inversion step, as the synthetic data waveforms may appear visuallyquite similar to the real data waveforms while still being generated in a flawed model (e.g., Pratt,2008). Also, each shot gather shows only a small proportion of the information in the full dataset.In order to develop a more complete view of the model quality, it is beneficial to consider boththe frequency-domain data fit and the proportion of the model illuminated by each frequencygroup. Figure 4.10 shows the amplitude and phase at 12 Hz (i.e., middle bandwidth) for thepreconditioned real data and the synthetic data generated in the final 2.5D FWI model (i.e., Figure4.4 e and f). The sensitivity can be inferred from Figure 4.6, which shows that the deepest parts87?30?12?18?24ln Amplitude?pi+pi0Phase1350699Source1350699Source1350699SourceTraveltime (ms)1 500 1000 1500 2000 23621350699SourceField Data Amp.Synth Data Amp.Field Data PhaseSynth Data Phasea)b)c)d)Figure 4.10: Frequency-domain panels showing amplitude and phase for each source/receiver pair at12 Hz. The amplitude is shown for each datum in the a) AVO-corrected field data and b) syntheticdata from the final 2.5D result. Amplitudes are displayed on a logarithmic scale. The phase is shownin radians for each datum in the c) field data and d) synthetic data.of the model are not well-constrained by the high-frequency parts of the dataset, whereas theshallowest parts of the model are interrogated by waveforms with a broader bandwidth.To assess the relative efficacy of the 2D and 2.5D waveform tomography methods, the reduc-tion in each respective objective function is presented in Table 4.1. The objective function ineach case was based on the log l2 norm and values are tabulated for several sets of iterations.Since the 2D method was only stable at lower frequencies, the last stage of inversion including88Iter 01?10 Iter 11?15 Iter 16?20 Iter 21?282D Result 29.6% 12.2% 8.0%2.5D Result 24.3% 7.7% 6.3% 19.9%Table 4.1: Relative reduction of objective function for the 2D and 2.5D inversion stages. As the erroris measured under the log l2 norm, the numerical size of the (dimensionless) misfit reduction is smallwhen compared to errors calculated using the l2 norm.bandwidth from 8 to 16 Hz and allowing updates to the model of attenuation was only carriedout in the 2.5D waveform tomography workflow. In all cases, the percent reduction in the objec-tive function was higher in the case of 2D full-waveform inversion processing of static-correctedfield data. However, since the data were static-corrected before 2D FWI, the absolute values ofthe data residuals are not comparable between the two methods. While a comparison betweenthe absolute errors is not meaningful in this case (due to the differences in the field data exposedto each workflow), a numerical measure of the relative model quality is desirable. The velocitymodel result after 20-iterations was recorded for each approach (Figure 4.4 c for 2D; Figure 4.4 dfor 2.5D). Synthetic data were calculated after forward modeling in 2.5D for each case (i.e., usingidentical numerical methods) and compared to the uncorrected field data. The model from 2Dprocessing produced a data residual 12.1% higher than the model from 2.5D processing, underidentical conditions. In the final 8 iterations that were not included in the 2D workflow, the 2.5Dwaveform tomography workflow further reduced the misfit by 19.9%.Figure 4.11 shows the data residual at various frequencies and iterations in more detail. Figure4.11 a summarizes the character of the reduction in data fit in each group of iterations. For thefirst 20 iterations of 2.5D FWI, only the phase residual was considered (and hence, the misfit inthe amplitude portion of the logarithmic l2 norm is not significant in the model updates). In thefinal 8 iterations, the norm is the linear combination of the (unweighted) phase and amplituderesiduals. Figure 4.11 a shows that while the phase residual increases modestly in the final 8iterations, the fit to the data amplitude is significantly improved and the corresponding amplitudeand total residuals are reduced. This trade-off is not unexpected, and a comparison betweenFigure 4.4 (d) and (e) shows that the changes to the velocity model are not significant. Notethat the misfit reduction shown in Figure 4.11 (a) is representative of the reduction within each89A B C DE2.5e60.0e6Phase ResidualIteration Number010201528Frequency (Hz)16131198DE2.5e60.0e6Amp. ResidualIteration Number010201528Frequency (Hz)16131198ED3.0e60.0e6Total ResidualIteration Number010201528Frequency (Hz)161311980.00e71.00e71.00e72.00e72.15e72.15e70.00e78.00e7Combined Data ResidualB 10?15A 00?10 C 15?20 D 20?28PhaseAmplitudeTotala) b)c) d)Figure 4.11: The data residual is shown for all iterations (0?28) and all frequencies (8-16 Hz). Thecombined residual (a) is composed of the sum of all the residuals for the frequencies processed ata given iteration. Each individual band of iterations (in a) comprises a different number of discretefrequencies (see regions in b for a representation of the number of summed frequencies), and the itera-tions with wider bandwidth therefore minimize a larger total objective function. Phase (b), amplitude(c) and total (d) residuals are illustrated. Five separate regions are indicated, labelled A through E.Region A corresponds to iterations 1?10 and comprises 5 frequencies from 8?9 Hz. Region B cor-responds to iterations 11?15 and comprises 13 frequencies from 8?11 Hz. Region C corresponds toiterations 16?20 and comprises 21 frequencies from 8?13 Hz. Regions A through C (iterations 1?20)consider only the phase residual. Region D corresponds to iterations 21?28 and comprises 33 fre-quencies from 8?16 Hz. Iterations 21?28 update the model based on the magnitude of the combinedamplitude and phase residuals (i.e., the total residual). Region E is presented to illustrate the effectsof low-frequency inversion stages on the data fit at higher frequencies, and the effects of phase-onlyinversion on the amplitude fit; however, no information from Region E is considered by FWI.90band of frequencies, but that the first iteration in each band produced a large reduction that isnot shown in the jump between bands (Table 4.1). Figure 4.11 (b) shows the phase residualcalculated at each frequency from 8?16 Hz in 0.25 Hz increments, for each iteration. However,only the regions labelled A?D are significant, and region E is presented only to provide furtherinsight into the behaviour of the wide-band misfit function. The misfit is mainly decreased inthe band of frequencies considered at each iteration, although some improvement is seen withinperhaps 1 Hz above the top inversion frequency (Figure 4.11 b, regions A?C). The amplituderesidual (Figure 4.11 c) is only considered in the final 8 iterations of inversion, and the decrease issignificant within region D. The initial 20 iterations produce no significant update in the amplituderesidual. The total residual (Figure 4.11 d) is actively considered only in the final 8 iterations ofinversion (it comprises the contributions of both phase and amplitude, weighted equally). Hence,the behaviour in region E is due only to the reduction of the phase residual in the first 20 iterationsof inversion. The front edge of region D at iteration 28 in Figure 4.11 (d) gives a rough guideto the data fit with respect to frequency in the final result. The contribution from the highestfrequency group (13?16 Hz) is minor, and almost nonexistent in the phase residual (Figure 4.11b) although the increased bandwidth helps in constraining the amplitude information needed inthe reconstruction of the attenuation model. A slightly elevated residual in the lowest ~1 Hzrange (8?9 Hz) is expected and is consistent with the lower signal to noise ratio in that part of thedataset.4.4.5 InterpretationFigure 4.12 shows the velocity and attenuation models resulting from 2.5D waveform tomogra-phy (reproduced from Figure 4.4 e and f with a 3? vertical exaggeration), with our interpretationoverlaid. The labels indicate the basins and faulted blocks that may exist in the region as a resultof Eocene extension and earlier compression. The attenuation model from 2.5D full-waveforminversion (Figure 4.12 b) provides information that is not available from previous geophysicalresearch in this study area. A number of other geophysical methods have been applied alongthe same profile as this work. Comparison with the prestack depth migration image of West-911500700051673333Velocity (m/s)0123Depth (km)0 5 10 15 20 25 30 35 40 45Projected line location (km)00.0670.0440.022Inverse Q0123 VE 3:1Figure 4.12: Velocity and attenuation results, reproduced from Figure 4.4 (d) and (e), here with 3:1vertical exaggeration. Interpreted faults and unit contacts are indicated, along with suggested rocktypes. A pull-apart basin is interpreted, filled by approximately 1.5 km of Cretaceous sedimentaryrocks, and bounded to the east and west by Hazelton group rocks of the Stikine Terrane. The shal-lowest 0.5?1.5 km appears to be dominated by Eocene volcanic and sedimentary rocks, which showlow velocity and attenuation in comparison to the older units.ernGeCo MDIC (2010) and the magnetotelluric inversion model of Spratt and Craven (2010) isinformative (see also Smithyman and Clowes (2012)), but our interpretation is based predomi-nantly on the geology and the models resulting from this study. Recent petrophysical analysessummarized by Kushnir et al. (2012) provide information about the expected physical propertiesof the rocks in the region, although we are cognizant that the in-situ samples may vary widelyfrom laboratory results.The shallowest 1 km of the model shows seismic velocities ranging between 3000 and 4000m/s and relatively low apparent attenuation (Figure 4.12). Based on these comparatively lowseismic velocities, this region most likely represents a mixture of Eocene volcanics and probablysome sedimentary rocks (e.g., the Ootsa Lake group), consistent with the results of Kushnir et al.92(2012). Although Chilcotin group basaltic rocks overlie much of the study area (Figure 4.2), otherworkers (Bordet et al., 2011) have indicated that these are comparatively thin, on the order oftens of meters. Accordingly, we consider it unlikely that the Chilcotin group rocks contributemeasurably to the seismic response.A sharp contrast in seismic velocity and attenuation is observed between the Eocene volcanicsand the underlying blocks that show rugged paleotopography (Figure 4.12). Below the volcanics,the seismic velocity increases to greater than 5000 m/s, which is diagnostic of either intact sedi-mentary rocks (e.g., the Cretaceous Skeena group or Sustut group rocks) or underlying basementrocks. The velocity gradient is particularly high from about 25?32 km along the line (Figure 4.12),an area which also corresponds closely with the presence of volcanic and volcaniclastic Hazeltongroup rocks of the Stikine Terrane that outcrop north of the seismic line (Figure 4.2). We inferthat Hazelton group rocks most likely extend below the seismic line at a depth of 1 km; however,it is also possible that the seismic response is directly due to these rocks north of the line, sincesome model non-uniqueness exists in the off-line direction.Based on the velocity and attenuation models, several features can be interpreted as faults.For example, the model characteristics at 12 km are interpreted as an extensionally reactivatedthrust fault (arrows at 12 km on Figure 4.12). In most cases, these features show up as low-attenuation anomalies (i.e., high Q), with very high wavenumber content; the features in themodel of attenuation sometimes correspond with areas of low-velocity, but the anomalies inthe velocity model appear less distinct. The coincident change of these properties may indicateunmodeled variability in the model of rock density, which is determined from velocity using theheuristic Gardner?s Relation. From about 13?24 km along the line, the velocities at depths from1 km to 2.5 km are moderate (4800?5500 m/s), values which are likely diagnostic of intact basinsediments. Thus we interpret this area on the cross section as a pull-apart sedimentary basin(Figure 4.12). It appears to be composed of several blocks separated by westward-dipping faultsand is substantially deeper in the western part of this region. To the east, the interpreted basinis bounded by the strong features interpreted to represent rocks of the Hazelton group, as notedpreviously. To the west, the basin is bounded by shallower high-velocity rocks for which the93velocities are not as well constrained as those to the east (due to the limitations of the seismicarray sensitivity in the data subset we used). However, examination of the surficial geology map(Figure 4.2) shows that outcrops of Stikinia are present northwest of the end of our line. Thus, wespeculate that this high-velocity material may also represent Hazelton group rocks.4.5 DiscussionThe use of a 2.5D assumption introduces some simplifications in comparison to 3D, which mustbe taken into account at the processing stage (and ideally at the survey-design stage). The y-coordinate of geometry is handled by phase-shifting individual plane-wave components in theFourier composition that generates the 3D wavefield from component 2D wavefields. However,the assumed support of each (sinusoidal) component is infinite in the y-direction; i.e., the 3Dpoint sources constructed at a given location are only valid in the case when the y-coordinate isparallel to geological strike. Without this criteria being fulfilled, the plane-wave decompositionof the 3D wave equation is invalid and the waveform forward modeling and inversion will notbe accurate. In practical terms, however, an active-source seismic survey only interrogates Earthstructures at some maximum distance from the source and receiver. In the Nechako Basin casestudy, the maximum offset is 14.4 km, the moveout velocity for the main part of the offset rangeis 4500 m/s and the end mute for the data waveforms is 800 ms after the first arrival. Hence,the maximum path length for a scattered wave arrival to be included in the windowed data ison the order of 18 km, and the maximum off-line distance interrogated is perhaps 5.4 km at themidpoint. While this calculation simplifies the geometry and the physics considerably, it is a rea-sonable estimate for the order of magnitude of the effects involved. Hence, the requirement thatthe true Earth be homogeneous in the y-direction can be relaxed to a requirement that the Earthbe smoothly varying in that direction in the locality of the survey, such that the bulk propertiesare close to those in the 2D model and there are no reflections generated. It is perhaps worthnoting that the same limitations apply in 2D modeling and inversion studies?which use onlythe zero-wavenumber plane-wave component?but are seldom mentioned. The most straightfor-94ward means of ensuring these conditions is to select survey configurations that cross geologicalstrike wherever possible. The coordinate system of the 2.5D approach may also be rotated toaccommodate strike at an oblique angle to the cross-line direction, with appropriate choice ofnumerical parameters to ensure that the synthetic wavefields are not aliased in the projected x?zplane. This also illuminates a practical limit to the degree of line curvature that can be accommo-dated in 2.5D processing: switchbacks in the acquisition array may result in wavepaths that arenormal to the finite-difference grid and synthetic wavefields with wavenumbers too high to modelpractically. However, elimination of some data may prove sufficient to mitigate these problems.In this study, the line is oriented perpendicular to some geological features (e.g., outcroppingHazelton-group rocks), but does not necessarily cross geological strike at depth.Whereas 2.5D waveform tomography has been demonstrated to yield improvements over 2Dprocessing, the increased computational cost may still prove prohibitive for cases when rapidturn-around of results is necessary. In offshore studies, when the seismic array may be deployedin a straight line configuration, the array geometry can usually be parameterized in 2D withoutdifficulties; other researchers have shown compelling successes in 2D FWI processing of datacollected in these circumstances. A 2.5D approach also provides improvements in modeling ac-curately the amplitude and phase of a 3D wavefield by simulating point sources; this remains areason to consider 2.5D processing as an improvement over the 2D equivalent. Corrections maybe applied to partially accommodate the problems of 3D geometry and maintain the computa-tional advantages of a 2D processing workflow: e.g., the static correction approach discussedpreviously, or a reprojection of sources and receivers on a per-trace basis to simulate 3D offset.In this case, proper treatment of any unresolved errors is of paramount importance, and may ne-cessitate processing with an alternative to the l2 norm that is typically used. Full 3D waveformtomography processing (with appropriate model constraints to accommodate the data null space)will almost certainly yield results with improved accuracy, especially when the array configura-tion allows for cross-line information to interrogate the 3D earth. However, as noted previously,the computational expense of 3D processing for a crooked-line 2D survey may not be justified.954.6 ConclusionsA method is developed for applying waveform tomography to early-arriving refraction waveformsfrom seismic reflection datasets that possess deviation from an ideal 2D survey configuration; e.g.,a typical crustal seismic profile acquired along existing and crooked roads. Thus, this newmethodis particularly applicable to datasets from academia or small-scale industry surveys that do notmake use of 3D seismic acquisition. The use of 2.5D full-waveform inversion, in combinationwith 3D traveltime inversion, enables the method to account for 3D geometry in the source andreceiver arrays. This is of greatest benefit when the seismic array configuration lends itself to a 2Dinterpretation stage, but the 3D geometry makes conventional 2D full-waveform inversion pro-cessing intractable. In addition, the 2.5D implementation of the full-waveform forward modelingand inversion steps results in improved fidelity in modeling the amplitude and phase effects ofgeometric spreading when compared to a 2D wave-equation method. The primary benefits aregreatly improved velocity models relative to traveltime inversion and generation of attenuationmodels within the region interrogated by the seismic waveforms. These improvements result inmoderate computational overhead, which in practice is typically on the order of 40?50 times thecost of 2D inversion for each iteration but significantly less than application of a full 3D algorithm.A case study from a crustal reflection profile in the Nechako Basin in south-central BritishColumbia, Canada exemplifies and validates the newly developed approach. Models of velocityand viscoacoustic attenuation provide valuable information that can be useful for interpretation,including identification of rock types and structural features. The high resolution exhibited bythese models provides significant improvements with respect to fault delineation and definitionof structural boundaries when compared to conventional tomography techniques. The use of low-frequency refracted waveforms limits the resolution when compared to an industry-style prestackdepth migration workflow; however, the high-resolution model of seismic velocity makes theplacement of features comparatively accurate. As well, the inversion of refraction waveformsmakes it possible to produce a detailed image including near-surface structures that may be in-visible to a conventional stacking or migration workflow.96Chapter 5New Geophysical Models for SubsurfaceVelocity Structure in theNechako?Chilcotin Plateau from 2.5DWaveform Tomography5.1 IntroductionSeismic inversion is applied to generate physical property models (P-wave velocity and numericalattenuation) for several profiles in the Nechako?Chilcotin plateau region of south-central BritishColumbia, Canada. These models are useful for the characterization of rock types in terms oftheir positions and thicknesses, which may be used in conjunction with geological ground truthto infer the extent of lithostratigraphic units in the subsurface. The velocity structure also maybe used for future reprocessing of the seismic reflection data to derive improved images basedon the better near-surface velocity models. These velocity and attenuation results are mainlyrelevant within the top 2?3 km of the subsurface because they are limited by the configuration ofthe seismic array used to acquire the 2008 Geoscience BC seismic data and the strong velocitycontrasts between volcanic and sedimentary rocks in the region.Interest in the Mesozoic (mainly Cretaceous) sedimentary units of the Nechako?Chilcotinplateau has been economically motivated by the potential for hydrocarbon development (Hayes97et al., 2004;Hannigan et al., 1994). Previous works most commonly refer to this study area as theNechako Basin, though this term is defined ambiguously (e.g., Riddell, 2011). As discussed by oth-ers (Calvert et al., 2009; Hayward and Calvert, 2009; Riddell, 2011; Calvert et al., 2011; Talingaand Calvert, 2012), a primary goal of recent (2004?2013) geological and geophysical work inthe area of the Nechako?Chilcotin plateau has been the delineation of the extent of the Ceno-zoic volcanic units. These units have historically presented significant challenges to geophysicalinvestigation of the Mesozoic sedimentary basins in the area because the seismic response isdominated by the strongly heterogeneous Eocene volcanic and Neogene basaltic rocks (Calvertet al., 2009). Figure 5.1 is modified from Calvert et al. (2011) and shows compiled geologicalinformation from Massey et al. (2005) as well as the location of relevant seismic lines and explo-ration wells. The surficial exposure of underlying Cretaceous sedimentary units is limited, andtheir characterisation is complicated by the presence of sedimentary rocks that were depositedcontemporaneously with Eocene extension and volcanism. The latter are 1) extensive, locally inexcess of 500 m thick (see Riddell, 2011, well log interpretations); 2) potentially indistinguishablefrom the underlying Cretaceous sediments (e.g., Kushnir et al., 2012); and 3) interbedded withvolcanic and volcaniclastic rocks that may exhibit widely varying seismic responses.Interpretation of well sonic logs and additional work in rock physics yield information aboutthe high-frequency seismic response of different rock types that can be correlated to changes inlithology (Mwenifumbo and Mwenifumbo, 2010; Riddell, 2011; Kushnir et al., 2012). However,geophysical investigations that make use of direct and refracted P-wave signals (Calvert et al.,2011; Hayward and Calvert, 2011; Smithyman and Clowes, 2012, 2013) have indicated thatthe in-situ response to low-frequency seismic waves cannot be reliably compared to the high-frequency sonic logs. In particular, the Neogene Chilcotin group basaltic volcanic rocks in thenear-surface region have shown seismic velocities that are 30?50% of velocities that are expectedfrom sonic logs (Kushnir et al., 2012). Comparisons between geophysical results and detailedgeological mapping (this paper; E. Bordet, personal communication) indicate that andesite unitswithin the Eocene volcanic sequences also exhibit surprisingly low P-wave velocities, on the sameorder as those observed for the Chilcotin basalt in situ.98Line 10Line 15Line 05Line 06 Line 11Line 12Exploration WellFigure 5.1: The surficial geology of the study area is depicted, modified from Calvert et al. (2011). Rel-evant seismic lines from the 2008 Geoscience BC seismic survey are depicted in red. The geologicalinformation database was synthesized from earlier results by Massey et al. (2005).Since rock types at depthmay not be distinguished unambiguously using their P-wave velocity,interpretation of the structures in the region must necessarily be done by careful extrapolationfrom known geological information. However, detailed information about the seismic (P-wave)velocities of the rocks in the upper few kilometres of the subsurface may still be applied to improvethe seismic image at depth, e.g., through the application of seismic migration techniques.995.1.1 Technical backgroundSeismic investigations are termed ?two-dimensional? if the data are acquired and processed toform a 2D image of the subsurface. The most common configuration for this is a nominally-straight seismic acquisition line that can be processed to form a cross-sectional image of thesubsurface. Three-dimensional seismic surveys are also commonly used (particularly in the re-source exploration industries) to characterize complex 3D subsurface structures by investigationusing a seismic array laid out in a plane or sheet. A class of methods exists, known as 2.5D(?two-and-one-half-dimensional?), which simulate 3D physics in a manner that makes specialconstraints on the geophysical model. These constraints make it possible to simulate (or invert)3D waveforms using a 2D model, and typically result in a substantial reduction in computationalcost when compared to 3D processing (Bleistein, 1986).Waveform tomography processing is applied to multi-channel seismic (MCS) data collectedin summer 2008 for Geoscience BC (Calvert et al., 2009). This method yields high-resolutionimages of in-situ P-wave velocity and numerical attenuation through nonlinear inversion of theseismic data waveforms. Synthetic waveforms are simulated using a viscoacoustic method in2.5D (Song and Williamson, 1995; Pratt, 1999; Smithyman and Clowes, 2013), which enablesdirect comparison with the field data waveforms that propagate in three dimensions. The seismicinversion grid is optimally aligned to cross geological strike if a dominant strike direction is knowna priori, which reduces significantly the approximation errors inherent in the 2.5D method.Waveform tomography makes use of eikonal-equation-based seismic tomography (Vidale,1990; Hole, 1992; Zelt and Barton, 1998) to develop a smooth model of subsurface P-wavevelocity based on inversion of first-arrival traveltime data. These data are interpreted manually(see discussion of picking strategy by Smithyman and Clowes, 2012) from the direct and refracteddata waveforms recorded in the 2008 Nechako Basin seismic dataset (Calvert et al., 2009). Thefirst-arrival traveltimes are inverted in 3D (using FAST; Zelt and Barton, 1998) and a representative2D model is extracted to be used for later full-waveform inversion processing.Full-waveform inversion (FWI) is applied (using a modified version of FULLWV; Pratt, 1999;100Smithyman and Clowes, 2013) to update the P-wave velocity model from traveltime inversion.FWI fits the phase of the data waveforms within a selected frequency band, which yields modelsof seismic velocity with improved resolution (Smithyman and Clowes, 2013). Additionally, FWIis applied on the data waveforms including amplitude information, which enables inversion fornumerical attenuation (parameterized as the inverse of the seismic Q factor at a centre frequencyof 10 Hz). However, the numerical model of attenuation that is recovered by FWI may notbe directly relatable to intrinsic Q (a physical property of rocks), because the FWI method maybe unable to distinguish between the disparate possible causes of reduced amplitude in seismicwaveforms. Scattering Q from P- to S-wave mode conversions and coda Q from structures smallerthan the resolution of the seismic model are both examples of phenomena that are not modelledby the 2.5D viscoacoustic FWI approach applied in this paper. Heuristic techniques may alsobe applied to accommodate amplitude variations between field seismic records and the syntheticrecords produced during FWI (e.g., log-linear AVO correction; see Smithyman and Clowes, 2013).These techniques may affect the recovered attenuation structures. Some ambiguities also exist inthe simultaneous recovery of P-wave velocity and numerical attenuation (e.g., Hak and Mulder,2011). However, models of numerical attenuation may be useful in distinguishing lithologicalboundaries based on the scattering potential of the rock units.5.1.2 Seismic lines and related worksThe seismic data discussed herein were acquired for Geoscience BC in the Nechako Basin (nowreferred to as the Nechako?Chilcotin plateau; Riddell, 2011) in 2008 by CGG Veritas (Calvertet al., 2009). Seven lines were acquired in total. However the discussion in this paper is focussedon Lines 05, 10 and 15, for which FWI results have been obtained, and Line 06 for which resultsfrom 3D traveltime tomography only have been derived. Each of the seven lines was processedby CGG Veritas to produce seismic stacks and post-stack time migrated (PostSTM) sections forinterpretation. These are discussed in detail by Calvert et al. (2009, 2011). Lines 05, 10 and15 were also processed byWesternGeCo MDIC (2010) using multiparameter joint inversion andsubsequent prestack depth migration (PreSDM) of the waveforms. Talinga and Calvert (2012)101interpret models from 3D traveltime tomography inversion of first-arrival data for all seven lines.The geometry of relevant seismic lines is shown in map view (Figures 5.1 and 5.2). The 2.5DFWI method characterizes each seismic profile with a cross-sectional (x, z) model and a cross-line direction (y). The cross-line direction is assumed to be the direction of geological strike, andshould be chosen based on a priori knowledge of the large-scale geology of the site in question(where possible; see Smithyman and Clowes, 2013). Figure 5.2 indicates the in-line and cross-linedirections chosen for each seismic profile, which define the direction of assumed homogeneity.Several multidisciplinary efforts have been undertaken that contribute to the geological infor-mation available along these profiles. Spratt and Craven (2010) and Calvert et al. (2011) processand interpret magnetotelluric data collected along eight profiles, including profiles along seismiclines 05, 06, 10 and 15. The positions of four exploration wells are indicated in Figures 5.1 and5.2; the well logs are interpreted by Mwenifumbo and Mwenifumbo (2010) (both a-4-L and b-16-J) and Riddell (2011) (all four wells); their interpretations provide information about the rockunits and lithology at depth. Surficial geology information is available from the BC-wide datasetofMassey et al. (2005), and detailed mapping along several of the profiles in the region has beencarried out (Bordet et al., 2011; personal communications, E. Bordet). Rock physics results, whichallow physical properties from geophysical inversion to be related to lithostratigraphic units, aretabulated by Kushnir et al. (2012).5.2 Waveform tomography resultsThe waveform tomography processing methodology employed in this study closely follows thatdiscussed by Smithyman and Clowes (2013), from which the results for Line 10 are included.Please see the previous study (my Chapter 4) for additional technical detail regarding the 2.5Dfull-waveform inversion methodology, which is summarized in our Section 5.1.1.Models of P-wave velocity and numerical attenuation are derived for three seismic lines, Lines5 and 15 and part of Line 10, from the 2008 Geoscience BC Nechako Basin seismic survey. Thesemodels fit first-arriving waveforms with bandwidth from 8 Hz to 13 Hz, with the exception of102yxyxyxLine 06Line 11Line 12Line 13LinLinLinLinUTM Northing (km)380 500400 420 440 460 480582059205840586058805900UTM Easting (km)Acquisition LineCDP Slalom LineFWI LineExploration WellLine 05Line 10Line 15AA?D D?CC?B?Bd-96-E a-4-L b-16-Jb-22-KFigure 5.2: The 2008 Geoscience BC Nechako Basin seismic survey comprised seven seismic acqui-sition lines. Three of these are processed to produce models from waveform tomography (Lines 05,10 and 15). The original acquisition line comprises the station positions occupied by geophones andvibrators. The CDP slalom line indicates the smoothed cross-sectional profile to which the CGG Ver-itas PostSTM image is projected (see text; Calvert et al., 2009). The orientations of the FWI in-line (x)and cross-line (y) axes are presented for each line discussed in this paper.Line 10, which includes bandwidth from 8 Hz to 16 Hz. The velocity models are derived usinga nonlinear multi-scale approach, with low data frequencies being inverted prior to high datafrequencies. In all cases, phase-only inversion is carried out until the last group of iterations,when logarithmic full-waveform inversion is carried out on the data phase and amplitude (Shinand Min, 2006). The model of attenuation is not modified until the last stages of full-waveforminversion using both amplitude and phase.103Due to ongoing research developments, the processing methodology used with Lines 05, 06and 15 differs slightly from that described by Smithyman and Clowes (2013):? For Lines 06 and 15, the traveltime tomography updates are calculated in a 3D model and a 2Dmodel image is extracted using ray-density information. This differs from the 2.5D traveltimetomography method used for Lines 05 and 10.? For Lines 05 and 15, log-linear correction of field data amplitudes is not applied, which leadsto scaling differences in the recovered numerical attenuation structure (discussed by Smithymanand Clowes, 2013; cf. Figure 5.10d with Figures 5.3d and 5.14d).? For Lines 05 and 15, the highest inverted data frequency is 13 Hz, since useful model updateswere not significant above 13 Hz in previous results (Smithyman and Clowes, 2013).? The vertical size of each model is increased, extending from 2000 m above sea level to ?3000 m.The velocity and attenuation models produced from waveform tomography processing of theseismic data collected in 2008 are critical tools for characterization of Eocene volcanic units andthe underlying strata. However, the physical properties that are investigated by these methods arenot unique to a particular unit or formation, since they correlate best with broad lithological types(see Kushnir et al., 2012). The identification of rock units must be carried out in relation to existinggeological information from surficial mapping and exploration wells. The models of numericalattenuation are generally less well-constrained than the models of velocity, and are mainly usefulas a diagnostic tool that relates information about localized energy loss, since several un-modelledphysical effects may contribute to the recovered values (see Section 5.1.1 and Smithyman andClowes, 2013). In particular, loss of energy at interfaces due to P-to-S-conversion may be mappedto low Q, resulting in high apparent attenuation at lithostratigraphic boundaries. Additionally,smoothing (a low-pass filter applied to the FWI gradient with length scale on the order of 1?2 km,dependent on line characteristics) was applied during the inversion for attenuation, resulting in aless detailed image than the model of P-wave velocity. It has been observed that the models of theinverse seismic Q-factor generally indicate high apparent attenuation in Cretaceous sedimentaryrocks (Smithyman and Clowes, 2013), which is confirmed from outcrop and well log comparisons104on Line 05 (this chapter, Section 5.2.1). This is consistent with previous results from Line 10 (viz.,Smithyman and Clowes, 2013). Previous work also suggests that the coherent volcanic rocks ofthe Eocene and Neogene generally exhibit low attenuation where they can be reliably identified.5.2.1 Line 05: southern NazkoLine 05 was acquired along the southern extent of the Nazko River valley, which trends northwest-southeast (see Figures 5.1 and 5.2). The presumed geological strike from regional-scale tectonicfeatures is approximately along-line, which makes cross-strike alignment of the inversion griddifficult or impossible with the line acquisition geometry. The +x-direction for FWI was set to158? to best fit the line of acquisition, with a priori assumed strike direction of 68?. However,recent mapping work (personal communications, E. Bordet) indicates that local geological strikeis eastward (approximately orthogonal to the line of acquisition) in some outcropping Cretaceousstrata near the Indian Head Promontory (approx. UTM 470.0 km E, 5847.5 km N; Figure 5.2; lineposition x ? 15 km). Figure 5.3 shows model results from waveform tomography for Line 05.In order to assess the quality of the inversion result, forward-modelled seismic data werecompared with the field seismograms in both time and frequency domains. Figure 5.4 shows threesource gathers (sources 100, 400 and 700) comparing field data and synthetic data. The phaseerror at 10 Hz is also shown for the initial model and the FWI model after 50 iterations of phase-only inversion, prior to updating the model of attenuation in iterations 51?60. These comparisonsindicate that the data fit is most optimal in the southern portion of the section, which is mainlyunderlain by Cretaceous sedimentary rocks. The data fit in the northern part of the section wasless optimal, which is visible on the gather for source 700 (Figure 5.4; A) where a package ofreflected arrivals is not matched properly leading to a cycle skip at middle offsets. However, theshallower part of the section is still resolved reliably in these areas. The reduction in phase errorat 10 Hz (the dominant data frequency) is visible in Figure 5.4, between the initial model (seelabel B and the model at the end of phase-only waveform inversion (see label C). The final 10iterations of phase + amplitude inversion resulted in a moderate increase in the phase residual,but a significant decrease in the amplitude residual.105Elevation (km)+2+10?1?2?3 0 5 10 15 20 25 30Position (km)Elevation (km)+2+10?1?2?3 0 5 10 15 20 25 30Position (km)Elevation (km)+2+10?1?2?3 0 5 10 15 20 25 30Position (km)Elevation (km)+2+10?1?2?3 0 5 10 15 20 25 30Position (km)1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0Velocity (km/s)0.00 0.05 0.10 0.15 0.20Attenuation (1/Q)?1.0 ?0.5 +0.5 +1.0? Velocity (km/s)0.0a)b)c)d)N SA BCFigure 5.3: Line 05; starting (a) and final (b) velocity models are shown. The difference betweenstarting and final models is shown (c) and the attenuation is presented (d). In (a) and (b), grey andblack contour lines show velocities as indicated on the velocity scale; in (d) the contour lines showattenuation values as indicated on that scale. In (c), bold letters correspond to those in Figures 5.5and 5.6 and are discussed in the text.1060.01.0Time (s)0.01.0Time (s)0.01.0Time (s)Source 100Source 400Source 700106011 2118 1 2118Receiver ReceiverSourcePhase error @ 10 Hz in initial model Phase error @ 10 Hz after 50 iter.100400700?pi+piPhase Error0AAB C~4.7 km N ~14.4 km S~4.7 km N ~14.4 km S~4.7 km N ~14.4 km S10 km15 km20 km25 km30 km15 km 20 km25 km 30 km5 km 5 km 10 km10 km15 km20 km25 km30 km15 km 20 km25 km 30 km5 km 5 km 10 kmFigure 5.4: Line 05; field data (black) and synthetic data (green) are shown for three shot gathers in theupper three panels. Source gathers are time-reduced, which flattens the arrivals to a reduction velocityof 4.5 km/s. The fourth set of panels reproduces phase residuals at 10 Hz for the initial model fromtraveltime inversion and the model after 50 iterations of phase-only FWI. For a perfect fit to noise-freedata, the phase error would be zero. Black lines identify the location of each shot gather on the panelsof phase residual information. Bold letters identify features that are discussed in the text.Figure 5.5 presents the final FWI velocity model from Figure 5.3b projected to the CDP slalomline A? A? (see Figure 5.2). Figure 5.6 presents the model of numerical attenuation from Figure5.3d. Figures 5.5 and 5.6 also overlay the PostSTM image fromCalvert et al. (2009), which enablescomparison between the two results. The reflectivity structure is well-resolved by PostSTM onLine 05. Specific features are identified using bold letters in text and images. Exploration well a-4-L is located in the Nazko River valley alongside Line 05 (see Figure 5.2), and has been interpretedbyMwenifumbo and Mwenifumbo (2010) and Riddell (2011). Figure 5.8 reproduces these results107including sonic log velocities fromMwenifumbo and Mwenifumbo (2010) and overlays them onsegments extracted from Figures 5.5 and 5.6 for direct comparison.10805001000150020002500Time (ms)13816North SouthCDPA A?2.50 4.25 6.00Velocity (km/s)1 km3 kmVE3:1Ground Surfacea-4-L0 5 10 15 20 25 30Position (km)ABCDENb Ef Csv NbFigure 5.5: Line 05; the PostSTM image from Calvert et al. (2009) is overlaid on the model of P-wave velocity from waveform tomographyprocessing, which is converted to time using its own velocities, a replacement velocity of 3 km/s and a datum of 1500 m elevation. Thewaveform tomography result is projected from the FWI profile to the CDP slalom line A-A?, located on Figure 5.2. Positions in km areapproximate, due to the projection from the FWI profile. The geological units at surface are indicated at the top of the image; Nb: Neogenebasalts; Ea: Eocene andesite; Ef: Eocene felsic rocks; Csv: Cretaceous sediments and volcanic rocks; Jev: Jurassic or earlier volcanics andvolcaniclastic rocks. Bold letters refer to features that are discussed in the text. Label a-4-L indicates the projected location of the explorationwell (Figure 5.2) for which well-log information is shown in Figure 5.8.10905001000150020002500Time (ms)13816North SouthCDPA A?0.100.00 0.20Attenuation (1/Q)1 km3 kmVE3:1Ground Surface0 5 10 15 20 25 30Position (km)a-4-LABCDENb Ef Csv NbFigure 5.6: Line 05; the PostSTM image from Calvert et al. (2009) is overlaid on the model of numerical attenuation (1/Q) from waveformtomography processing, which is converted to time as in Figure 5.5. The waveform tomography result is projected from the FWI profile tothe CDP slalom line A-A?, located on Figure 5.2. Positions in km are approximate, due to the projection from the FWI profile. The geologicalunits at surface are indicated at the top of the image; Nb: Neogene basalts; Ea: Eocene andesite; Ef: Eocene felsic rocks; Csv: Cretaceoussediments and volcanic rocks; Jev: Jurassic or earlier volcanics and volcaniclastic rocks. Bold letters refer to features that are discussed inthe text. Label a-4-L indicates the projected location of the exploration well (Figure 5.2) for which well-log information is shown in Figure5.8.1100 2 4 6 8Velocity (km/s)Well a-4-LLine 05x = 267500 2 4 6 8Velocity (km/s)Well d-96-ELine 05x = 267500 2 4 6 8Velocity (km/s)Well b-16-JLine 15x = 8750Sonic log T T velocity FWI velocity+10000?1000?2000?3000Elevation above MSL (m)Figure 5.7: Well sonic logs from exploration wells d-96-E, a-4-L and b-16-J are compared to 1Dvelocity profiles extracted from TT and FWI velocity models for Lines 05 and 15. The locations of thesewells are indicated on Figures 5.1 and 5.2. Wells d-96-E and a-4-L are adjacent, and are comparedwith the 1D velocity function from Line 05 (x = 26750 m); note that d-96-E is west of the line andbegins at a higher elevation. Well b-16-J is located approximately 21 km south of Line 15, and iscompared with the velocity profile from x = 8750 m. Gray background indicates portions aboveground level for the model position indicated.Figure 5.7 shows a direct comparison between the TT and FWI velocity models (Figure 5.3)with sonic logs from wells d-96-E and a-4-L (well b-16-J is compared with Line 15; section 5.2.4).Both the initial and final models from waveform tomography processing fit the well logs reason-ably well with the top 1500 m of the model. Well d-96-E shows generally slower velocities fromabout?500m elevation downwards, with a mismatch of approximately 1 km/s. However, the rel-1110720105013501600211023802500Well a-4-LShaleSandstoneShaleVolcanicsDepth (m)AlbianSedimentsCampanianSedimentsJurassicVolcanicsJurassicPluton32303311Riddell(2011)Mwenifumbo(2010)2.5 6.0SonicLog Vel.(km/s)0500100015002000Time (ms)2.5 6.0Velocity(km/s)0.0 0.2Attenuation(1/Q)Figure 5.8: Well-log interpretations are reproduced from Mwenifumbo and Mwenifumbo (2010) andRiddell (2011) for exploration well a-4-L, located on Figure 5.2. Riddell (2011) interprets the core to3311 m depth and provides unit descriptions that may be used to relate well logs and surficial geologyduring interpretation. Mwenifumbo and Mwenifumbo (2010) provides useful information about thesub-unit lithology of the shallowest 2500 m from the exploration well, which may be correlated tophysical property contrasts in the full-waveform inversion images (this work) and coherent reflectorsin the PostSTM image (Calvert et al., 2009) for seismic Line 05. Right-side panels are extracted fromFigures 5.5 and 5.6, which overlay reflectivity with velocity and attenuation, respectively.ative velocity profile of well d-96-E is similar below ?1000 m elevation. Well a-4-L shows largermismatches at depth, with a maximum 1.5 km/s difference from ?1000 to ?1200 m elevation.The deeper section of the waveform tomography model highly dependent on the initial modelfrom traveltime tomography, which may result in overestimation of velocities below about 1500m depth. Well-log interpretations from other workers (Figure 5.8) indicate that a succession of112Cretaceous sediments dominates the lithostratigraphy in the southern portion of the line. Theseextend to about 2400 m depth in the vicinity of well a-4-L (line position x ? 27 km, Figure 5.3),and are laterally continuous at surface at least as far north as x ? 15 km. Figure 5.3 shows P-wavevelocities between about 4 and 5 km/s that can be matched locally to this unit, and which agreewith the expected values from Kushnir et al. (2012). We infer that the Albian sediments are locallycontinuous at least between about x = 18 km and x = 32 km, based on consistent physical prop-erties visible in the velocity and attenuation models (Figure 5.3). These are underlain by Jurassicvolcanic rocks, possibly of the Stikine or (less likely) Cache Creek terranes (seeMwenifumbo andMwenifumbo, 2010, Riddell, 2011 and Calvert et al., 2011; well log and seismic interpretation).Jurassic rocks are interpreted in well a-4-L at approximately 2400 m depth (elevation z ? ?1.5km), but extrapolation from the waveform tomography velocity model suggests that these unitsmay be found 500?1000 m deeper to the north (see Figure 5.3b), especially from approximatelyx = 5 km to x = 20 km. Surface geological mapping efforts (personal communications, E. Bordet)observe southward-dipping strata at surface x = 15 km, on the Indian Head promontory. The ve-locity and attenuation models from FWI express similar structures, and combined interpretationof these in conjunction with the PostSTM image from Calvert et al. (2009) may allow extractionof in-plane dip information (see A, Figure 5.5). Structures in the model of P-wave velocity (Figure5.3b) indicate near-surface apparent dips close to 45? and dips at 2 km depth (x = 16 km) ofapproximately 20?.The PostSTM image possesses coherent reflectors that are well correlated with changes inseismic velocity from waveform tomography, particularly throughout the 500?1500 ms time win-dow (Figure 5.5), and the combination of these may be used to inform structural interpretation.A series of sub-horizontal features (dipping 10? to 15?) are visible on Figure 5.5 and 5.6 cross-ing well a-4-L at ? 1000 ms, labelled B. Figure 5.8 overlays two interpretations of the geologyfrom well a-4-L on excerpts from Figures 5.5 and 5.6. Seismic reflectivity (Calvert et al., 2009)and FWI velocity (Figure 5.5) correlate strongly with the lithological interpretation and sonic logsfrom Mwenifumbo and Mwenifumbo (2010). The sandstone/conglomerate encountered by wella-4-L is expected to possess a higher seismic velocity than the shales surrounding it, which is in113agreement with the model of P-wave velocity from waveform tomography. The occurrence ofa laterally-continuous, dipping low-velocity zone corresponds exactly with the interpreted shalelayering between 1050 m and 1350 m depth (B). This low-velocity zone may be extrapolatedfrom examination of the velocity from waveform tomography and the PostSTM image over atleast 4 km from x = 27 km to x = 31 km. The same structure may extend northward to x ? 20km, based on examination of the velocity difference image in Figure 5.3c, labelled C; however,this extension is speculative. The correlation between the velocity model and sonic logs is lessobvious below about 1200 ms, though a reduction in P-wave velocity is observed in proximity tothe presumed extent of Campanian sediments from the well (more obvious in Figure 5.3c). Cam-panian sediments that are identified in exploration well a-4-L (Mwenifumbo and Mwenifumbo,2010; Riddell, 2011) do not appear to outcrop at surface (personal communications; E. Bordet).The resolution of full-waveform inversion is substantially reduced at these depths, and whereasthe coarse velocity structure is compatible with the exploration well logs, we do not expect toidentify detailed structures in relation to the well below ? 1500 ms.Surface mapping efforts (personal communications, E. Bordet) indicate that the surface lithol-ogy is dominated by volcanic rocks to the north of the Indian Head promontory. These are mainlyEocene in age, but are overlain by Chilcotin Group basaltic rocks north of about UTM 5850.0 kmN (Figure 5.2). Dipping structures in the PostSTM image (Figures 5.5 and 5.6) appear to indicatecomplex paleostratigraphy, and comparison with the velocity and attenuation images suggeststhat two paleo-basin structures are present (north and south of a high-velocity/high-attenuationridge at D). Alternatively, these could be the expression of a single north?south basin to the eastof the line, whereas the feature at D could be representative of an adjacent structure that appearsin contrast due to the line curvature (Figure 5.2). Since Chilcotin group rocks are not expected tobe thicker than 100?200 m (and are usually < 50 m thick) in this region, the low-velocity zoneat surface (0 < x < 10 km, see Figure 5.3) is most likely not composed of exclusively Neogenevolcanics. One possible explanation for the low seismic velocities would be in-fill of incoherentEocene clastic material. In this case, we infer that the low-velocity/low-attenuation material at. 500m depth (Figure 5.3) may be Chilcotin volcanics and that the low-velocity/high-attenuation114materials in the next ? 1500 m may represent incoherent clastic sediments. The low P-wave ve-locities near surface at the southern end of Line 05 (E) are not well-constrained by the seismicacquisition, but earlier mapping efforts (summarized by Massey et al., 2005) indicate that thisarea is covered by Neogene basaltic rocks, which are observed to exhibit low velocities. Dippingstrata in the PostSTM image support the interpretation of a basin structure of ? 500 m depth inthis area.5.2.2 Line 06: northern BaezaekoLine 06 was acquired along logging roads west of the community of Nazko, crossing near theNazko cone and onto outcropping Jurassic bedrock (Figures 5.1 and 5.2; B-B?). Traveltime inver-sion was carried out on hand-picked first arrivals in 3D using FAST (Zelt and Barton, 1998). Inorder to assess the continuity of features between lines, Line 06 was added to the scope of thisstudy. However, due to time and budget constraints, full-waveform inversion was not carried outfor this line. The combination of the coarse (traveltime tomography) velocity model and preexist-ing PostSTM reflectivity image (Calvert et al., 2011) still enable a simplified interpretation of thestructures along Line 06. However, the quality of the PostSTM image for Line 06 is somewhatlow and significant reverberations are present near coherent events. This makes precise locationof stratigraphic boundaries difficult. First arrivals were picked by hand for all traces in the dataset,and every 4th shot gather was inverted to generate a 3D model. The final RMS traveltime misfitwas 31 ms. A 2D profile was extracted based on weighting from ray-density information andprojected to the CDP slalom line of Calvert et al. (2011). Figure 5.9a shows the model of P-wavevelocity projected to the UTM E?W (perpendicular to strike). Figure 5.9b presents the velocitymodel converted to time with overlaid PostSTM reflectivity.Jurassic bedrock outcrops at surface under the northwestern half of Line 06, resulting in astrong shallow high-velocity structure (A) visible in Figure 5.9. This structure appears to transitionabruptly to a low-velocity region at B. Mapping indicates mixed Eocene volcanics that dominatethe surficial geology in the southeastern half of the line. Reflectors in the area of B are disruptedadjacent to the high-velocity body at A, and it is likely that the younger material was emplaced in115a pull-apart basin during Eocene extension. The Nazko Cone is visible at C as a high-velocity ( 5km/s) anomaly, due to the presence of basaltic lavas that are somewhat faster than surroundingEocene rocks. Low velocities observed in the Nazko River valley (D) are most likely a continua-tion of Neogene basalts of the Chilcotin group that are observed at the northern end of Line 05.Dipping reflectors at depth (E and F) are piecewise-continuous from about 800 ms to 2000 ms inthe seismic image. These appear to correlate with an increase in velocity from about 4 km/s to5 km/s, which is consistent with either Cretaceous material (e.g., basin conglomerates) or oldervolcanic rocks. Some of the material at these depths may be related to Neogene volcanism, butif so it is indistinguishable from surrounding rocks by P-wave velocity at the resolutions availablefrom traveltime inversion.5.2.3 Line 10: southern BaezaekoLine 10 was acquired along logging roads west of the community of Nazko, with a total line lengthof approximately 77 km (Figures 5.1 and 5.2). However, we process and discuss a subset of theline comprising 699 vibroseis source points and 2362 distinct receiver locations, resulting in aprofile ? 45 km long. The direction of geological strike was taken to be 17?, which is reasonablegiven the trend of exposed basement and Eocene-age rocks at surface. However, this choice wasmade initially to optimally fit the orientation of the local section of the acquisition line for 2D FWIprocessing (Smithyman and Clowes, 2012) prior to the development of the 2.5Dwaveform tomog-raphy methodology discussed by Smithyman and Clowes (2013). The +x direction was alignedto 107? (see Figure 5.2). The model of P-wave velocity from traveltime inversion of the Line 10dataset is presented in Figure 5.10a, and the updated model from subsequent full-waveform in-version is presented in Figure 5.10b. The difference between these two results is shown in Figure5.10c, and the model of (numerical) viscoacoustic attenuation is shown in Figure 5.10d. Theseresults are reproduced from Smithyman and Clowes (2013), which also presents an expanded dis-cussion of the technical aspects of the processing. The traveltime tomography method used onLine 10 is 2.5D; i.e., updates are calculated in 3D and smoothed to enforce two-dimensionality.The full-waveform inversion step improves the interpretable resolution by perturbing the initial116model (see difference image; Figure 5.10c). Note, however, that restrictions to the data extentand imposed smoothing and muting cause the results to be mainly useful between x = 10 kmand x = 35 km. This is particularly noticeable in the model of attenuation (Figure 5.10d), sinceupdates were only permitted in the centre of the model during the final stages of inversion.Figure 5.11 presents three time-reduced source gathers (sources 200, 350 and 600) comparingfield data and synthetic data; these are equivalent to the comparisons by Smithyman and Clowes,2013, their Figures 7e, 8e and 9e, respectively). The phase error at 10 Hz is also shown for theinitial model and the FWI model after 20 iterations of phase-only inversion, prior to updating themodel of attenuation in iterations 21?28. The final 10 iterations of phase + amplitude inversionresulted in a moderate increase in the phase residual, but a significant decrease in the amplituderesidual, as discussed by Smithyman and Clowes (2013). The data fit is generally quite good alongthis line, with one notable exception being a region that generates a localized cycle skip (Figure5.11, see A), identified on the gather for source 600 and discussed by Smithyman and Clowes(2013)).Figure 5.12 reproduces the velocity model of Figure 5.10b. Figure 5.13 reproduces the atten-uation model. Both figures overlay PostSTM reflectivity from Calvert et al. (2009) for comparison.The quality of the PostSTM image for Line 10 is only moderate, and precise identification of fea-tures from the reflectivity is difficult. Features of interest are identified using bold labels in thefigures and text.11705001000150020002500Time (ms)13346Northwest SoutheastCDPB B?2.50 4.25 6.00Velocity (km/s)1 km3 kmVE3:1Ground Surface435 440 445 450Position (km)455 460b)A CB DE FElevation (km)+2+10?1?2?3 435 440 445 450 455 460Position (km)1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0Velocity (km/s)W Ea)N b Jev CsvEf N b Ef Ef N bEf N bFigure 5.9: Line 06; (a) The model of P-wave velocity from traveltime tomography is shown. Thedimmed region is not sampled by rays. Velocity contours (gray and black) are indicated on the colourbar. (b) The PostSTM image from Calvert et al. (2009) is overlaid on the model of P-wave velocity fromtraveltime tomography processing, which is converted to time using its own velocities, a replacementvelocity of 3 km/s and a datum of 1500 m elevation. The traveltime tomography result is projectedto the CDP slalom line B-B?, located on Figure 5.2. Positions in km are approximate, due to theprojection from the FWI profile. The geological units at surface are indicated at the top of the image;Nb: Neogene basalts; Ea: Eocene andesite; Ef: Eocene felsic rocks; Csv: Cretaceous sediments andvolcanic rocks; Jev: Jurassic or earlier volcanics and volcaniclastic rocks. Bold letters refer to featuresthat are discussed in the text.118Elevation (km)+2+10?1?20 5 10 15 20 25 30Position (km)35 40 45Elevation (km)+2+10?1?20 5 10 15 20 25 30Position (km)35 40 45Elevation (km)+2+10?1?20 5 10 15 20 25 30Position (km)35 40 45Elevation (km)+2+10?1?20 5 10 15 20 25 30Position (km)35 40 45a)1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0Velocity (km/s)0.0E?02 1.5E?02 3.0E?02 4.5E?02 6.0E?02Attenuation (1/Q)?1.0 ?0.5 +0.5 +1.0? Velocity (km/s)0.0b)c)d)W EFigure 5.10: Line 10; starting (a) and final (b) velocity models are shown. The difference betweenstarting and final models is shown (c) and the attenuation is presented (d). In (a) and (b), grey andblack contour lines show velocities as indicated on the velocity scale; in (d) the contour lines showattenuation values as indicated on that scale.1190.01.0Time (s)0.01.0Time (s)0.01.0Time (s)69911 2362 1 2362Receiver ReceiverSourcePhase error @ 10 Hz in initial model Phase error @ 10 Hz after 20 iter.200350600?pi+piPhase Error0Source 200Source 350Source 600~14.4 km W ~4.7 km E~14.4 km W ~4.7 km E~14.4 km W ~4.7 km EAA10 km15 km20 km25 km30 km15 km 20 km25 km5 km5 km10 km10 km15 km20 km25 km30 km15 km 20 km25 km5 km5 km10 kmFigure 5.11: Line 10; field data (black) and synthetic data (green) are shown for three shot gathersin the upper three panels. Source gathers are time-reduced, which flattens the arrivals to a reductionvelocity of 4.5 km/s. These panels are equivalent to those shown by Smithyman and Clowes (2013),their Figures 7e, 8e and 9e. The fourth set of panels reproduces phase residuals at 10 Hz for the initialmodel from traveltime inversion and the model after 20 iterations of phase-only FWI. For a perfect fitto noise-free data, the phase error would be zero. Black lines identify the location of each shot gatheron the panels of phase residual information. Bold letters identify features that are discussed in the text.12005001000150020002500Time (ms)20005500West EastCDPC C?2.50 4.25 6.00Velocity (km/s)1 km3 kmVE3:1Ground SurfaceACDEFBGH5 10 15 20 25 30Position (km)35Nb Ef Nb EfNb EfFigure 5.12: Line 10; the PostSTM image from Calvert et al. (2009) is overlaid on the model of P-wave velocity from waveform tomographyprocessing (Smithyman and Clowes, 2013), which is converted to time using its own velocities, a replacement velocity of 3 km/s and adatum of 1500 m elevation. The waveform tomography result is projected from the FWI profile to the CDP slalom line C-C?, located onFigure 5.2. Positions in km are approximate, due to the projection from the FWI profile. The geological units at surface are indicated at thetop of the image; Nb: Neogene basalts; Ea: Eocene andesite; Ef: Eocene felsic rocks; Csv: Cretaceous sediments and volcanic rocks; Jev:Jurassic or earlier volcanics and volcaniclastic rocks. Bold letters refer to features that are discussed in the text.12105001000150020002500Time (ms)20005500West EastCDPC C?0.00 0.04 0.08Attenuation (1/Q)1 km3 kmVE3:1Ground SurfaceACDEFBGH5 10 15 20 25 30Position (km)35Nb Ef Nb Ef Ef CsvFigure 5.13: Line 10; the PostSTM image from Calvert et al. (2009) is overlaid on the model of numerical attenuation (1/Q) from waveformtomography processing (Smithyman and Clowes, 2013), which is converted to time as in Figure 5.16. The waveform tomography result isprojected from the FWI profile to the CDP slalom line C-C?, located on Figure 5.2. Positions in km are approximate, due to the projectionfrom the FWI profile. The geological units at surface are indicated at the top of the image; Nb: Neogene basalts; Ea: Eocene andesite; Ef:Eocene felsic rocks; Csv: Cretaceous sediments and volcanic rocks; Jev: Jurassic or earlier volcanics and volcaniclastic rocks. Bold lettersrefer to features that are discussed in the text.122The central portion of the model from about UTM 425.0 km E eastward (Figure 5.2) is overlainby dacite of variable thicknesses (e.g., A) with some Quaternary cover (personal communications,E. Bordet), previously mapped as Ootsa Lake volcanic rocks (Massey et al., 2005). West of aboutx = 10 km the cover changes, and the geology of the portion of the line discussed is dominatedby Neogene-age Chilcotin Group basalts. These likely result in the low seismic velocities under3 km/s observed in the synclinal structure at B, to a maximum thickness of a few hundred metres.However, the Neogene basalts are usually seen to be less than 50 m thick in outcrop (personalcommunications, E. Bordet), which is insufficient to explain the full depth of the synclinal featureat B, and therefore it is likely in-filled by other material at depth (possibly Eocene-age volcaniclas-tics or sediments). Examination of the P-wave velocity structure and PostSTM reflectivity (Figure5.12) in conjunction with the velocity model in depth (Figure 5.3b) suggests that the depth toJurassic basement rocks in this area is variable from about 1500 m to less than 1000 m at x = 10km and approximately 500 m at x = 12 km (C). High P-wave velocities near surface indicatea probable southward extension of Hazelton Group volcanic rocks of the Stikine Terrane thatoutcrop north of the line. The depth to basement rocks is interpreted to be somewhat greater inthe centre of the line.Between 12 . x . 26 km, dacite rocks of a few hundred metres thickness are interpreted withP-wave velocity between about 3 km/s and 4.5 km/s and a low attenuation (A). A thin surfacelayer with anomalously low seismic velocities at D is believed to be due to Quaternary cover.The region of high seismic attenuation and moderate P-wave velocities (E) between 4.5 km/s and5.0 km/s is believed to represent a faulted Cretaceous sedimentary succession, reactivated andextended contemporaneously with Eocene volcanism. Layering is preserved locally (see PostSTMimage; Figures 5.12 and 5.13), but dips at 25? east where visible. The depth to Jurassic basementis 2000?2500 m in this area (F). A prominent high-velocity region (G) is overlain by volcanicrocks less than 500 m thick (likely dacite) between 25 . x . 33 km. The material at G exhibitsP-wave velocities in excess of 5 km/s, and is believed to represent a southward extension of theHazelton Group volcanic rocks of the Stikine Terrane that outcrop 2?3 km north of the line.Eocene volcanic rocks are interpreted to overlie Cretaceous sediments at H, with extensional123basin features that were likely contemporaneous with deposition. However, shallow materialswith P-wave velocities under about 3 km/s are most likely Chilcotin basalt, which are mappednorth of the line and extend under it towards the east.5.2.4 Line 15: TibblesLine 15 (Figures 5.1 and 5.2) was acquired along the east-west trending Tibbles Road. We exam-ined surficial geology, models from inversion of gravity data and magnetic susceptibility informa-tion (Bordet et al., 2011) in order to estimate a local geological strike direction of approximately23?. The +x grid direction for FWI was set to 113? (see Figure 5.2). Large-scale synthesis mapsindicate that the surficial geology of Line 15 is dominated by mixed Eocene volcanic rocks, butrecent detailed mapping efforts (personal communications; E. Bordet) identify distinct lithologiesthat may possess differing seismic responses.Figure 5.14 presents velocity models from traveltime inversion (a) and full-waveform inver-sion (b) as well as the model of numerical attenuation (d) from FWI. The difference between thetomography and FWI velocity models (c) is particularly useful in identifying isolated regions ofhigh and low velocity that are not identified by traveltime tomography. Due to the small acqui-sition footprint (Figure 5.2) with only 378 source stations and 755 receiver stations over a linelength of ? 14 km, full seismic fold (survey density) is only reached in the centre of the model.Consequently, the model is constrained by raytracing and traveltime inversion to about 1000 mdepth. Themodel of FWI is able to constrain the model somewhat better, but is only reconstructedreliably in the first ? 2000 m of the section, which is illuminated by the short-offset seismograms.Figure 5.7 compares 1D velocity profiles from TT and FWI (Figure 5.14) with well sonic logs fromexploration well b-16-J, located ? 21 km to the south. The fit between the extracted velocitymodels is relatively close, suggesting similarities between the structures found in well b-16-J andthose present along Line 15. A low-velocity zone in b-16-J (?300 m elevation) is not observed inthe results from Line 15, but may represent a local feature.Comparison of field and synthetic data and the frequency-domain phase residuals (Figure5.15) indicates that portions of the model illuminated by short-offset source/receiver pairs (to124Elevation (km)+2+10?1?2?3 0 4 8 12Position (km)Elevation (km)+2+10?1?2?3 0 4 8 12Position (km)Elevation (km)+2+10?1?2?3 0 4 8 12Position (km)Elevation (km)+2+10?1?2?3 0 4 8 12Position (km)1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0Velocity (km/s)0.00 0.05 0.10 0.15 0.20Attenuation (1/Q)?1.0 ?0.5 +0.5 +1.0? Velocity (km/s)0.0a) b)c) d)W E W EFigure 5.14: Line 15; starting (a) and final (b) velocity models are shown. The difference betweenstarting and final models is shown (c) and the attenuation is presented (d). In (a) and (b), grey andblack contour lines show velocities as indicated on the velocity scale; in (d) the contour lines showattenuation values as indicated on that scale.about 4500 m offset) fit the field seismograms accurately (A). However, the long-offset seismicdata outside about 4500 m offset are not reproduced accurately (B). Examination of the sourcegathers in Figure 5.15 shows that the synthetic data contain a wide-angle reflection at C thatcontinues to become a first-arrival in the synthetic data. The reflection event is present in thefield seismograms but does not possess sufficient amplitude to be visible in the first arrivals ofthe field data. This results in a mismatch between the field data and the model from traveltimeinversion?particularly at depth?that makes full-waveform inversion of the long-offset data very125difficult. The noise level is also quite high in some of the long-offset traces (e.g.,D), owing to poorcoupling of the acquisition equipment in sands on the acquisition line (Calvert et al., 2009). Thequestionable fit to data phase at long offset means that the deeper model?below about 1000?1500 m depth in the centre of the spread?may not be reliable. We therefore mainly constrainour discussion to the upper portion of the model.Figure 5.16 presents the final FWI velocity model from Figure 5.14b projected to the CDPslalom line D?D? (see Figure 5.2. Figure 5.17 presents the model of numerical attenuation fromFigure 5.14d. Both figures include overlaid seismic reflectivity from PostSTM processing (Calvertet al., 2009). The PostSTM image from Line 15 is relatively high-quality, and structures can becorrelated easily to the models of velocity and attenuation. Features of interest are identified usingbold letters on figures and in the text.Changes in seismic velocities near surface (above 500 ms; Figures 5.16 and 5.17) correlatestrongly with changes in the mapped surface lithology from personal communications (E. Bordet).Based on the surface geological mapping, low seismic velocities at A and B are believed to reflectthe presence of Eocene andesite rocks, with in-situ P-wave velocities less than about 3 km/s. Therocks in this region also possess anomalously high numerical attenuation (Figure 5.17), whichmay indicate that they are less intact than surrounding volcanics. Eocene dacite (C) is mappedunderlying the line of acquisition at 5 . x . 7.5 km and in close proximity to the line eastof x ' 5 km, especially near x ' 3.5 km. The dacite rocks appear to exhibit higher in situ P-wave velocities in the range of 3.5?4.5 km/s when compared to the andesite rocks. The PostSTMreflectivity possesses coherent reflectors in the time range< 500ms in the presence of interpreteddacite rocks, whereas the reflectivity does not show coherent structure in the shallow subsurface inregions A and B. This reinforces the interpretation that the rocks of andesite affinity are incoherent.Eocene rhyolite rocks are mapped west of about UTM 480.5 km E (x ' 2.5 km) and undividedEocene rocks are mapped east of about UTM 490.0 km E (x ' 11 km), but these are outside thereliable acquisition footprint and consequently in situ P-wave velocities are not known reliably.A strong reflector is identified in the PostSTM image (Calvert et al., 2009) at D that is laterallycontinuous for at least 6 km (from 6 . x . 12 km) and correlates with a localized high-velocity1260.01.0Time (s)0.01.0Time (s)0.01.0Time (s)3780011 755 1 755Receiver ReceiverSourceSource 100Source 200Source 300Phase error @ 10 Hz in initial model Phase error @ 10 Hz after 80 iter.100200300?pi ?pi/2 +pi/2 +piPhase Error02.7 km E10.0 km W5.4 km E7.3 km W9.1 km E3.6 km WA AABBBBAAAABBDDCCCDD6 km 8 km 10 km2 km 4 km6 km8 km10 km 2 km4 km6 km 8 km 10 km2 km 4 km6 km8 km10 km 2 km4 kmFigure 5.15: Line 15; field data (black) and synthetic data (green) are shown for three shot gathersin the upper three panels. Source gathers are time-reduced, which flattens the arrivals to a reductionvelocity of 4.5 km/s. The fourth set of panels reproduces phase residuals at 10 Hz for the initial modelfrom traveltime inversion and the model after 80 iterations of phase-only FWI. For a perfect fit to noise-free data, the phase error would be zero. Black lines identify the location of each shot gather on thepanels of phase residual information. Bold letters identify features that are discussed in the text.127Figure 5.16: Line 15; the PostSTM image from Calvert et al. (2009) is overlaid on the model of P-wavevelocity from waveform tomography processing, which is converted to time using its own velocities,a replacement velocity of 3 km/s and a datum of 1500 m elevation. The waveform tomography resultis projected from the FWI profile to the CDP slalom line D-D?, located on Figure 5.2. Positions inkm are approximate, due to the projection from the FWI profile. The geological units at surface areindicated at the top of the image; Nb: Neogene basalts; Ea: Eocene andesite; Ef: Eocene felsic rocks;
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Developments in waveform tomography of land seismic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Developments in waveform tomography of land seismic data with applications in south-central British Columbia Smithyman, Brendan Robert 2013
pdf
Page Metadata
Item Metadata
Title | Developments in waveform tomography of land seismic data with applications in south-central British Columbia |
Creator |
Smithyman, Brendan Robert |
Publisher | University of British Columbia |
Date Issued | 2013 |
Description | Waveform tomography (WT) is an advanced class of seismic imaging methods that may be applied to generate detailed models of subsurface velocity and attenuation by analysis of recorded field seismograms. WT comprises two steps: 1) initial velocity model building by traveltime inversion (raytracing) and 2) full-waveform inversion to update the model of velocity based on analysis of refracted data waveforms. The nature of this method requires that the initial model and simulated survey parameters reproduce closely those used for the field data acquisition. The application of WT to on-land seismic data is challenged by the geometry irregularities inherent in crooked-line two-dimensional (2D) seismic data acquisition. This thesis develops two methodologies that permit the inversion of on-land crooked-line vibroseis multi-channel seismic (MCS) reflection data by traveltime tomography (TT) and subsequent viscoacoustic full-waveform inversion (FWI) to generate 2D cross sections of P-wave velocity. The first method uses 2.5D TT and applies a TT-derived static correction to enable the subsequent application of 2D FWI. The second method uses 2.5D or three-dimensional (3D) TT followed by 2.5D full-waveform inversion that models 3D survey geometry (a new development). These technical developments are motivated and tested by a multi-part case study, which analyses data from the 2008 Geoscience BC Nechako Basin seismic survey. Both the static-correction method and 2.5D FWI method are applied to a single acquisition line. The results are analysed and compared with the conclusion that 2.5D WT provides a superior result, but at additional computational cost. This cost is still significantly lower than that associated with 3D FWI. The 2.5D WT method is applied to three datasets from the 2008 survey, and TT alone is applied to a fourth. The resulting models constrain the positions and lithology of Eocene volcanic rock units, Cretaceous strata and Jurassic or earlier basement rocks in the Nechako-Chilcotin Plateau region of south-central British Columbia. These results are synthesized and interpreted in conjunction with colocated surficial-geology maps, exploration well logs and geophysical results by other workers. This interpretation enhances the specific knowledge of lithostratigraphy along these acquisition profiles and informs general geological and geophysical analyses of the region. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2013-10-08 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0103357 |
URI | http://hdl.handle.net/2429/45219 |
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 | 2013-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2013_fall_smithyman_brendan.pdf [ 25.92MB ]
- Metadata
- JSON: 24-1.0103357.json
- JSON-LD: 24-1.0103357-ld.json
- RDF/XML (Pretty): 24-1.0103357-rdf.xml
- RDF/JSON: 24-1.0103357-rdf.json
- Turtle: 24-1.0103357-turtle.txt
- N-Triples: 24-1.0103357-rdf-ntriples.txt
- Original Record: 24-1.0103357-source.json
- Full Text
- 24-1.0103357-fulltext.txt
- Citation
- 24-1.0103357.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-0103357/manifest