UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Deep 3-D seismic reflection imaging of Precambrian sills in the crystalline crust of Alberta, Canada Welford, Joanna Kim 2004

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2005-995690.pdf [ 44.83MB ]
Metadata
JSON: 831-1.0052372.json
JSON-LD: 831-1.0052372-ld.json
RDF/XML (Pretty): 831-1.0052372-rdf.xml
RDF/JSON: 831-1.0052372-rdf.json
Turtle: 831-1.0052372-turtle.txt
N-Triples: 831-1.0052372-rdf-ntriples.txt
Original Record: 831-1.0052372-source.json
Full Text
831-1.0052372-fulltext.txt
Citation
831-1.0052372.ris

Full Text

D E E P 3-D S E I S M I C R E F L E C T I O N I M A G I N G O F P R E C A M B R I A N S I L L S IN T H E C R Y S T A L L I N E C R U S T O F A L B E R T A , C A N A D A By JOANNA KIM WELFORD B. Sc., McGill University, 1997 M . Sc., University of British Columbia, 2000  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF EARTH AND OCEAN SCIENCES We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA October 2004 © Joanna Kim Welford, 2004  Abstract  Using deep 3-D seismic reflection datasets collected by the Canadian petroleum exploration industry in southwestern and northwestern Alberta, the Head-Smashed-In and Winagami Precambrian sill complexes within the crystalline upper crust, previously identified on Lithoprobe 2-D multichannel reflection lines, are investigated to determine their 3-D geometries and reflective characteristics. During seismic processing of the dataset in southwestern Alberta, a recently developed wavelet-based method, Physical Wavelet Frame Denoising, is applied and shown to successfully suppress ground roll contamination while preserving low frequency signals from deeper structures. A new 3-D empirical trace interpolation scheme, DSInt, is developed to address the problem of spatial aliasing associated with 3-D data acquisition. Results from applying the algorithm to both datasets are comparable to available interpolation codes while allowing for.greater flexibility in the handling of irregular acquisition geometries and interpolated trace headers. Evidence of the Head-Smashed-In reflector in southwestern Alberta is obtained using a dataset acquired to 8 s T W T T (approx. 24 km depth). From locally coherent, discontinuous pockets of basement reflectivity, the dataset appears to image the tapering western edge of the deep reflections imaged by Lithoprobe. A statistical approach of tracking reflectivity is developed and applied to obtain the spatial and temporal distribution of reflections. Simple 1-D forward modelling results reveal that the brightest reflections likely arise from a 50 to 150 m thick body of high density/high velocity material although variations in the amplitudes and lateral distribution of the reflections indicate that the thickness of the sills is laterally variable. Thus, the results are consistent with imaging the tapering edge of the sill complex. Clear evidence of the Winagami reflection sequence in northwestern Alberta, emerges from  ii  the second dataset acquired to 5.1 s T W T T (approx. 15 km depth). Data sections outline a 3-D reflective sheet dipping to the southeast. From polarity comparisons from within the sedimentary sequence, the reflective signature of the deep body is inferred to result from higher density and higher velocity material than the surrounding host rocks, thus reinforcing the previous interpretation of the deep reflections as resulting from dolerite sills intruded into gneissic crystalline basement. Simple 1-D forward modelling results reveal that the thickness of the sheet is between 50 and 100 m throughout the survey region. From 3-D Kirchhoff forward modelling, the reflective sheet is detected between 11 and 16 km depth along the northwestern edge of the survey area. The absence of the Winagami reflections, and basement reflectivity in general, in the southeast of the survey region coincides with a positive aeromagnetic anomaly inferred to be caused by magmatic rocks. The presence of the magmatic rocks may have either influenced the geometry of the intrusion of the sills or overprinted their reflective signature, dependent upon the relative timing of emplacement of the two features. Both sill complexes appear to have been intruded horizontally into the crust from multiple injections of magma during a period of tectonic compression. The emplacement of these sills may have strengthened the crust and provided the rheological contrasts needed to initiate the formation of Paleozoic cratonic arches like the Peace River Arch of northwestern Alberta and the Montania Arch of southern Alberta. The results from this thesis represent the first opportunity in North America to examine upper-middle crustal structure to depths of approximately 20 km using industry-style 3-D seismic reflection techniques.  iii  Table of Contents  Abstract  ii  List of Tables  ix  List of Figures  x  Acknowledgements  xv  Dedication  xvii  1 Introduction  1  1.1  Upper crustal reflectors  2  1.2  Upper crustal reflectors in Canada  4  1.3  Investigation of the Head-Smashed-In reflector  5  1.4  Investigation of the Winagami sequence  8  1.5  Thesis goals  11  1.6  Thesis outline  12  2 The 3-D Surveys 2.1  2.2  13  The Pincher Creek (PCI) survey in southwestern Alberta  13  2.1.1  Tectonic setting of Pincher Creek, Alberta  13  2.1.2  General sedimentary stratigraphy in the Pincher Creek region  15  2.1.3  The Pincher Creek (PCI) 3-D dataset  19  The Pine Creek (PC2) survey in northwestern Alberta  20  2.2.1  20  Tectonic setting of the Pine Creek region, Alberta iv  2.3  2.2.2  General sedimentary stratigraphy of the Pine Creek region, Alberta . .  23  2.2.3  The Pine Creek (PC2) 3-D dataset  27  Chapter summary  29  3 Seismic Processing  30  3.1  Initial data preparation  32  3.2  Geometry  33  3.2.1  Acquisition of the PCI 3-D survey  34  3.2.2  Acquisition of the PC2 3-D survey  .  36  3.3  Trace editing and geometry quality control  37  3.4  Refraction statics  39  3.4.1  3-D model parameterization  40  3.4.2  Ray-tracing in 3-D  42  3.4.3  Starting models  42  3.4.4  Generalized linear inversion  45  3.4.5  Calculation of refraction statics  46  3.5  Ground roll removal  53  3.5.1  Conventional ground roll removal techniques  53  3.5.2  Physical Wavelet Frame Denoising  56  3.5.3  Application of PWFD  58  3.6  Geometrical spreading correction  64  3.7  Compensating for frequency attenuation  70  3.8  Velocity analysis  71  3.8.1  Velocity analysis in 3-D  73  3.8.2  Semblance analysis  74  3.9  Residual statics  77  v  3.10 Spiking deconvolution  79  3.10.1 Spiking deconvolution in practice  80  3.10.2 Spiking deconvolution for the PCI and PC2 datasets  82  3.11 Trace balance and stack  83  3.12 Migration  83  3.12.1 Challenges related to the migration of deep reflection data  86  3.12.2 Migration in practice  87  3.12.3 Migration of the PCI and PC2 datasets  88  3.13 f-x deconvolution  89  3.14 Seismic resolution  90  3.15 Chapter summary  91  4 Trace interpolation to mitigate spatial aliasing  92  4.1  What is spatial aliasing?  93  4.2  Conventional methods for trace interpolation  96  4.3  Practical issues related to interpolating real 3-D datasets  98  4.3.1  Interpolation codes available for this study  99  4.3.2  Specifying header information for interpolated traces  100  4.3.3  Rationale for a new practical 3-D interpolation scheme  102  4.4  DSInt 3-D interpolation  102  4.4.1  Copy distance as an error proxy  106  4.4.2  Interpolating missing receiver lines: the PCI example  106  4.4.3  Interpolating missing traces along receiver lines: the PC2 example . .  Ill  4.5  Conclusion  114  4.6  Chapter summary  114  vi  5 The Head-Smashed-In Reflector  115  5.1  The SALT'95 data  115  5.2  Pre-stack evidence for deep reflectors in the PCI dataset  118  5.3  Post-stack evidence for deep reflectors in the PCI dataset  118  5.4  Statistical analysis of the deep reflections in the PCI data volume  122  5.5  Waveform modelling  127  5.6  Multiples?  130  5.7  Do the PCI deep reflections represent further evidence of the HSI reflector?  5.8  Chapter summary  .  134 134  6 The Winagami Reflection Sequence  136  6.1  The CAT'92 and PRAISE'94 data  137  6.2  Pre-stack evidence for deep reflectors in the Pine Creek region  139  6.3  Post-stack evidence for deep reflectors in the Pine Creek region  140  6.3.1  149  Migration tests  6.4  Polarity constraints  154  6.5  Waveform modelling  155  6.6  Modelling the Winagami reflection sequence in 3-D  161  6.7  Lack of reflectivity in the southeast  167  6.8  Do the PC2 deep reflections represent further evidence of the Winagami reflec-  6.9  tion sequence?  172  Chapter summary  172  7 Discussion  174  7.1  Role of sills in craton stabilization  174  7.2  Sill emplacement mechanisms  175  7.3  Emplacement of the HSI reflector  177 vii  7.3.1 7.4  178  Emplacement of the Winagami reflection sequence  181  7.4.1  183  Regional tectonic activity during emplacement  7.5  Constraints on the mode of emplacement  186  7.6  Role of sills in cratonic arch formation  188  7.6.1  Rheological considerations  189  7.6.2  Uplift of cratonic arches  194  7.7 8  Regional tectonic activity during emplacement  Chapter summary  196  Conclusions  197  Bibliography  201  viii  List of Tables  2.1  Acquisition parameters for the PCI 3-D dataset in southwestern Alberta. . . .  19  2.2  Acquisition parameters for the PC2 3-D dataset in northwestern Alberta. . . .  27  3.1  Initial 1-D velocity model used in gli3d for the PCI dataset  44  3.2  Initial 1-D velocity model used in gli3d for the PC2 dataset  45  3.3  Statistical distributions of the individual refraction static corrections applied to the PCI and PC2 datasets  50  4.1  Threshold frequency of spatial aliasing for the PCI dataset  96  4.2  Threshold frequency of spatial aliasing for the PC2 dataset  96  7.1  Uniaxial compressive strengths of various rock types at near-surface conditions 190  ix  List of Figures  1.1  Locations of several bright crustal reflections imaged around the world . . . .  3  1.2  Canadian occurrences of bright upper crustal reflections  5  1.3  Unrmgrated stack of the HSI reflector on the SALT seismic lines  6  1.4  Location of the PCI 3-D survey relative to the HSI reflector  7  1.5  The Winagami reflection sequence on the CAT and PRAISE seismic lines  1.6  Winagami reflection sequence cross-cutting the Ksituan-Chinchaga basement  . .  9  domain boundary  10  1.7  Location of the PC2 3-D survey relative to the Winagami reflection sequence .  11  2.1  Basement domain map of the Pincher Creek (PCI) region  14  2.2  Schematic stratigraphic column for the Pincher Creek (PCI) survey  16  2.3  Geological assemblage map for the Pincher Creek (PCI) survey area  18  2.4  Schematic diagram of the Pincher Creek triangle zone  18  2.5  Survey geometry for the PCI 3-D dataset  20  2.6  Source-receiver distributions for all shots from the PCI 3-D dataset  21  2.7  Basement domain map of the Pine Creek (PC2) region  24  2.8  Schematic stratigraphic column for the Pine Creek region  25  2.9  Spatial relationship between the Pine Creek region and the Peace River Arch .  26  2.10 Survey geometry for the PC2 3-D dataset  28  2.11 Source-receiver distributions for all shots from the PC2 3-D dataset  28  3.1  Processing flows applied to the PCI and PC2 datasets  31  3.2  Simple schematic diagrams of shots and receivers for a 3-D shot and a 3-D CMP 34 x  3.3  Midpoint cluster in a bin from the PCI dataset  35  3.4  Stacking fold map of the PCI 3-D dataset  35  3.5  Scatter of midpoints in a 35 m by 280 m bin from the PC2 dataset  37  3.6  Stacking fold maps of the PC2 3-D dataset for different bin sizes  38  3.7  Model parameterization in gli3d  41  3.8  A simple diagram of gli3d's approach to 3-D ray-tracing  43  3.9  Diagram of the parameters used in the refraction statics calculations  47  3.10 Surface topography for the PCI and PC2 surveys  49  3.11 Diagram illustrating the components of the delay time calculation of short wavelength statics  50  3.12 Effect of refraction statics on inline gathers from the PCI dataset  51  3.13 Effect of refraction statics on inline gathers from the PC2 dataset  52  3.14 Shot gather contaminated with ground roll with filtering results  54  3.15 A selection of 2-D physical wavelets used in ground roll suppression  57  3.16 Outline of the PWFD processing sequence  59  3.17 Results of ground roll suppression using PWFD  61  3.18 Further results of ground roll suppression using PWFD  63  3.19 Comparison of geometrical spreading corrections for a far offset trace  . . . .  67  3.20 Comparison of geometrical spreading corrections for a near offset trace . . . .  68  3.21 Effect of i  1 , 5  geometrical spreading correction on inline gathers from the PC2  dataset  69  3.22 N M O explanatory diagram  72  3.23 A synthetic CMP gather and its semblance spectrum  75  3.24 A CMP gather from the PCI dataset and its semblance spectrum  76  3.25 A CMP gather from the PC2 dataset and its semblance spectrum  77  3.26 Effect of residual statics on a CMP gather from the PCI dataset  79  xi  3.27 Simple diagram of the subsurface beneath both surveys  82  3.28 Diagrams showing the effect of migration on a diffraction and on a dipping reflection  84  3.29 Quantitative analysis of the migration process  85  4.1  f-k spectrum for a dipping reflector  94  4.2  Spatial aliasing in the f-k domain  95  4.3  Steps involved in dummy trace seeding of the 3-D datasets  101  4.4  Actual and dummy receiver locations for DSInt interpolation  103  4.5  Schematic diagram of DSInt trace copying  103  4.6  Bin configuration (15 degrees by 240 m) for a shot at the center of the PCI survey 104  4.7  Increase in live traces from DSInt for 15 degree by 120 m bins for shot 200 from the PCI survey  105  4.8  Adapted Megabin geometry of the PCI survey  107  4.9  Live trace concentrations before and after DSInt interpolation for additional receiver inlines  108  4.10 Comparison of DSInt and suinterp results for interpolated receiver inlines  . .  109  4.11 Live trace concentrations before and after DSInt interpolation for missing receivers along inlines  Ill  4.12 Comparison of DSInt and suinterp results for interpolating missing inline receiver traces  113  5.1  Detailed Archean basement domain map of the PCI region  116  5.2  Unmigrated stacked section of line 30 from SALT'95  117  5.3  Inline sections through the stacked PCI data volume  119  5.4  Crossline sections through the stacked PCI data volume  121  5.5  One-to-one comparison between deep reflectors from SALT'95 and PCI . . .  122  xii  5.6  Distribution of common reflector picks as a function of T W T T  5.7  Spatial distribution of the number of reflector picks common to both inline and  124  crossline sections  124  5.8  Four interpreted HSI common reflector surfaces  126  5.9  1-D Modelling results for the PCI dataset  128  5.10 Impact of multiples on deep reflection detection  132  5.11 Impact of sill thickness on reflection detection in the presence of multiples . .  133  6.1  Detailed basement domain map of the PC2 region  138  6.2  Schematic diagram of the anastamosing shape of the Winagami bands of reflectivity  6.3  138  Spatial relationship between the topmost Winagami reflective band on Lithoprobe lines to the north and the PC2 survey  139  6.4  Inline sections through the unmigrated PC2 data volume  143  6.5  Crossline sections through the unmigrated PC2 data volume  144  6.6  Time slices through the unmigrated PC2 data volume  145  6.7  Chair-cut views through the unmigrated PC2 data volume  146  6.8  Spatial tie between time slice (b) and the chair-cut view through the unmigrated PC2 data volume  6.9  147  Spatial tie between time slice (e) and the chair-cut view through the unmigrated PC2 data volume  148  6.10 Inline 2-D post-stack Kirchhoff time migration results for a northern inline section of the PC2 data volume using different velocities  151  6.11 Inline 2-D post-stack Kirchhoff time migration results for a southern inline section of the PC2 data volume using different velocities 6.12 Comparison of basement diffractions from the Lithoprobe and PC2 results . .  xiii  152 153  6.13 Well logs from the overlying sedimentary sequence used to obtain polarity constraints for the WRS  156  6.14 Polarity comparison between synthetic seismograms and the WRS reflections  157  6.15 1-D Modelling results for the PC2 dataset  158  6.16 Comparison of 1-D modelling results with coherent reflections from the PC2 data volume  160  6.17 3-D model of the top of the Winagami sequence from 3-D Kirchhoff modelling 163 6.18 3-D Kirchhoff modelling results for selected stacked inlines  165  6.19 3-D Kirchhoff modelling results for selected stacked crosslines  166  6.20 Spatial link between the lack of reflectivity in the southeast and the aeromagnetic signal  169  6.21 Interpreted aeromagnetic map of northern Alberta  171  6.22 Spatial link between the Winagami reflection sequence and the deep reflections in the PC2 data volume 7.1  173  Schematic representation of the present tectonic configuration in southern A l berta and northern Montana  7.2  178  Representation of the tectonic evolution of the Medicine Hat Block and its surroundings  179  7.3  Temporal constraints on the emplacement of the HSI sills  180  7.4  Diagram of dike and sill propagation in response to localized stresses  182  7.5  Tectonic evolution of northwestern Laurentia  185  7.6  Surface deflections resulting from the juxtaposition of strong and weak crust .  190  7.7  Inferred uniaxial compressive strength of dolerite and gneiss at depth  192  7.8  Schematic diagram of the folding of sills of different thicknesses  7.9  Suggested spatial connection between the Peace River Arch and the Ksituan arc 195  xiv  .  193  Acknowledgements  First and foremost, I would like to thank my supervisor, Professor Ron Clowes, for the wonderful advice he has given me throughout my graduate school career and for allowing me the freedom to pave my own way in research. His trust in my abilities gave me the confidence to pursue whatever research tangents struck me as interesting. I would also like to thank the members of my doctoral supervisory committee, Professors Lori Kennedy, Tad Ulrych and Andrew Calvert. Their unwavering faith and constant encouragement provided the foundation on which this thesis was built.  Without the help of Jonathan Ravens from the Institute of Geological and Nuclear Sciences in New Zealand, this thesis would have been lost long ago. Jonathan's willingness to answer my questions and work through my processing crises at all times of day and night were invaluable to my research. Derek Woodward and the other GNS workers also deserve great credit for patiently allowing me to overload their inboxes while still offering helpful advice.  Special thanks to EnCana Corporation in Calgary for donating the PCI dataset to UBC, for providing part of the financial support for my Natural Sciences and Engineering Research Council (NSERC) Industrial Post-graduate Scholarship and for allowing me to work on site for almost six months in early 2002. My company supervisor at EnCana, Janet Porter-Chaudhry, deserves particular recognition for her dedication to helping me with my thesis both at EnCana and via email and for her efforts to make me feel at home in Calgary.  I am very grateful to BP Canada, particularly to Dr. Tom Koleszar, for donating the PC2  xv  dataset to U B C . Their generous bequest provided the necessary balance to this thesis and allowed me to hone my 3-D seismic processing skills. Thanks also to Kelman Archives in Calgary who ultimately provided the data to U B C . LogTech Canada Limited graciously donated the well logs for the PC2 survey area.  I am thankful to Professor David Eaton from the University of Western Ontario for providing the 3-D modelling code used in Chapter 6.  I am truly blessed to have been a member of the Geophysics group in the Department of Earth and Ocean Sciences at U B C for so many years. Throughout my tenure in the department, I have had the privilege of crossing paths with many spectacular people. Of those I hold dearest to my heart, I would particularly like to acknowledge Phil Hammer, Gwenn Flowers, Kris Innanen, Megan Sheffer, Kate Bone, Andrew Gorman, John Amor, Stephane Rondenay, Camille Li, Charly Bank and Gaby Fernandez-Viejo. The entire Rock Dogs women's ice hockey team, past and present, are also to be thanked for teaching me how to work as part of a team and for helping me to focus my thesis frustrations into winning hockey games.  Many thanks to my family back in Montreal for all their encouragement over the years and to my best friend, Jennifer Taylor (nee Davies), for always being my biggest fan.  Lastly, thanks to my wonderful husband, Colin Farquharson, for his constant love and support. He has stayed by my side through two graduate degrees and has resisted the urge to run away from all my complaining. I promise to be a lot easier to live with from now on.  xvi  To Colin  xvii  C  h  a  p  t  e  r1  I n t r o d u c t i o n  As a blank canvas, the Earth's crust materialized from the differentiation of our planet over four billion years ago. Since then, tectonic processes such as collision, accretion and extension, have transfigured its surface with features such as mountains, volcanoes and ocean basins. While many tectonic processes continue to actively shape many parts of the planet today, similar and perhaps different processes were active in the Earth's tectonic past. The superposition of many of these tectonic influences is recorded within our present-day crust. The North American landmass, with its extensive temporal and spatial record of a vast array of tectonic forcings, represents an ideal research environment for unraveling the complex evolution of our planet. From the high peaks of the upthrust Rocky Mountains to the exposed bedrock of the ancient Canadian Shield, areal geological mapping with its associated dating and chemical studies can provide great insight into many of these evolutionary processes. Still, lacking the depth dimension, our understanding remains incomplete. While surface outcrops of deeper crustal structures, drillhole information and xenoliths that have sampled the crust and mantle often can provide a glimpse into the third dimension, deeper crustal deformation generally remains unknown. Seismic reflection techniques provide one means by which deep structures can be probed. The use of seismic techniques in the imaging of deep crustal structure has been actively utilized for over fifty years; however, due to both acquisition costs and computer limitations, these studies have generally only involved the acquisition and interpretation of 2-D seismic lines. While 3-D near-vertical incidence surveys have been in use in the oil industry for the  1  Chapter 1. Introduction  2  last thirty years, 3-D experiments to delineate crustal structure have been undertaken for only the last fifteen years. Most 3-D crustal studies to date have primarily focused on wide-angle refraction techniques (Milkereit et al., 1986; Kanasewich et al., 1987) with broadside shooting into 2-D receiver lines of varying orientation. Attempts at near-vertical incidence 3-D experiments for crustal studies have been undertaken to a far lesser extent. In fact, these efforts can be better described as quasi-3-D experiments as they have only involved broadside shooting into intersecting 2-D profiles (Cook and Coflin, 1990; Kanasewich et al., 1995). To date, there has been only one industry-style 3-D experiment collected for crustal studies, the German Integrated Seismics Oberpfalz (ISO 89) experiment around the K T B borehole (Stiller, 1991). This experiment was the world's first high resolution industry-style 3-D survey intended to image the crystalline upper crust. The lack of deep 3-D surveys such as ISO 89 can be attributed partly to the significantly higher acquisition costs of 3-D and partly to the fact that unlike sedimentary basins, crystalline crust generally shows lower reflectivity. Nonetheless, there are areas around the world where crustal 2-D seismic reflection studies have identified reflectors within the crystalline crust that could be good targets for 3-D investigations.  1.1  U p p e r crustal reflectors  From large-scale crustal investigations such as Canada's Lithoprobe, Britain's BIRPS, the U.S.'s COCORP and other similar projects, a wealth of deep crustal images of different tectonic environments from all over the world now exists. A curious feature that has appeared on many of these deep survey results is the presence of bright upper crustal reflections between 5 and 25 km in otherwise homogeneous crust (Fig. 1.1) (Hyndman, 1988; Korsch et al., 1992; Litak and Hauser, 1992; Jarchow etal., 1993; Pratt etal., 1993; Brown etal., 1996; Law and Snyder, 1997; Mandler and Clowes, 1997; Papasikas and Juhlin, 1997; Ross and Eaton, 1997; Mandler and  Chapter 1. Introduction  3  Clowes, 1998). While these reflections were initially thought to be quite anomalous, an increasing number of deep seismic reflection datasets have shown that they are more common than originally thought. In fact, some geoscientists believe that these reflectors exist everywhere in the crust and that their detection has been limited solely by acquisition constraints (D. Eaton, personal communication, 2001).  Figure 1.1: Locations of several bright crustal reflections imaged around the world are shown as dots.  Several causes have been proposed to explain bright upper crustal reflectivity. These include intrusive mafic sills (Litak and Hauser, 1992; Pratt et al., 1993; Mandler and Clowes, 1997; Papasikas and Juhlin, 1997; Ross and Eaton, 1997; Mandler and Clowes, 1998), mylonitic shear zones (Korsch et al., 1992; Law and Snyder, 1997), melt bodies (Jarchow et al., 1993; Brown et al., 1996) and fluid-filled porous zones (Hyndman, 1988). While some of these causes have been confirmed through drilling (Juhlin, 1990) and others inferred from nearby surface outcrops (Mandler and Clowes, 1997), most interpretations of upper crustal reflectivity are based on detailed seismic analyses combined with results from lab studies and other geological and geophysical methods. The character of the seismic response from a deep reflector contains information about the contrast in physical properties that generated the reflection. Consequently, examination of the seismic signature can provide lithological constraints on the reflector itself. Since both melt  Chapter 1. Introduction  4  bodies andfluid-filledporous zones are generally less dense and have lower velocities than their hosts, their seismic response will have an opposite polarity to reflections generated by high velocity and high density mafic sills. Meanwhile, since mylonitic zones can produce seismic responses similar to both melts and sills (Law and Snyder, 1997), shear zone interpretations for deep reflections are generally argued based on geometrical relationships and regional tectonics.  1.2  U p p e r crustal reflectors i n C a n a d a  Anomalous upper crustal reflectors have been detected on a number of 2-D Lithoprobe multichannel reflection lines (Fig. 1.2). While the Wollaston Lake reflector of northern Saskatchewan has been imaged on one seismic line (Mandler and Clowes, 1997), both the Winagami reflection sequence of northwestern Alberta (Ross and Eaton, 1997) and the Head-Smashed-In (HSI) reflector in southwestern Alberta (Mandler and Clowes, 1998) have been recorded along multiple seismic lines. Polarity constraints for both the Wollaston Lake and the HSI reflectors indicate that the reflections are originating from high velocity and high density reflective bodies. While such results can be consistent with both mafic sills and mylonitic shear zones, surface outcrops of sills near the Wollaston Lake reflector support its interpretation as a sill complex. Based on similarities in reflection geometries and characteristics among all three reflectors, previous investigators have extended the mafic sill interpretation to all three of these examples of upper crustal reflectivity in western Canada. Since sill emplacement is believed to be a fundamental step in the process of cratonization (Nelson, 1991), studying the nature of these inferred sills and determining the ancient stress environment within which they were emplaced can provide insight into the evolution of the Canadian craton. Two of these sequences were imaged within the crystalline basement below the Western Canada Sedimentary Basin in Alberta, a locale of extensive seismic exploration for oil and gas. The present focus on 3-D seismic exploration by the petroleum industry provides opportunities  Chapter 1. Introduction  5  Figure 1.2: Proposed areal extent of bright upper crustal reflections in Canada as constrained by overlain 2-D Lithoprobe multichannel reflection lines. to further probe these unusual upper crustal features if the surveys were recorded over them and to sufficiently long recording times. Fortunately, two such datasets have been obtained from industry by U B C ; hence my work focuses on achieving this task using both a dataset acquired near the HSI reflector and another larger dataset that overlies the proposed extent of the Winagami. My efforts represent the first opportunity in North America to examine upper-middle crustal structure to depths of approximately 20 km using industry-style 3-D seismic reflection techniques.  1.3 Investigation of the Head-Smashed-In reflector The Head-Smashed-In (HSI) reflector, so named due to its proximity to the Head-Smashed-In Buffalo Jump UNESCO World Heritage site, was imaged along three of the 2-D multichannel  Chapter 1. Introduction  6  crustal seismic reflection lines collected by Lithoprobe as part of the Southern Alberta Lithospheric Transect (SALT) in 1995 (Fig. 1.3). This sequence is characterized by between one and four bands of bright reflections with slightly varying dips confined within a depth range of 6 to 21 km (2-7 seconds two-way-traveltime (TWTT)). Based on seismic coverage to the north and to the east, the reflections were attributed to a large reflective body with an areal extent estimated to be at least 6000 km (Mandler and Clowes, 1998). 2  SALT line 32  Figure 1.3:  Unmigrated stack of the HSI reflector on the S A L T seismic lines. On all seismic plots, approximate  depths are based on an average velocity of 6 km/s. Since velocities in the sediments are slower, the base of the sediments (which occurs between 1.5 to 3 s T W T T on the plots) isn't as deep as shown. The displayed seismic profiles are highlighted in black on the location plot.  The character of the HSI reflections varies from monocyclic to multicyclic across the SALT  Chapter 1.  Introduction  1  seismic lines. Mandler and Clowes (1998) interpret this to indicate that the reflective body is laterally heterogeneous with afinitethickness. By correlating the polarity of thefirstamplitudes from the HSI with that of reflections of known impedance contrast within the overlying sedimentary sequence, they conclude that the HSI reflections result from a high-velocity/density structure. Due to the inconsistent variability in dips across the HSI, this high-velocity/density structure is attributed to the presence of mafic sills and not to an organized mylonitic shear zone (Mandler and Clowes, 1998). With the HSFs intriguing reflective character and location within the upper crust, a unique opportunity presented itself in 1998 when it was learned that EnCana Corporation and its partners Anadarko Petroleum Corporation and Paramount Resources Limited would be acquiring a 3D seismic reflection dataset 30 km to the southwest of the location where Lithoprobe recorded the HSI reflections (Fig. 1.4). Negotiations were conducted between U B C and the industry partners. Consequently, in January 1999, the data were recorded to 8 s TWTT, almost 5 s longer than the planned acquisition, to enable academic studies to be conducted on the deep dataset. This Pincher Creek dataset (referred to hereafter as PCI for clarity), which was subsequently donated to U B C and formed the starting point for my doctoral project, spanned an area of 6 km by 9 km. The main goal for this dataset is to detect the HSI and if present, attempt to extract information about its 3-D geometry and characteristics.  Figure 1.4:  Location of the PC  1 3-D survey relative to the HSI  constrain the areal extent of the HSI are overlain for reference.  reflector. The numbered S A L T seismic lines that  Chapter 1. Introduction  1.4  8  Investigation of the W i n a g a m i sequence  The Winagami reflector sequence of northwestern Alberta, named for a nearby town, was first imaged along four seismic reflection profiles acquired during Lithoprobe's Central Alberta Transect (CAT) experiment of 1992 (Fig. 1.5). In 1994, the sequence was seen again along eight of the profiles acquired during Lithoprobe's Peace River Arch Industry Seismic Experiment (PRAISE) in northwestern Alberta (Fig. 1.5). Along the different seismic lines, the sequence is characterized by between one and five bands of bright reflections within a depth range of 9 to 24 km (3-8 s (TWTT)). Locally the reflections tend to be horizontal to sub-horizontal but together they appear to form a large-scale anastomosing fan that converges to the southeast (Ross and Eaton, 1997). From the seismic coverage, this extensive reflective body appears to cover an area of more than 120,000 km (Ross and Eaton, 1997). 2  The Winagami reflections are consistently characterized by a simple waveform with a dominant frequency of 20 Hz (Ross and Eaton, 1997). For the uppermost reflection, this waveform appears to be a composite waveform resulting from the superposition of reflections from both the top and bottom of the reflective layer. From forward seismic modelling using a horizontal diabase sill intruded into a granitic background, Ross and Eaton (1997) estimate that this uppermost reflection arises from a horizontal sill of 70 ± 10 m thickness. Their sill interpretation is based on evidence that the Winagami reflections crosscut and overprint dipping reflectivity within the granitic host rocks (Fig. 1.6).  Chapter 1.  Introduction  PRAISE line 14 0  - , - - .....  > XJ  6  3  x  CD T3  18 -T •> 10 km 24  PRAISE line 20 S NE  N  3  SW Winagami reflection sequence  10 km  *  CAT line 1  N  S  V.E.2:1 Figure 1.5: Unmigrated stack of the Winagami reflection sequence on the CAT and PRAISEseismic lines. On all plots, approximate depths are based on an average velocity of 6 km/s. Since velocities inthe sediments are as deep as shown, slower, the base of the sediments (which occurs between 1.8 and 2.2 s TWTT on the plots) isn't displayed seismic The seismic lines on which the Winagami is seen are overlain in white on the location plot. The profiles are highlighted in black on the location plot.  Chapter 1.  10  Introduction  PRAISE line 12  Winagami reflection sequence  Figure 1.6: Migrated stack of the Winagami reflection sequence (WRS) cross-cutting the Ksituan-Chinchaga basement domain boundary on PRAISE line 12. On all seismic plots, approximate depths are based on an average velocity of 6 km/s. Since velocities in the sediments are slower, the base of the sediments (which occurs between 1.8 and 2.2 s TWTT on the plots) isn't as deep as shown. The PRAISE and CAT lines on which the Winagami is seen are overlain in white on the location plot. The displayed seismic profile is highlighted in black on the location plot.  The procurement of an industry 3-D dataset that overlapped the Winagami sequence and probed to sufficient depth was achieved following a meeting with BP Canada in May of 2002. In this instance, the area of interest was presented and appropriate datasets were isolated. As a result, the Pine Creek 3-D dataset (referred to hereafter as PC2 for clarity) was released by BP Canada in May 2003 and was ultimately provided by Kelman Archives in June 2003. Since this dataset was acquired strictly for exploration purposes, recording times were not as long as the PCI dataset. Still, with 5.1 s T W T T of recording, it was hoped that at least the top half of the Winagami sequence could be imaged and mapped in 3-D. This dataset was significantly larger than the PCI dataset both in terms of the number of shots and receivers but also with  Chapter 1. Introduction  11  respect to areal extent (18 km by 22 km). Goals for this dataset are much like those for PCI but the location of the survey immediately above the proposed extent of the Winagami (Fig. 1.7) suggests greater chances of detecting and analyzing the sequence.  Winagami . reflection ^.. \ s e q u e n c e  Figure 1.7: Location of the PC2 3-D survey relative to the Winagami reflection sequence. The PRAISE and CAT seismic lines that constrain the areal extent of the Winagami are overlain for reference.  1.5  Thesis goals  The fundamental goal of this thesis is to evaluate the applicability of industry-style 3-D seismic reflection datasets to upper crustal studies. To this end, I am employing two 3-D datasets collected by the Canadian petroleum industry to investigate inferred sill complexes within the crystalline crust of Alberta, Canada. These sill complexes were previously identified on 2-D multichannel reflection profiles collected by Lithoprobe. Provided the sill complexes are detectable on the 3-D data, further goals of this thesis are:  • to map the sill geometries in 3-D • to exploit reflection polarity relationships within the sedimentary sequence to constrain the possible causes of the deep reflections • to obtain constraints on sill thicknesses from 1-D forward modelling of their seismic signatures  Chapter 1. Introduction  12  • to further constrain the sill geometries through forward 3-D modelling • to examine the role of the interpreted sill complexes in the tectonic history of their respective regions • to investigate possible sill emplacement mechanisms based on the imaged geometries of the sills. As mentioned, my efforts represent the first opportunity in North America to examine uppermiddle crustal structure to depths of approximately 20 km using industry-style 3-D seismic reflection techniques.  1.6  Thesis outline  Chapter 2 of this thesis introduces the two industry 3-D datasets in greater detail and outlines the geological and tectonic settings in which they were acquired. Chapter 3 outlines the steps involved in the processing of both datasets. In addition, the issue of ground roll removal from the PCI dataset is addressed. Chapter 4 presents a new technique for data interpolation of 3-D data and contrasts the two acquisition geometries dealt with in this thesis. In Chapters 5 and 6 respectively, the results from the PCI and PC2 datasets are presented. The 3-D results from this thesis along with earlier Lithoprobe results are used to discuss the topics of sill emplacement, craton stabilization and cratonic arch formation in Chapter 7. Finally, Chapter 8 presents the final conclusions from this thesis.  Chapter 2 The 3-D Surveys  Acquisition of 3-D seismic reflection surveys is a costly undertaking. To cover a comparable area to a series of 2-D surveys, a larger number of sources and receivers is required and the sources used are often explosive and therefore drilled into place rather than surficial like vibrating sources. While these extra expenses prove to be money well spent for the oil industry, extending 3-D surveys for academic crustal studies, which require many more sources and receivers to provide adequate aperture for deeper investigations, becomes prohibitively expensive. As such and for now, 3-D academic studies must rely on the availability of datasets from the oil industry in regions of interest and these will be limited to investigations of upper crustal structures. Given these constraints, probing anomalous bright reflectors within the upper crystalline crust of Alberta, a locale of extensive seismic exploration for oil and gas, provides a unique opportunity for using 3-D data for upper crustal investigations.  2.1 The Pincher Creek (PCI) survey in southwestern Alberta 2.1.1 Tectonic setting of Pincher Creek, Alberta Pincher Creek is located in southwestern Alberta within the foothills of the Rocky Mountains. As the region lies beneath the thickened southwest portion of the Western Canada Sedimentary Basin, it has been the focus of oil exploration since the early 20th century (Newson, 2001). While the PCI 3-D dataset was initially intended to look solely at the sedimentary sequence, EnCana and its partners agreed to collect extra data to allow academic studies to be conducted  13  Chapter 2. The 3-D Surveys  14  at UBC which would focus on the crystalline basement and its interactions with the overlying sediments.  •113°  -112°  -110°  Archean crust Loverna Block  Vulcan structure  Medicine Hat Block  Figure 2.1: Basement domain map of the region surrounding the Pincher Creek 3-D survey area. The locations of the SALT'95 2-D reflection lines collected by Lithoprobe are also shown for reference. Note that the basement domains have been delineated based on my own interpretation of the aeromagnetic data.  Beneath the Western Canada Sedimentary Basin lie the Precambrian basement domains that together form the westward extension of the exposed Canadian Shield (Fig. 2.1). Since few of these domains are exposed at the surface, their boundaries have been determined on the basis of aeromagnetic data and K-Ar and U-Pb dating of a limited number of drill cores that intersect basement (Burwash et al., 1962; Ross et al., 1991; Villeneuve et al., 1993; Pilkington et al., 2000). Most of the domains making up the Precambrian basement in southern Alberta are Archean blocks whose amalgamation history into the Archean Hearne province remains  Chapter 2. The 3-D Surveys  15  an area of active research. While Hoffman (1988) suggests that the blocks were sutured together along Proterozoic orogenic belts by 1.78 Ga, more recent studies (Clowes et al., 2002; Gorman et al., 2002) propose that their assembly into western Laurentia was complete by the end of the Archean. The PCI 3-D survey region lies directly above one of these Archean domains, the Medicine Hat Block. This block, which is composed of granitic gneisses, has been dated at 2.6-3.3 Ga (Villeneuve et al., 1993). From deep imaging on seismic reflection data from Lithoprobe's SALT'95 experiment, Lemieux et al. (2000) interpret that the Medicine Hat Block formed as a result of an Archean collisional orogen involving two Archean domains. This collisional orogen accommodated an indeterminate amount of eastward and northeastward shortening. To date, most of the evidence for the HSI reflector indicates that it is confined within the Medicine Hat Block. The nature of the collision between the Medicine Hat Block and the Archean Loverna Block to the north remains a controversial issue with some suggesting the collision involved a southdipping subduction zone (Eaton et al., 1999a) and others a north-dipping subduction zone (Clowes et al., 2002; Gorman et al., 2002), both corresponding to the Vulcan Structure. From dating of basement drill core samples, this collision is thought to have occurred between 2.1 and 1.8 Ga (Burwash et al., 1962; Villeneuve et al., 1993). Due to lack of evidence for the HSI crossing into the Vulcan Structure, the sills are interpreted to have intruded into the Medicine Hat Block prior to 1.8 Ga (Mandler and Clowes, 1998).  2.1.2 General sedimentary stratigraphy in the Pincher Creek region At depth, the boundary between the Cambrian sequence and the Precambrian basement of the Medicine Hat Block corresponds roughly with a strong coherent reflector at 2.4 s T W T T (~ 5.5  Chapter 2. The 3-D Surveys  16  km depth for an average velocity of 4.5 km/s in the sedimentary sequence). This depth to basement is consistent with the basement reflector at 2.6 s T W T T seen on Line 30 of the Lithoprobe Southern Alberta Lithospheric Transect (SALT'95) data, 30 km to the north (Lemieux et al., 2000).  Upper Cretaceous  -500 -1000 97.5  CO  £ -1500 E  c -2000 o o  I  Lower Cretaceous  -2500  144 Middle Jurassic  -3000 -3500  208 320  Mlsslsslpplan 360 Middle and Upper Devonian 387 505 545  -4000 -4500  PrecomDrtan casement  2600  Figure 2.2: General schematic stratigraphic column for the Pincher Creek (PCI) survey area showing the dominant lithologies overlying the Precambrian Medicine Hat Block for each geological period. The depths to the stratigraphic boundaries are extrapolated from the seismic reflection results of Lemieux (1999) to the north and are adjusted relative to the basement reflector depth seen on the PCI 3-D data (assuming an average velocity of 4.5 km/s for the sedimentary sequence). The stratigraphic column is adapted from Lemieux's (1999) presentation of a stratigraphic correlation chart from Core Laboratories in Calgary.  -5000 Shales  I  Sandstones/si Its tones H M  ~\J^  r  >  Carbonates  Unconformity  A general schematic stratigraphic column for the PCI survey area is presented in Figure 2.2; it shows the dominant lithologies overlying the Precambrian Medicine Hat Block for each geological period. Three distinct pulses of sedimentation can be interpreted from the two unconformities in the stratigraphic record. Lemieux (1999) summarizes these three pulses as corresponding to 1) superposed rift and passive margin sequences from the late Precambrian to the  Chapter 2. The 3-D Surveys  17  Devonian, 2) subsidence from the Devonian to the Triassic and 3) formation of a foreland basin from the Jurassic to the Early Tertiary in response to the development of the mountainous Foreland Belt of the Canadian Cordillera. The initial pulse of sedimentation in the Pincher Creek area was halted by the uplift of Montania, a paleotopographic high which existed from the Upper Cambrian until the Upper Devonian (Deiss, 1941). Montania roughly corresponds geographically with the Medicine Hat Block (Lemieux et al., 2000). Due to the hiatus in sedimentation during Montania's uplift, there is a complete lack of Ordovician and Silurian sediments south of 51° N in the foreland basin of the Canadian Rockies and within the state of Montana (Deiss, 1941). The subsidence of the Pincher Creek region during the Devonian caused the re-submergence of Montania and as such, a resumption in sedimentation. This pulse, which is evidenced by Mississippian sediments, slowly ceased with the demise of passive margin sedimentation during the Permian (Kent, 1994). The final pulse of sedimentation, initiated in the Jurassic, resulted from the formation of the Foreland Belt to the west and represents the thickest package of sediments within the stratigraphic column. This predominantly clastic Cretaceous package contains some of the youngest stratigraphic producers of hydrocarbons such as the sandstone Cardium Formation of the Upper Cretaceous (Newson, 2001). In terms of surface geology, the PCI 3-D survey straddles Mesozoic and Cenozoic sedimentary strata (Fig. 2.3). Due to shortening from the Cordilleran orogen, these strata have been upthrust northeastward along low-angle thrust faults and are exposed as imbricated Jurassic to Cretaceous rocks (Hiebert and Spratt, 1996).  18  Chapter 2. The 3-D Surveys  Figure 2.3: Geological assemblage map of the region surrounding the Pincher Creek (PCI) 3-D survey area. The surface geology is dominated by clastic Mesozoic and Cenozoic sedimentary rocks of the Western Canada Sedimentary Basin.  Extent of Pincher Creek Survey  10 km  Figure 2.4: Schematic diagram of the Pincher Creek triangle zone adapted from Begin et al. (1996), from Spratt and Lawton (1996) and from Hiebert and Spratt (1996). Dark lines delineate thrust faults whiledashed lines show stratigraphic boundaries. The relative east-west extent of the PCI 3-D survey is shown for reference.  Chapter 2. The 3-D Surveys  19  Structurally, the Pincher Creek region is located above the triangle zone that marks the eastern edge of the Foothills thrust belt of the southern Canadian Cordillera (Jones, 1996). A triangle zone is characterized by foreland directed thrust faults interacting with backthrusts. A schematic diagram of the Pincher Creek triangle zone is shown in Figure 2.4. As a consequence, the Pincher Creek region is one of complex subsurface structures.  2.1.3  T h e P i n c h e r C r e e k ( P C I ) 3-D dataset  The Pincher Creek 3-D dataset was collected by Venture Seismic Limited in January 1999 for EnCana Corporation and its partners. The acquisition parameters for this survey are presented in Table 2.1.  Total number of shots 452 2-4 kg Dynamite Shot size/source 4-18 m Shot depth 240 m Shot interval 120 m Receiver group interval 240-480 m Receiver line spacing 23 Maximum number of live receiver lines 910 Maximum number of live channels 6 Geophones per receiver station 4m Geophone spacing 2 msec Sample rate 8 sec Record length 8064 m Maximum source-receiver offset 3046 m Average source-receiver offset 66° Azimuth of source/receiver lines Table 2.1: Acquisition parameters for the PCI 3-D dataset in southwestern Alberta.  The geometry for the acquisition of the PCI 3-D dataset, which is shown in Figure 2.5, is a slightly altered version of the Megabin geometry pioneered by EnCana. Shots were recorded by up to 23 receiver lines that were spaced 240 to 480 meters apart. Along each receiver line, 19 to 55 receiver groups were connected at a spacing of 120 m. Shots at the NNW and SSE ends  Chapter 2. The 3-D Surveys  20  of the survey area were recorded by 12 receiver lines while more central shots were recorded by the entire spread. As such, each shot was recorded by between 460 and 910 receiver groups. While the maximum source-receiver offset recorded was 8064 m, the average offset was 3046 m as shown in the source-receiver offset distribution plot (Fig. 2.6a). The azimuthal coverage from the survey was equally distributed although certain azimuths were better sampled than others (Fig. 2.6b).  Figure 2.5: Survey geometry for the PCI 3-D dataset.  2.2 2.2.1  T h e Pine Creek (PC2) survey i n northwestern A l b e r t a Tectonic setting of the Pine Creek region, A l b e r t a  The Pine Creek region of northwestern Alberta lies above the Western Canada Sedimentary Basin, east of the Cordilleran deformation front. As in the south, the area has been the focus of  21  Chapter 2. The 3-D Surveys  N  o  Offset (m)  Number of source-receiver pairs  Figure 2.6: Source-receiver distributions for all shots from the PCI 3-D dataset; a) histogram of the source-receiver offsets, b) angle histogram of the source-receiver azimuths. oil exploration for several decades and as such, has significant seismic coverage. While most of these surveys are too shallow for crustal studies, the PC2 3-D dataset collected by BP Canada was recorded to 5.1 s T W T T (~ 15 km depth) making it appropriate for upper crustal investigations. The basement domains that lie beneath the Western Canada Sedimentary Basin in the Pine Creek region are more numerous and somewhat younger than those in southwestern Alberta (Fig. 2.1). Similar to the south, few of these domains are exposed at the surface, so their boundaries have been determined on the basis of aeromagnetic data and K-Ar and U-Pb dating of drill cores that intersect basement (Burwash et al., 1962; Ross et al., 1991; Villeneuve et al., 1993; Pilkington et al., 2000). Most of the domains in the region are Proterozoic with those most relevant to the PC2 survey, the Wabamun and the Chinchaga, ranging in age from 2.0 to 2.4 Ga (Ross, 2002). From overlapping crystallization ages of 2.19 to 2.08 Ga, the Chinchaga and Buffalo Head  Chapter 2. The 3-D Surveys  22  domains are thought to have acted as a distinct entity during the accretionary history of northern Alberta (Ross and Eaton, 2002). As both their Nd isotopic signatures show an Archean component (Theriault and Ross, 1991), they are interpreted to have originated as rifted slivers of the Rae province to the northeast (Ross and Eaton, 2002). By 1.98 Ga, the Chinchaga - Buffalo Head unit was flanked on both sides by outward-dipping subduction zones corresponding to the Ksituan magmatic zone in the west (1.90-1.98 Ga) and the Thelon-Taltson magmatic zone in the east (1.92-1.98 Ga) (Ross, 2002; Ross and Eaton, 2002). The closing of these corresponding ocean basins stitched these domains to the Rae province and contributed to the growth of western Laurentia in northern Alberta. The Snowbird Tectonic Zone (STZ), which extends to the northeast in Figure 2.7, is an enigmatic crustal feature that spans from Hudson Bay to the foothills of the Cordillera. Where it is exposed as high grade mylonitic rocks that outline crustal-scale augen, the zone marks the division between the Rae and Hearne Archean provinces (Ross et al., 1991). In Alberta, the zone bifurcates into a northern zone that marks the northern edge of the Wabamun domain and a southern zone that is associated with the Thorsby domain. The northern zone truncates prominent north-south aeromagnetic patterns seen in the Chinchaga and Buffalo Head domains (Ross et al., 1991). While the Wabamun represents one of these crustal-scale augen, the Thorsby domain is interpreted to have developed as an extensional basin from 1920 to 1850 Ma in response to the combined influences of the termination of strike-slip faults and the counterclockwise rotation of the Hearne province due to plate consumption in the Trans-Hudson Orogen to the southeast (Ross et al., 2000; Ross, 2002). The subsequent closing of this basin is evidenced by the Rimbey magmatic belt (1.79 to 1.85 Ga) and the accompanying Lacombe metasediments which resulted from southeastward subduction of the Thorsby beneath the Hearne province.  Chapter 2. The 3-D Surveys  23  From the Lithoprobe deep seismic results, the Winagami reflection sequence appears to extend through several Proterozoic crustal domains. From cross-cutting relationships with dipping crustal fabric and reflection offsets from aged brittle events, the emplacement of the Winagami reflection sequence is thought to have occurred between 1890 and 1760 Ma (Ross and Eaton, 1997). As such, it post-dates the assembly of most of the domains in the study region and was perhaps coeval with tectonism associated with the opening and closing of the Thorsby basin to the southeast. This may explain the truncation of Winagami-type reflections further to the southeast. The PC2 3-D survey straddles the STZ between the Wabamun and the Chinchaga domains. While no Lithoprobe lines sample this particular boundary elsewhere, the PRAISE lines to the east image the transition from the Buffalo Head domain to the Wabamun across the STZ as a south-dipping detachment (Ross and Eaton, 2002). Coincidentally, this boundary corresponds with the northward truncation of the inferred Winagami reflections in the area. As to whether the STZ beneath the PC2 dataset will reveal a similar detachment is uncertain. However, evidence for the Winagami in both the Chinchaga and Wabamun domains suggests that these reflections may cut through the 3-D survey region.  2.2.2  G e n e r a l sedimentary stratigraphy of the Pine Creek region, A l b e r t a  A general schematic stratigraphic column for the Pine Creek area is presented showing the dominant lithologies overlying the Precambrian basement for each geological period (Fig. 2.8). As in the south, three main pulses of sedimentation, separated by unconformities, can be deduced. These include 1) superposed rift and passive margin sequences from the late Precambrian to the Cambrian, 2) platformal sequences deposited within an inland sea from the Devonian to the Triassic and 3) formation of a foreland basin from the Cretaceous to the Early Tertiary in response to the development of the mountainous Foreland Belt of the Canadian Cordillera to the west.  24  Chapter 2. The 3-D Surveys  -120"  -118'  Proterozoic crust  | Buffalo Head ~| Chinchaga  •  Wabamun Thorsby  -116" Proterozoic magmatic belts  -114'  -112*  -110"  Proterozoic metasediments  Lacombe  Ksituan Thelon-Taltson Rimbey  Archean crust  Hearne  LSi!,  Rae  Figure 2.7: Basement domain map of the region surrounding the PC2 3-D survey area. The locations of the CAT and PRAISE 2-D reflection lines collected by Lithoprobe are also shown for reference. The Snowbird Tectonic Zone and its inferred bifurcation are plotted in the northeast of the map. Note that the basement domains have been delineated based on my own interpretation of the aeromagnetic data.  25  Chapter 2. The 3-D Surveys  Age (m.y.)  Dominant Lithology  58  1000 500  Tertiary and Upper Cretaceous  0  84  -500 Middle Cretaceous  -1000  104  -1500  Lower Cretaceous  j"u?assic .y-~-s,  -/*""~\  ./"'~\.  ./  Misslsslpplan  -2000 -2500  Upper Devonian •  Middle Devonian Lower to Middle Devonian  *~ \^ ^/ ^\^^ f \_- y r  -3000  t  V. :  119 350 360 369 373  Figure 2.8: General schematic stratigraphic column for the Pine Creek region showing the dominant lithologies overlying the Precambrian basement for each geological period. The depths to the stratigraphic boundaries were obtained from the Atlas of the Western Canada Sedimentary Basin (Wright et al., 1994). A stratigraphic correlation chart from Core Laboratories in Calgary was used to obtain the approximate ages of the boundaries.  407  y *"v__>505 545  -3500  Precambrian basement  ~~\ Shales  \ I r  \ y  Sandstones/slltstones s  '  2080  Carbonates  Unconformity  As in the south, the first pulse of sedimentation was halted by the uplift of a cratonic block. In northwestern Alberta, this block corresponds to the Peace River Arch (PRA), a long-lived structure (outlined in Fig. 2.9) that acted as a topographic high from the Late Proterozoic to the Early Paleozoic and subsequently collapsed from the Carboniferous to the Triassic (O'Connell, 1994; Eaton et al., 1999b). With its varied history and diverse depositional environments, it has been the focus of oil exploration for many decades. The uplift of the PRA along the western passive margin of Laurentia has been attributed to elastic deflection due to intraplate stresses. These stresses are thought to have resulted from the extension onto land of oceanic transform faults running perpendicular to the rifted margin (O'Connell, 1994). Eaton et al. (1999b) suggest that the presence of the Winagami sequence and its contribution to strengthening the crust also influenced the stress buildup. Along with the Western Alberta Arch to the southwest, the PRA cut off the Pine Creek region  26  Chapter 2. The 3-D Surveys -120"  -118"  -116"  -114"  Figure 2.9: The spatial relationship between the Pine Creek region and the Peace River Arch (rough outline plotted as dashed line), a potential influence on sedimentary deposition. The Cordilleran deformation front (DF) and the shaded outline of the Winagami reflection sequence are shown for reference.  -120'  -118'  -116'  -114'  from the open sea to the west. It was only in the Devonian that sedimentation resumed when an epeiric sea formed in the region. Throughout the Devonian, the PRA remained exposed as evidenced by the lack of Devonian sediments on the structure itself and a gradual thickening of sediments away from the arch. The collapse of the arch and the coeval demise of passive margin sedimentation by the end of the Carboniferous halted the second pulse of sedimentation. Localized subsidence related to the enigmatic Dawson Creek Graben Complex is thought to have produced the thin layer of Jurassic sediments in the Pine Creek area (O'Connell, 1994). As in the south, the final pulse of sedimentation was initiated by the formation of the Cordilleran Foreland Belt to the west and the corresponding loading and subsidence of the cratonic margin. Presently, the surface geology in the Pine Creek region uniformly consists of Cenozoic sedimentary strata. In that both the Pincher Creek and Pine Creek seismic survey regions overlie the Western Canada Sedimentary Basin, the stratigraphic units below the two survey areas are similar. They do however differ greatly in that the Pincher Creek (PCI) region lies west of the Cordilleran deformation front and the Pine Creek (PC2) region lies to the east. As such, the sedimentary sequence in the Pine Creek region was subjected to far less deformation following deposition  Chapter 2. The 3-D Surveys  27  compared to the Pincher Creek region to the south and so is less structurally complex.  2.2.3  The Pine Creek (PC2) 3-D dataset  The PC2 3-D dataset was collected by Veritas Energy Services Inc. in January 2001 for BP Canada. The acquisition parameters for this survey are presented in Table 2.2.  Total number of shots Shot size/source Shot depth Shot interval Source line spacing Receiver group interval Receiver line spacing Maximum number of live receiver lines Maximum number of live channels Geophones per receiver station Geophone spacing Sample rate Record length Maximum source-receiver offset Average source-receiver offset Azimuth of receiver lines Azimuth of source lines  3724 4 kg Dynamite 18 m 140 m 630 m 70 m 560 m 11 2410 6 4m 2 msec 5.108 sec 16,190 m 2904 m 136° 46°  Table 2.2: Acquisition parameters for the PC2 3-D dataset in northwestern Alberta. The acquisition geometry of the PC2 dataset is shown in Figure 2.10. For this much larger survey, 3724 shots were recorded by patches of 5 to 11 receiver lines each consisting of 12 to 254 receiver groups. These receiver lines were spaced 560 m apart and the receiver groups were connected at a spacing of 70 m. Due to this variable acquisition, each shot was recorded by between 240 and 2410 receiver groups. While the maximum source-receiver offset recorded was 16,190 m, the average offset was 2904 m as shown in the source-receiver offset distribution plot (Fig. 2.11a). The source-receiver azimuthal coverage from the survey was equally distributed (Fig. 2.11b).  Chapter 2. The 3-D Surveys  Figure 2.10: Survey geometry for the P C 2 3-D dataset.  a)  3 5  b)  xio5  N  0  Offset (m) Figure 2.11:  ®  source-receiver pairs  Source-receiver distributions for all shots from the P C 2 3-D dataset; a) histogram of  source-receiver offsets, b) angle histogram of the source-receiver azimuths.  Chapter 2. The 3-D Surveys  2.3  29  Chapter summary  In this chapter, the tectonic environments and stratigraphic sequences beneath the two 3-D surveys were presented along with their acquisition details. The next chapter outlines the steps involved in the processing of these datasets.  Chapter 3  Seismic Processing  Following the sinking of the Titanic in 1912, an invention based on the use of acoustic waves for iceberg detection led, in 1917, to the first U.S. patent suggesting the use of seismic waves for exploration (Sheriff and Geldart, 1999). Since that early start and with countless subsequent successes in locating economic targets, particularly in the petroleum industry, seismic methods have dominated applied geophysics. Over the years, as acquisition, data processing and interpretation techniques have evolved and improved, seismic data have yielded more and more important geological information. With the significant increase in the size of acquired datasets and the need to image finer and finer structures, seismic data processing arguably has seen the greatest number of developments. Consequently, standard seismic processing flows now involve the execution of an elaborate sequence of sophisticated techniques designed to create images as close to the reality underfoot as possible. In this chapter, the general steps involved in the processing of both 3-D datasets are outlined (Fig. 3.1). Dataset-specific details are added as needed. Prior to the donation of the PCI 3-D dataset to U B C , the shallow portion of the dataset corresponding to the sedimentary sequence (down to 3.6 s TWTT) was processed by Rob Tilson of Kelman Seismic Processing in Calgary, Alberta. His velocity analyses were used as the starting point to the velocity analyses undertaken for the PCI dataset in this thesis. Despite this earlier industry processing, the work presented here represents the first attempt to process the PCI 3-D dataset down to 8 s TWTT.  30  Chapter 3. Seismic Processing  (  Geometry (3.2) Trace editing  31  Q Geometry (3.2)  )  (3.3)J)  Trace editing (3.3)J Refraction statics (3.4) )  Refraction statics (3.4) )  ^ . ^ ^ g K roM^OJfflai;^^!^  )  ( Geometrical spreading correction  Geometrical spreading correction (3.6) )  '(3^6)^)  ( CDP sort)  ( CDP sort)  Velocity analysis (3.8) )  x2 Velocity analysis (3.8) )  Residual statics (3.9)  )  x2 c  Residual statics  ( Spiking deconvolution  Spiking deconvolution  (3.10)~)  Trace balance and stack (3.11))  (V-JC deconvolution (3.13))  (3.10) )  Q Trace balance and stack  ( Migration  (3.H/T)  (3.12) )  /-A: deconvolution (3.13))  Figure 3.1: Processingflowsapplied to the PCI and PC2 datasets. Bold numbers indicate the section in which the processing step is discussed in the thesis. Note that an extra ground roll removal step was applied to the PCI dataset while migration was not performed. Two passes of velocity analysis and residual statics were applied to each dataset.  Chapter 3. Seismic Processing  32  The PC2 3-D dataset donated by BP Canada was processed by two separate contractors in Calgary prior to this study but none of their results were made available. As with the PCI dataset, the deep data have not been previously investigated.  Software Most of the processing results presented in this thesis were obtained using the Globe Claritas interactive seismic processing software developed by the Institute of Geological and Nuclear Sciences Limited in Wellington, New Zealand. Other processing packages used include the Colorado School of Mines' Center for Wave Phenomena's Seismic Unix software and the FreeUSP package from Amoco Production Company's Tulsa Research Center. The computation of refraction statics was done using Hampson Russell's gli3d software.  3.1  Initial data preparation  Raw S E G Y field records from the PCI survey were provided by EnCana to U B C in the Fall of 2000. These data were initially examined using Landmark ITA's INSIGHT processing package (with limited 3-D capability) until the Globe Claritas system became available in April of 2001. The initial raw records contained a number of shorter traces corresponding to tests (e.g. instrument noise, impulse, distortion) undertaken by the acquisition team. These were removed prior to reading the S E G Y data into Claritas which required constant trace lengths. The next step in the initial data preparation involved the inclusion of geometry information in the trace headers. This information was obtained from the survey and observer notes provided with the data. This step proceeded efficiently as the observer notes provided by the recording team were complete and accurate. The PC2 data arrived from Kelman Archives as raw SEGD records in June of 2003. This dataset was significantly larger than the PCI dataset (~ 7 times larger) and originally arrived  Chapter 3. Seismic Processing  33  with no geometry information provided. In September 2003, some recording information was located and using this, I proceeded to work out the proper geometry over the next two months. By November 2003, this task was complete and processing was begun. The proper specification of geometry information within the trace headers of seismic data is absolutely crucial to most procedures in data processing. As such, great care and attention was devoted to this preparatory step.  3.2  Geometry  Once the seismic traces from each shot can be assigned to their specific receivers and so be located spatially, the midpoints between all receivers and shots can be determined (Fig. 3.2a). For flat lying layers, the midpoints on the surface lie directly above the reflecting points on those layers. By re-sorting the traces according to common-midpoints (CMPs) (Fig. 3.2b), correcting for their travel paths (Section 3.8) and stacking them together, a new trace can be obtained that characterizes the reflections at that midpoint location. Since this summed trace results from the exploitation of data redundancy, it is generally less noisy than pre-stack traces. The number of traces used in this summation corresponds to the CMP fold. Since the Earth's reflective layers are generally not perfectly flat, C M P locations and the common reflecting points at depth, known as common-depth points (CDPs), do not generally occur at the same horizontal x,y location. Because of this, the concept of stacking over a commondepth point is generalized to stacking over a common-depth horizontal 2-D area. In the case of a perfectly straight 2-D seismic profile, this horizontal subsurface 2-D area collapses to a line segment parallel to the seismic profile and the traces which correspond to reflections along that specific line segment are stacked together. The length of the segment nominally corresponds to half the receiver spacing and the segment itself is centered on the CMP location. For 3-D  Chapter 3. Seismic Processing  34  Figure 3.2: Schematic diagram of a) the 3-D geometry for a simple shot with four receivers and b) the 3-D geometry for a common-midpoint (CMP) made up of two pairs of shots and receivers. Dark gray ellipses indicate the surface locations of the source-receiver midpoints in (a) and the CMP location in (b). Light gray ellipses indicate the locations of the source-receiver midpoints in (a) and the CMP location in (b) on the reflecting surface at depth. For a flat reflector (as in this diagram), the CMPs and common-depth points (CDPs) occur at the same horizontal x,y location. surveys, the horizontal subsurface 2-D area corresponds to a C M P bin with bin dimensions depending on the receiver spacing along (inline) and across (crossline) receiver lines. Crooked 2-D seismic profiles generally involve CMP bins with crossline dimensions arbitrarily chosen by the processor. Despite the fact that CMPs and CDPs rarely coincide at the same horizontal location in practical situations, the terms are still used interchangeably in the petroleum industry and so will also be used interchangeably within this thesis.  3.2.1  Acquisition of the P C I 3-D survey  The PCI survey was acquired using an adapted Megabin geometry with co-linear receiver lines and source lines. This type of acquisition was pioneered and patented by EnCana Corporation. Constraining all sources and receivers to parallel lines has the cost-saving effect of reducing the number of lines that need surveying and consequently reducing both the number of land permits required and the environmental impact of the survey. In this type of acquisition, sources and receivers are co-deployed along straight lines with the source spacing generally double that of  Chapter 3. Seismic Processing  35  the receivers. Since this geometry is so regular, the midpoints between shots and receivers tend to all coincide with the center of the CMP bins (Fig. 3.3), assuming that the underlying reflectors are flat.  Figure 3.3: Midpoint cluster in a bin from the center of the PCI dataset.  Figure 3.4: Stacking fold map of the PCI 3 - D dataset. All offsets are included.  Chapter 3. Seismic Processing  36  The PCI 3-D survey geometry provides for an optimal stacking bin size of 60 m by 120 m. After stack, this binning geometry yields a data volume of 109 traces in the WSW-ENE direction by 73 traces in the NNW-SSE direction. The maximum fold for all the bins is 211 and out of the total 7957 CDP bins, 1894 bins did not contain any traces due to the irregular survey shape. As many of the inline receiver lines are separated by 480 meters, the resulting stacking fold map using all offsets has a distinctly stripy character (Fig. 3.4). Since this particular Megabin geometry results in denser receiver coverage in the inline direction than in the crossline direction, the survey was strategically oriented so that the crossline direction would run along the strike of the Cambrian to Mississippian sequence which doesn't vary greatly.  3.2.2  Acquisition of the PC2 3-D survey  The PC2 survey was acquired using the more common 3-D approach of having receiver lines and shot lines perpendicular to one another. These source and receiver lines generally follow more tortuous paths than those in Megabin surveys as they avoid cultural features. Consequently, the source-receiver midpoints are far more scattered within CMP bins (Fig. 3.5). By halving the receiver group spacing and the receiver line spacing to obtain the binning geometry, the starting bin size used for the PC2 dataset was set to 35 m by 280 m. Using this bin size, the scatter of midpoints in Figure 3.5 suggests that halving the crossline bin lengths to 140 m would actually reduce the degree of smearing within bins and consequently improve the stack. This ability to further optimize bin sizes is one advantage of having crooked source and receiver lines. If left with 35 m by 280 m bins and using all source-receiver offsets, the stacked PC2 data volume consists of 528 traces in the NW-SE direction by 79 traces in the SW-NE direction with a maximum fold of 201 (Fig. 3.6A). With the reduced bin size of 35 m by 140 m, the number of NW-SE inline traces remains the same but the number of SW-NE crossline traces is increased  Chapter 3. Seismic Processing  37  Figure 3.5: Scatter of midpoints in a 35 m by 280 m bin from the center of the PC2 dataset. to 157. Consequently, the resulting maximum fold is reduced to 107 (Fig. 3.6B). Due to the diamond-shape of the survey, fold drops off in the corners regardless of the bin sizes. It should be noted that the displayed PC2 fold maps were constructed after 702 shot records had been eliminated from the dataset due to excessive noise. In some cases, this noise was cultural while in others, weak shots resulted in poor energy propagation. Despite the removal of these records, the stack fold throughout the survey remains consistently high with both bin sizes. Although the ideal binning geometry for this standard orthogonal acquisition should correspond to a bin size of 35 m by 70 m, the bins were kept at 35 m by 140 m to increase the stacking fold and so help bring out deep reflections. Smearing effects within the bins were deemed negligible given the horizontal to sub-horizontal geometry of the reflective targets.  3.3  Trace editing and geometry quality control  Editing of individual noisy traces on both datasets was achieved using the interactive sqc program within Claritas. Traces were visually deemed too noisy and were killed if no useful signal could be extracted from any part of the trace. These noisy traces were often caused by receiver malfunctions (as recorded in the observer notes) or excessive cultural noise in proximity to the receiver location. Trace editing for the PC2 dataset was done more liberally than for the PCI  Chapter 3. Seismic Processing  38  Figure 3.6: Stacking fold maps of the PC2 3-D dataset for (A) 35 m by 280 m bins and (B) 35 m by 140 m bins. Note that with a smaller bin size, coverage remains uniform but with lower folds. All offsets were used for these maps. dataset since data coverage was significantly better. Trace editing for both datasets was conducted on swaths of traces corresponding to individual receiver lines for each shot. Since these swaths were defined based on information in the trace headers, any geometry errors were readily apparent. On both datasets, some of these geometry errors resulted from typos in the observer notes. However, a problem for several of the shots from the PC2 dataset involved five groups of up to five receivers which were either connected in the reverse order or out of sequence along the line. In addition, polarity reversals were observed on a limited number of traces from the PC2 dataset. These polarity flips were found to correspond to 14 specific receivers. These problems were all diagnosed and corrected prior to further processing.  Chapter 3. Seismic Processing  3.4  39  Refraction statics  Since seismic surveys are rarely acquired over perfectly flat topography and a weathering layer of varying thickness generally underlies this topography, seismic data can suffer from distortions that misalign the traces within a CMP gather and so deteriorate the quality of the stack. To combat these distortions, refraction statics must be calculated. These statics are trace-bytrace shifts which correct for near-surface irregularities and rearrange the data into what they would have looked like if they have been acquired over a flat plane with no weathering layer present. Computing these shifts is an essential step in the processing of land seismic data. While topographical information for seismic surveys is generally provided by the surveyors staking out the acquisition geometry in the field, details about the underlying weathering layer are not obvious from the surface. Fortunately, this information can be obtained from the first breaks of the acquired seismic data. A first break on a seismic trace is the first arrival of the seismic signal at that receiver location. This arrival corresponds to the fastest traveltime for seismic energy refracting in the shallow subsurface layers and so provides information about the velocity structure of those layers. From the redundant coverage in a seismic reflection survey, a wealth of first break picks is available to help construct a velocity model of the near surface. Using this model, the appropriate static shifts can be computed for each shot and receiver location. First break arrivals from both datasets were picked interactively using a trace-to-trace crosscorrelation routine within Claritas' sqc program. To align the first breaks closely enough to allow the cross-correlation routine to function properly, the first arrivals wereflattenedusing a linear moveout velocity of 3000 m/s, representative of the weathered sediments through which the refracted arrivals travelled. Once the first breaks were aligned, a pick was made on the first trace of a swath of traces corresponding to a complete receiver line and the routine automatically determined where the picks should be made on allthe other traces in the swath by comparing  Chapter 3. Seismic Processing  40  adjacent traces from left to right. While the method was generally quite accurate on the first attempt, all first breaks were individually checked visually to ensure that the proper part of the trace had been picked. While this extra quality control step expended some of the extra time gained from using the automatic routine, it greatly increased the confidence level in the picks. Refraction statics were computed for both the PCI and PC2 datasets using the generalized linear inversion program gli3d from Hampson Russell (Hampson and Russell, 1984). The program iteratively converges to a 2-D or 3-D velocity model of the near-surface by performing three tasks during each iteration. Firstly, ray tracing is run through a starting model and traveltimes for the traced rays are computed. Secondly, model perturbations to reduce the error between the computed and observed traveltimes are calculated. Lastly, the starting model is updated using the calculated model perturbations.  3.4.1  3-D model parameterization  Near-surface models in gli3d consist of an arbitrary number of low velocity layers, N, overlying a half-space. Model parameters are defined using a series of continuous functions in x and y, which describe the surface topography E(x, y), the elevations at the base of each layer Di(x, y) where i = 1,N and the velocities at the top of each layer Vi(x, y) where i = 1, N + 1 which includes the velocities at the top of the half-space (Fig. 3.7). With this parameterization, a model consisting of one low velocity layer is described using four functions while a two layer model requires six functions. Constructed models can be laterally complex since both the layer thicknesses and velocities can vary spatially. Conversely, velocities within a layer must remain constant in the vertical direction. Since the topmost weathering layer in the model is generally very thin and poorly sampled, its velocity is assumed to be known and must be assigned. One of the main differences between 2-D and 3-D surveys is that, for 3-D surveys, traces from a given shot do not generally correspond to a uniform range of offsets. Trace locations can  Chapter 3. Seismic Processing  41  Figure 3.7: Model parameterization in gli3d. For this simple model with two low velocity layers, the continuous functions E(x, y), Di (x, y) and Vi(x, y) describe the surface topography, the elevation at the base of the i layer and the velocity at the top of the i layer. th  th  cluster at certain offsets depending on the acquisition geometry. As such, for undersampled offsets, refractions from the base of the weathering layer may not be recorded and no information about the thickness of that layer at that location will be available. Furthermore, offset gaps can result in areas of the model being undersampled by rays and consequently velocities in those areas will be poorly resolved. To circumvent these problems, gli3d parameterizes the 3-D velocity model layer boundaries using cubic polynomial patches that are smoothly joined together with the number of patches dictating the smoothness of the model layer (i.e., the fewer number of patches, the smoother the layer boundary). These patches are forced to be continuous and have continuous first derivatives across their boundaries. Consequently, the whole model can be self-consistently updated after each inversion step even in areas with little to no ray coverage. Furthermore, since this parameterization allows both layer thicknesses and velocities to be determined at any point in the model using interpolation, ray-tracing can be performed in any direction.  Chapter 3. Seismic Processing  3.4.2  42  R a y - t r a c i n g i n 3-D  Ray-tracing for each shot in gli3d is performed within vertical 2-D planes that contain both the shot and receiver locations (Fig. 3.8). Since these planes can be oriented in any arbitrary horizontal direction within the 3-D starting models, the algorithm is described as a full 3-D raytracing scheme. For any given plane, a computationally efficient ray-tracing algorithm satisfying Snell's law is run which considers the raypath from the shot to the receiver to be made up of three parts where 1. the raypath from the shot to the base of a layer is roughly vertical 2. the raypath in the underlying layer is straight and roughly horizontal 3. the raypath from the base of a layer to the receiver is roughly vertical. The further assumption that the layers beneath the shot and the receiver are locally horizontal speeds up the algorithm by allowing the traveltimes for the first and third parts of the raypath to be calculated only once for each shot and receiver location respectively. The dip of the middle raypath is computed by determining the elevation difference between the base of the overlying layer at the source and the receiver locations. Similarly, the velocity for the horizontal ray path is taken from the average of the velocities beneath the source and the receiver locations within the layer that the ray is travelling. For a shot within the topmost layer, a direct arrival is also calculated.  3.4.3  S t a r t i n g models  Simple 3-D starting models for the gli3d algorithm were constructed by first examining the slopes of the first break arrivals as a function of offset for representative shot gathers. Since the slope of the first break trend for energy refracting in a layer of uniform velocity will be constant, layers of sharply differing velocities will produce trends of different slopes separated by  Chapter 3. Seismic Processing  43  i  Figure 3.8: A simple diagram of gli3cTs approach to 3-D ray-tracing. On the 2-D vertical slice, only the refracted arrival at the base of the second layer is shown. kinks. For layers with vertical velocity gradients that don't result in sharp velocity contrasts, the first break trend will vary smoothly as a function of offset without distinct kinks. Such trends can be modelled using many thin constant velocity layers. For a representative shot gather from the center of the PCI survey, the average first break trend consisted of three relatively flat segments of slightly varying dip connected by faint kinks. As such, a 1-D starting model of two weathering layers overlying a half-space was chosen. The velocity of the shallowest weathering layer was set to 762 m/s (2500 ft/s), a standard industry parameter used in southwestern Alberta (R. Tilson, personal communication, 2001). The velocities for the second layer and the top of the basement were obtained by matching the slopes of the individual trend segments with those computed for the 1-D starting model. Once the slopes matched and the velocities were set, the thicknesses for both layers 1 and 2 (and correspondingly, the traveltime intercepts for layer 2 and the top of the half-space) were adjusted until the kinks were lined up and the best visual fit was achieved. The resulting 1-D starting model parameters are listed in Table 3.1. The starting model was extended to a 3-D model by copying the 1-D model to all other shot point locations. These individual 1-D models were hung from  Chapter 3. Seismic Processing  44  the surface topography creating a 3-D model with layer boundaries that mimicked the surface topography and with layers of constant thickness. For a total of 306,248 observed picks, the maximum root-mean-squared (RMS) error computed for the traveltimes from this 3-D starting model was 35 ms. While this approach to obtaining a starting model is non-unique, the closer any starting model is to the true solution, the faster the inversion will converge.  Layer 1 2 3  Velocity (m/s) 762 2800 3700  Thickness (m) 30 44  -  Intercept (ms) 0 56 78  Table 3.1: Initial 1-D velocity model used to construct the 3-D starting velocity model for use in gU3d for the PCI dataset.  As with the PCI dataset, the average first break trend for a representative shot gather from the center of the PC2 dataset was relatively smooth with approximately three flat segments of slightly varying dip. Consequently, a 1-D starting model of two weathering layers overlying a half-space was again chosen. Since up-hole times were not provided for the PC2 dataset and for lack of an industry standard for northwestern Alberta, the velocity of the shallowest weathering layer was arbitrarily set to 1000 m/s, slightly faster than the value used for southern Alberta. Velocities for both layer 2 and the top of the basement and layer thicknesses for both layers 1 and 2 (and correspondingly traveltime intercepts for layer 2 and the top of the half-space) were obtained in the same way as for the PCI dataset. The resulting 1-D starting model parameters are listed in Table 3.2. Again, this 1-D model was copied to all other shot locations to construct a 3-D starting model with layer boundaries that mimicked the surface topography and with layers of constant thickness. For a total of 3,001,352 observed picks, the maximum RMS error computed for the traveltimes from this 3-D starting model was 62 ms.  Chapter 3. Seismic Processing  Layer 1 2 3  45  Velocity (m/s) 1000 1500 3200  Thickness (m) 46 69  -  Intercept (ms) 0 51 151  Table 3.2: Initial 1-D velocity model used to construct the 3-D starting velocity model for use in gU3d for the PC2 dataset.  3.4.4 Generalized linear inversion After Hampson and Russell (1984), the forward ray-tracing problem is posed as a simple matrix equation  t = Am,  (3.1)  where the vector m consists of all the thickness and velocity model parameters, the vector t comprises the calculated traveltimes for those model parameters and the matrix A describes the ray-tracing operator. Once a starting model is denned, tis calculated. If an error vector At is constructed where At is the traveltime error between the calculated and observed picks and the matrix B is known and represents the sensitivities of specific traveltime perturbations to specific model perturbations, then a Am vector of model perturbations can be computed by inverting the linear system of equations  A i = BAm.  (3.2)  Note that while the model perturbations are calculated relative to the updated model obtained in the previous iteration, which is the starting model for the first iteration, the traveltime perturbations are always computed relative to the observed picked traveltimes. Given that seismic surveys usually result in much larger quantities of first break picks than  Chapter 3. Seismic Processing  46  model parameters, the inverse problem is overdetermined with many more knowns than unknowns. As such, a least-squares fit,  A m = (B B)- B M, T  1  T  (3.3)  is sought (Lines and Treitel, 1984). Non-random picking errors are reduced using an automatic pick editing step after the first iteration. Editing is performed by determining the RMS error for all the observed and calculated picks and then removing picks with individual errors greater than twice this RMS value. For the PCI dataset, this resulted in the removal of 12,867 first break picks (4.2% of all picks) after the first iteration. For PC2, 66,870 first break picks (2.2% of the total picks) were eliminated after the first iteration. Since the inversion of B B T  in Equation 3.3 at every iteration is computationally expens-  ive, Hampson and Russell (1984) substitute the inversion with an iterative Conjugate-Gradient approximation. In the end, the final near-surface velocity model obtained for the PCI dataset yielded a residual RMS error of 18 ms compared with 35 ms for the starting model. For the PC2 dataset, the final residual error was 13 ms compared with 62 ms for the starting model.  3.4.5  Calculation of refraction statics  The calculation of refraction statics using gli3d involves four basic corrections that account for the depth of the buried shot, the low velocity weathering layer(s), the topography and any residual short wavelength distortions. Together, these static corrections transform the dataset from one acquired over a variable topography h (x, y) with shots buried at depth d(x,y) within an a  underlying low velocity weathering layer to the corresponding dataset that would have been acquired over aflathorizontal surface hj, with no weathering layer present and with the sources located on the flat surface. In all these corrections, the assumption is made that the seismic energy travels vertically in the weathering layers, a reasonable assumption for strong velocity contrasts  Chapter 3. Seismic Processing  47  between the weathering layers and the underlying half-space. The parameters involved in the calculation of these refraction statics are shown in Figure 3.9.  hs(x,y)  2 ^ = shot  Figure 3.9:  Diagram of the parameters used in the refraction statics calculations; h, (x, y) and hd correspond to  the elevations of the actual topography and the flat horizontal datum to which the data are being moved respectively; d(x, y) is the depth to the shot at a particular spatial location; Vk{x, y) and tk{x, y) correspond to the velocity and thickness of the k  th  layer.  Since explosive sources are generally buried within the topmost low velocity weathering layer, the first correction S  up  that can be determined from the near-surface velocity model is the  one-way up-hole correction which brings shots to the surface according to  S p{x, y) — U  d(x,y)  (3.4)  Since shot depths for the PCI dataset varied between 4 and 18 m, different up-hole corrections resulted for each shot. Meanwhile, for the PC2 dataset which consistently involved shots at 18 m depth, no up-hole times were provided by the observer and so a constant up-hole correction was applied for all shots. The signing convention employed in this thesis translates positive corrections into downward shifts (later times) and negative corrections into upward shifts (earlier times). This is the same convention as that adopted by the Society of Exploration Geophysicists (SEG). Since moving the shots to the surface increases the traveltimes for the seismic energy, this correction is always positive.  Chapter 3. Seismic  Processing  48  Correcting for the existence of low velocity weathering layers beneath a seismic survey involves determining the one-way difference in traveltimes that arises if the weathered material is replaced by material with a faster velocity V . For both the PCI and PC2 datasets, V was chor  r  sen to be 3200 m/s. For PCI, this value was selected based on the earlier industry processing of the dataset for which good results in the sediments were obtained. For PC2, the replacement velocity was based on the gli3d results. For the two-layered near-surface models used for both surveys, the weathering correction S  at a specific location (x, y) is  w  (3.5) Since the goal of the weathering correction is to replace low velocity layers with material of higher velocity V , this correction is always negative. r  The elevation correction S (x,y)  accounts for the one-way difference in traveltimes be-  e  tween rays travelling to the surface h (x, y) and those travelling to the reference datum hd aca  cording to  (3.6) Since the velocity of the medium at this stage is already set to the replacement velocity V , the T  weathering correction must already have been applied otherwise significant error will be incorporated into this correction. For the scenario in Figure 3.9 where h is consistently deeper than d  h (x, y), the correction will be negative. Intuitively this makes sense since the seismic energy s  will have less material to travel through and as a consequence will arrive earlier. For the PCI dataset, the datum was set to 1200 m while the topography at the receivers varied between 1214 and 1364 m (Fig. 3.10a). As such, the elevation corrections were all negative. Conversely, the PC2 datum was set to 1185 m while the topography at the receivers varied between 958 and 1412 m (Fig. 3.10b). As such, elevation corrections for the PC2 dataset had both positive and  Chapter 3. Seismic Processing  49  negative values.  —i  800  i  i  i  i  i  i  900  1000  1100  1200  1300  1400  Figure 3.10: Surface topography for the PCI (a) and the PC2 (b) surveys. The same scale, perspective and orientation were used for both plots. Receiver lines are superimposed on both plots for reference.  The final correction involved in the computation of refraction statics is the short wavelength correction Ssw- Unlike most methods for calculating refraction statics, gli3d is able to compute both long wavelength statics associated with the computed velocity model as well as short wavelength statics needed to compensate for errors introduced by the application of smoothers to the near-surface velocity model at each iteration during the inversion. These smoothers eliminate high frequency variations in the velocity model and so produce an incomplete solution. Computation of surface consistent short wavelength statics is done using the delay time method. This method calculates the difference in traveltimes for a critically refracted ray and for a ray travelling along the projection of the refracted path along the refractor itself (Fig. 3.11). This delay  Chapter 3. Seismic Processing  50  time can be divided into a source delay time and a receiver delay time. For a horizontal refractor geometry, these delay times will be identical.  Figure 3.11: Diagram illustrating the components of the delay time calculation of short wavelength statics. The overall delay time can be split into shot and receiver components.  PCI  c  Minimum Maximum Median Mean Standard deviation  PC2  Ssw -36 48 0 -0.4 8.0  E  -29 -1 -15 -13.7 7.0  18 18 18 18  -85 -29 -51 -52.9  S -81 71 -7 -7.2  Ssw -43 140 0 0.5  0  11.9  28.6  11.0  up  u  Minimum Maximum Median Mean Standard deviation  -51 -5 -31 -30.2 8.1  S  5 24 20 19.7 2.4  S  W  e  ST  -80 25 -41 -39.3 14.1 ST  -181 183 -50 -54.3 41.6  Table 3.3: Statistical distributions of the individual refraction static corrections applied to the PC 1 and PC2 datasets. A l l corrections are in milliseconds.  The total refraction static correction ST applied to each shot and receiver location is the summation of the up-hole, weathering, elevation and short wavelength corrections,  ST(X, y) = S (x, up  y) + S (x, y) + S (x, y) + S w{x, w  e  S  y).  (3.7)  Chapter 3. Seismic Processing  51  The statistical distributions of the individual and total corrections for the two datasets are presented in Table 3.3. The longer offsets and the greater variation in surface topography for the PC2 dataset (Fig. 3.10b) result in a wider range of elevation corrections and consequently a wider range of total corrections than for the PCI dataset. Thus, relative to the PCI data, the PC2 data benefit more from the application of refraction statics. The gathers from the PCI dataset (Fig. 3.12) show only minor improvements in the alignment of the reflection hyperbolas while the gathers from the PC2 dataset (Fig. 3.13) show significant improvement with many fewer distortions. Any remaining distortions are corrected at a later stage through the application of residual statics (Section 3.9).  Figure 3.12: Effect of refraction statics on inline gathers from the PCI dataset with the before and after results shown in (A) and (B) respectively. The data have had a 500 ms AGC applied to help bring out the reflections.  Chapter 3. Seismic Processing  52  Figure 3.13: Effect of refraction statics on inline gathers from the PC2 dataset with the before and after results shown in (A) and (B) respectively. The data have had a 500 ms AGC applied to help bring out the reflections.  Chapter 3. Seismic Processing  3.5  53  Ground roll removal  1  The generation of dispersive surface Rayleigh waves on land seismic surveys is a common phenomenon caused by the coupling of P-waves and vertically polarized shear waves (SV). This ground roll propagates along the free surface at low group velocities and when not accounted for in the acquisition design, often dominates seismic gathers as high amplitude, low frequency coherent noise. In many instances, this noise can completely obscure important reflections. For the PCI dataset, ground roll contamination was significant and most deep reflections in the data were obscured, threatening to render the dataset useless for the study of deep sills (Fig. 3.14a). In order to make the dataset viable, ground roll suppression was required. While ground roll was also present on the PC2 data, the larger offsets from that survey meant that far fewer traces were contaminated and so only a conservative high-pass filter (12-14 Hz) was applied to deal with the ground roll problem.  3.5.1  Conventional ground roll removal techniques  As land seismic datasets are often contaminated with high amplitude, low frequency dispersive ground roll, many methods have been developed to suppress this coherent noise. Most of the more conventional methods rely on the ability to discriminate between ground roll and reflections based on temporal frequency and apparent velocity (Sheriff and Geldart, 1999). Of these, the commonly used Fourier-based approaches, which include temporal frequency filtering and multichannel f-k filtering and offset/time variant f-k filtering, assume stationarity of the seismic wavelet. Due to preferential attenuation of high frequencies with depth, shallow seismic A version of this section has been published. Welford, J.K. and Zhang, R. (2004) Ground roll suppression from deep crustal seismic reflection data using a wavelet-based approach: a case study from western Canada, Geophysics, 69:877-884. The work reflects a collaboration between myself and Rongfeng Zhang of the Consortium for the Development of Specialized Seismic Techniques at UBC. Rongfeng's contribution was the development of the Physical Wavelet Frame Denoising algorithm and the adaptation of the code for use with my dataset. My contribution was to apply this ground roll removal scheme to deep seismic data. 1  54  Chapter 3. Seismic Processing  ioH  ,  1  Figure 3.14: Shot gather (a) contaminated with ground roll and (c) after f-k filtering. The low frequency portion of the f-k spectra for plots (a) and (c) are shown in (b) and (d) respectively. Due to the spatial aliasing of the ground roll, evidenced by the wrapping of the spectral lines between 4 and 10 Hz in the enlarged plot in (e) (shown schematically in plot (f)), f-kfilteringof this gather provides no better result than applying a high-passfilter.Only a time gain has been applied to the plots. On plots (a) and (c), the label B denotes the approximate location of the sediment-basement contact.  Chapter 3. Seismic Processing  55  investigations generally maintain distinct frequency bandwidths for reflections and ground roll while deep seismic data often have overlapping low frequency bandwidths for both reflections and ground roll. As such, for deep reflection data, frequencyfilteringcannot remove ground roll without also greatly affecting reflections. While overlapping bandwidths do not preclude the use of f-k filtering which can separate events based on dip, for ground roll with very low group velocities and hence, very steep dips, spatial aliasing causes the f-k spectrum of the ground roll to wrap around, complicating all muting attempts. In southwestern Alberta, the ground roll travels through a surface weathering layer with a depth extent ranging from 0 m to 50 m (Fenton et al., 1994) and velocities that vary between 610 m/s and 762 m/s. As such, spatial aliasing can be severe. This is evident in Figure 3.14b, which shows the f-k spectrum up to 20 Hz for the gather in Figure 3.14a. Due to the conflicting dips and the wrapping of the spectral lines of the aliased ground roll (Fig. 3.14e), muting in the f-k domain provides no better result than high-pass filtering (Fig. 3.14d). As a result, for this dataset, conventional Fourier-based ground roll suppression techniques were ineffective at protecting low frequency reflections. In addition to the standard Fourier-based methods, various alternate approaches exist for suppressing ground roll. The use of noise canceling algorithms (Linville and Meek, 1995), adaptive lattice filters (Saatcilar and Canitez, 1994), adaptive phase-matched filters (Herrmann and Russell, 1990; Saatcilar and Canitez, 1988), the hybrid Radon transform (Trad et al., 2001), the Karhunen-Loeve transform (Liu, 1999) and depth filtering (McMechan and Sun, 1991) are but a few. Recently, the wavelet transform (WT) has been investigated as a promising new means of removing ground roll from seismic data (Deighan and Watts, 1997; Yu et al., 2002; Corso et al., 2003). The 1-D wavelet transform examines the local character of the seismic trace by decomposing it into a 2-D time-frequency plane where dispersive noise can be more easily isolated and filtered. Thus, this decomposition does not assume stationarity of the wavelet and provides a more appropriate tool for dealing with deep seismic data. Deighan and Watts (1997)  Chapter 3. Seismic Processing  56  and Yu et al. (2002) have both demonstrated that ground roll can be attenuated successfully by performing offset-dependent trace by trace muting in the wavelet domain. Taking these wavelet successes to a higher dimension, Physical Wavelet Frame Denoising (Zhang and Ulrych, 2003) provides a new way of removing ground roll from reflection data using a 2-D wavelet that exploits the shape differences between linear coherent noise like ground roll and hyperbolic reflections from deep structures.  3.5.2  Physical Wavelet F r a m e Denoising  Physical Wavelet Frame Denoising (PWFD), introduced by Zhang and Ulrych (2003), is a technique based on mathematical frame theory. Briefly, a frame can be defined as a family of basis functions that characterize the inner products between a finite signal and the basis functions and for which the energy of these inner products is bounded from above and below. Using the values of these projections and the basis functions themselves, a signal decomposed using the basis functions can be reconstructed (Mallat, 1999). Filtering of the signal can be achieved by thresholding the values of the inner products prior to the reconstruction, with the quality of the filtering depending on how well the basis functions resemble the character of the signal. By ensuring that the basis functions resemble the signal, the data in the transform domain will be more compact and the noise and signal can be more easily separated. For the purpose of seismic applications, the signal is considered to be a 2-D common shot gather and the basis functions to be scaled and translated versions of some 2-D wavelet. Conventional 2-D wavelet techniques often use wavelets that are separable in both time and space (i.e., just as for ib(x, t) = ib(x) x ib(t), ib is separable). While this separability simplifies the reconstruction, the related 2-D wavelets are not ideal for use with pre-stack seismic data as they do not resemble the signals of interest. PWFD was designed specifically to deal with reflection shot, receiver and CDP gathers and utilizes a 2-D wavelet that combines both the time  Chapter 3. Seismic Processing  -2000  0  57  2000  -2000  Offset (m)  0  2000  -2000  Offset (m)  0  2000  Offset (m)  Figure 3.15: A selection from the family of physical 2-D wavelets used in Physical Wavelet Frame Denoising (PWFD) as a function of offset and two-way-traveltime (TWTT). By varying wavelet coefficients, the 2-D wavelets can shift in offset and apex time while also changing the curve of the hyperbola and the sharpness of the temporal wavelet. From plot to plot in thisfigure,s increases from left to right and s increases from top to bottom. On all plots, both b and b are keptfixedat 3 s and 0 m respectively. t  t  x  x  element (the seismic wavelet) and its spatial distribution (the hyperbola). The basic 2-D seismic wavelet used in PWFD is called a physical wavelet (Fig. 3.15) and is given by the equation:  2(2o) \a*»t 3TT  t  1  2-1  t ~ b  t  V  x  —  b-,  ) J  Chapter 3. Seismic Processing  xexp  58  llt-bt  2  (  + H  V  2  — a(x  —  b)  ,  2  x  (3.8) where H is a constant usually set to the maximum time of the data, a is a positive real number that quantifies the energy attenuation (as a increases, attenuation increases), s and s are the x  t  scaling parameters in space and time respectively, and b and b are the translation parameters t  x  in space and time respectively. As the scaling parameters for the physical wavelet are adjusted, each new basis wavelet changes in both temporal wavelet shape and spatial curvature (Fig. 3.15). As such, at each shot gather, a range of hyperbolas that could arise given the spatial extent of the gather and its finite frequencies, is computed. Generally, 3 to 5 values of s and 10 to 18 values of s are sufficient t  x  to capture most scaling variations. In addition to the scaling, the physical wavelets are translated in both space and time across the shot gather with the translation parameter increments corresponding to the time and space sampling intervals of the data. This choice of parameter increments has been found to ensure the stability of the reconstruction (Zhang and Ulrych, 2003). Once a frame of 2-D physical wavelet bases has been constructed, the shot gather is projected onto the 2-D wavelets and coefficients Ck are computed for each projection. The seismic signal is then reconstructed from the coefficients and the frame basis functions using the algorithm presented by Zhang and Ulrych (2003). Since the physical wavelet basis functions are not orthogonal, the reconstruction is more efficiently performed in the Fourier domain.  3.5.3  Application of PWFD  PWFD was applied to 2-D common shot gathers using a two-step hybrid approach (Fig. 3.16). Firstly, a 1-D wavelet transform (WT) filter was applied as a preprocessing step to the data  Chapter 3. Seismic Processing  59  to partly suppress noise. This filter was adapted from the one designed by Trad and Travassos (2000) for use infilteringnoise spikes from magnetotelluric data. Use of this preprocessing filter renders obsolete the need for surgical muting of noisy traces and ensures that PWFD will not be affected by noise spikes. The 1-Dfilterwas applied to all traces in each gather at once by concatenating all traces end to end into a super trace. Since both random noise and reflections along the super trace approximate white noise, their effect influences coefficients at all scales and translations in the wavelet domain. Conversely, high amplitude coherent noise or energy bursts are localized in time and so will only have strong influences on specific coefficients. By statistically computing scale-dependent threshold values above which coefficients are deemed to be due to such noise, the coefficients can be scaled down, resulting in the suppression of high amplitude coherent noise and energy bursts. Further details of this procedure can be found in Zhang and Trad (2002). To avoid artifacts from the concatenation of the traces, a taper was applied to the ends of the traces prior to construction of the super trace (Fig. 3.16). (Sort 3-D gathers into inline 2-D gathers)  <^  (Apply t-scaling) (Apply linear taper to each trace) (~ (Concatenate all traces in a 2-D gather^.ft into one supertrace (AppK 1-D WT filtering)  °"|  OH.,  Re-sort concatenated supertrace into a 2-D gather T.-r,^ "(Apply 2-D P W F D X £ y f , t ^ , „  ^  (Remove t-scaling) (Sort inline 2-D gathers into 3-D gathers) Figure 3.16: Outline of the PWFD processing sequence. The highlighted portion demonstrates the steps involved in the Hybrid Physical Wavelet Frame Denoising (PWFD) approach.  Chapter 3. Seismic Processing  60  Once noise has been suppressed using the 1-D WT thresholding algorithm, the traces are sorted back into shot gathers and the PWFD method is applied. As the 2-D shot gather is projected onto different scaled and translated versions of the 2-D physical wavelet, portions of the gather that are similar to the 2-D physical wavelet (hyperbolas) will correspond to high coefficients while linear events like ground roll will have low coefficients. By zeroing out the low coefficients, a 2-D shot gather can be reconstructed with ground roll further suppressed. As an added benefit, thresholding of low coefficients also successfully removes first breaks (Fig. 3.17d), dispensing of the need for independent first break muting. The cut-off at which the thresholding occurs is determined empirically for the whole dataset by testing the results from different cutoffs on a number of representative gathers. Since PWFD is a 2-D-wavelet-based method and the PCI dataset is 3-D, each 2-D receiver line for each shot in the survey was processed through the 1-D W T and 2-D PWFD filtering schemes as an independent common shot gather. While this caused the algorithm to only enhance hyperbolas in the inline direction, the inline direction provides the highest shot fold and thus the most statistically reliable results. Due to the sparse 3-D acquisition geometry of this survey, low folds precluded the use of PWFD on CDP gathers. Prior to application of the filters, a ^-scaling was applied to the data to partially correct for geometrical spreading while not over-gaining the ground roll. Following the application of PWFD, this scaling was removed (Fig. 3.16). Refraction statics had already been applied to the data to ensure that the departures between reflections and smooth hyperbolas were minimal. Each receiver line common shot gather from the survey contained 19 to 55 traces of 2000 samples each. As such, computer memory limitations required that each gather be split into a few super traces for the 1-D W T filtering step. While this may have minimized the contrast between the white trace and noise bursts, the effect is not deemed to be significant. The results of putting the data from Figure 3.17a through the 1-D W T preprocessing step are shown in Figure 3.17c.  Chapter 3. Seismic  Processing  61  Figure 3.17: Shot gather (a) contaminated with ground roll, (b) after application of a 13-15 Hz high-pass filter, (c) after 1-D wavelet transform (WT) preprocessing, and (d) after 2-D Physical Wavelet Frame Denoising (PWFD) filtering. Plots (e) and (f) show enlarged versions of portions of plots (b) and (d) respectively (outlined by dashed gray boxes). Note on plots (e) and (f) how reflection moveouts have been preserved. On all plots, the label B denotes the approximate location of the sediment-basement contact. In (d), note the zone of enhanced deep reflectivity (outlinedby label R) compared to the results in (b) which were obtained using a more conventional approach. Plots (b), (c), (d), (e) and (f) have been scaled by constants to make the maximum amplitudes in the solid gray boxes (between 1300 and 1700 m and 2.08 and 2.24 s) on all plots the same. A time gain has also been applied to the plots.  Chapter 3. Seismic  Processing  62  While the 1-D WT filter did not succeed in removing the ground roll, its amplitudes are reduced relative to the surrounding reflections. These data were subsequently run through the 2-D PWFD algorithm. For this step, the spatial scaling parameter s ranged from 1.5 to 5.5 (normalx  ized by trace spacing) and the temporal scaling parameter s ranged from 6.125 to 25 (normalt  ized by sampling rate). The translation parameters b and b ranged over all offsets and times. x  t  The final results of the application of both the 1-D W T thresholding and the 2-D PWFD filtering of the data from Figure 3.17a are shown in Figure 3.17d. These results demonstrate that the 1-D W T and 2-D PWFD filtering successfully suppressed the ground roll energy, thereby enhancing deep reflections. Comparison of CDP gathers constructed from shot gathers that have been conservatively high-pass filtered (Fig. 3.18a) versus those that have been put through PWFD (Fig. 3.18d) shows that PWFD enhances deep reflections more effectively. As a consequence, the plots of semblance velocity spectra that result after PWFD (Fig. 3.18e) allow for more detailed velocity analyses than those from high-pass filtering (Fig. 3.18b). Since PWFD indiscriminately recovers hyperbolic signals, it should be noted that hyperbolic back scattered ground roll, diffractions from singular scattering points, and multiples will be enhanced by the method. The case for multiples is evidenced by high semblance values at low velocities in the semblance plot from Figure 3.18e. These multiples can be further muted using the Radon transform (Foster and Mosher, 1992) or other techniques for multiple removal from land data. Back scattered energy and diffractions are harder to isolate but may benefit from a second application of PWFD to common receiver gathers. Since this approach would have provided minimal improvement to the PCI results, PWFD was applied only to shot gathers for this thesis. From a comparison of brute stacks resulting from gathers that have been high-pass filtered (Fig. 3.18c) versus those processed with PWFD (Fig. 3.18f), it can be seen that PWFD resolves a greater amount of coherent deep reflectivity, often despite low stack fold (Fig. 3.18g). For the  Chapter 3. Seismic  1000  Processing  Offset (m) 2000 3000  63  Velocity (m/s) 3000 5000  M M 1.9 « 4 . 5  >  E 3  7.5?  I-I  x  > cn  3'  J= 4 >. co  10  O  13  5  CD  a a. CD "O  c  •1  6  £  M6  iUMf H 9  7\  22 (a)  (b)  (C)  o 40  Figure 3.18: Comparison between (a) a CDP gather that has been through a 13-15 Hz high-pass filter, (b) its semblance plot, and (c) the resulting brute stack, and (d) a CDP gather that has been processed through Physical Wavelet Frame Denoising (PWFD), (e) its semblance plot, and (f) its resulting brute stack; (g) shows the change in fold across the stacks (c) and (f). Arrows on (c) and (f) indicate the location of the stacked trace resulting from the CDP in (a) and (d). Variable offset plotting has been used for plots (a) and (d). The label B denotes the approximate location of the sediment-basement contact on plots (c) and (f). In (f), the label R outlines the zone of enhanced deep reflectivity.  Chapter 3. Seismic  Processing  64  PCI dataset, this enhanced reflectivity has a similar character and occurs at similar times to the HSI reflector recorded to the north, perhaps suggesting that both datasets are sampling the same target. Since PWFD uses hyperbolic 2-D wavelets to decompose seismic gathers, the question of how the algorithm will handle non-hyperbolic signals is an obvious one. While Zhang and Ulrych (2003) admit that non-hyperbolic signals will adversely affect the results, they show that smooth deviations from hyperbolic signals have almost no effect. Still, application of statics before PWFD is advisable to minimize deviations. Since the deep reflectors of interest in the PCI region should appear horizontal based on the limited spatial extent of the 3-D survey and on their character on nearby seismic lines (Mandler and Clowes, 1998), the hyperbolic signal assumption inherent in PWFD is thought to be appropriate. Also, since PWFD can handle apex-shifted hyperbolas, offline reflections will still be preserved. For the PCI dataset, the retention of low frequency signals at the expense of ground roll has resolved coherent reflectivity at depths where even conservative high-pass filters either leave behind significant ground roll energy or lose both ground roll and reflection energy. As such, the intended study of deep mafic sills in the area can be undertaken, a task that would have been severely compromised without the application of PWFD and the resulting ground roll suppression.  3.6  Geometrical spreading correction  When a point source embedded in a homogeneous and isotropic Earth is detonated, the resulting wavefront is spherical as it moves away from the source. As this wavefront expands, amplitudes over the wavefront surface decrease as a function of 1/r where r is the distance from the source. This decrease in amplitude, termed geometrical spreading, must be corrected prior to processing to enable late reflections to maintain a closer approximation to their true reflection strength at  Chapter 3. Seismic Processing  65  depth. For a heterogeneous and layered Earth, the wavefront expansion diverges from its spherical shape depending on velocity variations in the subsurface and on amplitude loss from reflections and refractions at layer boundaries. As such, the process of correcting for amplitude divergence is more complicated than just applying a scaling on the order of r, the radius. Common forms of geometrical spreading corrections include the V^ t scaling from NewMS  man (1973) and the t scaling suggested by Claerbout (1985). The t scaling approach was used 2  2  in the preliminary processing of the PCI dataset down to 3.6 seconds TWTT. The advantage of this particular approach is that it is velocity independent and that it approximately recovers both amplitude divergence and amplitude attenuation. The main disadvantage to this approach for deep seismic data is that the correction essentially overcorrects at late times. In order to determine an appropriate scaling for the PCI dataset, the trace envelope for each ungained trace was computed. The trace envelope is obtained by smoothly connecting either the troughs or the peaks along the seismic trace to produce a curve. The resulting curve acts as a reliable proxy of the amplitude decay rate (Yilmaz, 1987). Once the amplitude decay rate for a given trace was known, an appropriate gain function to represent that trace was calculated by determining the maximum value along the decay curve, and for each sample along the trace dividing that maximum by the sample amplitude. The gain function was then applied to the trace by multiplying the sample amplitude by the gain value for that trace sample. As such, sample amplitudes equivalent to the maximum were assigned a gain value of 1.0 and remained unchanged. For smaller sample amplitudes, the gain value was greater than 1.0 and so these were appropriately amplified. To ensure that the traces were not being scaled solely relative to the offset-decreasing amplitudes of the dominant refracted arrivals, the gains were computed once the traces had been muted down to 400 ms past their first break picks as this seemed to correspond to where the refracted arrivals resided.  Chapter 3. Seismic Processing  66  The degree of envelope smoothing in the trace envelope gain method will have an equivalent smoothing effect on the gain function. For instance, if the envelope is smoothed too much, large amplitudes from strong reflections will be smeared with respect to time and the resulting maximum of the trace envelope will be less than the actual amplitude of the trace at that sample point. As such, the gain at that sample point will be less than 1.0 and the output trace will have some of its high amplitudes diminished. Similarly, if the envelope is not smoothed enough, all amplitudes along the trace will be scaled relative to what may be the excessively high amplitude of one distinct reflection. To examine the impact of the smoothing window length, envelope gains were computed for various smoothing windows and these were compared to the results for t scaling, t  15  scaling  and t scaling of the original ungained trace. The results for a trace recorded 4536 meters away 2  from the shot are shown in Figure 3.19. To investigate the offset dependence of the results, Figure 3.20 shows the equivalent plots for a trace recorded 433 meters from a shot. Examination of Figures 3.19 and 3.20 reveals that regardless of smoothing window lengths, trace envelope scaling of the far offset trace is equivalent to t scaling as a first approximation. Conversely, trace envelope scaling of the near offset trace is roughly equivalent to t scaling at early times. 2  At late times in Figure 3.20, the high values for the trace envelope scaling are likely due to a lack of deep reflections along the trace.  Chapter 3. Seismic Processing  70 60 50 o o « 40 c 'co 30 O  20  67  50 ms smoothing window  200 ms smoothing window  7  7  gain from smoothed —~. trace envelope \  gain from smoothed trace e n v e l o p e ' s /  / /  10 aST  0 70  1000 ms smoothing window  60 50 o o 40 45 c '(0 30 O 20 i  gain from smoothed trace envelope  10 3  4  5  6  7  8  0  1  2  3  4  5  TWTT (s) , TWTT (s) Figure 3.19: Comparison of geometrical spreading corrections for different smoothing window lengths (50,200, 500 and 1000 ms) for a trace recorded 4536 meters from the shot. Trace envelope gains are plotted against different types of t scaling. a  While trace envelope scaling provides an efficient means of determining a customized gain function for each trace, its use depends on a uniform distribution of reflections along the trace which is rarely the case for deep seismic data. Further, the use of trace envelope scaling can alter amplitude versus offset (AVO) information which can provide insight into the nature of the reflectors themselves. While this alteration can be useful when trying to emphasize reflectors where amplitude decreases as a function of offset, for the inverse situation where amplitude increases with offset, the technique may cancel out and perhaps reverse the AVO trend. While AVO analysis cannot be used to determine the nature of the deep reflectors investigated in this thesis due to the limited offsets, the trace envelope gain technique was nonetheless judged to be too harmful to AVO effects and an offset-independent compromise of t  1 5  scaling was  Chapter 3. Seismic Processing  68  TWTT (s) TWTT (s) Figure 3.20: Comparison of geometrical spreading corrections for different smoothing window lengths (50,200, 500 and 1000 ms) for a trace recorded 433 meters from the shot. Trace envelope gains are plotted against different types of t scaling. a  deemed to be the best method of correcting for geometrical spreading while conserving A V O trends. For consistency, the same geometrical spreading correction was applied to the PC2 dataset (Fig. 3.21). While a more complete approach to correcting for geometrical spreading is to compute corrections that are both time and offset dependent using the regional velocity field, the t - correction was deemed sufficient for this thesis. 1  5  Chapter 3. Seismic Processing  69  Chapter 3. Seismic  3.7  Processing  70  Compensating for frequency attenuation  In addition to amplitude decay due to geometrical spreading, seismic waves preferentially lose their higher frequencies as they travel deeper into the Earth. This attenuation of frequencies is attributed to conversion of seismic energy to heat and may be related to pore fluid saturation (Yilmaz, 1987). Since seismic resolution deteriorates as the dominant frequency of the seismic wave decreases (i.e. high frequencies are lost), both vertical and lateral resolution will diminish with depth unless the attenuated frequencies are compensated. Similarly, the compensation for frequency attenuation and dispersion helps ensure the stationarity of the seismic wavelet (i.e. it does not change as it travels through the Earth). Frequency attenuation is quantified in terms of a quality factor, Q. For very small values of Q, attenuation is large while as Q approaches infinity, attenuation is effectively zero. The process by which Q is compensated is often referred to as inverse Q filtering. While many favor the simplicity of the constant Q approach (Kjartannson, 1979), a great deal of effort has recently been focused on obtaining Q models based on real observations. Most researchers working to obtain Q models from seismic data have focused on the use of Vertical Seismic Profiling (VSP) data (Spencer et al., 1982). VSP data are an obvious candidate for Q estimation due to the location of the receivers within boreholes at the points where Q is being determined. Still, given the relatively small database of publicly available Q estimates from VSP surveys, the need for methods to determine Q directly from surface seismic surveys seems pressing. While a number of studies have attempted to obtain Q estimates from surface seismic data (Dasgupta and Clark, 1998), few have had particularly convincing results. Given the lack of feasible methods to estimate frequency attenuation in either the PCI or PC2 dataset, processing approaches used within this thesis have been undertaken using either time variant methods or wavelets whenever possible. Both of these approaches do not assume  Chapter 3. Seismic Processing  71  stationarity of the seismic wavelet and so faithfully deal with the time varying frequency content.  3.8  Velocity analysis  Prior to stacking traces within a CMP gather to form the final representative stacked trace, the different travel paths of the reflected arrivals need to be compensated using a normal-moveout (NMO) correction. This correction is based on the assumption that, for small offsets, reflections on CMP gathers will be hyperbolic and centered at zero offset (Fig. 3.22). From this assumption, for a single uniform layer of velocity  VINT  that exists at the Earth's surface and is underlain by  a half-space, the reflection traveltime t at any given offset x can be obtained by simply knowing the zero-offset traveltime 1 according to the equation 0  From this relationship, a correction can be obtained to remove the effect of offset and align all reflections as if they had all been acquired at zero-offset. This N M O correction can be expressed as  (3.10) While for the simple example in Figure 3.22  VNMO  is equal to the interval velocity  VINT  of that  layer, for models with multiple horizontal layers, V^MO approximates the root-mean-squared velocity  VRMS  down to the reflecting layer. If reflective layers are dipping, the stacking veloc-  ity will actually be larger than the RMS velocity. As such, when converting from N M O velocities back to interval velocities for migration purposes, the effect of the dipping layers must be removed to avoid erroneous conversions. One direct approach to dealing with these dipping layers is to compute a dip-moveout (DMO) correction which will migrate the affected arrivals to  Chapter 3. Seismic Processing  72  their equivalent zero-offset CMPs. Due to computational constraints, the lack of a 3-D D M O algorithm and the horizontal geometry of the deep targets, D M O was not performed for this thesis.  CMP Gather Figure 3.22:  N M O explanatory diagram.  From Equation 3.9, two important characteristics of hyperbolic reflections become apparent. Firstly, for a given velocity, the moveout of a reflection increases with offset. Secondly, as velocity increases, moveout at large offsets will decrease. Consequently, since velocities in the Earth generally increase with depth, stacking velocities for deep reflections will be difficult to estimate if the reflections are not sampled over large enough offsets. Consequently, the resolution with which N M O velocities can be determined generally decreases with depth. This resolution is further degraded at large offsets where the small-spread hyperbolic assumption breaks down. The resulting difference between the stacking velocity and VJVMO is termed the spread length bias (Yilmaz, 1987). In contrast with static corrections which shift entire traces up or down, N M O corrections actually stretch the individual traces in a CMP gather with the extent of stretching depending on the offset and the N M O velocity. From Equation 3.10, the largest corrections and corresponding  Chapter 3. Seismic Processing  73  stretching occur at shallow times and at large offsets. If the original sampling points along the trace are used following the stretching, the resulting sampling rate will vary along the trace. To avoid incorporating parts of the trace that have experienced extreme stretching into later processing steps, a stretch mute can be applied which zeroes portions of the trace where the sampling rate has been altered by an amount greater than some percentage of the initial sampling rate. For both the PCI and PC2 datasets which were resampled to 4 ms, a 40% stretch mute implies that only sample rates between 2.4 ms and 5.6 ms will be conserved. As the percentage mute is reduced, the corresponding muting will be more severe.  3.8.1  Velocity analysis i n 3-D  In contrast to 2-D CMP gathers, 3-D CMP gathers generally have broad azimuthal coverage. Consequently, a potential advantage of performing velocity analysis on 3-D CMP gathers is that azimuthal variations in the offset dependent moveout can be examined. Such variations can provide direct evidence of seismic anisotropy if the anisotropy has at least one axis of symmetry that is horizontal (e.g. horizontal transverse isotropy or orthorhombic). In sedimentary areas of the crust, azimuthal velocity studies using three-component data have proven useful in the hydrocarbon industry where seismic anisotropy has been used to map slow zones such as fractures, which contribute both to the hydrological system and to hydrocarbon reservoirs (Rabbel and Mooney, 1996). Unfortunately, since azimuthal variations in P-wave velocities rarely exceed a few percent (Rabbel and Mooney, 1996), datasets like the ones under investigation in this thesis which were recorded using only vertical-component instruments, are not good candidates for azimuthal velocity studies. Furthermore, very few of the 3-D CMP gathers in either the PCI or PC2 datasets have sufficient folds, offsets and azimuthal coverage to warrant such an attempt. As such, only the offset dependent moveout information was used for the velocity analyses undertaken in this thesis.  Chapter 3. Seismic  3.8.2  Processing  74  Semblance analysis  Stacking velocities can be obtained using a variety of velocity analysis techniques. For both the PCI and PC2 datasets, I opted for using semblance velocity spectra. Semblance analysis for a given C M P gather involves obtaining NMO-corrected and stacked traces using a range of constant velocities. The stacked traces are then plotted as a function of velocity versus twoway zero-offset time and the highest stacked amplitudes are picked as representing the optimal stacking velocities (Fig. 3.23). If multiples are present in the CMP gather, they will show up as high amplitudes at approximately the same velocity as their shallower (generally slower) source reflections. Consequently, automatic semblance picking routines may mistakenly pick these velocities unless multiple suppression is performed first. Since multiples on land data are generally harder to identify and remove than on marine data, I manually picked all the velocities from the semblance spectra to avoid this pitfall. For the PCI dataset, semblance plots were computed for every 4th CMP gather using velocities that ranged from 2500 to 6500 m/s with an increment of 25 ms. To increase the signal to noise ratio, three CMP gathers on either side (in the inline direction) of the CMP gather in question were incorporated into the semblance calculation. These nearby gathers had different weights applied, decreasing from 0.9 to 0.6 to 0.3 with increasing offset. The autocorrelation window for the semblance calculation was set to 40 ms and a 250 ms A G C was applied prior to the semblance calculation. When no clear velocity trend was apparent on a semblance spectra plot due to low C M P fold, no velocities were picked. Figure 3.24 shows an example of the velocity analysis results for a CMP gather from the PCI dataset.  Chapter 3. Seismic Processing  -6000  Figure 3.23:  -3000  75  Offset (m) 0  3000  6000  Velocity (m/s) 30004000 5000 6000 7000  A synthetic CMP gather and its semblance spectrum. Random noise with a signal to noise ratio  of 20 has been added to the gather to stabilize the spectrum calculation. The interval velocity model used to construct the synthetic gather is overlain on the spectrum as a dashed gray line and the dotted circle outlines the high amplitudes associated with a multiple at 3.4 s.  Due to the much larger size of the PC2 dataset and the simpler subsurface geology, only every 10th C M P gather was used for semblance analysis. As for the PCI dataset, three CMP gathers on either side (in the inline direction) of the C M P gather in question were incorporated into the semblance calculation to increase the signal to noise ratio. These nearby gathers again had different weights applied, decreasing from 0.9 to 0.6 to 0.3 with increasing offset. The velocities used for computing the semblance plots ranged from 2500 to 6500 m/s with an increment of 25 ms. The autocorrelation window for the semblance calculation was again set to 40 ms and an A G C of 250 ms was applied prior to the calculation. As with the PCI dataset, no stacking velocities were picked from CMPs with low folds. Figure 3.25 shows an example of the velocity analysis results for a CMP gather from the PC2 dataset.  Chapter 3. Seismic Processing  76  1000  Absolute offset (m) 2000 3000 4000  5000  Velocity (m/s) 3000 5000  1000  Absolute offset (m) 2000 3000 4000  5000  1000  Absolute offset (m) 2000 3000 4000  5000  Velocity (m/s) 3000 5000  1000  Absolute offset (m) 2000 3000 4000  5000  Figure 3.24: Velocity analysis example for the PCI dataset showing (a) a CMP gather, (b) its semblance velocity spectrum and (c) the CMP gather following NMO correction using the velocity profile shown by the dashed light gray line in (b). Plots (d), (e) and (f) show enlarged portions of plots (a), (b) and (c) respectively. The enlarged regions are outlined by the dashed dark gray boxes in (a), (b) and (c). A 1000 ms AGC gain has been applied to all the seismic data plots.  Chapter 3. Seismic  77  Processing  Absolute offset (m) 1000 2000 3000 4000 5000  Velocity (m/s) 3000 5000  Absolute offset (m) 1000 2000 3000 4000 5000  Figure 3.25: Velocity analysis example for the PC2 dataset showing (a) a CMP gather, (b) its semblance velocity spectrum and (c) the CMP gather following NMO correction using the velocity profile shown by the dashed light gray line in (b). A 1000 ms AGC gain has been applied to all the seismic data plots.  3.9  Residual statics  Following velocity analysis, the NMO-corrected CMP gathers were used to calculate a residual statics correction for each shot and receiver location. These surface consistent residual statics are another set of trace-to-trace shifts which further align reflections within a small time window for each C M P gather. These shifts are designed to optimize the quality of the stack.  Chapter 3. Seismic  Processing  78  Calculation of residual statics for both the PCI and PC2 datasets was done using the stackpower method (Ronen and Claerbout, 1985). This approach requires both the pre-stack N M O corrected data and an unnormalized stack made from that same NMO-corrected data. The first step in the calculation involves subtracting one shot gather from the stacked volume and then cross-correlating that shot gather with the portion of the stacked volume from where it came. The residual static shift for that shot is picked as corresponding to the maximum cross-correlation value. This shift is applied to all the traces in the shot gather and then the gather is added back to the stacked volume. This is repeated for all shot gathers. Next, the input traces are re-sorted to receiver gathers and the subtracting/cross-correlating/shifting procedure is repeated for each receiver gather. In so doing, the stacked volume is constantly being updated as the shifts are calculated. The stack-power method is iteratively performed until the shot and receiver shifts fall below some acceptable level. For the PCI dataset, residual statics were computed over a 1500 to 2400 ms window that corresponded to bright reflections near the base of the sedimentary sequence. Similarly, a 1700 to 2200 ms window was used for the PC2 data. The stack-power method was repeated twenty times and for each iteration, only shifts up to 8 ms were permitted for any given shot or receiver position. Thus, in theory, this iterative procedure could lead to shifts of up to 160 ms. After every eighth iteration, the cumulative shifts werefilteredto remove any linear trends and the reference stack was rebuilt to reflect this filtering. Figure 3.26 shows an example of the significant improvement provided by applying residual statics to a CMP gather from the PCI dataset. Since residual statics may alter the values of the optimal stacking velocities, a second pass at velocity analysis was undertaken once residual statics had been computed and applied. Prior to this second attempt at semblance analysis, thefirstpass N M O corrections were removed. To further enhance the stacks for the PCI and PC2 datasets, a second pass of residual statics was computed using the refined velocity picks.  Chapter 3. Seismic Processing  79  Figure 3.26: A CMP gather from the PCI dataset before (a) and after (b) the application of residual statics. A 500 ms AGC gain has been applied to both plots. Only the first 3 s TWTT from each gather is shown for clarity.  3.10  Spiking deconvolution  The convolutional forward model of a seismic trace x(t) states that for coincident source and receiver positions over a horizontally stratified Earth, a seismic trace can be constructed directly from the seismic wavelet w(t), the Earth's reflectivity series e(t) and noise n{t) according to  x(t) = w(t)®e(t)  + n(t),  (3.11)  where the symbol ® denotes convolution (Yilmaz, 1987). For simplicity, w(t) is assumed to encompass the effects of reverberation, attenuation, source signature, coupling between source and near-surface, coupling between receiver and near-surface and the recording instrument response. Using this convolutional forward model, the goal of deconvolution is to use the seismic trace x(i) to extract the Earth's reflectivity series e(i). In so doing, deconvolution increases temporal resolution by compressing the seismic wavelet. Since the seismic trace is usually the only known entity in Equation 3.11, assumptions must be made in order for an inverse solution to be constructed. The necessary assumptions are:  Chapter 3. Seismic Processing  80  1. the Earth's reflectivity is random (white) 2. the Earth is made up of horizontal layers of constant velocity 3. the noise component is zero 4. the source wavelet does not change shape as it travels (stationarity assumption) 5. the source wavelet is minimum phase. If the above conditions hold true, the only remaining unknown in equation 3.11 is the reflectivity series e(t) which can be obtained from  e(t) = a(t) <g> x(t)  (3.12)  where a(t) is the inverse filter of the wavelet such that a(t) (g> w(t) produces a delta function. This collapsing of wavelets to delta functions along the seismic trace is the essence of spiking deconvolution.  3.10.1 Spiking deconvolution in practice Since seismic reflection experiments generally involve propagating an unknown source wavelet into complex geological layers in the presence of noise, deconvolution would appear to be a hopelessly flawed pursuit from the start. Nonetheless, over its many years of applications in both industry and academia, deconvolution has proven to be a robust tool even when many of the enumerated assumptions are incorrect. While the validity of assumptions #1 and #2 relies on the chosen survey location and the depth of investigation, noise can be minimized to an extent during acquisition or subsequently filtered during processing so as to better satisfy assumption #3. As explained in Section 3.7, assumption #4 also can be satisfied by compensating for lost frequencies or by performing deconvolution in a time-windowed manner. While an approximation of the seismic source wavelet can be directly recorded in the field, this is rarely done. Instead, the wavelet is usually extracted from the data itself. Yilmaz (1987)  Chapter 3. Seismic Processing  81  demonstrates that for a seismogram derived from a white reflectivity series, the amplitude spectrum of the source wavelet approximates the smoothed amplitude spectrum of the seismic trace. Alternately, he shows that the autocorrelation of the source wavelet is equivalent to the truncated autocorrelation of the seismic trace. As such, for a white reflectivity series, the amplitude information of the source wavelet can be extracted directly from the seismic trace using its amplitude spectrum or its autocorrelation. In practice, the wavelet length is most often chosen to correspond to the first transient zone in the autocorrelation of the seismic trace (Yilmaz, 1987). Since a suite of different wavelets can produce the same amplitude spectrum or autocorrelation function, the phase of the source wavelet must also be considered. The two main properties of a wavelet are that 1) it is one-sided (i.e. all values of the wavelet before its origin time are zero) and 2) it has finite energy (Robinson and Treitel, 1980). This is particularly intuitive for the seismic case with an explosive shot as it means that the wavelet begins when the shot goes off and that it has a finite length as the explosion is of finite duration. The phase of the resulting wavelet describes the energy distribution within that wavelet. If the wavelet has most of its energy concentrated near its start (i.e. front loaded) it is termed minimum phase. The minimum phase wavelet is the optimal wavelet to use in the deconvolution process as it has a stable inverse. While deconvolution may still work if the source wavelet deviates from minimum phase, the closer the wavelet is to minimum phase, the better the deconvolution results. As explosive sources generally produce source wavelets that are close to minimum phase, the minimum phase wavelet with the amplitude spectrum determined from the seismic trace is considered to be the source wavelet, thus assumption #5 can be satisfied. By collapsing the seismic wavelet to a spike, spiking deconvolution essentially flattens the amplitude spectrum by boosting high frequencies. The negative consequence of this is that the deconvolution can enhance unwanted noise. For this reason, applying spiking deconvolution prior to stack can be advantageous since the added noise can be cancelled out during the stacking process.  Chapter 3. Seismic Processing  3.10.2  82  Spiking deconvolution for the PCI and PC2 datasets  In simplistic terms, the subsurface beneath both surveys investigated in this thesis essentially consists of a block of sedimentary layers overlying a uniform basement block containing the odd reflector (Fig. 3.27). Consequently, the validity of the deconvolution assumptions will vary from one block to the other. While the white reflectivity assumption is generally appropriate for the sedimentary block, the sparseness of basement reflectors means that this assumption fails for the basement block. This lack of basement reflectivity can further aggravate the problem by limiting the ability to estimate the source wavelet at depth and so make time-windowed spiking deconvolution inappropriate for the PCI and PC2 datasets.  Figure 3.27: Simple diagram of the subsurface beneath both surveys. Frequency attenuation with depth is demonstrated using the broadening wavelet on the left.  Surface-consistent deconvolution as performed by Globe Claritas estimates the source wavelet for the entire trace based on only one time window per trace. Given that the random reflectivity assumption would be violated by using a time window from the basement arrivals, the same windows that were used for residual statics were used to estimate the source wavelet. Thus, for  Chapter 3. Seismic Processing  83  the PCI dataset, the source wavelet was estimated from the amplitude spectrum of a 1500 to 2400 ms window and for the PC2 dataset, the amplitude spectrum of a 1700 to 2200 ms window was used. The spiking deconvolution operator length was set to 160 ms for both datasets and 0.1% prewhitening was applied. Due to computational constraints related to the large data file sizes, spiking deconvolution for both datasets could only be applied in a shot-consistent manner.  3.11  Trace balance and stack  Once refraction statics, a geometrical spreading correction, ground roll suppression, N M O corrections, residual statics and deconvolution had been applied to each dataset, the traces were sorted into CDP gathers, balanced and then stacked. The trace balancing was applied in the spatial sense according to absolute amplitudes. This step ensured that the average amplitude for all the traces in a CDP gather was constant. The algebraic sum or stack of the resulting CDP traces was normalized by the square root of the number of stacked traces to reduce the influence of noise.  3.12  Migration  While a stacked data volume provides a glimpse into the subsurface, the view is somewhat distorted since the reflections along a trace are all assumed to have come from directly beneath the C M P location. For horizontal reflections, this assumption is correct. However, dipping reflections and diffraction hyperbolas in the stack are improperly imaged. The process by which these diffraction hyperbolas are collapsed and the dipping reflections are restored to their proper length and subsurface position is known as migration. Figure 3.28 illustrates the basic principles of migration. For diffractions, the amplitudes along the hyperbola are summed and placed at the apex of the hyperbola thus imaging the point  Chapter 3. Seismic Processing  84  Unmigrated stack  Unmigrated stack  Migrated stack  Migrated stack  Figure 3.28: Diagrams showing the effect of migration on a) a diffraction and b) on a dipping reflection. The top of each column shows a simple model containing the pattern of normal-incidence rays (black arrows) from a) a scattering point (black dot) and b) a dipping reflector (black line). The boxes along the tops of the models represent coincident source and receiver locations. The unmigrated stacked sections corresponding to the models in the top row are plotted in the middle row. The shapes of the stacked reflections are highlighted in dark gray on both the top and middle row plots. The bottom of each column shows the migrated version of the stack in the middle row. The shapes of the actual reflectors are underlain as the black dot and the black line on the bottom row plots. The light gray arrows in the top row plots outline the displacement directions from the stacked reflections to their migrated counterparts.  Chapter 3. Seismic Processing  85  scatterer that generated the reflections (Fig. 3.28a). Meanwhile, dipping planar reflections in the stack are shortened and moved up-dip by migration (Fig. 3.28b). From Yilmaz (1987), the horizontal d and time d displacements of the center of a dipping migrated reflection can be x  t  quantified according to  (v itan 6 ) 2  t  (3.13)  and  (u tan 6) 2  d = t  2  t  t\l  (3.14)  where v is the medium velocity, t is the traveltime and tan 0 is equal to At/ Ax, the dip of the t  reflection on the unmigrated stack. Following migration, the new steeper dip of the migrated reflection tan 9 can be determined according to t  tan#*  t an Ot i c9  (3.15)  2  t  The spatial relationship between these displacements and dips is shown in Figure 3.29.  tane.  Figure 3.29: Quantitative analysis of the migration process adapted from Yilmaz (1987). The horizontal and time displacements of the reflection are represented by d and d respectively. After migration, the dip of the reflection tan 6 has increased to tan 6 . x  t  t  t  Chapter 3. Seismic Processing  86  A number of general conclusions about migration can be drawn from Equation 3.13. Firstly, as medium velocities increase, the horizontal displacements of reflections due to migration will increase. Secondly, within a medium of constant velocity, horizontal displacements will increase with increasing traveltime. Lastly, the effect of migration becomes more pronounced with increasing reflection dip.  3.12.1 Challenges related to the migration of deep reflection data Migration of deep reflection data differs greatly from migration of reflections from shallower sedimentary sequences (Warner, 1987). For one, many of the arrivals needed to migrate the deep reflections such as diffractions from discontinuous edges, can suffer from greater attenuation as they travel the longer distances back to the surface and can also be distorted by lateral near-surface velocity variations. A further challenge related to migrating deep reflection data relates to migration aperture. Since velocities generally increase with depth and the curvature of a diffraction hyperbola depends on the velocity function (Yilmaz, 1987), diffraction hyperbolas at late times in a stacked section will be broader and thus will span more traces than shallower diffractions. Consequently, migration algorithms which involve the summation of amplitudes along hyperbolas will require a larger number of traces in order for the summation to properly collapse the diffractions. The number of traces involved in the summation is referred to as the migration aperture. Yilmaz (1987) defines the migration aperture necessary for successful migration as  2NX+1  where NX is the horizontal displacement caused by migration d divided by the stacked trace x  spacing. Since horizontal displacements increase with increasing medium velocity and traveltime, deep reflections can require more traces in the summation than were actually recorded. Consequently, migration of deep data may be dip-limited purely based on the number of stacked traces available.  Chapter 3. Seismic Processing  87  The steepest dips that can be properly migrated given the number of stacked traces in the PC2 volume can be assessed using Equation 3.13. For 2-D migration, using all 528 traces in the inline direction as the migration aperture, NX corresponds to approximately 263 traces. Since the stacked traces are spaced 35 m apart, the maximum horizontal displacement possible for a dipping reflection in the middle of the inline section corresponds to 9205 m. Thus, for a traveltime of 3.5 s and a medium velocity of 4500 m/s, the maximum dip that can be migrated given the data coverage is 18 ms per trace which corresponds approximately to a reflector dip of 50° from the horizontal. Results for the PCI data volume are far more restrictive. For 2-D migration using all 109 traces in the inline direction as the migration aperture, NX corresponds to just 54 traces. With a stacked trace spacing of 60 m, the maximum horizontal displacement possible for a dipping reflection in the middle of the inline section is only 3240 m. Consequently, using the same parameters as above, the maximum dip that can be migrated given the data coverage is 11 ms per trace which corresponds approximately to a reflector dip of 22° from the horizontal.  3.12.2  M i g r a t i o n i n practice  A large family of migration algorithms presently exists. Three common algorithms which solve for the scalar wave equation are Kirchhoff migration, finite difference migration and f-k migration. Each of these types of migration can be implemented in a number of different ways; in time or in depth, before or after stack, in 2-D or 3-D. Similarly, each type of migration has its own limitations with respect to velocities and/or dips. The most general Kirchhoff migration algorithms, which include lateral velocity variations, can successfully migrate steeply dipping reflections. However, certain implementations require a constant migration velocity and so can be limiting in areas with variable velocities. On the other hand, finite difference migration can handle both lateral and vertical velocity gradients but may have dip limitations depending on the type of approximation used. Lastly, f-k migration works for all dips but only allows vertical  Chapter 3. Seismic Processing  88  velocity gradients. The choice of which migration algorithm to use and how to implement the chosen algorithm depends on a number of factors including the complexity of the subsurface, the knowledge of subsurface velocities, the availability of migration codes and computer limitations. As the importance of individual factors will vary from dataset to dataset, so will the approach used for migration. Given that both datasets examined in this thesis are 3-D, a 3-D migration code would seem the obvious choice. However, while 3-D migration codes (particularly post-stack codes) are available in most commercial processing packages, none presently exists within the Globe Claritas package. Instead, several two-pass 2-D post-stack migration algorithms are available. These algorithms first perform a 2-D migration along the inline direction followed by a second 2-D migration along the crossline direction, an approach first suggested by Gibson et al. (1983). In addition to significantly reducing the required computation times, these algorithms do not explicitly require the same trace spacing in both directions and so trace interpolation (Chapter 4) need not necessarily be performed prior to migration. While the two-pass 2-D migration approach is equivalent to the one-pass 3-D approach for a constant-velocity medium (Jakubowicz and Levin, 1983), differences between the approaches increase as dips become steeper and lateral velocity variations are introduced. Nonetheless, given the horizontal geometry of the targeted sills and their relatively uniform host rocks, the two-pass approach used in this thesis is thought to be appropriate.  3.12.3  M i g r a t i o n of the P C I a n d P C 2 datasets  As mentioned in Chapter 2, the PCI survey was acquired within the foothills of the Canadian Rocky Mountains where strata of varying velocity have been upthrust northeastward along lowangle thrust faults. Given the orientation of the survey, the two-pass 2-D migration approach  Chapter 3. Seismic Processing  89  remains appropriate for the sedimentary sequence as most variations are oriented in the inline direction. Had the focus of this thesis been the sedimentary sequence, the strong lateral velocity variations would have required that depth migration be used. Furthermore, conflicting dips would have required pre-stack migration. However, since the deep reflections recorded by the PCI survey were essentially horizontal (as will be shown in Chapter 5), the benefits of applying migration were deemed negligible and so no migration was done to the PCI stacked results. In contrast with the PCI survey, the PC2 survey was acquired within the undeformedportion of the Western Canada Sedimentary Basin. As such, the sedimentary sequence overlying the basement is less complex and layer velocities are more laterally continuous. Given these conditions, simple post-stack time migration is sufficient to properly reposition reflections. Since the stacked results for the PC2 dataset show a southeastward dipping reflective surface (as will be shown in Chapter 6), inline 2-D post-stack time migration was applied to gain more insight into the geometry of the reflective feature. Since the migration results tie in directly with the description of the feature, these results are presented in Chapter 6.  3.13  f-x deconvolution  In contrast to reflections from sedimentary layers, deep reflectivity generally appears more discontinuous and consequently can be more difficult to interpret. As such, for display and interpretation purposes, f-x deconvolution was applied to both final data volumes in order to extract coherent linear reflections. This step involves transforming each trace to the frequency-space domain and then performing complex Weiner deconvolution along a particular spatial direction for each frequency. The results are then transformed back to the time domain. For both 3-D data volumes, f-x deconvolution was applied in the inline direction followed by the crossline direction.  Chapter 3. Seismic Processing  3.14  90  Seismic resolution  The steps involved in seismic data processing ultimately aim to provide the best in focus view of the structures at depth. Nonetheless, due to acquisition limitations and shallow geological influences, deeper structures will always remain somewhat blurred. As a consequence, when interpreting seismic reflection data profiles or volumes, resolution limitations must be kept in mind so that only resolvable structures are interpreted. Temporal or vertical resolution refers to the minimum vertical separation needed between two reflecting surfaces for them to be distinguishable as two separate reflections. This separation is defined by Sheriff and Geldart (1999) as corresponding to one quarter of the dominant wavelength A. Since A = v/f  where v is velocity and / is frequency, the dominant wave-  length will increase with depth as both velocities increase and high frequencies get attenuated. As such, vertical resolution will become poorer with depth. A seismic signal of approximately 30 Hz travelling through a 6000 m/s medium will have a dominant wavelength of 200 m. Consequently, reflections will have to be 50 m apart to be imaged as two distinct reflections. In reality, this need not necessarily be the case as closer reflections can constructively interfere and so exhibit tuning effects that can be used to obtain structural constraints (Widess, 1973). Nonetheless, the A/4 rule of thumb provides a good estimate of vertical resolution limits. Lateral resolution refers to the minimum horizontal separation needed between two reflecting points for them to be distinguishable as being separate. This resolution is dictated by the radius of the Fresnel zone,  (3.16) where t is the two-way traveltime of the reflection. For the same seismic signal of approxi0  mately 30 Hz travelling through a 6000 m/s medium with a reflection at 5 s, reflecting points will have to be 2450 m apart (twice the Fresnel zone radius) to be imaged as distinct. As with  Chapter 3. Seismic Processing  91  temporal resolution, lateral resolution becomes poorer with depth. One way to improve upon this lateral resolution is to apply migration which approximately collapses the Fresnel zone to the dominant wavelength A (Stolt and Benson, 1986).  3.15  Chapter summary  In this chapter, the steps involved in the processing of both the PCI and PC2 datasets are presented. Resolution constraints are also outlined. These set the stage for the interpretations presented in Chapters 5 and 6.  Chapter 4 Trace interpolation to mitigate spatial aliasing  The common acquisition geometry used for land 3-D seismic reflection surveys involves laying out parallel lines of receivers with perpendicular lines of sources. These define the inline and crossline directions respectively. The PC2 survey is a typical example of this type of recording geometry. Megabin, the 3-D acquisition geometry used for the PCI survey, differs from standard acquisitions in that source and receiver lines are co-deployed and parallel. For both geometries, data acquisition is performed by having a patch of receiver lines record a given shot. This common 3-D recording technique generally results in CMP bins with broad azimuthal coverage. Given the high costs of acquiring 3-D data, it is common practice in the petroleum industry to undertake 3-D surveys with a smaller receiver spacing along inlines than along crosslines. While horizontal reflectors in the crossline direction will be unaffected by this cost-saving approach, reflectors with a component of dip in the crossline direction will be undersampled and can suffer from spatial aliasing, causing migration problems. To circumvent these problems while keeping acquisition costs low, trace interpolation in the crossline direction has become a standard practice in 3-D data processing. This step is necessary for both standard and Megabin acquisition geometries. Based on a review of the literature, many seismic trace interpolation schemes presently exist. These range from simple routines to find and connect coherent events of arbitrary dip to sophisticated methods based on prediction error filtering. All existing approaches have specific  92  Chapter 4. Trace interpolation to mitigate spatial aliasing  93  input requirements and acquisition geometries for which they work best. Consequently, finding the best technique for a given dataset isn't always straightforward. Since the Precambrian targets of this research project are mostly horizontal to sub-horizontal reflectors, spatial aliasing shouldn't be a significant problem and migration results should not be adversely affected. Nonetheless, this thesis provides a special opportunity to test potential interpolation schemes on two different real 3-D datasets. It should be stressed however that the work presented in this chapter was undertaken to address the practical issues related to interpolating the PCI and PC2 3-D datasets using available software and so does not represent an exhaustive investigation of trace interpolation techniques.  4.1  What is spatial aliasing?  Dipping events on a stacked seismic section can be quantified on the basis of their step-out, which is measured by following the peak of the reflection from trace to trace. This step-out is expressed as Ax/At where A x is the gap between two traces and A i is the time gap between the two peaks from one trace to the other. This step-out is normally specified in units of km/s. Yilmaz (2001) demonstrates that this step-out is equal to the ratio of frequency (units of cycles/s) to wavenumber (units of cycles/km) such that:  /  Ax  k = AT  ( 4 1 )  The significance of this relationship is that for a reflector with a constant dip, higher frequencies will correspond to higher wavenumbers. Since a sharp reflection on a seismic section is composed of a broad range of frequencies, it follows that when a 2-D Fourier transformation is applied to a seismic section, dipping reflections in time will appear as dipping lines in f-k space while flat reflections will be confined to the frequency axis (Fig. 4.1). As a convention in  Chapter 4. Trace interpolation to mitigate spatial aliasing  94  this thesis, reflections dipping to the right are plotted as positive wavenumbers while reflections dipping to the left are plotted as negative wavenumbers. Offset (m)  Wavenumber  Figure 4.1: Dipping reflector (a) on a stacked section and its f-k spectrum (b). The seismic section contains 101 traces with a trace spacing of 10 m. If reflections were dipping to the left, they would plot as negative wavenumbers on the f-k spectrum.  The spatial sampling of a stacked section determines the Nyquist wavenumber above which spatial aliasing occurs. Above the Nyquist, spatial aliasing causes the dipping lines on the f-k spectrum to wrap around into the opposite half of the plot as shown in Figure 4.2. Since migration moves energy up-dip, the aliased frequencies appear to have the wrong dips and their energy is dispersed. While interpolation does not create any new data, it is a cost-effective way of unwrapping the f-k spectrum so that the reflections are mapped to the proper location in f-k space (Yilmaz, 2001). The maximum threshold frequency ( /  m a  x ) above which spatial aliasing occurs for a given  trace spacing is presented by Yilmaz (2001) as  Chapter 4. Trace interpolation to mitigate spatial aliasing  200  Offset (m) 400 600  800  1000  -0.4  95  -0.2  Wavenumber 0  0.2  0.4  (b)  Figure 4.2: (a) Seismic section that involves the same reflector as Fig. 4.1 but this time the trace spacing has been increased to 20 m; (b) its f-k spectrum. Note the wrapping around of lines in f-k space above the Nyquist wavenumber due to undersampling of the reflection.  where Ax is the trace spacing, 6 is the dip of the reflector and v is the medium velocity. Using the gaps between stacked traces for the PCI and PC2 datasets, a medium velocity of 4.5 km/s and dip angles of 10 to 60 degrees, tables of threshold frequencies are presented (Tables 4.1 and 4.2). Since the dominant frequency range for the data investigated in this thesis ranges from 20 to 30 Hz, it becomes apparent that for both surveys, reflectors along the crossline direction with dips greater than 20 degrees will be spatially aliased in the f-k domain. The situation is particularly bad for the traces spaced 240 m apart at the NNW and SSE ends of the PCI survey and for the crossline traces of the PC2 survey spaced 280 m apart. Table 4.2 demonstrates that even a reduced crossline spacing of 140 m for the PC2 dataset may still leave moderately dipping reflections spatially aliased. Consequently, for both acquisition geometries presented in this thesis and for reflections with dips greater than 20 degrees, trace interpolation in the crossline direction should be performed prior to migration. Note that deeper reflections will be less affected  Chapter 4. Trace interpolation to mitigate spatial aliasing  96  by spatial aliasing due to the higher medium velocities.  Dip angle (deg) 10 20 30 40 50 60  60 m 108 Hz 55 Hz 38 Hz 29 Hz 25 Hz 22 Hz  120 m 54 Hz 27 Hz 19 Hz 15 Hz 12 Hz 11 Hz  240 m 27 Hz 14 Hz 9 Hz 7 Hz 6 Hz 5 Hz  Table 4.1: Threshold frequency of spatial aliasing for different stacked trace spacings for the PC 1 dataset. These values correspond to a medium velocity of 4.5 km/s.  Dip angle (deg) 10 20 30 40 50 60  35 m 185 Hz 94 Hz 64 Hz 50 Hz 42 Hz 37 Hz  70 93 Hz 47 Hz 32 Hz 25 Hz 21 Hz 19 Hz  140 m 46 Hz 23 Hz 16 Hz 13 Hz 10 Hz 9 Hz  280 m 23 Hz 12 Hz 8 Hz 6 Hz 5 Hz 5 Hz  Table 4.2: Threshold frequency of spatial aliasing for different stacked trace spacings for the PC2 dataset. These values correspond to a medium velocity of 4.5 km/s.  4.2 Conventional methods for trace interpolation With increased acquisition of sparse 3-D surveys for geophysical exploration over the last twentyfive years, trace interpolation has seen its footprint in the geophysical literature broaden significantly. During this time, new and increasingly complex methods for interpolation have been introduced. The earliest unpublished trace interpolation routines were based on simple trace copying and/orflexbinning algorithms. Flex binning, which involves expanding the size of empty CMP  Chapter 4. Trace interpolation to mitigate spatial aliasing  97  bins until live traces are encountered or a prescribed maximum bin size is reached, was a popular tool for the interpolation of marine seismic data for many years. Despite its prevalent use, flex binning and variations on the technique remain poorly represented in the literature. The first published methods for trace interpolation were based on isolating linear coherent dipping events and interpolating the amplitudes along those dips (Larner et al., 1981; Bardan, 1987). More recently, such methods have been adapted to allow for reflections whose dips change with time (Martinson and Hopper, 1992). The use of prediction error filters for interpolation in x-t (Claerbout, 1992), f-x (Spitz, 1991; Porsani, 1999) and f-x-y (Wang, 2002) space has seen much success. However, since f-x based methods rely on the principle that a prediction operator estimated at a given frequency can be used to predict data at a higher frequency but smaller trace spacing, the techniques do not work well for large gaps in data coverage and/or irregularly spatially sampled traces. Additional assumptions that reflecting events are linear and that amplitudes don't vary along the reflections also limit the practicality of such methods unless they are applied to very small data windows. Another approach to interpolation involves transforming the x-t data to a new domain and then using the data's signature in that domain to construct the equivalent signature for perfectly spatially sampled data. This fabricated signature is then transformed back to the x-t domain, essentially producing the interpolated data. Successful applications of this approach have been undertaken using the Fourier transform (Gulunay, 2003), the Radon transform (Trad et al., 2002), and semblance velocity spectra (Sacchi and Ulrych, 1995). Alternatively, linear inversion routines exist to recover missing traces by using wave-equation based migration operators (Ronen, 1987) or by estimating and correcting for the Fourier coefficients from irregularly sampled data (Hindriks and Duijndam, 2000; Wang, 2003).  Chapter 4. Trace interpolation to mitigate spatial aliasing  98  4.3 Practical issues related to interpolating real 3-D datasets In the acquisition of 3-D datasets on land, a variety of circumstances can lead to deviations from the ideal planned geometry. Cultural features can force the receiver and/or source lines to be crooked, local obstructions can cause receivers and/or sources to be skipped or relocated (skidded) and land permitting problems can leave large areas unsampled. The elimination of noisy traces further decimates the dataset. All these influences can result in large gaps in the data and in traces that are irregularly sampled in space. The combination of these effects with the intentional undersampling of data in the crossline direction provides for a potentially complicated interpolation problem. Prior to stack, the multiplicity of available data means that interpolation can be undertaken in a wide variety of ways using a wide variety of gathers, such as common-midpoint, commonreceiver, common-shot, common-offset or common-azimuth. This flexibility is advantageous in that geometrical acquisition relationships can be exploited in the interpolation routine and the influence of the interpolated traces can be better incorporated into the stack. The increased number of traces will also help to improve the signal-to-noise ratio of the stack. The downside to interpolating pre-stack data is that computation times and data storage requirements will be greater since many more traces need to be processed. The very act of sorting traces into CMP bins and stacking can alleviate many of the problems that arise due to irregular data sampling. Since the CMP bins are regularly spaced and hopefully contain at least one trace, the stacked volume will often heal itself of gaps (provided the gaps are not extremely large) and the stacked trace spacing will be regular. While the crossline trace spacing typically remains coarser than the inline spacing, f-x type interpolation schemes can usually adequately fill in for the missing crossline traces in the stacked volume. Most of the interpolation schemes introduced in the previous section are shown in the literature to work successfully on both synthetic and simple real data examples. However, examples  Chapter 4. Trace interpolation to mitigate spatial aliasing  99  of the implementation of these techniques on large 3-D datasets, such as the ones investigated in this thesis, have not been published. Codes for only a few of the simpler methods are available in commercial processing packages and generally these can only be run on 2-D data. The main challenge to running these algorithms on large 3-D datasets is to implement them such that they are computationally efficient and autonomous (i.e. requiring minimal user input).  4.3.1  Interpolation codes available for this study  The one interpolation module available within Claritas is trinterp. This algorithm takes the log of the fast Fourier transform of an input trace and then uses a polynomial fit (linear by default) to interpolate each log frequency across a given number of traces for each time sample. In so doing, the amplitudes and phases of each event are interpolated separately and all dips are interpolated at the same time. As the algorithm only inputs interpolated traces at the locations of dead traces, these must first be seeded into the dataset prior to interpolation. While potentially more time consuming in the setup stage, trace seeding allows for greater flexibility in that the interpolation can be designed to only work on localized acquisition gaps. Seismic Unix (SU) offers one interpolation scheme suinterp based on automatic event picking. This algorithm inserts and interpolates a fixed number of traces between existing traces by estimating the dip of the main event for a given time sample using the low-pass filtered existing traces. The drawbacks to this method are that only one dip can be interpolated per time sample, the filtered input traces must not be spatially aliased and the interpolation can't be so easily targeted to specific gaps. Also, since the algorithm incorrectly assigns header information (as admitted in the documentation), the method cannot be easily applied to large 3-D datasets without requiring significant book keeping. The FreeUSP processing package offers a windowed f-k interpolation scheme, fkterp, based  Chapter 4. Trace interpolation to mitigate spatial aliasing  100  on the work of Gulunay (2003). As with SU's suinterp code, fkterp does not easily adapt to irregular acquisition gaps. Lastly, FreeUSP contains a number of regridding algorithms for use with both irregularly sampled 2-D and 3-D data. These algorithms can employ either linear or cubic spline interpolation but further details are unavailable due to the limited documentation.  4.3.2 Specifying header information for interpolated traces Regardless of the interpolation technique used, one of the main practical challenges to interpolating real 3-D datasets arises from the fact that seismic traces, both real and interpolated, reside at specific locations in space. Since most processing schemes require some knowledge of this spatial information, interpolation schemes need to be able to insert new traces into the dataset with the proper header information as needed or else dummy traces with the proper header information need to be stuffed into the dataset and then written over by the interpolation scheme. The latter arguably is the most appropriate approach for large 3-D datasets prior to stack as it ensures that the locations of the new traces are accurately known and their header information is properly assigned. Seeding 3-D datasets with dummy traces having the correct geometry information in their headers prior to interpolation is essential, particularly for pre-stack traces. Generally, these new dummy trace headers require not only the new receiver's spatial position and that of its source, but also its corresponding CDP bin. Since the goal of interpolation is usually to decrease bin sizes and increase their number, the new binning geometry needs to be designed prior to the interpolation so that the bin numbers and locations can be properly assigned. Since individual receiver and shot locations are usually assigned some kind of index number, these must also be adjusted to allow for new traces to be inserted. For this thesis, pre-stack dummy traces were inserted into each individual shot gather in three steps. First, dummy receivers were inserted along the inlines to fill in for acquisition gaps  Chapter 4. Trace interpolation to mitigate spatial aliasing  101  x£> ® < ©  0  0  0  Q)  0  g  »  Existing trace  Dummy trace Linear interpolation — direction for dummy trace locations O  9  0  Figure 4.3: Steps involved in dummy trace seeding of the 3-D datasets: A) insertion of dummy inline receivers, B) renumbering of inline receiver indices, C) insertion of dummy crossline receivers. The dashed lines show the linear spatial interpolation lines for new dummy receiver locations. In this mock acquisition geometry, the first and second digits of the receiver indices refer to the inline and crossline numbers respectively. Dummy receiver indices are printed in italics.  (Fig. 4.3A). For crooked inlines, new receiver locations were inserted by linearly interpolating between existing receiver locations. Next, the indices of the individual receivers were renumbered to allow for new inlines of traces to be added (Fig. 4.3B). Lastly, dummy receivers were inserted along the crosslines by linearly interpolating between the existing crossline receiver locations (Fig. 4.3C). At each dummy receiver location, a zero trace was added to each shot gather.  Chapter 4. Trace interpolation to mitigate spatial aliasing  4.3.3  102  Rationale for a new practical 3-D interpolation scheme  Given the limited number of interpolation codes available to me and the logistical complications involved in transferring my entire pre-stack 3-D datasets to alternate processing packages, I decided to design my own interpolation scheme within Claritas. My goals were to develop a simple, robust pre-stack algorithm that would require minimal user interaction and that would exploit geometrical acquisition relationships using individual shot gathers. Shot gathers were chosen for the interpolation as they generally have the highest trace folds compared to receiver and CMP gathers, and they contain traces that have all been influenced by the same source signature and so should look relatively similar. The resulting interpolation algorithm DSInt, which is presented in the next section, is reminiscent of the trace copying andflexiblebinning methods that were initially used in seismic processing.  4.4  DSInt 3-D interpolation  DSInt, which stands for Dead Simple Interpolation, is an empirical 3-D pre-stack trace copying algorithm that relies on the assumption that traces recorded within a limited offset and azimuthal range should be relatively similar as they've sampled roughly the same part of the subsurface. By limiting the azimuthal extent over which the trace copying is undertaken, anisotropic information in the data can be preserved. To implement the DSInt algorithm, dummy traces are first inserted into the dataset at locations where data would ideally have been acquired using the procedure outlined in section 4.3.2. These ideal locations for the PCI dataset are plotted along with the actual receiver locations in Figure 4.4.  Chapter 4. Trace interpolation to mitigate spatial aliasing  103  Figure 4.4: Locations of actual receivers (black boxes) and extra receivers required for an ideal 120 m by 120 m acquisition geometry (black crosses) for the PCI 3-D survey. Prior to interpolation, dead zero traces were inserted at the locations of all the crosses for each shot of the dataset.  /  ^Cshot 0  live trace  | trace copying O dead trace  Figure 4.5: Schematic diagram of DSInt trace copying.  Once the dataset has been populated with dummy traces, both the live and dead traces are sorted into bins based on their offsets and azimuths. The interpolation is then performed for each bin that contains at least one live trace by simply searching around a dummy/dead trace  Chapter 4. Trace interpolation to mitigate spatial aliasing  104  for the nearest live trace and copying the values from the live trace into the dummy trace position (Fig. 4.5). To limit smearing effects, the bin configuration can be oriented such that the center of the first azimuthal bin is aligned with the geological strike in the area (i.e., it is aligned with the strike of the survey). This reorientation was performed for the PCI survey binning and can be seen in Figure 4.6 where the original receiver locations are plotted on top of bins with an azimuthal range of 15 degrees and an offset range of 240 m. 0°  180' Figure 4.6: Bin configuration (15 degrees by 240 m) for a shot at the center of the PC 1 survey. Receiver positions are overlain as small black circles. Interpolation is only performed in bins that contain at least one live trace.  The choice of bin sizes for the DSInt algorithm involves a tradeoff between offset range and azimuthal range. Boxes with limited azimuthal range ensure that information about anisotropy won't be lost while boxes with limited offset range minimize the effect of moveout on the traces. Still, if these ranges are too restrictive, fewer dummy traces will fall in boxes with live traces and as such, little interpolation will actually be done. To obtain appropriate ranges for a given survey, a reasonable starting tactic is to use small bins and increase their size until the desired number of new traces have been generated. Since the optimal PCI inline receiver spacing was 120 m, this distance was deemed to be the  105  Chapter 4. Trace interpolation to mitigate spatial aliasing  most appropriate choice of minimum offset range for the PCI dataset. With no other potential indicators of anisotropy in the region apart from the geological strike, 15 degrees was chosen arbitrarily as the minimum azimuthal range required to conserve the finest details of anisotropy. A graphical representation of this bin configuration for a shot at the center of the survey along with the number of live traces in each bin before and after interpolation is shown in Figure 4.7.  2  4  6  8  10  12  14  number of live traces  Figure 4.7: Original live trace concentration within bins and increased live trace concentration within bins resulting from interpolation using a 15 degree by 120 m bin geometry for shot 200 at the center of the PC 1 survey (shown by star). The dots overlain on each plot represent the original acquisition geometry and the idealized full acquisition geometry respectively. Dummy traces that fall outside of the bins remain dummy traces after the interpolation is complete.  The patchy distribution of filled bins in Figure 4.7 represents 375 bins containing at least one live trace and within these, 1311 new traces can be generated. For all the shots in the survey, this configuration results in 138,002 bins containing at least one live trace and within these, 567,588 new traces can be generated. While this increase in the number of live traces seems  Chapter 4. Trace interpolation to mitigate spatial aliasing  106  large, many of the ideal receiver locations remain dummy traces following the interpolation. To improve upon this coverage, the offset and/or the azimuthal ranges can be increased. One drawback to increasing the offset range for DSInt is that moveout becomes increasingly important. To circumvent this problem, an N M O correction using a rough stacking velocity model can be applied prior to the interpolation and then removed immediately afterwards. This approach allows for larger offset ranges to be used and results in interpolated results that look less blocky. Any shifts that result from the use of an inappropriate velocity model can be fixed later using residual statics.  4.4.1  Copy distance as an error proxy  Since even small azimuthal ranges result in large bins at large offsets, some measure is needed to gauge the reliability of the trace copying algorithm. The distance over which the trace was copied provides one such measure. By tracking this distance and storing it within the trace headers, some record is kept of from where the trace values actually originated. While the nearest live trace will always be chosen as the most appropriate trace to copy, a secondary filter can be implemented to prevent trace copying over distances that exceed some prescribed distance within the bin. While copy distance was tracked in the tests of DSInt, no such secondary filter was applied in this thesis.  4.4.2  Interpolating missing receiver lines: the PCI example  Since the main goal of interpolating 3-D datasets is to generate extra traces in the crossline direction (i.e., new receiver lines between the recorded ones), that is the focus of the first test example of DSInt. The PCI dataset was chosen for this test since the survey design lent itself particularly well to such a test. As mentioned in Section 3.2.1, the PCI dataset was acquired using an adapted Megabin geometry. This adaptation involved adding four extra shot/receiver  Chapter 4. Trace interpolation to mitigate spatial aliasing  107  lines to the original acquisition geometry so that the stacking bin size could be reduced from 120 m by 480 m to 120 m by 240 m. These extra lines are highlighted in Figure 4.8. By removing these lines from the dataset and attempting to recreate them using DSInt, a direct comparison between the recorded data and the interpolated results can be made.  Figure 4.8:  Adapted Megabin geometry of the PCI survey with extra inlines indicated by the arrows.  The shot chosen for this example was from the northwestern part of the PCI survey and was recorded by three of the extra receiver lines (1106, 1114 and 1122). To start, the traces along these receiver lines were set to dummy zero traces. Next, both the live and dummy traces were sorted into bins using the 90 degree by 240 m bin geometry shown in Figure 4.9. The azimuthal range was made large in this example so as to ensure that the maximum number of dummy traces would fall into bins. Despite the large bin sizes at large offsets, the maximum distance over which any trace was copied was 338 m. To avoid biasing the interpolation results with any influence from the 2-D PWFD ground roll suppression algorithm which was applied to the inline receiver gathers (Section 3.5), the data  Chapter 4. Trace interpolation to mitigate spatial aliasing  2  4  6  108  8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 number of live traces  Figure 4.9: Live trace concentrations before and after DSInt interpolation for additional receiver inlines. The interpolation was done using a 90 degree by 240 m bin geometry for a shot in the northwestern part of the PCI survey (shown by star). The dots overlain on each plot represent the locations of live traces before and after the interpolation. Note that in addition to the three extra inlines of receivers (indicated by the arrows), gaps along the original receiver lines have been filled. used for this example were not subjected to PWFD. Instead, only refraction statics and a high pass (12-14 Hz) filter were applied. Furthermore, to prevent any blockiness from arising in the final results due to the use of an offset range greater than the inline receiver spacing, an N M O correction was applied before the trace copying and removed thereafter. The original recorded inline gathers and those predicted using DSInt are shown for all three extra inlines in the first two columns of Figure 4.10. The data are shown only to 3 seconds T W T T to allow for a clearer comparison. To provide a further comparison for the DSInt results, a second interpolation was done using the suinterp algorithm from Seismic Unix. Since suinterp automatically generates new traces during the interpolation, the data used for this interpolation were not seeded with dummy traces  Chapter 4. Trace interpolation to mitigate spatial aliasing  109  j Copy distance (m)  Recorded Figure 4.10:  DSInt  suinterp  Comparison of recorded inline data with DSInt and suinterp results for interpolated receiver inlines.  From left to right, the panels show the original recorded receiver inline gather, the same receiver inline gather predicted using DSInt and the same receiver inline gather predicted using suinterp for inlines 1106 (top row), 1114 (middle row) and 1122 (bottom row). The predicted inline gathers were obtained from all the shot traces except those recorded along inlines 1106,1114and 1122. Shown above the DSInt results are graphs of the distances over which each new trace was copied. The DSInt results shown are for a 90 degree by 240 m bin geometry. A 250 ms A G C gain was applied to all the data sections.  Chapter 4. Trace interpolation to mitigate spatial aliasing  110  prior to the interpolation. Instead, the extra inline receiver gathers were physically removed from the dataset and the remaining live traces were sorted into crossline receiver gathers. Next, suinterp was run and the algorithm proceeded to insert one new interpolated trace between each pair of existing receivers using the automatic event picking described previously. Since the resulting traces had incorrect headers, these were reassigned prior to resorting into the inline receiver gathers shown in the third column of Figure 4.10. From the comparison of the DSInt and suinterp results presented in Figure 4.10, DSInt arguably does an equivalent if not better job of generating receiver inlines with reflective features that most resemble those recorded in the field. For example, the reflection at 1.9 s T W T T on inline 1114 which is barely resolved by suinterp, is clearly reproduced by DSInt. While the smoothness of the reflections produced by DSInt varies from inline to inline depending on the line's orientation relative to the shot, the local character of the reflections generally agrees well with the recorded reflections. In contrast, suinterp produces smoother results but fails to resolve many of the recorded reflections. One obvious feature of the results in Figure 4.10 is that DSInt produces more coherent events than are evident in the real data. These phantom reflections can arise if more than two dummy traces are replaced by the same live trace within a bin. Ultimately, since these extra events aren't likely to be generated at the same times for the other shots, they will be cancelled out once the traces are sorted into CMP bins and stacked. The simple implementation of the DSInt algorithm, its ability to take advantage of arbitrary acquisition geometries and its success at recovering the main reflection events for missing inlines make it a reasonable and accessible tool for interpolating large 3-D datasets. When the 2-D limitations, the strict input requirements and the header problems associated with suinterp are further taken into account, the practicality of DSInt becomes even more apparent.  Chapter 4. Trace interpolation to mitigate spatial aliasing  4.4.3  111  Interpolating missing traces along receiver lines: the P C 2 example  To assess the ability of DSInt to simply interpolate missing traces along receiver lines, the PC2 dataset was used. The shot chosen for this example was from the northwestern part of the PC2 survey. To start, the dataset was decimated by replacing every second live trace with a dummy trace. Next, both the live and dummy traces were sorted into bins using the 40 degree by 70 m bin geometry shown in Figure 4.11. The offset range was set to half the inline receiver spacing of the decimated dataset to force the algorithm to search further than the neighbouring inline trace location for a trace to copy. By limiting the offset range, the need to apply an N M O correction before the trace copying was negated.  2  4  6  8  10  number of live traces Figure 4.11: Live trace concentrations before and after DSInt interpolation for missing receivers along inlines. The interpolation was done using a 40 degree by 70 m bin geometry for a shot in the northwestern part of the PC2 survey (shown by star). The dots overlain on each plot represent the locations of live traces before and after the interpolation. The arrows indicate the receiver inline shown in Figure 4.12.  From the comparison of the recorded and the DSInt interpolated inline traces shown along the top row in Figure 4.12, DSInt does a good job of reconstructing the missing inline traces.  Chapter 4. Trace interpolation to mitigate spatial aliasing  112  The shallow reflections occurring between 1.6 and 2.3 s T W T T are restored as is a strong deep reflection at 4.6 s TWTT. While the difference plot appears to show considerable disagreement between the recorded and DSInt interpolated results, the main features are recovered. Ultimately, the high frequency jitter in the DSInt results can be corrected at a later stage using residual statics.  ,  To again provide a further comparison for the DSInt results, a second interpolation was done using the suinterp algorithm from Seismic Unix. Since suinterp automatically generates new traces during the interpolation, the dataset was decimated by physically removing every second live trace. Next, suinterp was run and the algorithm proceeded to insert one new interpolated trace between each pair of existing receivers using the automatic event picking routine. While no sorting of the traces was required prior to the application of suinterp, the headers on the resulting interpolated traces were again incorrect. These would require reassignment prior to further processing. The interpolation results from suinterp are shown in the middle of the bottom row in Figure 4.12. While suinterp also successfully managed to recover both the shallow and deep reflections from the recorded data, the reconstructed gather appears more wormy. This worminess arises due to the steep dips of both the refracted arrival and the ground roll which violate the spatial aliasing condition discussed in Section 4.3.1. Despite these artifacts, the differences between the recorded and the suinterp results (plotted in the bottom right panel of Fig. 4.12) are comparable to the differences for the DSInt results. While DSInt and suinterp equally manage to reconstruct missing traces along this 2-D inline, DSInt exploited the 3-D information available from the survey and so could have dealt with non-uniform gaps in the trace spacing. This added flexibility and the previously mentioned limitations of suinterp make DSInt a more practical tool for interpolating 3-D datasets.  113  Chapter 4. Trace interpolation to mitigate spatial aliasing  Copy distance (m)  10001  Copy distance (m)  Recorded-DSirit] I''  Recorded suinterp]  j"  |'l  1  I'  III ii  Figure 4.12: Comparison of DSInt and suinterp results for interpolating missing inline receiver traces. The recorded and the decimated inline traces are shown in the top left and bottom left panels respectively. The DSInt interpolated results are shown in the top middle panel along with a graph of the distances over which each new trace was copied. The top right panel shows the difference between the recorded and the DSInt interpolated inline traces. The copy distance graph is again shown for reference. The suinterp interpolated results are shown in the bottom middle panel. The bottom right panel shows the difference between the recorded and the suinterp interpolated inline traces. The DSInt results shown are for a 40 degree by 70 m bin geometry. A t scaling has been applied to all six panels and the data are all plotted using the same amplitude scale. 1-5  Chapter 4. Trace interpolation to mitigate spatial aliasing  4.5  114  C o n c l u s i o n  Overall, my DSInt scheme provides a simple interpolation method that is easily applied to 3-D pre-stack shot gathers. The use of interpolation bins serves to restrict the area over which the interpolation is deemed reliable. By applying the algorithm to shot gathers, the new traces can contribute to a greater number of larger CDP gathers. During stacking, these larger CDP gathers help to cancel out more noise and the resulting stack can be migrated without problems due to spatial aliasing.  4.6  C h a p t e rs u m m a r y  In this chapter, spatial aliasing was explained and traditional interpolation schemes were outlined. My own empirical 3-D trace copying algorithm, DSInt, was introduced and examples from both the PCI and PC2 datasets were used to show the method's effectiveness.  Chapter 5  The Head-Smashed-In Reflector  1  The Head-Smashed-In reflector, first imaged in 1995 using Lithoprobe 2-D multichannel seismic reflection profiles, lies within the Archean Medicine Hat Block of southwestern Alberta. With its horizontal to sub-horizontal bands of bright reflections confined within a depth range of 6 to 21 km (2-7 s TWTT), the HSI reflector has been interpreted as a deep mafic sill complex with an areal extent of at least 6000 km (Mandler and Clowes, 1998). The HSI reflector's loc2  ation within the upper crust of southern Alberta, a region of significant oil and gas exploration, was the initial motivation for using deep industry-acquired 3-D seismic reflection data to undertake a first-of-its-kind localized 3-D seismic investigation of sills within the upper crystalline crust. In this chapter, the Lithoprobe evidence for the HSI is briefly reviewed and the results from the PCI dataset are presented within this context. The regional implications of the results for the HSI reflector are further discussed in Chapter 7.  5.1  T h e S A L T ' 9 5 data  Lithoprobe's Southern Alberta Lithospheric Transect (SALT'95) experiment consisted of a total of nine 2-D multichannel seismic reflection lines extending over an area of over 75 000 km  2  (approximately 250 km east-west by 400 km north-south). The seismic lines in proximity to the A version of part of this chapter has been published. Welford, J.K. and Clowes, R.M. (2004) Deep 3-D seismic reflection imaging of Precambrian sills in southwestern Alberta, Canada, Tectonophysics, 388(1-4): 161-172. I undertook the work independently under Ron's supervision. The figures presented in the chapter reflect more recent results than those in the paper. 1  115  Chapter 5. The Head-Smashed-In  Reflector  116  HSI reflector are shown in Figure 5.1. Line 30 images the HSI reflector best (Fig. 5.2) although lines 29 and 32 also show truncated expressions of the reflectors; these truncations have been used to delimit the edges of the HSI to the north and east. Due to lack of seismic coverage, the southern extent of the HSI is unknown. As mentioned in Chapter 1, the character of the HSI reflections varies from monocyclic to multicyclic across the SALT seismic lines suggesting that the reflective body is laterally heterogeneous with a finite thickness. Polarity constraints based on comparisons with the overlying sedimentary sequence indicate that the HSI reflections result from a high-velocity/high-density structure interpreted as mafic sills (Mandler and Clowes, 1998).  -112°  -  Lithoprobe seismic lines HSI outline — Provincial border — Cordilleran deformation front SAREX refraction line with shots  -114°  -113°  -112°  •111°  -  -  -111°  -110°  Figure 5.1: Detailed Archean basement domain map of the region surrounding the PC 1 3-D survey area adapted from Figure 2.1. The locations of the numbered SALT95 2-D reflection lines and the SAREX seismic refraction seismic line with corresponding shots collected by Lithoprobe are also shown for reference. Note that the Archean basement domains have been delineated based on my own interpretation of the aeromagnetic data.  The HSI reflector appears as packages of one to four horizontal and sub-horizontal bands  Chapter 5. The Head-Smashed-In  Reflector  117  of reflected energy. Two bands of reflectivity that intersect along lines 29 and 30 provide true dip estimates of 13° to the NNW for the shallower reflection and 15° to the N W for the deeper one (Mandler and Clowes, 1998). While the reflections are bright along all three lines, the HSI reflector dips westward and tapers out at the western half of line 30 toward the limit of the Medicine Hat Block (Fig. 5.2). Unfortunately, the location of the PCI 3-D survey, carried out for exploration purposes, is closest to the region of the HSI reflector that is the most poorly delineated. Nonetheless, with 8 seconds T W T T of recorded data, the depth extent over which the HSI reflector may exist is well sampled.  SALT line 30  HSI reflector Figure 5.2: Unmigrated stacked section of line 30 from SALT'95 showing a) the complete line down to 8 s TWTT (~ 24 km depth) and b) an enlarged section of the tapering western edge of the HSI reflector. Outlined box in plot b) shows the portion of the line used in Fig. 5.5. No vertical exaggeration has been applied to the plots.  Chapter 5. The Head-Smashed-In  5.2  Reflector  118  Pre-stack evidence for deep reflectors in the PCI dataset  Upon receipt of the PCI dataset in 2000, preliminary examination of the raw shot gathers revealed that deep reflections between 3-7 s TWTT were obscured due to the significant ground roll contamination discussed in Section 3.5. Consequently, brute stacks were generated at that time to assess whether evidence of the HSI reflector had been recorded by the PC 1 survey. These stacks, while very noisy, did show some evidence for coherent sub-horizontal reflectivity between 2.5 and 5.5 s TWTT. Hence, the processing sequence shown in Figure 3.1, which included PWFD, was undertaken (Section 3.5.3). While the application of PWFD helped enhance deep reflections on the shot gathers as demonstrated in Figure 3.17, stacked results, the primary output from the processing sequence, are necessary for making any meaningful interpretation.  5.3 Post-stack evidence for deep reflectors in the PCI dataset Evidence for deep reflectors is observed on both inline and crossline sections throughout the processed PCI data volume (Figs. 5.3 and 5.4). Due to confidentiality issues, the data are only shown below 3 s T W T T (below the base of the sediments). However, one exception is the easternmost crossline section, where the base of the sediments is shown for comparison with the deep reflectors (Fig. 5.4).  Chapter 5. The Head-Smashed-In  Reflector  119  Figure 5.3: Inline sections through the stacked PC 1 data volume plotted between 3 and 6 s TWTT (approximate depth of 9 to 18 km). The locations of the sections with respect to the data cube are shown at the top and are colour coded with the light gray corresponding to the southernmost slice and the darkest gray corresponding to the northernmost slice. The sections have been horizontally exaggerated 2:1. The black arrows on the sections point to examples of coherent horizontal to sub-horizontal reflections. The gray box overlain on the northernmost slice outlines the region shown in Figure 5.9. Apart from the geometrical spreading correction, no gain has been applied to the data sections. On all plots, approximate depths are based on an average velocity of 6 km/s. Since velocities in the sediments are slower than 6 km/s and the base of the sediments occurs between 2 to 2.5 s TWTT, reflections in the upper crystalline crust are not as deep as indicated by the approximate depth scale.  Chapter 5. The Head-Smashed-In Reflector  120  As explained in Section 3.12.3 and further illustrated by the stacked results, the deep reflections recorded by the PCI survey are essentially horizontal to sub-horizontal. Consequently, migration would leave them unaffected and as such, all subsequent analyses and interpretations of the PCI results are based purely on the stacked volume. On most of the inline sections from the PCI data volume, more than one distinct coherent reflection is imaged as (indicated by the multiple arrows in Fig. 5.3). Similarly, the crossline sections (Fig. 5.4) show multiple coherent reflections. While locally continuous, none of the reflections spans the entire data volume. Instead, pockets of deep reflectivity are located in the southwest, southeast and the northeast of the data volume. Given this irregular geometry, the reflections are not likely to be originating from a single continuous reflective interface. While the PCI data volume does not reveal continuously coherent bright reflections like the ones imaged along most of the SALT'95 lines, as illustrated in Fig. 5.5, the geometry and the traveltimes of the reflections recorded by the PCI survey are similar to the local geometry and the traveltimes of the HSI reflector to the north. Still, their continuity and prominence are less pronounced. This is perhaps to be expected since the western edge of the HSI, as imaged on Lithoprobe SALT line 30 directly to the north, diminishes to the west as noted earlier (Fig. 5.2).  Chapter 5. The Head-Smashed-In  NNW  Reflector  121  SSE  Figure 5.4: Crossline sections through the stacked PCI data volume plotted between 3 and 6 s TWTT (approximate depth of 9 to 18 km). The easternmost crossline section is shown from 2 to 6 s TWTT (approximate depth of 6 to 18 km) to contrast the geometry of the deep reflectors with shallower ones in the sedimentary package. The locations of the sections with respect to the data cube are shown on the right and are colour coded with the light gray corresponding to the westernmost section and the darkest gray corresponding to the easternmost section. The sections have been horizontally exaggerated 2 to 1. The black arrows on the sections point to examples of coherent horizontal to sub-horizontal reflections. Apart from the geometrical spreading correction, no gain has been applied to the data sections. On all plots, approximate depths are based on an average velocity of 6 km/s. Since velocities in the sediments are slower, the base of the sediments (which occurs between 2 to 2.5 s TWTT on the bottom plot) and reflections in the upper crystalline crust are not as deep as indicated by the approximate depth scale. Note that the trace polarities appear reversed as the seismic plots have been nipped to allow NNW to be on the left.  Chapter 5. The Head-Smashed-In Reflector  122  W  >  "D —^ r:.;'"!;::';^'lin™hl;r:;'"S O 12 >< ,  ;r  D 15  fi- *  , —  T ^  » f  18 1 km  3  H.E. 2:1  Figure 5.5: One-to-one comparison between (A) deep reflectors at the western end of line 30 from SALT'95 and (B) the northernmost inline slice from Fig. 5.3. As the data along line 30 had a 1000 ms AGC applied, the same gain has been applied to the data in B. The data sections have been horizontally exaggerated 2:1.  5.4  Statistical analysis of the deep reflections in the PCI data volume  In contrast with the sedimentary targets of shallow 3-D seismic investigations, reflectors within the crystalline basement are generally less coherent and rarely span the entire recorded data volume. Consequently, seismic interpretation packages designed for shallow studies are not ideal tools. For instance, horizon tracking algorithms are not designed to find the edges of reflectors but rather try to extend the horizons to the edges of the data volume. In order to investigate laterally finite pockets of reflectivity, alternate methods must be devised. While several of the reflectors seen on the PCI dataset extend over multiple inlines and crosslines, time slices are unable to provide a consistent enough image of their 3-D geometry to allow for an adequate presentation or interpretation. Similarly, the variations in reflector  Chapter 5. The Head-Smashed-In Reflector  123  length and reflection character make picking of horizons a fairly subjective undertaking. To circumvent this problem and reduce the picking bias in the interpretation, a statistical approach to mapping the reflectors in 3-D is followed. Indeed, and as an historical note, some of the earliest evidence for deep crustal reflections in 2-D surveys derived from a statistical analysis of events recorded on 1950s vintage industry data in Germany that extended to adequate travel times (Dohr and Fuchs, 1967). Given that 3-D planar reflective bodies will appear as linear reflectors along both inline and crossline sections, a spatial coherency filter (Kong et al., 1985) is applied to all the inline and crossline sections in order to automatically extract all linear events that span at least 3 km (the approximate Fresnel zone at 4.5 s for a 10 Hz signal in a medium velocity of 5000 m/s). Once these events are isolated, only linear reflectors from the inlines and crosslines that intersect each other (within ± 10 ms) are selected. Subsequent elimination of common picks that are not bordered by at least one other common pick (within ± 20 ms) imposes a minimum areal constraint on the picked common reflectors, which further reduces the population of common reflector picks. Using this approach, a self-consistent 3-D distribution of reflectors is obtained without any picking bias. For the PCI dataset between 3 and 8 s TWTT, a total of 127,608 coherent picks were obtained for all the inline sections and 44,260 picks for all the crossline sections. From these, a total of 10,915 common picks were found to intersect within 10 ms of one another. With the added areal constraint of requiring common picks to be bordered by at least one other pick, the number of common picks was reduced to 7891. From the distribution of common reflector picks as a function of T W T T for the PCI dataset (Fig. 5.6), most of the reflectors appear to lie between 3.2 and 5.5 s TWTT, with the greatest accumulation occurring between 4 and 5 s. These values agree well with the times at which the HSI was imaged on the SALT'95 lines (Figs. 5.2 and 5.5). A plot of the number of common picks per stacked trace shows that traces with several picks are not randomly distributed in space  Chapter 5. The Head-Smashed-In Reflector  124  but are focused in three main areas (Fig. 5.7), perhaps indicating the presence of discontinuous layered reflectors. From examination of the easternmost crossline section where more than one reflector is imaged (Fig. 5.4), this would appear to be a reasonable conclusion.  number of common picks Figure 5.6: Distribution of reflector picks common to both inline and crossline sections as a function of TWTT. The coloured boxes on theright-handside of the plot show the location of the interpreted reflector surfaces displayed in Fig. 5.8.  Figure 5.7: Spatial distribution of the number of reflector picks common to both inline and crossline sections per stacked trace. The PCI acquisition geometry is shown for reference. Three main regions with multiple picks, possibly showing areas with layered reflectors, are outlined by dashed ellipses.  Chapter 5. The Head-Smashed-In  Reflector  125  Further investigation of Figure 5.6 shows that specific time ranges show the highest numbers of common reflector picks, perhaps indicating the presence of distinct reflection surfaces. Using the four highest peaks in the pick distribution as a guideline, the common reflector picks are subdivided into reflector 1 (from 4180 to 4260 ms, indicated by the purple box on Fig. 5.6), reflector 2 (from 4300 to 4380 ms, indicated by the blue box on Fig. 5.6), reflector 3 (from 4720 to 4800 ms, indicated by the green box on Fig. 5.6), and reflector 4 (from 4830 to 4920 ms, indicated by the orange box on Fig. 5.6). The time spans selected for each reflector were based on the time range of the histogram partition and were only expanded if neighbouring partitions also had a large number of common picks. By enlarging the time ranges, reflector dips can be incorporated as can unwanted neighbouring reflectors. Since these time ranges can be customized, they allow for the topography of the reflectors to be captured and so are easier to interpret and more informative than time slices through the data volume. The four interpreted reflector surfaces show very different geometries (Fig. 5.8). Reflector 1 spans most of the survey with clusters of common picks located in the southwest, southcenter and north. Reflector 2 appears as a more organized version of reflector 1 with distinct pockets of coherent reflectivity. In contrast to the two upper reflectors, reflector 3 is mostly concentrated in the southwestern part of the survey. Meanwhile reflector 4 shows the greatest degree of coherence in the north with disconnected reflector picks in the south. While sparsely distributed and discontinuous, these interpreted reflector surfaces are generally consistent with the distribution of coherent reflections observed along the data volume sections (Figs. 5.3 and 5.4). Nonetheless, while the reflector geometries in Figure 5.8 do not preclude the notion that these reflectors form part of the southwestern extension of the HSI reflector sequence, they do point to a complex intrusive pattern if they do in fact represent mafic sills.  Chapter 5. The Head-Smashed-In Reflector  126  4200 4300 4400 •-4500 4600  ^  —I 3  4700 4800 4900  Figure 5 . 8 : Four interpreted common reflector surfaces determined from Fig. 5.6. The four plots show common reflector picks (1) between 4180 and 4260 ms (purple box on Fig. 5.6), (2) between 4300 and 4380 ms (blue box on Fig. 5.6), (3) between 4720 and 4800 ms (green box on Fig. 5.6), and (4) between 4830 and 4910 ms (orange box on Fig. 5.6), respectively. The PCI acquisition geometry is shown on all plots for reference.  Chapter 5. The Head-Smashed-In Reflector  5.5  111  Waveform modelling  For the three instances of bright upper crustal reflectivity attributed to sills in western Canada, waveform modelling was undertaken to gain insight into the reflector configurations needed to generate the recorded seismic signals. Mandler and Clowes (1997) performed 1-D reflectivity modelling using a layer of sill material intruded into granitic country rock to conclude that the inferred sills making up the Wollaston Lake reflector ranged in thickness from 50 to 150 m. Modelling results from Ross and Eaton (1997) using slightly different velocity and density values constrained the top Winagami reflector to a thickness of 70 ± 10 m. In an attempt to reproduce the HSI reflector's seismic response in 2-D, Mandler and Clowes (1998) stitched together the computed synthetic seismograms from a sequence of 1-D reflectivity models. Their modelling results, while non-unique, showed that lateral variations in the HSI's reflectivity over a few kilometers could be reproduced by laterally varying the number of reflective sills (between 2 and 9) and their individual thicknesses (from 10 to 150 m). Given the discontinuous character of the deep reflections recorded by the PCI dataset, attempting to model the reflective structures in 2-D or 3-D would be exceedingly time-consuming and the results would be far more non-unique than the results from Mandler and Clowes (1998) given the poorer quality of the PCI results. Nonetheless, some insight can be drawn from performing simple 1-D modelling. Using an assumed minimum phase wavelet set to 120 ms in length and estimated from the amplitude spectrum of a 1500 to 2400 ms window of the PCI stacked data as the source wavelet, simple 1-D synthetic seismograms were constructed using a range of different three-layer 1-D models of a sill intruded into a host rock (Fig. 5.9). The algorithm for generating the synthetic traces involved computing the reflection coefficients at the appropriate zero-offset traveltimes and convolving them with the source wavelet. Multiples were not included in the generation of the synthetic traces. As with the earlier modelling results of Ross and Eaton (1997), the sill and  Chapter 5. The Head-Smashed-In Reflector  25 50 100 150 Sill thickness (m)  128  •  200  Granite Dolerite  Figure 5.9: 1-D Modelling results for the PCI dataset; a) the minimum phase wavelet, estimated from a 1500 to 2400 ms window of the data, b) computed seismograms using the extracted wavelet and the simple three layer models of different thicknesses of dolerite sills emplaced in granite shown in (c). In d), a portion of the northernmost inline slice from Figure 5.3 is shown for comparison against the computed seismograms. Apart from the geometrical spreading correction, no gain has been applied to the recorded data section. The synthetic and real seismograms are plotted with the same time scale. The models were constructed using velocity and density values of 6000 m/s and 2650 kg/m for the granite, and 6500 m/s and 3000 kg/m for the dolerite, respectively (Juhlin, 1990). These layer properties correspond to normal incidence reflection coefficients of 0.1 and -0.1 at the top and bottom of the sill layer respectively. Transmission losses at the top of the sill were not incorporated into the modelling. 3  3  Chapter 5. The Head-Smashed-In  Reflector  129  the host rock material properties used were those from dolerite sills emplaced into unfractured granite from the Baltic shield (Juhlin, 1990). Thus, the models were constructed using velocity and density values of 6000 m/s and 2650 kg/m for the granite, and 6500 m/s and 3000 kg/m 3  3  for the dolerite, respectively. These values are appropriate given that the basement rocks of the Medicine Hat Block are granitic gneisses. Sill thicknesses of 5 m, 10 m, 25 m, 50 m, 100 m, 150 m and 200 m were used to span the range of thicknesses deduced from the earlier studies. A portion of the northernmost slice from Figure 5.3 was extracted for comparison (Fig. 5.9d). The dominant frequency of the wavelet in Figure 5.9 is approximately 12.5 Hz. From the discussion of temporal resolution in Section 3.14, the dominant wavelength for a medium of 6000 m/s and a frequency of 12.5 Hz corresponds to 480 m. Therefore, according to Sheriff and Geldart (1999), reflective boundaries must be at least 120 m apart to be distinguished as separate. From the results from Figure 5.9, this estimate seems quite good. For a sill with a thickness of 100 m, the waveform looks much like the source wavelet. However, for the 150 and 200 m sills, the dominant negative peak develops an extra notch signaling that two reflective boundaries are present. Nonetheless, as argued by Widess (1973), tuning effects due to the interference of the top and bottom reflections can be observed from the amplitude and the length of the recorded waveform and thus, the thickness of the reflector body can still be estimated for values less than 120 m. From the comparison in Figure 5.9, the computed waveforms for sill thicknesses between 50 and 150 m seem to best correspond with the recorded waveforms. Considering that the data used for the comparison contains some of the brightest basement reflectivity recorded by the PCI survey, these values of sill thickness likely represent the maximum possible thickness in the survey area. Furthermore, the modelling results show that as the sill becomes thinner and thinner, destructive interference between the reflections from the top and bottom of the layer cause the resulting waveform to weaken such that it would be easily lost within the noise. From these modelling results, the deep reflections (or lack thereof) from the PC 1 dataset are consistent  Chapter 5. The Head-Smashed-In  Reflector  130  with a highly heterogeneous sill complex where sill thicknesses vary significantly over short distances.  5.6 Multiples? One concern with investigating deep reflections beneath sedimentary sequences is whether the presence of free surface or intrabed multiples from shallower reflections will influence or be mistaken for deep reflectors. Since multiples on land data are generally less regular than on marine data and thus more difficult to detect and isolate, 1-D forward modelling provides one method of estimating their influence. Using the reflectivity method of Fuchs and Miiller (1971), synthetic seismograms were calculated with and without multiples for 1-D subsurface models with and without a 150 m sill (Fig. 5.10). The minimum phase source wavelet in Figure 5.9a was again used. The sedimentary sequence for the models was obtained by converting the stratigraphic chart from Figure 2.2 into depth layers with velocities, densities and Q values appropriate to the given lithologies. The velocities were obtained from those used at EnCana for the 2-D pre-stack depth migration of the seismic data within the sedimentary sequence. Average densities to correspond with those velocities were obtained from Barton (1986). Seismic attenuation, Q, values were obtained for sandstones and granites from Lay and Wallace (1995) while the value for carbonates was provided by a carbonate geologist at EnCana (J. Porter-Chaudhry, personal communication, 2004). The properties assigned to the sill were the same as used for the modelling in Figure 5.9. As demonstrated in Figure 5.10, the seismograms are dominated by two bright reflections at 1.8 and 2.3 s TWTT. The first arises from a strong velocity and density contrast within the shales and sandstones (corresponding to the transition from Upper Cretaceous to Lower Cretaceous sediments) and the second is due to the presence of higher velocity and higher density  Chapter 5. The Head-Smashed-In  Reflector  131  carbonates. The boundary between the base of the sediments which corresponds to highly compacted shales and sandstones and the granitic basement rocks results in only a weak pulse. In Figure 5.10b, the presence of a 150 m sill results in a relatively strong negative pulse at approximately 4.6 s TWTT. Based on these modelling results, reflections from a 150 m sill should be detectable even in the presence of multiples. While the depth at which the sill exists will influence the traveltime of the reflection and perhaps cause its arrival to coincide with that of a multiple, the traveltimes of these primary and secondary arrivals will vary with offset thus allowing the two to be separated during velocity analysis. Provided that the stacking velocity corresponding to the sill reflection is chosen, the multiple should be at least partially cancelled out during stacking. To assess the influence of sill thickness on its resolvability in the presence of multiples, a series of seismograms like the one including multiples shown in Figure 5.10b, were calculated for a sill with thicknesses of 10 m, 50 m, 100 m and 150 m (Fig. 5.11). These results demonstrate that as the sill thins to 50 m and less, its reflective signal falls below that of the multiples and thus it cannot be readily detected. Such a conclusion is consistent with the simple 1-D forward modelling results shown in Figure 5.9. If the deep reflections recorded by the PCI dataset are sampling a sill complex with sills of laterally varying thickness, the results in Figure 5.11 potentially provide one explanation for their discontinuous character.  Chapter 5. The Head-Smashed-In Reflector  2H  400r76000 Velocity (m/s)  b)  4000 6000  6  s—* £  a 8  a> Q  No multiples  With multiples  No multiples  With multiples  J_P_  1 I  ?  132  —r—i— 2.5 3.0 Density (g/cc) 2 5 3.0 r  —i r 100 200  Q 100 200  m •'•'•5;  10 12 14 16  i !  I Shales and Sandstones  Carbonates  ffiffiffi  Granite  Dolerite  Figure 5.10: Impact of multiples on deep reflection detection for a 1-D subsurface model a) without a sill and b) with a 150 m sill. Shown from left to right for both parts a) and b) are the 1-D model, the velocity, density and Q profiles and the computed seismograms with and without multiples. Apart from a geometrical spreading correction, the seismograms (which are each repeated six times for clarity) have not been gained. The grey arrows relate the jumps in the property profiles with their corresponding reflections.  Chapter 5. The Head-Smashed-In Reflector  a)  0  10 m  50 m  133  100 m  150 m  Figure 5.11: Impact of sill thickness on reflection detection in the presence of multiples. Using the property profiles from Figure 5.10b and including multiples, the plots from left to right in a) show the complete seismograms for sills with thicknesses of 10 m, 50 m, 100 m and 150 m and in b) show enlarged portions of the seismograms (indicated by the dashed gray box in (a)) focusing on the time at which the sill reflection arrives. Apart from a geometrical spreading correction, the seismograms (which are each repeated six times for clarity) have not been gained. Black arrows indicate the location of the sill reflection among the multiples.  While the possibility remains that some of the deep reflections recorded on the PCI dataset are simply multiples from shallower reflections, a few lines of evidence suggest that this is not the case. Firstly, the stacking velocities for the deeper reflections as determined from semblance analysis tend to be higher than for the shallower reflections (Fig. 3.24). The opposite would be true for multiples. As these higher stacking velocities are used to stack the deeper  Chapter 5. The Head-Smashed-In  Reflector  134  reflections, multiples from shallower sedimentary reflectors, if present, are at least partly cancelled out. Secondly, the SALT'95 data, which were collected in the same tectonic environment, do not exhibit significant evidence of multiples (Fig. 5.2) so neither should the PCI dataset, although multiple suppression should be better on the SALT data given the longer offsets. Lastly, the pattern of the deep reflections does not follow that of the overlying geological trends as they would if the multiples were to be generated in any kind of consistent manner (see easternmost crossline section in Fig. 5.4). Given these points, it is likely that the reflections are indeed being generated from within the Precambrian basement.  5.7  D o the P C I deep reflections represent further evidence of the H S I reflector?  Although there are no obvious ways to directly relate the basement reflections seen on the PCI data with the HSI reflector imaged to the northeast by Lithoprobe, a number of general consistencies are encouraging. Firstly, the reflectors appear at about the same times in both areas. Secondly, the lower amplitudes of the basement reflectors in Pincher Creek, while inconsistent with those from the main HSI body reflectors, do resemble those directly to the north along the western end of SALT line 30 (Figs. 5.2 and 5.5). Similarly, the discontinuity of the reflectors themselves also agrees with the discontinuous reflectors along the western end of line 30. Lastly, the 1-D forward modelling results suggest that the brightest reflections could be arising from a sill between 50 and 150 m which agrees with the 1-D modelling results presented by Mandler and Clowes (1998) for a section of the HSI reflector to the northeast.  5.8  Chapter summary  In this chapter, the evidence for deep reflections in the PCI data volume was presented. A statistical approach to interpreting 3-D data volumes for deep reflectors was also introduced and 1-D forward modelling of the deep reflections was undertaken. The PCI results and those from  Chapter 5. The Head-Smashed-In  Reflector  135  both the PC2 dataset (Chapter 6) and the earlier Lithoprobe surveys are used to discuss the topics of sill emplacement, craton stabilization and cratonic arch formation in Chapter 7.  Chapter 6  The Winagami Reflection Sequence  The existence of the Winagami Reflection Sequence (WRS) in central Alberta was revealed from 2-D multichannel seismic reflection lines collected by Lithoprobe in 1992. Subsequent Lithoprobe seismic investigations in northwestern Alberta in 1994 recorded further evidence of the WRS. With its horizontal to sub-horizontal bands of bright reflections confined within a depth range of 9 to 24 km (3-8 s TWTT), the Winagami reflection sequence, like the HSI reflector discussed in Chapter 5, has been interpreted as a deep mafic sill complex (Ross and Eaton, 1997). If this interpretation is correct, with an inferred areal extent of more than 120,000 km , 2  the Winagami sequence is comparable in size to many of the world's large Phanerozoic igneous provinces (Coffin and Eldholm, 1994). The considerable lateral extent of the WRS and its location within the upper crystalline crust underlying the Western Canada Sedimentary Basin provided the impetus for seeking an industry-acquired 3-D dataset that sampled the sequence. The PC2 3-D dataset, which was acquired above the proposed extent of the WRS and was recorded to 5.1 s TWTT, provided a suitable candidate for the second localized 3-D seismic investigation of sills presented in this thesis. In this chapter, the Lithoprobe evidence for the WRS is briefly reviewed and the results from the PC2 dataset are presented within this context. The regional implications of the results for the WRS are further discussed in Chapter 7.  136  Chapter 6. The Winagami Reflection Sequence  6.1  137  The C A T ' 9 2 and PRAISE'94 data  In total, the Winagami reflection sequence is observed on four of the ten multichannel lines from the Central Alberta Transect (CAT) experiment and on eight of the ten lines from the Peace River Arch Industry Seismic Experiment (PRAISE) (Fig. 6.1). As introduced in Chapter 1, the sequence is characterized by between one and five bands of bright reflections which are locally horizontal to sub-horizontal but together form a large-scale anastomosing fan that converges to the southeast (Ross and Eaton, 1997) (Fig. 6.2). Individual reflections are consistently characterized by a simple waveform with a dominant frequency of 20 Hz. From forward seismic modelling of the uppermost reflective band, Ross and Eaton (1997) infer that the top reflection arises from a horizontal body with a thickness of 70 ± 10 m. The sill interpretation for the WRS is based on seismic evidence from PRAISE line 12 that shows the Winagami reflections cross-cutting a dipping suture zone within the granitic host rocks without offset (Fig. 1.6). Due to the southeastward dipping geometry of the Winagami reflections as imaged by Lithoprobe, a dataset recorded to 5.1 s T W T T over the northern extent of the WRS will sample a greater number of the Winagami reflections than an identical dataset to the south. Nonetheless, assuming that the bands of reflectivity imaged by Lithoprobe maintain the same apparent dip of approximately 1° as they continue southward, the topmost band of reflections should still be recorded at earlier times than 5.1 s T W T T in the PC2 survey region (Fig. 6.3). Furthermore, as the reflective bands appear to converge as they extend southeastward, the possibility of recording two bands of reflectivity remains.  Chapter 6. The Winagami Reflection Sequence  -120°  -118° I  -116° I  CAT'92 seismic lines PRAISE'94 seismic lines — — Winagami outline  ; n  -120°  138  -114° !  -112°  Provincial b o r d e r - - - Cordilleran ^ deformation front  1  1  1  -118°  -116°  -114°  rc  -112°  Figure 6.1: Detailed basement domain map of the region surrounding the PC2 3-D survey area adapted from Figure 2.7. The locations of the CAT'92 and PRAISE'94 2-D seismic reflection lines collected by Lithoprobe (with line numbers shown in white circles) are also shown. The inferred bifurcation of the Snowbird Tectonic Zone is plotted in the northeast of the map. Note that the basement domains have been delineated based on my own interpretation of the aeromagnetic data.  reflective  r  bands V  host rock  Figure 6.2: Schematic diagram of the anastamosing shape of the Winagami bands of reflectivity.  Chapter 6. The Winagami Reflection Sequence  139  Winagami reflection sequence 3D Survey  Figure 6.3:  t N  Spatial relationship between the topmost Winagami reflective band o n L i t h o p r o b e lines to the north  and where it should appear on the PC2 data assuming the feature's dip remains constant.  6.2  Pre-stack evidence for deep reflectors i n the Pine C r e e k region  To assess whether the PC2 dataset contained any evidence of the Winagami reflection sequence to justify a reprocessing of the data, ten shot gathers collected in the northeastern part of the survey were sent to U B C in the summer of 2003. These raw gathers revealed evidence of at least one band of bright reflections at the appropriate times for the Winagami sequence. Based on these encouraging preliminary results, the entire PC2 3-D dataset was obtained. Bright Winagami-like reflections occurring between 4 and 5.1s T W T T are evident throughout most of the PC2 dataset on both shot (Fig. 3.21b) and C M P (Fig. 3.25a) gathers. Correspondingly, velocity semblance spectra plots point to high stacking velocities for the deep reflections (Fig. 3.25b) which suggest that they originate from within the crystalline basement.  Chapter 6. The Winagami Reflection Sequence  6.3  140  Post-stack evidence f o r deep reflectors i n the P i n e C r e e k region  Deep reflections are easily observed between 4 and 5.1 s T W T T on inline (Fig. 6.4), crossline (Fig. 6.5) and time slices (Fig. 6.6) throughout the PC2 data volume. Due to the consistency in the arrival times and the relatively high amplitudes of the deep reflections observed on most of the CDP gathers from the PC2 dataset, brute stacks of the data are able to provide an interpretable 3-D image of the deep reflector. With further processing, the reflector's coherence is improved. Figure 6.4 shows nine inline sections through the unmigrated PC2 data volume. Progressing from the southwest to the northeast of the data cube (sections a to i), these sections consistently show one and often times two Winagami-like reflections between 4.5 and 5.1 s TWTT. For sections a and c, the horizontal to sub-horizontal reflections are swamped to the southeast by the limbs of shallower diffractions. Meanwhile in b, the Winagami-like reflection crosscuts the dipping diffraction limbs. Moving northeast through the volume, sections d through i provide a clearer image of a band of reflectivity dipping to the southeast, at its deepest disappearing beyond the limit of recorded time. These dipping events often exhibit their own smaller diffractions suggesting topography on the reflector. Evidence of a second reflective band is observed at the limits of recorded time on inline sections a, b, c and f. The crossline sections shown in Figure 6.5 show a coherent reflection at roughly 4.5 s T W T T on the northwesternmost sections (a through d). Moving southeastward (sections e to f), this reflection begins to dip southwestward. Given the southeastward dip of the reflections along the inline sections (Fig. 6.4), crossline sections further southeast than section f show no Winagamilike reflections. Of note, there is evidence for a second band of reflections at the limits of recorded time on both crossline sections a and b. Time slices extracted from the data volume between 4.2 and 5.1s T W T T provide an areal view of the distribution of the reflections at different traveltimes (Fig. 6.6). Slices a through c  Chapter 6. The Winagami Reflection Sequence  141  demonstrate the concentration of reflections in the northwest of the volume at the earliest times at which the inferred Winagami reflections occur. At later times, a sheet outlined by the inline and crossline dipping reflections is hinted at from the quasi-linear pattern of reflections in slices d, e and f (highlighted in slices d and e by the arrows). In slice g, linear evidence of a dipping reflection remains in the northeast while the overall trend has mostly disappeared. Finally, near the limits of recorded time, slices h and i show evidence for a second reflective event in the northwest. While sections through a data volume are a useful interpretation tool, volume diagrams provide the spatial link between individual sections. In Figure 6.7, three different chair-cut views through the data volume provide a link between instances of reflectivity on the inline, crossline and time sections. In plot a, the crossline section (WSW to ENE) is dominated by a horizontal reflection. The intersection of that reflection with the inline section (WNW to ESE) is suggestive of a sheet dipping to the southeast. The time slice illustrates where the sheet cuts through the bottom of the diagram. Similar patterns are evidenced in plots b and c. To better illustrate the connection between the time slices in Figure 6.6 and the volume diagrams in Figure 6.7, Figures 6.8 and 6.9 show the spatial locations where time slices b and e were extracted. Time slice b cuts through the top of the sub-horizontal reflection on the crossline section of Figure 6.7a while time slice e cuts through the reflective sheet as it dips southeastward. Both the Winagami-like reflections and the shallower diffractions observed throughout the data volume have geometries which greatly differ from the horizontal reflections from the overlying sediments (arriving before 2.5 s TWTT) on all of the vertical data plots. Consequently, the later arrivals are likely not multiples but are reflections that originate from within the crystalline basement. This interpretation is supported by the velocity analyses which show that the deep reflections required higher stacking velocities than the shallower sedimentary arrivals.  Chapter 6. The Winagami Reflection Sequence  142  Reflectivity throughout the data volume is generally quite pronounced although the southeastern portion of the volume is particularly lacking in reflections. Since the amplitudes of the sedimentary arrivals are generally uniform across all the vertical data sections, lateral variations in reflectivity at later times likely represent the true variations in reflectivity. In all, the Winagami-like reflections in the PC2 data volume outline a 3-D reflective surface that is shallowest in the northwest and dips to the southeast, eventually disappearing beneath the processed data volume. Along the northwestern edge of the data volume, some evidence for a second reflective band is observed although only at the limits of the recorded time. Based on the post-stack evidence, a number of general observations can be made about the PC2 data volume:  1. Winagami-like reflections are consistently observed between 4 and 5.1s T W T T throughout the data volume 2. these reflections, particularly in the north, appear to dip southeastward with an estimated unmigrated dip of 2° 3. shallower diffractions with apexes between 3 and 4 s T W T T whose eastern limbs crosscut Winagami-like reflections are regularly observed in the volume, particularly in the southwest 4. reflectivity is least pronounced in the southeastern portion of the data volume 5. this lack of reflectivity applies to both the Winagami-like reflections and shallower diffractions 6. the Winagami-like deep reflections are not likely to be multiples based on velocity and geometrical observations.  Chapter 6. The Winagami Reflection Sequence  143  Figure 6.4: Inline sections through the unmigrated PC2 data volume plotted between 2 and 5.108 s TWTT (approximate depth of 6 to 15 km). The locations of the sections with respect to the data cube are shown at the top using the lettered labels. The sections have not been horizontally or vertically exaggerated. The locations of the reflections attributed to the Winagami sequence are shown with the black arrows. White arrows indicate shallower diffractions. Apart from the geometrical spreading correction, no gain has been applied to the data sections. On all plots, approximate depths are based on an average velocity of 6 km/s. Since velocities in the sediments are slower, the base of the sediments (which occurs between 2 to 2.5 s TWTT on the plots) isn't as deep as shown.  Chapter 6. The Winagami Reflection Sequence  144  Figure 6.5: Crossline sections through the unmigrated PC2 data volume plotted between 2 and 5.108 s TWTT (approximate depth of 6 to 15 km). The locations of the sections with respect to the data cube are shown at the top using the lettered labels. The sections have not been horizontally or vertically exaggerated. The locations of the reflections attributed to the Winagami sequence are shown with the black arrows. White arrows indicate shallower diffractions. Apart from the geometrical spreading correction, no gain has been applied to the data sections. On all plots, approximate depths are based on an average velocity of 6 km/s. Since velocities in the sediments are slower, the base of the sediments (which occurs between 2 to 2.5 s TWTT on the plots) isn't as deep as shown.  Chapter 6. The Winagami Reflection Sequence  145  Figure 6.6: Time slices through the unmigrated P C 2 data volume. The time at which each slice was extracted is shown on the TWTT scale on the right using the lettered slice labels. The arrows highlight the quasi-linear pattern of reflections attributed to a dipping sheet.  Chapter 6. The Winagami Reflection Sequence  146  Figure 6.7: Chair-cut views through the unmigrated P C 2 data volume. Arrows indicate regions of coherent reflectivity. Apart from the geometrical spreading correction, no gain has been applied to the data.  Chapter 6. The Winagami Reflection Sequence  147  Figure 6.8: Spatial tie between time slice (b) of Figure 6.6 and the chair-cut view through the unmigrated PC2 data volume.  Chapter 6. The Winagami Reflection  Sequence  148  Figure 6.9: Spatial tie between time slice (e) of Figure 6.6 and the chair-cut view through the unmigrated PC2 data volume. The arrows highlight the quasi-linear pattern of reflections attributed to a dipping sheet.  Chapter 6. The Winagami Reflection Sequence  6.3.1  149  Migration tests  Migrating the deep reflections of interest in the PC2 data volume provides three challenges. Firstly, as explained in Section 3.12.1, many of the arrivals needed to properly migrate deep reflections such as diffractions from discontinuous edges, can be lost due to the greater amount of attenuation they experience as they travel the longer distances back to the surface and, if present, can also be distorted by near-surface lateral velocity variations. The second challenge specific to the PC2 dataset is that the reflections of interest are at the limit of the recorded time. Consequently, much of the diffracted energy required for migration arrived too late to be recorded. Thirdly, for the late reflections, the migration aperture is too small to adequately collapse diffractions. To assess the impact of these influences on the migration results for the PC2 dataset, 2-D post-stack Kirchhoff time migration with a constant velocity function was performed on two inline sections through the data volume using a range of velocities. The original stacked data along with the migrated results are shown in Figures 6.10 and 6.11 for a northern and southern section respectively. While a migration velocity between 5000 and 5500 m/s is an appropriate choice for reflections within granitic host rocks overlain by thick sediments, the results look severely over-migrated (i.e. smiley). In fact, over-migration seems to be a problem for migration velocities greater than 3000 m/s. Further aggravating the situation are the ends of the stacked traces which are over-migrated into smiles that interfere with the deep reflections of interest. These migration artifacts were minimized to an extent by the application of a temporal taper to the ends of the traces prior to migration. Despite the presence of multiple migration artifacts, the deep Winagami-like reflections of interest remain coherent regardless of the migration velocity used. While the positions of the sub-horizontal reflections in the southwestern part of the data volume are unaffected by migration (Fig. 6.11), the dipping reflection at 4.5 s T W T T in Figure 6.10 is moved up-dip in the plane  Chapter 6. The Winagami Reflection Sequence  150  of the inline section for velocities of 3000 m/s and greater. Meanwhile, the lack of recorded data beyond 5.1s T W T T results in the southeastward dipping edge of the migrated reflection appearing stunted. While these 2-D post-stack time migration results show some re-positioning of the deep reflections of interest, the improvements are slight compared to the harmful incorporation of migration artifacts into the data volume. Furthermore, given that the dip of the reflective sheet as evidenced in the unmigrated stacked volume is not aligned with either the inline or the crossline direction, 3-D migration of this reflective feature would appear to be essential. Lacking such an algorithm, a 3-D forward modelling approach to determining the geometry of the reflective sheet is undertaken in Section 6.6. While the 2-D migration tests on the deep reflections associated with the Winagami reflection sequence did not aid in its interpretation, the shallower diffractions (shown by white arrows in Figs. 6.10 and 6.11) have benefited to an extent. Using migration velocities from 3000 to 6000 m/s, these diffractions have collapsed to truncated horizontal to sub-horizontal reflections. Due to their stunted reflective characters and their earlier traveltimes, these are interpreted as disconnected reflections arising from inhomogeneities in the gneissic basement that are not associated with the Winagami sequence. This interpretation is supported by the results from line 20 of Lithoprobe's PRAISE experiment which show similar evidence of upper crustal diffractions separate from the coherent Winagami reflections (Fig. 6.12). It should be stressed however that the shallower diffractions could be originating offline from within the sedimentary sequence or the shallow basement. Thus, 3-D migration is required to properly collapse and locate the sources of the diffractions.  Chapter 6. The Winagami Reflection Sequence  151  Figure 6.10: Inline 2-D post-stack Kirchhoff time migration results for a northern inline section of the PC2 data volume using different velocities (shown in the top left corner of each plot). The location of the reflection attributed to the Winagami sequence is shown with the black arrow on the unmigrated stacked section at the top left. A shallower diffraction is indicated by the white arrow. All sections have had a 1000 ms AGC gain applied. On all plots, approximate depths are based on an average velocity of 6 km/s. Since velocities in the sediments are slower, the base of the sediments (which occurs between 2 to 2.3 s TWTT on the plots) isn't as deep as shown.  152  Chapter 6. The Winagami Reflection Sequence  stacked  BJUw  1000 m / s ^ ^ ^ ^ ^ s H XI XI  -~~ 3  3  9  J*  a a> xi  .  *f  j  V  m  15  >  XI XI  ~  3  M  BBifs- • • lliSI!lt*lP  '  mSr 5000 m / s ^ g ^ r ? '  -3T  3  •  li 15  Figure 6.11: Inline 2-D post-stack Kirchhoff time migration results for a southern inline section of the PC2 data volume using different velocities (shown in the top left corner of each plot). The location of the reflections attributed to the Winagami sequence are shown with the arrows on the unmigrated stacked section at the top left. A shallower diffraction is indicated by the white arrow. All sections have had a 1000 ms AGC gain applied. On all plots, approximate depths are based on an average velocity of 6 km/s. Since velocities in the sediments are slower, the base of the sediments (which occurs between 2 to 2.3 s TWTT on the plots) isn't as deep as shown.  Chapter 6. The Winagami Reflection  Sequence  153  Figure 6.12: Comparison of basement diffractions from the Lithoprobe and PC2 results (a) without and (b) with migration using a constant velocity of 6000 m/s. The same horizontal scale is used for both plots. The Winagami reflections are shown with the black arrows while shallower diffractions are indicated by the white arrows. On all plots, approximate depths are based on an average velocity of 6 km/s. Since velocities in the sediments are slower, the base of the sediments (which occurs between 2 to 2.3 s TWTT on the plots) isn't as deep as shown.  Chapter 6. The Winagami Reflection Sequence  6.4  154  P o l a r i t y constraints  Due to the significant depth at which the WRS resides, no direct physical sampling of the Winagami complex is possible. Nonetheless, further information about the nature of the deep reflectors can be inferred from relationships between the shallow reflections and nearby well logs. Specifically, such comparisons can provide polarity constraints which can help characterize the nature of the impedance contrast from which the reflections have originated. The PC2 survey area has been the focus of oil and gas exploration for over twenty years. Consequently, the area has been densely sampled using well logs. Density and sonic logs for a well at the center of the survey are shown in Figure 6.13b. The density and sonic logs are multiplied to produce the acoustic impedance log. Furthermore, the density and sonic logs are used to calculate a log of vertical incidence reflection coefficients according to  ; PlOLl + P2&2  (6.1)  where p\,a.\, p and a are the density and P-wave velocity values above and below the reflec2  2  tive boundary respectively. Once calculated, the log of reflection coefficients is convolved with a source wavelet to produce a synthetic trace. Using the Kingdom Suite software from Seismic Micro-Technology Inc., the source wavelet in Figure 6.13a was estimated from the amplitude spectrum of the near-most seismic trace to the well location. Only arrivals from the sedimentary sequence corresponding to the sampled depth of the well log were used in the estimation. Since the estimated wavelet was output as zero phase by the software and the resulting synthetic seismogram showed little resemblance to the recorded data, constant phase shifts were applied to the source wavelet until the highest correlation coefficient between the synthetic and recorded trace was obtained. As such, the source wavelet shown had a 154 degree constant phase shift applied, which is why it is dominated by a negative peak at zero time. Both the resulting synthetic trace and the recorded trace correspond well with a correlation coefficient of 0.727 and  Chapter 6. The Winagami Reflection Sequence  155  most reflections can be aligned. The Nordegg horizon highlighted in Figure 6.13 marks the transition from overlying low velocity, low density shales and sandstones to higher velocity, higher density carbonates, as evidenced by the sharp increases on both the sonic and density logs. These increases result in a strong positive reflection coefficient. From the comparison between the synthetic seismogram and the recorded trace, the Nordegg positive reflection coefficient corresponds to a negative peak. This is not surprising given the dominant negative peak at zero time for the source wavelet (Fig. 6.13a). Assuming that this relationship holds for the entire data volume, positive reflection coefficients should correspond with dominant leading negative peaks throughout the PC2 data volume. The Nordegg horizon and its corresponding logs from Figure 6.13 are compared against a portion of an inline section through the stack in Figure 6.14. An enlarged segment of the stack (outlined in gray) demonstrates that the leading dominant peak for the deep reflection appears negative. This is reinforced by the horizontal sum of the highlighted traces which clearly shows a negative peak as the first dominant arrival (shown by dark gray arrow). If the relationship derived for the Nordegg horizon holds true, the polarity of the deep reflection suggests that it is arising from a positive reflection coefficient. Since doleritic material has both higher velocities and densities than the host granitic gneisses, these results support the interpretation that the reflections are arising from mafic sills.  6.5  Waveform modelling  The thickness of the topmost sill body from the Winagami reflection sequence was previously estimated at 70 ± 10 m (Ross and Eaton, 1997). With multiple data sections containing the Winagami-like reflections, the PC2 data volume offers an opportunity to examine changes in the reflector's thickness over the span of the survey using simple 1-D waveform modelling.  Chapter 6. The Winagami Reflection Sequence  156  Time (ms)  a) •3 0.53  =5.  E  0  i 11 I \)  41  \7  <-0.5 -  h U  -1 Sonic Density (m/s) (kg/m*3}  v J  1000 5000 2000 3000 Al  RC  Synth.  Real  1200"  1400-  1600-  1800-  Q. 2000Q 2200-t  2400-  2600-I  NORDEGG  1.6  Figure 6.13: a) Source wavelet extracted from the recorded trace closest to the well location, b) well logs from the overlying sedimentary sequence used to obtain polarity constraints for the WRS. From left to right, shown are the sonic log, the bulk density log, a log of acoustic impedance (AI), a log of calculated reflection coefficients (RC), a calculated zero-offset synthetic seismogram (repeated five times for clarity) and the near-most recorded data to the well location (repeated five times for clarity) for comparison. The highlighted Nordegg horizon, a high velocity/high density contrast with the overlying sediments, is used as the reference horizon and is characterized by a strong negative peak.  Chapter 6. The Winagami Reflection Sequence  157  Figure 6.14: Polarity comparison between synthetic seismograms and the WRS reflections. Top left plot shows a section extracted from Fig. 6.13 highlighting the reference Nordegg horizon. Bottom plot shows an unmigrated inline section through the stacked PC2 data volume. A portion of the section outlined by the gray box is enlarged for comparison with the synthetic waveforms in the top plot. The horizontal stack of these highlighted traces is also plotted for comparison both in the highlighted box and again at the top right of the figure (labeled WRS real) where it is repeatedfivetimes for clarity and plotted with the same scale as the traces in the top left. The similarity in the reflection character and polarity suggests that the WRS results from a high velocity/high density contrast.  Using the wavelet from Figure 6.13 as the source wavelet, simple 1-D synthetic seismograms were constructed using a range of different three-layer 1-D models of a sill intruded into a host rock (Fig. 6.15). As with the waveform modelling for the PCI dataset, the algorithm for generating the synthetic traces involved computing the reflection coefficients at the appropriate zero-offset traveltimes and convolving them with the source wavelet. Multiples were not included in the generation of the synthetic traces. As with the earlier modelling results, the sill and the host rock material properties were assigned velocity and density values of 6000 m/s and 2650 kg/m for the granite, and 6500 m/s and 3000 kg/m for the dolerite, respectively. 3  3  Chapter 6. The Winagami Reflection Sequence  158  1  •S 0.5  (  I  3  J Vr  =5. 0 E <-0.5  1  •1 -200  c \>  -100 0 100 Time (ms)  200  b)  c) F l  10  50 25 150 Sill thickness (m) Granite  200  Dolerite  Figure 6.15: 1-D Modelling results for the PC2 dataset; a) the wavelet from Figure 6.13, b) computed 1-D seismograms using the wavelet and the simple three layer models of different thicknesses of dolerite sills emplaced in granite shown in (c). The models were constructed using velocity and density values of 6000 m/s and 2650 kg/m for the granite, and 6500 m/s and 3000 kg/m for the dolerite, respectively (Juhlin, 1990). 3  3  Chapter 6. The Winagami Reflection Sequence  159  Sill thicknesses of 5 m, 10 m, 25 m, 50 m, 100 m, 150 m and 200 m were used to span the range of thicknesses deduced from the earlier studies. The dominant frequency of the wavelet in Figure 6.15a is approximately 20 Hz. From the discussion of temporal resolution in Section 3.14, the dominant wavelength for a medium of 6000 m/s and a frequency of 20 Hz corresponds to 300 m. Therefore, according to Sheriff and Geldart (1999), reflective boundaries must be at least 75 m apart to be distinguished as separate. From the results from Figure 6.15, this estimate seems slightly optimistic. While sills with thicknesses between 150 and 200 m clearly show distinct negative peaks for the top and bottom reflections, sill thicknesses of 100 m and less simply look like stretched versions of the source wavelet. Nonetheless, since the length of the wavelet can be indicative of the sill thickness, thickness estimates can still be obtained for values less than 150 m. Figure 6.16 shows the comparison between the 1-D modelling results (Fig. 6.15) and examples of coherent reflectivity from three sections through the data volume. The horizontal stack of the highlighted traces from the southwestern section (Fig. 6.16a) reveals a wavelet with two positive peaks flanking a dominant negative peak. From comparison with the 1-D modelling results, this wavelet most resembles the results for a sill with a thickness of 100 m both in terms of wavelet length and character. Similarly, both the middle and northeastern sections also produce horizontal stacked traces with wavelets of similar length and reflective character. From these combined results, the Winagami sill imaged by the PC2 dataset is inferred to have a constant thickness of approximately 100 m throughout the survey area.  Chapter 6. The Winagami Reflection Sequence  160  Figure 6.16: Comparison of 1-D modelling results from Figure 6.15 with coherent reflections from a) the southwest, b) the center and c) the northeast of the PC2 data volume. Shown for each section are: the location of the section (top left with shown data section highlighted in dark gray), the data section (top right), the modelling results (bottom left), the traces highlighted by the gray box in the data section (bottom right), and the horizontal stack of the highlighted traces (bottom center). The modelling results and the stacked data are plotted to the same scale.  Chapter 6. The Winagami Reflection Sequence  6.6  161  M o d e l l i n g the W i n a g a m i reflection sequence i n 3-D  Although the migration results for the deep reflections in the PC2 data volume provide no more interpretable information than the stacked results, further insight can be obtained by forward modelling of the sequence. To this end, the kmod3d code developed by David Eaton at the University of Western Ontario was used (Eaton and Clarke, 2000). This efficient program is based on the Kirchhoff integral method and calculates the scattered 3-D wavefield generated by a body of arbitrary shape in a homogeneous elastic medium. The parameters of the wavefield are determined for any point in the medium using the wavefield and its normal derivatives on the closed surface of the scattering body. A number of high frequency assumptions are made to simplify the solution. The parameterization for kmod3d involves approximating the surface of the scattering body using plane triangular surface elements (Eaton and Clarke, 2000). For each triangle, the wavefield is approximated as a plane wave. The influence of each of the triangle's vertices toward scattering is quantified by a weighting function and the triangle's total contribution is determined by numerical integration. For modelling the top reflection of the inferred Winagami sequence in the PC2 data volume, the scattering surface is simply specified as a layer that extends laterally to the bounds of the model space. While kmod3d can handle coarse sampling of the scattering surface on the order of the Fresnel zone radius (approximately 1.5 km for the deep reflections), the model was gridded using a 1 km by 1 km mesh to avoid computational instabilities. An initial model with variable sill topography was constructed based on geometrical variations observed from the stacked volume. To properly locate the resulting reflections in time, the velocity of the background medium was set to 5500 m/s, a depth-weighted average of the velocities of the sediments and the basement in the PC2 region. The contrast in velocities and densities across the scatterer's surface were assigned the same values as were used for the 1-D  Chapter 6. The Winagami Reflection Sequence  162  seismic modelling in Chapter 5 and in the previous section. Thus, the model was constructed using velocity values of 6000 m/s for the granite and 6500 m/s for the dolerite. For the densities, values of 2650 kg/m for the granite and 3000 kg/m for the dolerite were used. 3  3  The flexibility of kmod3d is particularly evident in its ability to handle any source and receiver configuration. These are placed along the top of the model and can be specified with topography. Synthetic traces are calculated for every shot and receiver pairing listed in a parameter file. To reproduce some of the stacked results from Figures 6.4 and 6.5, coincident source and receiver pairs were located above the sill model at locations corresponding to the CDP bins of the stacked traces. While multi-component data can be generated using kmod3d, only vertical component traces were computed. The goal of the 3-D modelling of the top of the reflective sheet was to construct a model whose geometry would reproduce the shapes and to some extent the character of selected inline and crossline sections through the stacked volume. By iterative fine adjustments of the geometry of the top of the sill body in the initial model, a 3-D model of the sill is generated (Fig. 6.17). Since a dataset recorded to 5.1 s T W T T will not contain reflections from reflectors deeper than 16 km, only the portion of the model surface above 16 km depth is shown. Figure 6.18 shows a comparison between inline sections through the stacked PC2 data volume and the computed traces for coincident sources and receivers along the same inlines (shown as labeled red lines in Fig. 6.17). Along inline A, the truncation in the reflection at 4.5 s T W T T is closely reproduced by the modelled traces. To achieve this fit, the eastward dip of the modelled surface was made to drop off suddenly. Further to the northeast along inline B, the southeastward dip of the reflector is made less steep to capture the dip of the reflections. As the southeastward dip of the model is made more gradual to the northeast, a good match is achieved between the sub-horizontal reflections along the stacked and modelled data sections (inlines C and D).  Chapter 6. The Winagami Reflection Sequence  163  sill depth (km) Figure 6.17: 3-D model of the top of the Winagami sequence from 3-D Kirchhoff modelling. Since a dataset recorded to 5.1 s TWTT will not contain reflections from reflectors deeper than 16 km, only the portion of the model surface above 16 km depth is shown. For the parameterization, the model surface extended across model space. The locations of the inline model sections from Figure 6.18 are shown in red and labeled to correspond with thatfigure.Similarly, the locations of the crossline model sections from Figure 6.19 are shown in green and labeled to correspond with thatfigure.Receivers lines from the PC2 survey are overlain for reference.  Chapter 6. The Winagami Reflection Sequence  164  Figure 6.19 shows the equivalent comparison between crossline sections through the stacked PC2 data volume and the computed traces for coincident sources and receivers along the same crosslines (shown as labeled green lines in Fig. 6.17). Along crossline A, an excellent match is achieved between the modelled and stacked data in terms of both arrival times and the geometry of the reflections. The loss of amplitude at either end of the recorded data section is due to a drop in the stack fold. The results from crossline B are less encouraging with none of the strong horizontal reflections appearing on the stacked data. However, the southwestward dip of the reflection from the stacked data does correspond in shape but not amplitude with the reflections arising from the abrupt change in dip orientation of the model geometry at that location. Crossline C shows excellent correspondence between the model results and the reflections recorded in the northeast. Meanwhile, the computed reflections in the southwest are not seen in the stacked section. However, if present, such reflections would likely be lost in the noise. The modelled reflections from Figure 6.19C illustrate the need for 3-D migration. As evidenced by Figure 6.17, crossline C does not overlie the model surface within 16 km. Consequently, the modelled and stacked arrivals in the northeast of the data sections are originating offline from further to the northwest. This makes complete sense considering that dipping features on stacked data will be migrated up-dip. Thus, had it been possible to properly migrate the PC2 data volume in 3-D, the Winagami reflections along crossline C would have been moved updip to the northwest and no deep reflections would be visible on the migrated version of the section.  Chapter 6. The Winagami Reflection  Sequence  165  Figure 6.18: 3-D Kirchhoff modelling results for several of the stacked inlines shown in Figure 6.4. The locations of the displayed sections with respect to the data cube are shown on the right using the lettered labels.  Chapter 6. The Winagami Reflection  wsw  Sequence  166  ENE  Figure 6.19: 3-D Kirchhoff modelling results for several of the stacked crosslines shown in Figure 6.5. The locations of the displayed sections with respect to the data cube are shown on the right using the lettered labels.  Chapter 6. The Winagami Reflection Sequence  167  While the kmod3d modelling results are non-unique, the final 3-D model does generally reproduce both the geometry and the reflective character of the deep reflections in the stacked volume. Also, in the absence of a 3-D migration code, the modelling results allow for an estimation of the proper location and geometry of the reflective sheet. From these results, the fortuitous location of the PC2 survey for modelling the Winagami sills is highlighted. Assuming that the reflections continue their southeastward dip as modelled, had the survey between acquired but 10 km to the southeast, no evidence for deep Winagami reflections would have been detected.  6.7 Lack of reflectivity in the southeast The 3-D forward modelling results reinforce the conclusion that the PC2 volume captures reflections from a southeastward dipping reflective sheet which dips below the limits of recorded time in the southeastern portion of the data volume. In addition to losing the Winagami-like reflections, the southeast portion of the volume beneath the base of the sediments appears devoid of basement reflectivity entirely. To quantify this lack of reflectivity and obtain a spatial view of the distribution of reflectivity within the PC2 data volume as a whole, the statistical approach to event picking developed for the PCI data in Chapter 5 was applied. Again, a spatial coherency filter (Kong et al., 1985) was applied to all the inline and crossline sections in order to automatically extract all linear events that spanned at least 3 km (the approximate Fresnel zone at 4.5 s for a 10 Hz signal in a medium velocity of 5000 m/s). Once these events were isolated, only linear reflectors from the inlines and crosslines that intersected each other (within ± 1 0 ms) were selected. Common picks that were not bordered by at least one other common pick (within ± 20 ms) were subsequently eliminated to impose a minimum areal constraint on the picked common reflectors.  Chapter 6. The Winagami Reflection Sequence  168  For the PC2 dataset between 3.2 and 5.108 s TWTT, a total of 763,266 coherent picks were obtained for all the inline sections and 141,709 picks for all the crossline sections. From these, a total of 90,552 common picks were found to intersect within 10 ms of one another. With the added areal constraint of requiring common picks to be bordered by at least one other pick, the number of common picks was reduced to 81,966. The spatial distribution of the number of common picks per stacked trace is shown in Figure 6.20A. The plot clearly reinforces qualitative observations of the seismic plots that the southeast of the survey is devoid of reflections. Curiously, the aeromagnetic pattern over the survey contains a strong positive anomaly exactly in the location where reflectivity is lost (Fig. 6.20B). Since the region's distance from the magnetic pole can cause magnetic anomalies to appear shifted laterally or distorted due to phase contributions from the non-vertical dipole (Blakely, 1996), the data were reduced to pole using a subroutine provided by Blakely (1996). This transformation readjusts the data to what they would look like if measured at the Earth's pole where the ambient magnetic field is vertical. The ambient magnetic field in the PC2 region was taken to have an inclination of 76.5° and a declination of 67.3° (Hope and Eaton, 2002). Since the inclination and declination of the remanent magnetic field in the region are unknown, these were assigned the same values as the ambient field. The results of the reduction to pole are shown in Figure 6.20C. From parts A and C from Figure 6.20, the lack of reflectivity in the southeast and a corresponding aeromagnetic high may somehow be linked. As introduced in Chapter 2, the PC2 3-D survey straddles a zone that has been interpreted as the boundary between the Chinchaga and the Wabamun domains across the northern splay of the Snowbird Tectonic Zone. This boundary is generally interpolated between strong aeromagnetic lows to the west and east of the PC2 survey region. However, immediately in the vicinity of the survey, the signal becomes more complicated with high aeromagnetic anomalies from the Wabamun domain (known as the Wabamun High) extending into the aeromagnetically low  Chapter 6. The Winagami Reflection Sequence  Figure 6.20:  169  Spatial link between lack of reflectivity in the southeast and the aeromagnetic signal. A) Spatial  distribution of the number of reflector picks common to both inline and crossline sections per stacked trace; B) aeromagnetic signal over the PC2 survey area; C) reduced to pole aeromagnetic signal over the PC2 survey area. The PC2 acquisition geometry is shown on all plots for reference.  Chapter 6. The Winagami Reflection Sequence  170  Chinchaga domain (known as the Chinchaga Low) (Fig 6.21). Since the Wabamun High is characterized by up to 20 km wide sub-circular positive anomalies bordered by narrow lows (Villeneuve et al., 1993), the aeromagnetic signature around the PC2 survey (Fig. 6.20C) would suggest that the dataset is entirely located within the Wabamun domain and that the Snowbird Tectonic Zone, if present, exists to the north. The aeromagnetic highs within the Wabamun have been interpreted as being caused by undeformed magmatic rocks (Villeneuve et al., 1993). Given the abrupt cessation of reflectivity corresponding to the high anomaly in Figure 6.20, two end-member emplacement scenarios can be envisaged. The emplacement of the magmatic rocks could have predated the intrusion of the Winagami sills which for some reason were later forced around the magmatic rocks. Conversely, the emplacement of the magmatic rocks could have postdated the Winagami sills and so overprinted their reflections. The latter scenario would help explain the abrupt southeastward truncation of Winagami-like reflections in the southwest of the data volume (Fig. 6.4). Unfortunately, there is no conclusive way to determine why the southeast of the PC2 data volume is lacking in reflectivity. While the zone without reflectivity corresponds presently with an aeromagnetic high, with only one basement drill core dated at 2.3 Ga from the southeastern portion of the Wabamun domain (Villeneuve et al., 1993), inferences made about the composition and age of the aeromagnetic highs in the region are purely speculative.  171  Chapter 6. The Winagami Reflection Sequence  -116"  -118-  -116"  -114-  -114-  -112-  -110*  -112'  Figure 6.21: Interpreted aeromagnetic map of northern Alberta. A gap in data coverage in the north is left blank. Basement domain boundaries inferred from the aeromagnetic signatures are overlain by white lines. The locations of the CAT and PRAISE 2-D reflection lines collected by Lithoprobe are also shown for reference. The Snowbird Tectonic Zone and its inferred bifurcation are plotted as both thick solid and thick dashed gray lines in the northeast of the map.  Chapter 6. The Winagami Reflection Sequence  6.8  111  Do the P C 2 deepreflectionsrepresent further evidence of the Winagami reflection sequence?  The traveltimes, dips and amplitudes of the deep reflections from the PC2 dataset suggest that they represent the southward extension of the top of the Winagami reflection sequence imaged by Lithoprobe. Assuming that the topmost reflection from the PRAISE experiment maintains its apparent dip southward, the spatial link between the PC2 and the Lithoprobe results can be made despite the offset between the surveys (Fig. 6.22). Note that both the base of the sediments and the top of the WRS from the Lithoprobe data project to the equivalent reflections on the PC2 data section. One curious aspect of this comparison is that the apparent dip of the deep reflection from the PC2 data section is approximately 2° and thus is steeper than the constant 1° apparent dip used to tie the surveys together. This discrepancy may simply reflect a localized change in the geometry of the reflective sheet, perhaps associated with the cause of the absence of basement reflectivity observed in the southeastern portion of the data volume.  6.9  Chapter summary  In this chapter, the evidence for deep reflections in the PC2 data volume was presented. These results plus those for the PCI dataset presented in Chapter 5 are used to discuss the topics of sill emplacement, craton stabilization and cratonic arch formation in Chapter 7.  Chapter 6. The Winagami Reflection Sequence  173  Figure 6.22: Spatial link between the Winagami reflection sequence and the deep reflections in the PC2 data volume, assuming the feature's dip remains constant. The base of the sediments can also be linked from the Lithoprobe results to the PC2 data volume.  Chapter 7  Discussion  Following from previous interpretations of the upper crustal reflection sequences in western Canada (Mandler and Clowes, 1997; Ross and Eaton, 1997; Mandler and Clowes, 1998) and based on the results from this thesis, the mafic sill interpretation for the HSI and Winagami reflections is supported. While polarity constraints from the Winagami reflections in the PC2 dataset would also be consistent with a mylonitic shear zone as the source of the reflections, the sill interpretation is preferred based on the sharpness of the reflections, the consistency of the waveform throughout the data volume and the locally variable geometry of the reflections. Based on the sill interpretation, the following sections present a discussion of the emplacement mechanisms, the timing of emplacement and the tectonic role of the HSI and Winagami sill complexes.  7.1  Role of sills in craton stabilization  Cratons are defined as stable portions of the Earth's crust which have not been affected by orogenic activity for over 1 billion years (Allaby and Allaby, 1991). In the process of cratonization, the emplacement of mafic sills within the continental crust is thought to be a fundamental process (Nelson, 1991). These sills are inferred to originate from the underplating of continental crust by basaltic magmas followed by melting of the lower crust. Through repeated underplating events, cratons become thicker and are stabilized. Investigating sill occurrences within the crust can provide insight into this cratonization process and can shed light on the tectonic conditions and stress regime during emplacement. 174  Chapter 7. Discussion  175  Nelson (1991) argued that the formation of stable cratons follows a unified scheme in which, • new basaltic material is added to the crust at island arcs • island arcs are amalgamated into larger continental blocks through collision • lithospheric delamination during collision causes upward advection of asthenospheric material which in turn causes thermal uplift of the orogen • melting of the lower crust following contact with the asthenospheric material results in mafic underplating of the crust and sill emplacement within the crust • orogenic collapse follows once collisional forces have ceased. Each of these events can influence the composition of the crust, its structure and consequently its seismic character.  7.2 Sill emplacement mechanisms The emplacement of basaltic rock into the continental upper crust results from the upward transport of magma from the base of the crust or the upper mantle along vertical conduits or dikes. To explain how these dikes are reoriented into horizontal sill intrusions, a variety of emplacement mechanisms have been proposed. From observations of dolerite sills intruded into sedimentary sequences, Bradley (1965) postulated that sill emplacement is gravitationally driven and that density contrasts between basaltic magma and young sediments dictate the depth of emplacement. According to this commonly accepted mechanism, magma is intruded along a contact such as a bedding surface and essentially floats the sill roof. As disputed by Gretener (1969), this mechanism would result in most sills being emplaced at or near the basement-sediment contact, a situation which is not observed to be the case. Instead, Gretener (1969) provides an alternative mechanism in which  Chapter 7. Discussion  176  sills are emplaced when the crust and/or sediments are in a state of horizontal compression. Provided that the least principal compressive stress is oriented vertically, magma will take the path of least resistance and intrude horizontally. Observations by Parsons et al. (1992) of dikes intruding along planes perpendicular to the least principal compressive stress support Gretener's emplacement mechanism. In zones of tectonic extension, these dikes will extend vertically all the way to the Earth's surface while in regions of confined tectonic compression where the least principal compressive stress is vertical, the dikes will be reoriented horizontally and will produce sheet-like intrusions. According to Zoback et al. (1989), most intraplate regions of the world are presently characterized by compressional stress regimes with extension only occurring in thermally uplifted regions. Assuming that such conditions also prevailed during the Proterozoic, occurrences of horizontal sills should be widespread in intraplate regions of Precambrian crust that have experienced magmatism. Since tectonic stresses typically are uniform over distances many times the thickness of the elastic part of the lithosphere (down to 10-35 km) (Zoback et al., 1989), the stress distribution within plate interiors is dominated by the net forces associated with the regional tectonics. Regional tectonics will therefore dictate large-scale intrusional patterns. Meanwhile, local stress variations which have wavelengths that are a fraction of the lithospheric thickness, will control the local variations in the intrusional pattern. At depth in the crystalline crust, few open fractures exist due to the high confining pressures. Consequently, intruding magma must wedge apart the host rocks to provide itself with a pathway. From similarities with shallower fracturing associated with highly pressurized water, the sill emplacement process is often likened to hydraulic fracturing (Philpotts, 1990).  Chapter 7. Discussion  7.3  177  Emplacement of the H S I reflector  The emplacement of the HSI reflector sequence as mafic sills has been previously interpreted as resulting from pressurized injection of magma into a Theologically layered upper-middle crust under long-term horizontal tectonic stress (Mandler and Clowes, 1997). Given the compositional uniformity of the Medicine Hat Block into which the HSI was emplaced and the significant inferred lateral extent of the HSI, large-scale tectonic horizontal compression likely drove the emplacement. Given that the Pincher Creek region and the SALT'95 seismic lines of interest both overlie the Medicine Hat Block, similar upper crustal conditions likely exist and a similar explanation for emplacement of the sills is justified. Pollard et al. (1975) argued on the basis of outcrops of intrusions in sedimentary regions that local inhomogeneities and stress field variations could lead to complex intrusional patterns. Since the western edge of the Medicine Hat Block likely experienced the greatest variations in local stress relative to its center where most of the HSI resides, these stress conditions could result in the reduced amplitudes and discontinuous character of the reflections imaged on both the western end of SALT line 30 and the PCI dataset. Consistent with this interpretation, the similar reflective characteristics along line 30 and in the PCI volume may simply be imaging the tapering edge of the reflector sequence itself, whatever the cause of its diminishing character. If Nelson's link between continental underplating and upper crustal sills is accepted, the presence of a 20 km thick Proterozoic underplate, inferred on the basis of Lithoprobe refraction interpretations to lie beneath the Medicine Hat Block (extending from SAREX shot 3 southward on Fig. 5.1; (Clowes et al., 2002; Gorman et al., 2002)), may provide an explanation for the HSI mafic intrusions. The underplated region does not extend north beneath the Vulcan structure but does thicken to the south. This could indicate that the HSI's emplacement is related to this Proterozoic underplate and that the HSI intrusions likely extend southward.  Chapter 7.  7.3.1  Discussion  178  R e g i o n a l tectonic activity d u r i n g emplacement  As mentioned in Chapter 2, from the lack of seismic evidence for the HSI reflector extending into the Vulcan Structure whose formation was complete by 1.8 Ga (Burwash et al., 1962; Villeneuve et al., 1993), emplacement of the HSI sills is inferred to have occurred prior to 1.8 Ga (Mandler and Clowes, 1998). Without any cross-cutting relationships to provide an estimate for the earliest time of emplacement, timing constraints must be inferred from the regional tectonics at the inferred time of emplacement. The present configuration of tectonic domains is shown in Figure 7.1 for reference.  -120°  -110°  Figure 7.1: Schematic representation of the present tectonic configuration adapted from Mueller et al. (2002). Tectonic domains and features: GFTZ, Great Falls Tectonic Zone; HP, Hearne Province); MHB, Medicine Hat Block; THO, Trans-Hudson Orogen; VS, Vulcan Structure; WP, Wyoming Province.  Following from Clowes et al. (2002), the tectonic evolution of southwest Laurentia spanned the Late Archean to the Paleoproterozoic (Fig. 7.2). By 2.6 Ga, the Medicine Hat Block had been formed from the collision of two Archean blocks along a westward dipping subduction zone (Lemieux et al., 2000). Later, the Archean Wyoming Province in the south collided with the assembled Medicine Hat Block along a northeastward subduction zone which ultimately  Chapter 7.  Discussion  179  oceanic Y  MHB-E  MHB-W  \ •  •  1  MHB •  A  •  •  MHB •  ;  • GFTZ  WP  WP 2.6 Ga  WP  2.5 Ga  2.1-1.8 Ga  H (LB) VS  GFTZ  GFTZ WP  WP  \  O DB  DB 1.85 Ga  1.80 Ga  n  Collisional magmatic arc  • Subduction zone Figure 7.2:  v  T  ~) HSI evidence ) Mafic underplate evidence  Representation of the tectonic evolution of the Medicine Hat Block and its surroundings from 2.6 to  1.8 Ga. The figure was inspired and adapted from Clowes et al. (2002). The detailed description of the figure is provided in the text. Tectonic domains and features: D B , Dakota Block; G F T Z , Great Falls Tectonic Zone; H(LB), Hearne (Loverna Block); M H B , Medicine Hat Block; VS, Vulcan Structure; WP, Wyoming Province.  developed into the Great Falls Tectonic Zone. The age of this collisional event remains a topic of controversy. While Clowes et al. (2002) argue that the collision took place in the Archean, Mueller et al. (2002) suggest a Paleoproterozoic collision age of 1.9 Ga based on dating of zircons. Depending on which scenario took place, the collision between the Medicine Hat Block and the Loverna Block of the Archean Hearne province to the north along a short-lived northward dipping subduction zone either followed or was coeval with the collision of the Wyoming Province and Medicine Hat Block. The formation of the Vulcan Structure, a crustal collisional  Chapter 7. Discussion  180  zone between the Medicine Hat Block and the Loverna Block, occurred along a southward (Eaton et al., 1999a) or northward (Clowes et al., 2002; Gorman et al., 2002) dipping subduction zone (with both possible directions shown in Fig. 7.2) between 2.1 and 1.8 Ga (Burwash et al., 1962; Villeneuve et al., 1993). Finally, by 1.85 Ga, in response to the western Trans-Hudson Orogen to the east, a westward dipping subduction zone developed between the Wyoming Province and the Archean Dakota Block to the southeast. The lack of a lower bound on the timing of emplacement for the HSI provided by Mandler and Clowes (1998) ignores the link between mafic underplating and sill emplacement within the crust (Fig. 7.3a). From the dating of lower crustal xenoliths in the Medicine Hat Block, the mafic underplate beneath southern Alberta and northern Montana is inferred to have developed between 1.85 and 1.72 Ga (Davis et al., 1995). If the emplacement of the HSI sequence was coeval with the formation of the underplate, the timing of sill emplacement can be further constrained to between 1.85 and 1.8 Ga (Fig. 7.3b). However, since the lower crustal heating event that reset the age of the xenoliths could have occurred some time after the emplacement of the underplate, the lower bound on emplacement may be slightly earlier than 1.85 Ga. 2.1  1  2.0  1  1.9  1  1.8  1  1.7  2.1  1  2.0  1.9  1.8  1.7  Ga  fia  •T- VULCAN STRUCTURE \ , ?  HSi REFLECTOR ;  a  ? HSI  Figure 7.3: Temporal constraints on the emplacement of the HSI sills a) as suggested by Mandler and Clowes (1998) and b) as suggested if the emplacement of the sills was coeval with the emplacement of a mafic underplate beneath southern Alberta and northern Montana.  Chapter 7. Discussion  181  7.4 Emplacement of the Winagami reflection sequence From cross-cutting relationships with dipping crustal fabric and reflection offsets from dated brittle events, the emplacement of the Winagami reflection sequence is argued to have occurred between 1890 and 1760 Ma (Ross and Eaton, 1997). As with the HSI reflector, the Winagami sequence is inferred to have been emplaced into a Theologically layered upper-middle crust under horizontal compression. From the timing constraints, this compression has been associated with nearby brittle indentation of the Slave province into the Rae province to the northeast (Ross and Eaton, 2002). This collision was accompanied by temporally sporadic but regionally extensive magmatism mostly evidenced as vertical dike intrusions (Ross and Eaton, 1997). One aspect of the Winagami reflection sequence that has not been previously investigated is its unusual anastamosing shape, evidenced by both the PC2 and Lithoprobe results (Fig. 6.2). Considering that the emplacement of a sill will locally alter the stress field such that the stresses perpendicular to the sill propagation direction will increase, nearby sills propagating parallel will be deflected away from each other as they each seek out an orientation that better favours the opening of cracks (Pollard, 1973). Based on this scenario, the Winagami sills must originate in the southeast and diverge as they propagate to the northwest. While physically plausible, such a propagation direction contradicts the indentation of the Slave province into the Rae province as the source of the magmatism. From the high upper crustal RMS velocities obtained to the northeast of the WRS from the refraction/wide-angle reflection (RWAR) Peace River Arch Seismic Experiment (PRASE) (Zelt and Ellis, 1989), Ross and Eaton (2002) proposed that the Winagami sills originated in the northeast. Based on this scenario, the sills would have to converge as they propagate southward. From mechanical modelling results, Pollard (1973) outlined three conditions under which sills could converge: 1. if the host rocks are strongly anisotropic and thus control the path of intrusion  Chapter 7. Discussion  182  2. if there is a time gap between the intrusion of individual sills such that the applied stress can relax between intrusions 3. if the sills encounter regional extension normal to the propagation direction which favours the opening of cracks.  Based on numerical simulations of vertical fluid-filled dikes propagating and converging towards the topographic load of a volcano (Dahm, 2000), a fourth condition for convergence can be envisaged where horizontal sills propagate towards an increased and localized maximum compressive stress (Fig. 7.4). In order for the southward propagation of the Winagami sills to be plausible, at least one of the four proposed conditions must be satisfied.  Figure 7.4: Diagram of dike and sill propagation in response to a) the topographic load from a volcano (based on results from Dahm (2000)) and b) a localized horizontal compressive force.  Since the basement rocks into which the WRS is intruded are generally gneissic, metamorphic segregation and banding of minerals may have imparted an anisotropy into the upper crust which controlled the direction of emplacement. From drillcore evidence from several of the basement domains, both weak and strong foliations are observed which could support this scenario (Villeneuve et al., 1993). Thus, condition #1 is plausible. Since the only way to discount  Chapter 7.  Discussion  183  condition #2 involves drilling to the sills and dating them, in the absence of contradictory evidence, it remains a viable contender to support a southward propagation of the sills. Meanwhile, tectonic candidates for satisfying conditions #3 and/or #4 can be investigated by examining the regional tectonic activity during emplacement.  7.4.1  Regional tectonic activity during emplacement  Following from Ross (2002), the tectonic assembly of northwestern Alberta took place during the Paleoproterozoic in three general regions separated by the Great Slave Lake shear zone in the north and the enigmatic Snowbird Tectonic Zone in the south (Fig. 7.5). In the central block, by 1980 Ma, the Buffalo Head and Chinchaga domains acted as one unit flanked to the east and west by outward dipping subduction zones (Ross, 2002; Ross and Eaton, 2002). The westward subduction zone resulted in the formation of the Ksituan arc along a rifted portion of the Archean Slave province known as the Nova domain. To the east, subduction resulted in the Thelon-Taltson magmatic arc along the western edge of the Archean Rae province. By 1940 Ma, the Buffalo Head-Chinchaga domain was sutured to the Thelon-Taltson belt while subduction continued in the west. To the north of the Great Slave Lake shear zone, the Slave province had collided with the Rae province and a passive margin overlain by the Coronation supergroup had developed along its western margin. Outboard to the west, the Hottah arc was initiated along a westward dipping subduction zone. Tectonic activity in the south was initiated by 1900 Ma in response to the counter-clockwise rotation of the Hearne province due to plate consumption in the southeast as part of the Trans-Hudson Orogen (Ross et al., 2000; Ross and Eaton, 2002). In the central block, subduction had been replaced by continental collision and further north, the Hottah terrane had collided with the western Slave province and the Great Bear magmatic arc had resulted. Within 50 million years, the Fort Simpson terrane collided with the amalgamated Hottah terrane and Slave province to complete the assembly of  Chapter 7. Discussion  184  western Laurentia north of the Great Slave Lake shear zone. Meanwhile, the central block was relatively quiescent apart from shifting of the Taltson arc southward along lateral shear zones. At 1850 Ma, a possible emplacement time for the Winagami sequence, the greatest amount of tectonic activity in the region was occurring south of the Snowbird Tectonic Zone. As the plate consumption in the Trans-Hudson Orogen continued, the Thorsby ocean basin was rifted open between the Wabamun domain and the Hearne province. With the cessation of Trans-Hudson plate consumption to the southeast by 1800 Ma, the Hearne province started moving northwestward again and consequently the Thorsby basin was partially subducted under the Hearne forming the Rimbey granites. While this tectonic framework does not provide an obvious candidate to satisfy condition #3, the opening of the Thorsby basin and the corresponding northwestward ridge push could have satisfied condition #4 by providing a localized and increased compressive stress towards which the sills could have propagated. If correct, this scenario could further constrain the timing of emplacement of the Winagami sequence to between 1890 and 1850 Ma, synchronous with the opening of the rift basin. Since no direct methods are available to confirm that the conditions necessary for sill convergence were satisfied for the emplacement of the Winagami sequence, southeastward sill convergence can't be ruled out. This possibility coupled with the higher upper crustal velocities to the north and the coeval magmatism associated with the indentation of the Slave province into the Rae province are suggestive of a source of magmatism to the northeast of the imaged extent of the intrusions. While a temporal and spatial link can be made between the HSI sills and a mafic underplate in the southwest of Alberta, no equivalent underplate has been recognized in the Winagami region and its magmatic influence appears limited to the upper crust (Ross, 2002). As argued by Ross and Eaton (2002) on the basis of results from Ernst and Buchan (2001), the lack of an associated underplate does not preclude its existence but may indicate that the magmas travelled a great distance from their source.  Chapter 7.  Discussion  185  WRS Emplacement 1890 to 1760 Ma Figure 7.5: Tectonic evolution of northwestern Laurentia from 1980 to 1800 Ma. Thefigurewas inspired and adapted from Ross (2002). The detailed description of thefigureis provided in the text. Tectonic domains and features: BH-C, Buffalo Head-Chinchaga domain; C, Coronation supergroup; FS, Fort Simpson terrane; G, Great Bear magmatic arc; GSLsz, Great Slave Lake shear zone; H, Hottah magmatic arc; K, Ksituan magmatic arc; N, Nova domain; Q, Queen Maud uplift; R, Rimbey granites; S, Slave province; STZ, Snowbird Tectonic Zone; T, Thelon magmatic arc; Ta, Taltson magmatic arc; Th, Thorsby domain; Wa, Wabamun domain.  186  Chapter 7. Discussion  7.5  C o n s t r a i n t s o n t h e m o d eo f e m p l a c e m e n t  Based on studies of exposed sill complexes worldwide, individual sills appear to range in thickness from millimeters to hundreds of meters and are capable of spanning tens of kilometers in length (Philpotts, 1990). For horizontal sills to propagate, the magma pressure must exceed the lithostatic load acting normal to the intrusion (Pollard, 1973). Since the Winagami sills were intruded prior to the deposition of the Western Canada Sedimentary Basin, their present-day depth distribution between 10 and 20 km does not represent the depths at which they were intruded. Compensating for the thickness of the sedimentary cover, Ross and Eaton (2002) suggest that sill emplacement occurred between 6.5 and 16.5 km depth. Assuming a density of 2650 kg/m  3  for the host rocks, these depths correspond to a range of lithostatic pressures from 169 to 429 MegaPascals (MPa). Consequently, magma pressures must have exceeded these values for intrusion to have taken place. Likening sill emplacement to hydraulic fracturing, Philpotts (1990) provides the following equation to relate the width, W, and length, L, of a sill with the physical properties of the host rock such that  — = ™*L where P  ex  (7.1)  is the excess pressure of the magma and pho,t and V represent the density and Pp  wave velocity of the host rock. From this equation, sills are expected to be thinner and extend over longer distances at depth where host rock velocities and densities are generally higher. With estimated thicknesses between 50 and 150 m, the HSI and Winagami sills are excessively thick for their depth if they were emplaced by one pulse of magmatism. According to equation 7.1, to emplace the 70 m thick top sill of the Winagami sequence over a distance of 600 km, with a density of 2650 kg/m and a velocity of 6000 m/s for the host rocks, the ex3  cess magma pressure would have to be 5 MPa (3% of the lithostatic pressure at 6.5 km depth).  Chapter 7. Discussion  187  Conversely, if the imaged sill body was the location of more than one injection of magma, a commonly observed phenomenon that arises when earlier intruded sheets have time to cool and shrink forming fractures that the later intrusions exploit (Philpotts, 1990), much lower excess magma pressures are required. For instance, a 2 m sill could be intruded over the same distance in the same host rocks with an excess magma pressure of only 0.14 MPa (0.1% of the lithostatic pressure at 6.5 km depth). With no way of knowing the excess magma pressure during emplacement of the sills, the lower values are preferred. As such, both the HSI and Winagami sills imaged by Lithoprobe are interpreted to have resulted from multiple injections of magma. Such an emplacement mechanism could explain the lateral variability in their geometries, thicknesses and reflective characters. Carslaw and Jaeger (1959) determined that the ratio of cooled temperature T to emplacement temperature T at a distance x from the center of a sill of arbitrary thickness 2a over a 0  given time t could be described by  (7.2) where  (7.3) and where the thermal diffusivity k of the magma and the host rock are assumed the same. Simplifying this equation to only consider the temperature at the center of the sill (x = 0), the equation becomes  (7.4) Since the boundary conditions for both solutions assume that the host rock has a temperature of 0°C, the emplacement temperature T must be adjusted accordingly. 0  Chapter 7.  Discussion  188  From Equation 7.4, assuming a host rock temperature of 200°C (Hyndman and Lewis, 1999), a magma temperature of 1200°C (adjusted to 1000°C) and a thermal diffusivity of 1 0  -6  m /s 2  (Philpotts, 1990), the time it takes for the center of a sill to cool to one tenth of its adjusted emplacement temperature (100°C) for a sill of thickness 2a is  (7.5) Thus, respectively, a 2 m sill, a 70 m sill and an 150 m sill will each take approximately 370 days, 1200 years and 5600 years to cool under the conditions above. Assuming that the 70 m thick Winagami sill was constructed from a stack of thirty-five individually cooled 2 m sills, the cooling time would correspond to approximately 35 years. All these results suggest that thick, laterally extensive sill complexes like the Winagami and the HSI are formed instantly on geological timescales.  7.6 Role of sills in cratonic arch formation For over 300 million years, from the Middle Ordovician to the Middle Jurassic, the cratonic platform underlying the present-day Western Canada Sedimentary Basin was subdivided into a network of uplifted arches and intervening basins (Porter et al., 1982). These features consequently had a strong influence on sedimentation patterns. The geometries of the arches and the duration of their uplift have since been estimated from the stratigraphic record. The cause of the uplift of the arches, known as epeirogenesis, over such long periods of time remains a topic of active research. Proposed driving forces for epeirogenesis include lithospheric contraction due to cooling (Sleep and Snell, 1976), mantle flow (Pysklywec and Mitrovica, 2000), lithospheric flexure due to surface loading (Price, 1973; Beaumont, 1981) and horizontal intraplate  Chapter 7. Discussion  189  stresses (Cloetingh, 1986; Heller et al., 1993). The uplift of the Peace River Arch has specifically been interpreted as resulting from the extension onto land of oceanic transform faults running perpendicular to the rifted margin (O'Connell, 1994). Based on the coincidental occurrence of the Winagami reflection sequence and the Peace River Arch, Eaton et al. (1999b) proposed a link between the existence of mafic sills within the crystalline crust and the development of cratonic arches. Citing the higher yield strength of mafic rocks versus quartz-bearing lithologies (Ord and Hobbs, 1989), they argued that the sills served to stiffen the crust. Given that the HSI reflector to the south is also intruded into a block which developed into its own cratonic arch, Montania, during the Paleozoic, this potential link deserves further investigation. From finite-element modelling, Heller et al. (1993) argue that weak intraplate stresses can lead to uplift when zones of high and low strength are juxtaposed (Fig. 7.6). By applying a compressional horizontal load on the order of 10 MPa to a 2-D model containing a deep block of strong material surrounded by weak material, they were able to generate surface deflections of a few meters over a width of approximately 100 m, broadly comparable with the dimensions of the Peace River Arch. As such, their results suggest that the formation of cratonic arches may be related to rheological variations in the crust and that local stress instabilities could be sufficient to initiate their uplift. This scenario is appropriate for the formation of Montania and the Peace River Arch on an otherwise quiescent passive margin platform.  7.6.1  Rheological considerations  The strength or rheology of the continental crust depends on a number of factors including the geothermal gradient, the presence of water, the stress rate and the rheology of the crust's individual components. All other factors being equal, lateral variations in rheology within the crust can lead to the juxtaposition of weak rock units and strong rock units.  Chapter 7.  Discussion  190  E  * I  S  _•-  &20-P  Q  cc -6 (/)  0-rn -50  -200  mi ii u m -25 0 25 HORIZONTAL POSITION (km)  -100 0 100 HORIZONTAL POSITION (km)  50  200  -50  -200  -25 0 25 HORIZONTAL POSITION (km)  -100 0 100 HORIZONTAL POSITION (km)  -M-20 50  -? 3  200  Figure 7.6: Surface deflections resulting from the juxtaposition of strong and weak crust as modelled by Heller et al. (1993). The shaded areas within the models represent strong crust at either the top (left) or the bottom (right) of the model. 7 denotes the Young's modulus ratio between the strong and weak blocks. The modelled surface deflections result from the application of a compressive horizontal load on the order of 10 MPa to the edges of the model. The strength of a rock unit, otherwise defined as its ability to support a differential stress (van der Pluijm and Marshak, 1997), depends on the orientation of the applied stress. From uniaxial compressive stress tests which involve applying a vertical load to a rock sample until it fails, certain rock types are observed to support compressive stress better than others. For both the PCI and PC2 study areas, the host rocks are essentially granitic gneisses (Ross et al., 1991) while the inferred intrusions are interpreted as doleritic sills. Experimental values for the uniaxial compressive strengths of different rock types determined at near-surface conditions for engineering applications (Table 7.1) indicate that doleritic material is considerably stronger when it comes to compressive stress than granite or gneiss.  Rock type Granite Gneiss Quartzite Dolerite  Min 153 159 200 227  Max 233 256 304 319  Mean 188.4 195.0 252.0 280.3  Table 7.1: Uniaxial compressive strengths of various rock types in units of MegaPascals (MPa) at near-surface conditions. Experimental values are taken from Bieniawski (1973). The brittle-ductile transition zone corresponds to approximately 15 km in the continental  Chapter 7.  Discussion  191  crust but can vary depending on the geothermal gradient and the amount of heat production within the crust (Sibson, 1982). Above this zone, confining pressure from the overlying rock mass can increase the differential stress that a rock can withstand by suppressing the rock's ability to fracture (van der Pluijm and Marshak, 1997). Meanwhile, below the brittle-ductile transition, higher temperatures generally reduce the strength of the rocks such that they deform by flowing instead of fracturing. Assuming that the Winagami sills were located between 6.5 and 16.5 km depth at the start of the Paleozoic, their presence was concentrated in the brittle portion of the crust. Using an average linear geothermal gradient of 20°C/km (van der Pluijm and Marshak, 1997), the temperatures and confining pressures at the top and bottom of the crust containing sills would have ranged from 130° to 330°C and from 169 to 429 MPa, respectively. If the crust was hotter in the Proterozoic such that a higher geothermal gradient of 30°C/km should be used, temperatures in the crust containing sills would range from 195° to 495°C while confining pressures would remain the same as for the 20°C/km thermal gradient. Based on these estimates, provided that the sills remained strong relative to the host rocks for temperatures up to 500°C and confining pressures on the order of 400 MPa (approximately 15 km depth for a host rock density of 2650 kg/m ), the existence of the sills could have influenced the rheology 3  of the upper crust. From high pressure and high temperature experimental results, gneiss is observed to withstand compressive stresses on the order of 600 MPa at confining pressures of 400 MPa and temperatures of 700°C (Kronenberg et al., 1990). Meanwhile, clinopyroxenite, which has a similar composition to dolerite, can withstand compressive stresses of up to 1000 MPa for approximately the same conditions (Kirby and Kronenberg, 1984). Based on these results, the crust containing sills would have remained strong at depths down to 16.5 km even under hotter conditions than the present. Furthermore, the experimental results indicate that the strength difference between gneissic host rocks and dolerite sills becomes more pronounced with depth. This trend is illustrated by the extrapolated linear compressive strength curves in Figure 7.7.  Chapter 7.  Discussion  192  100  Figure 7.7:  Uniaxial compressive strength (MPa) 300  500  700  900  Inferred uniaxial compressive strength of dolerite and gneiss at depth. T h e curves were generated by  connecting the near-surface mean uniaxial compressive strengths (Table 7.1) with the high temperature and high pressure experimental results from K i r b y and Kronenberg (1984) and Gottschalk et al. (1990). T h e inferred depth range o f emplacement for the Winagami intrusions is shown as the shaded gray box. T h e extrapolated projection of gneissic strength past 15 k m depth (the approximate confining pressure o f 400 M P a ) is shown as the dashed black line.  Given the contrast in compressive strength between dolerite and gneiss and assuming a buckling fold model (van der Pluijm and Marshak, 1997), when a gneissic block intruded with horizontal sills is subjected to a horizontal compressive stress, the stronger doleritic layers will produce folds while the weaker host rocks will deform themselves to accommodate the folds (Fig. 7.8). Since the wavelength of the folding of a layer depends in part on the contrast in strength between the layer and the host rock, the trends in Figure 7.7 indicate that the deepest sills in the gneissic block will be the strongest and will fold the least. Since the wavelength of the folding of a layer also depends on the thickness of the layer (with thicker layers folding less), when multiple layers of various thickness are present, the wavelength of the folding exhibited by the thinner layers will increase due to the influence of thicker layers in close proximity. As such, these layers will reinforce each other such that a gneissic block intruded with many horizontal sills of varying thickness will be stronger than a  Chapter 7. Discussion  193  neighbouring block without sills (Fig. 7.8d).  a)  before  after  d)  b)  c)  before  p  l  r  MM  diorite sill  — •  horizontal  h M-  after  l granitic host rock  compression  • host rock deformation  Figure 7.8: Schematic diagram of the folding of sills of different thicknesses in the presence of a uniform horizontal compressive force; a) thick sill, b) sill of intermediate thickness, c) thin sill. Plot d) shows the effect of having multiple sills of varying thickness in the presence of a horizontal compressive force. Diagram inspired by van der Pluijm and Marshak (1997).  Based on the general spatial correspondence between the northern extent of the WRS and a localized increase in upper crustal RMS velocities from 6.45 to 6.7 km/s as determined from the PRASE seismic refraction experiment (Zelt and Ellis, 1989), Ross and Eaton (2002) suggest that if the higher velocities are due solely to the presence of high velocity mafic sills then sills could constitute up to 50% of the upper crust at the northern extent of the Winagami complex. Given that the imaged sills make up less than 5% of the upper crust, many thin sills, too thin to be resolved seismically, would be required to generate the observed higher velocities. If such a distribution of thick and thin sills exists, the basement blocks in that region could potentially be much stronger than the surrounding sill-free blocks. Relating these speculations to Heller's results, given that the range in Young's moduli for dolerite and gneiss span from 8.19-9.94 and  Chapter 7. Discussion  194  from 1.04-1.91 N/m respectively (Carmichael, 1989), a composite block containing both ma2  terials will likely have a higher effective Young's modulus than gneiss alone. Consequently, the inferred rheological contrast provided by sills could be sufficient to cause uplift in the presence of weak intraplate stresses.  7.6.2  Uplift of cratonic arches  In reality, while intruded sill material potentially has the ability to strengthen crustal blocks, more conditions must be met before any uplift will take place. Firstly, since even sill-free gneiss can withstand compressive stresses on the order of 600 MPa to significant depth (Kronenberg et al., 1990), in the presence of weak intraplate stress on the order of 10 MPa (Heller et al., 1993), rheological variations imparted by sills will be incapable of fracturing the gneissic basement and allowing part of it to be uplifted. Consequently, the only way rheological variations due to the presence of sills can result in the uplift of crustal arches is if steep fractures or detachment surfaces already exist in the crust along which uplift can be accommodated. For the case of the Peace River Arch (PRA), a number of sutured basement domains were involved in its uplift, namely the Buffalo Head, the Chinchaga and the Ksituan. Given the similarities in shape between the Ksituan arc and the PRA (Fig. 7.9), perhaps the suture between the Ksituan arc and the Chinchaga domain provided the detachment surface along which the uplift occurred. While this scenario is contradicted by the results from line 12 of the PRAISE experiment which image the suture as a steep east-dipping zone of reflectivity that is crosscut by the Winagami sills without offset (Fig. 1.6), perhaps the uplift was focused along the southern suture boundary where the geometries of the Ksituan and the PRA are most similar. In southern Alberta, the Paleozoic Montania Arch corresponded roughly with the present day extent of the Medicine Hat Block (Lemieux et al., 2000). Thus, the Vulcan Structure, the  Chapter 7. Discussion  195  -120'  -118  -116°  -  -114" N  s •  /  /•  s  56'  \  \ \ \  Winagami * sequence  -120'  -118"  -116  -  -114"  Figure 7.9: Suggested spatial connection between the Peace River Arch (PRA, rough outline plotted as the white dashed line) and the Ksituan arc. The Cordilleran deformation front (DF) and the outline of the Winagami reflection sequence are shown for reference.  northern boundary of Montania (Price and Sears, 2000), may have acted as the detachment surface for the uplift. While it is possible that the proposed basement sutures acted as the detachment surfaces for the uplift of the Peace River Arch and Montania, these ideas cannot be directly tested and so remain pure speculation. Nonetheless, these musings point to the possibility that sills not only contribute to the thickening and stabilization of cratons but also play a role in their long term tectonic development.  Chapter 7. Discussion  7.7  196  Chapter summary  In this chapter, the results from the PCI and PC2 datasets along with the earlier Lithoprobe results were used to address the emplacement of the HSI and Winagami inferred sill complexes in terms of mechanics, tectonics and timing. Furthermore, the possible influence of these sill complexes on the formation of cratonic arches through their strengthening of the crust was also addressed.  Chapter 8  Conclusions  The results from this thesis represent the first opportunity in North America to examine uppermiddle crustal structure to depths of approximately 20 km using industry-style 3-D seismic reflection techniques. Using datasets collected by the Canadian petroleum industry in southwestern and northwestern Alberta, Precambrian sill complexes within the crystalline upper crust, which were first identified along Lithoprobe 2-D multichannel seismic reflection profiles, are investigated in 3-D. In so doing, new developments and applications of data processing procedures are implemented. Based on information derived from the seismic results, inferences about the emplacement mechanisms of the sills and associated tectonic development at the time of emplacement are developed. During seismic processing of the dataset in southwestern Alberta (PCI), a recently developed wavelet-based method, Physical Wavelet Frame Denoising (Zhang and Ulrych, 2003), is adapted and applied to suppress significant ground roll contamination while preserving low frequency signals from deeper structures. Through this approach, a severe problem of identification of coherent reflectivity at depths where even conservative high-passfilterseither left behind significant ground roll energy or diminished both ground roll and reflection energy is resolved, enabling study of deep sills in the area. To address the common problem of spatial aliasing associated with 3-D data acquisition, a new 3-D empirical trace interpolation scheme, DSInt, is developed. It produces results comparable to available interpolation codes while providing greater flexibility in handling irregular  197  Chapter 8.  Conclusions  198  acquisition geometries and interpolated trace headers. The method is tested on both 3-D datasets and successfully recovered both missing inline receivers and entire missing receiver lines. Further evidence of the Head-Smashed-In reflector, the inferred sill complex in southwestern Alberta (Mandler and Clowes, 1998), is obtained using a dataset acquired specifically to 8 s T W T T (approx. 24 km depth) by industry to enable this study. Discontinuous basement reflectivity is detected, although it differs in amplitude characteristics and coherence from the main HSI body. Images of the sills are comparable to the near-most Lithoprobe results, which appear to image the tapering western edge of the deep reflections. To ascertain the 3-D distribution of coherent reflectivity within the data volume, a new statistical approach of tracking reflectivity is developed and applied. Using this method, laterally consistent pockets of reflectivity are revealed at different times within the time range that the HSI arrivals are generally observed. Simple 1-D forward modelling results reveal that the brightest reflections could be arising from a 50 to 150 m thick body, consistent with earlier estimates of the reflector's thickness (Mandler and Clowes, 1998). Since modelling tests, which included multiples, reveal that deep reflections arising from bodies of less than 50 m in thickness could not be easily detected, sills with laterally varying thicknesses could have produced the deep reflections evidenced in the PCI data volume. As such, these results are consistent with imaging the tapering edge of the sill complex. Clear evidence of the Winagami reflection sequence (Ross and Eaton, 1997), the inferred sill complex in northwestern Alberta, emerges from the second dataset (PC2) acquired to 5.1 s T W T T (approx. 15 km depth). Deep reflections recorded along inline, crossline and time sections of the data reveal a 3-D reflective sheet dipping to the southeast. Polarity comparisons of known reflections from within the sedimentary sequence are used to approximate the sign of the impedance contrast from which the deep reflections are arising. These constraints indicate that a reflective body of higher density and higher velocity than the surrounding gneissic basement  Chapter 8. Conclusions  199  is responsible for the deep reflections. Such a finding is consistent with the earlier interpretation that the Winagami reflection sequence represents doleritic sills (Ross and Eaton, 1997). Simple 1-D forward modelling results reveal that the thickness of the sheet is between 50 and 100 m throughout the survey region. From 3-D Kirchhoff forward modelling, the geometry of the reflective sheet is constructed and found to exist between approximately 11 and 16 km depth beneath the survey area. Deeper evidence of the reflective sheet is unavailable due to the limited recording time. Comparison of the distribution of reflectivity within the PC2 data volume and the regional aeromagnetic pattern reveal that the absence of Winagami reflections, and basement reflectivity in general, coincides with a positive aeromagnetic anomaly, previously interpreted as corresponding to magmatic rocks. I suggest that the two features may be linked. The emplacement of the magmatic rocks could have preceded the emplacement of the Winagami sills and somehow influenced their intrusional geometry. Alternatively, emplacement could have postdated the Winagami sills and so overprinted their reflective signature. Because the latter scenario would be more likely to obliterate all reflectivity, it is preferred. However, without dating of the relevant basement drillcore, the question of timing remains unanswered. The tectonic environments related to the emplacement of both sill complexes are investigated to help further understand the temporal, rheological and thermal constraints on their emplacement. In particular, the anastamosing shape of the Winagami reflections is used to argue the case for two possible directions of emplacement. Since neither direction could be refuted based on available evidence, the magmatism source in the northeast and the southward propagation direction suggested by Ross and Eaton (2002) is preferred. Lastly, based on analogies with fold mechanics, a link between the existence of deep sill complexes and the eventual development of Paleozoic arches along quiescent cratonic platforms is suggested. I propose that the sills strengthen the crust and that reactivated basement sutures near the sill margins act as the loci of basement uplift.  Chapter 8. Conclusions  200  Generally, this thesis demonstrates that 3-D seismic reflection datasets collected by the petroleum industry can contain meaningful information for tectonic investigations of the upper crystalline crust if they are recorded to adequate times. Given the exceptionally high costs of acquiring such datasets solely for academic purposes, industry donations of data likely provide the only means, for now, of undertaking such studies. However, due to the significant computational and time demands of working with large 3-D datasets, the deep targets of interest should be well defined by other methods and also should be readily apparent on representative data gathers to ensure meaningful results. Communication between industry and academia is essential to isolating future deep targets of interest and to ensuring that data are acquired to sufficiently long listening times so as to allow academic studies of the crust below the sedimentary section.  References  Allaby, A., and Allaby, M . , 1991, Concise dictionary of earth sciences: Oxford University Press, first edition. Bardan, V., 1987, Trace interpolation in seismic data processing: Geophysical Prospecting, 35, 343-358. Barton, P. J., 1986, The relationship between seismic velocity and density in the continental crust - a useful constraint?: Geophysical Journal of the Royal Astronomical Society, 87,195— 208. Beaumont, C , 1981, Foreland basins: Geophysical Journal of the Royal Astronomical Society, 65,291-329. Begin, N., Lawton, D., and Spratt, D., 1996, Seismic interpretation of the Rocky Mountain thrust front near the Crowsnest deflection, southern Alberta: Bulletin of Canadian Petroleum Geology, 44, 1-13. Bieniawski, Z. T , 1973, Engineering classification of jointed rock masses: Transactions of the South African Institution of Civil Engineers, 15, 335-344. Blakely, R. J., 1996, Potential theory in gravity and magnetic applications: Cambridge University Press. Bradley, J., 1965, Intrusion of major dolerite sills: Transactions of the Royal Society of New Zealand, 3, 27-55.  201  References  202  Brown, L . D., Zhao, W., Nelson, K. D., Hauck, M . , Alsdork, D., Ross, A., Cogan, M . , Clark, M . , Liu, X., and Che, J., 1996, Bright spots, structure and magmatism in southern Tibet from INDEPTH seismic reflection profiling: Science, 274, 1688-1690. Burwash, R. A., Baadsgaard, H., and Peterman, Z. E., 1962, Precambrian K - A r dates from the Western Canada Sedimentary Basin: Journal of Geophysical Research, 67, 1617-1625. Carmichael, R. S., 1989, Practical handbook of physical properties of rocks and minerals: CRC Press, Inc. Carslaw, H. S., and Jaeger, J. C , 1959, Conduction of heat in solids: Oxford University Press, Inc., second edition. Claerbout, J. E , 1985, Imaging the Earth's interior: Blackwell Scientific Publications, Boston, first edition. Claerbout, J. F., 1992, Earth soundings analysis: processing versus inversion: Blackwell Scientific Publications, Boston. Cloetingh, S., 1986, Intraplate stresses: a new tectonic mechanism for fluctuations of relative sea level: Geology, 14, 617-620. Clowes, R. M . , Burianyk, M . J., Gorman, A. R., and Kanasewich, E. R., 2002, Crustal velocity structure from SAREX, the Southern Alberta Refraction Experiment: Canadian Journal of Earth Sciences, 39, 351-373. Coffin, M . F., and Eldholm, O., 1994, Large igneous provinces: crustal structure, dimensions, and external consequences: Reviews of Geophysics, 32, 1-36. Cook, F. A., and Coflin, K. C , 1990, Experimental three-dimensional imaging of crustal structure in the northwestern Canadian Arctic: Tectonophysics, 173, 43-52.  References  203  Corso, G., Kuhn, P. S., Lucena, L . S., and Thome, Z. D., 2003, Seismic ground roll timefrequency filtering using the Gaussian wavelet transform: Physica A: Statistical mechanics and its applications, 318, 551-561. Dahm, T., 2000, Numerical simulations of the propagation path and the arrest of fluid-filled fractures in the Earth: Geophysical Journal International, 141, 623-638. Dasgupta, R., and Clark, R. A., 1998, Estimation of Q from surface seismic reflection data: Geophysics, 63, 2120-2128. Davis, W. J., Berman, R., and Kjarsgaard, B., 1995, U-Pb geochronology and isotopic studies of crustal xenoliths from the Archean Medicine Hat Block, northern Montana and southern Alberta: Paleoproterozoic reworking of Archean crust in Ross, G., Ed., 1995 Alberta Basement Transect Workshop:: Lithoprobe Secretariat, The University of British Columbia, Lithoprobe Report 47, 329-334. Deighan, A. J., and Watts, D. R., 1997, Ground-roll suppression using the wavelet transform: Geophysics, 62, 1896-1903. Deiss, C , 1941, Cambrian geography and sedimentation in the central Cordilleran region: Bulletin of the Geological Society of America, 52, 1085-1116. Dohr, G., and Fuchs, K., 1967, Statistical evaluation of deep crustal reflections in Germany: Geophysics, 32, 951-967. Eaton, D. W , and Clarke, G., 2000, A Kirchhoff integral method to model 3-D elastic-wave scattering: 70th Annual International Meeting of the Society of Exploration Geophysicists Expanded Abstracts. Eaton, D. W , Ross, G. M . , and Clowes, R. M . , 1999a, Seismic-reflection and potential-field  References  204  studies of the Vulcan structure, western Canada: a Paleoproterozoic Pyrenees?: Journal of Geophysical Research, 104, 255-269. Eaton, D. W., Ross, G. M . , and Hope, J. 1999b, The rise and fall of a cratonic arch: A regional seismic perspective on the Peace River Arch, Alberta: Bulletin of Canadian Petroleum Geology, 4 7 , 346-361. Ernst, R. E., and Buchan, K. L., 2001, The use of mafic dyke swarms in identifying and locating mantle plumes in Ernst, R. E., and Buchan, K. L., Eds., Mantle plumes: their identification through time:: Geological Society of America, Special Paper 352, 247-265. Fenton, M . M . , Schreiner, B. T , Nielsen, E., and Pawlowicz, J. G., 1994, Quaternary geology of the Western Plains in Mossop, G. D., and Shetsen, I., Eds., Geological Atlas of the Western Canada Sedimentary Basin:: Canadian Society of Petroleum Geologists and Alberta Research Council, 413-420. Foster, D. J., and Mosher, C. C , 1992, Suppression of multiple reflections using the Radon transform: Geophysics, 5 7 , 386-395. Fuchs, K., and Miiller, G., 1971, Computation of synthetic seismograms with the reflectivity method and comparison with observations: Geophysical Journal of the Royal Astronomical Society, 23, 417^133. Gibson, B., Larner, K., and Levin, S., 1983, Efficient 3-D migration in two steps: Geophysical Prospecting, 33, 1-33. Gorman, A. R., Clowes, R. M . , Ellis, R. M . , Burianyk, M . J. A., Kanasewich, E. R., Asudeh, I., Hajnal, Z., Spence, G. D., Henstock, T. J., Levander, A. R., Keller, G. R., and Miller, K. C , 2002, Deep Probe - Imaging the roots of western North America: Canadian Journal of Earth Sciences, 39, 375-398.  References  205  Gretener, P. E., 1969, On the mechanics of the intrusion of sills: Canadian Journal of Earth Sciences, 6, 1415-1419. Gulunay, N., 2003, Seismic trace interpolation in the Fourier transform domain: Geophysics, 68,355-369. Hampson, D., and Russell, B., 1984, First-break interpretation using generalized linear inversion: Journal of the Canadian Society of Exploration Geophysicists, 20, 40-54. Heller, P. L., Beekman, E , Angevine, C. L., and Cloetingh, S. A. P. L., 1993, Cause of tectonic reactivation and subtle uplifts in the Rocky Mountain region and its effect on the stratigraphic record: Geology, 21, 1003-1006. Herrmann, R. B., and Russell, D. R., 1990, Ground roll: Rejection using adaptive phasematched filters: Geophysics, 55, 776-781. Hiebert, S. N., and Spratt, D. A., 1996, Geometry of the thrust front near Pincher Creek, Alberta: Bulletin of Canadian Petroleum Geology, 4 4 , 195-201. Hindriks, K., and Duijndam, A. J. W., 2000, Reconstruction of 3-D seismic signals irregularly sampled along two spatial coordinates: Geophysics, 65, 253-263. Hoffman, P. E , 1988, United plates of America, the birth of a craton: Early Proterozoic assembly and growth of Proto-Laurentia: Annual Reviews of Earth and Planetary Sciences, 16, 543-603. Hope, J., and Eaton, D. W., 2002, Crustal structure beneath the Western Canada Sedimentary Basin: constraints from gravity and magnetic modelling: Canadian Journal of Earth Sciences, 39,291-312.  References  206  Hyndman, R. D., and Lewis, T. J., 1999, Geophysical consequences of the Cordillera-Craton thermal transition in southwestern Canada: Tectonophysics, 306, 397-422. Hyndman, R. D., 1988, Dipping seismic reflectors, electrically conductive zones, and trapped water in the crust over a subducting plate: Journal of Geophysical Research, 93, 13391— 13405. Jakubowicz, H., and Levin, S., 1983, A simple exact method of 3-D migration - theory: Geophysical Prospecting, 31, 34-56. Jarchow, C. M . , Thompson, G. A., Catchings, R. D., and Mooney, W. D., 1993, Seismic evidence for active magmatic underplating beneath the Basin and Range province, western United States: Journal of Geophysical Research, 98, 22095-22108. Jones, P. B., 1996, Triangle zone geometry, terminology and kinematics: Bulletin of Canadian Petroleum Geology, 44, 139-152. Juhlin, C , 1990, Interpretation of the reflections in the Siljan Ring area based on results from the Gravberg-1 borehole: Tectonophysics, 173, 345-360. Kanasewich, E. R., Hajnal, Z., Green, A., Cumming, G. L., Mereu, R. F., Clowes, R. M . , MorelA-L'Huissier, P., Chiu, S., Macrides, C. G., Shahriar, M . , and Congram, A. M . , 1987, Seismic studies of the crust under the Williston Basin: Canadian Journal of Earth Sciences, 24,21602171. Kanasewich, E . R., Burianyk, M . J. A., Dubuc, G. P., Lemieux, J. E , and Kalantzis, E , 1995, Three-dimensional seismic reflection studies of the Alberta basement: Canadian Journal of Exploration Geophysics, 31,1-10. Kent, D. M . , 1994, Paleogeographic evolution of the cratonic platform: Cambrian to Triassic in  References  207  Mossop, G. D., and Shetsen, I., Eds., Geological Atlas of the Western Canada Sedimentary Basin:: Canadian Society of Petroleum Geologists and Alberta Research Council, 69-86. Kirby, S. H . , and Kronenberg, A. K., 1984, Deformation of clinopyroxenite: evidence for a transition in flow mechanisms and semibrittle behavior: Journal of Geophysical Research, 89,3177-3192. Kjartannson, E., 1979, Constant Q - wave propagation and attenuation: Journal of Geophysical Research, 84, 4737-4748. Kong, S. M . , Phinney, R. A., and Roy-Chowdhury, K., 1985, A nonlinear signal detector for enhancement of noisy seismic record sections: Geophysics, 50, 539-550. Korsch, R. J., Wake-Dyster, K. D., and Johnstone, D. W , 1992, Seismic imaging of Late Paleozoic-Early Mesozoic extensional and contractional structures in the Bowen and Surat basins, eastern Australia: Tectonophysics, 215, 273-294. Kronenberg, A. K., Russell, J. E., andHandin, J., 1990, Mechanical anisotropy of gneiss: failure criterion and textural sources of directional behavior: Journal of Geophysical Research, 95, 21613-21634. Larner, K., Gibson, B., and Rothman, D., 1981, Trace interpolation and the design of seismic survey: Geophysics, 46, 407. Law, A., and Snyder, D. B., 1997, Reflections from a mylonitized zone in central Sweden: Journal of Geophysical Research, 102, 8411-8425. Lay, T , and Wallace, T. C , 1995, Modern global seismology: Academic Press, Inc., first edition.  References  208  Lemieux, S., Ross, G. M . , and Cook, F. A., 2000, Crustal geometry and tectonic evolution of the Archean crystalline basement beneath the southern Alberta plains, from new seismic reflection and potential-field studies: Canadian Journal of Earth Sciences, 37, 1473-1491. Lemieux, S., 1999, Seismic reflection expression and tectonic significance of Late Cretaceous extensional faulting of the Western Canada Sedimentary Basin in southern Alberta: Bulletin of Canadian Petroleum Geology, 47, 375-390. Lines, L., and Treitel, S., 1984, A review of least-squares inversion and its application to geophysical problems: Geophysical Prospecting, 32, 159-186. Linville, A. F., and Meek, R. A., 1995, A procedure for optimally removing localized coherent noise: Geophysics, 60, 191-203. Litak, R. K., and Hauser, E. C , 1992, The Bagdad reflection sequence as tabular mafic intrusions: Evidence from seismic modeling of mapped exposures: Bulletin of the Geological Society of America, 104,1315-1325. Liu, X., 1999, Ground roll suppression using the Karhunen-Loeve transform: Geophysics, 64, 564-566. Mallat, S., 1999, A wavelet tour of signal processing: Academic Press, second edition. Mandler, H . A. E , and Clowes, R. M . , 1997, Evidence for extensive tabular intrusions in the Precambrian shield of western Canada: A 160-km-long sequence of bright reflections: Geology, 25, 271-274. Mandler, H . A. F., and Clowes, R. M . , 1998, The HSI bright reflector: further evidence for extensive magmatism in the Precambrian of western Canada: Tectonophysics, 288, 71-81.  References  209  Martinson, D. G., and Hopper, J. R., 1992, Nonlinear seismic trace interpolation: Geophysics, 57,136-145. McMechan, G. A., and Sun, R., 1991, Depthfilteringof first breaks and ground roll: Geophysics, 56, 390-396. Milkereit, B., Bittner, R., and Meissner, R., 1986, Off-line acquisition of crustal reflection and refraction data: Geophysical Research Letters, 13,1161-1164. Mueller, R A., Heatherington, A. L., Kelly, D. M . , Wooden, J. L., and Mogk, D. W , 2002, Paleoproterozoic crust within the Great Falls tectonic zone: Implications for the assembly of southern Laurentia: Geology, 30, 127-130. Nelson, K. D., 1991, A unified view of craton evolution motivated by recent deep seismic reflection and refraction results: Geophysical Journal International, 105, 25-35. Newman, P., 1973, Divergence effects in a layered Earth: Geophysics, 38,481^188. Newson, A. C., 2001, The future of natural gas exploration in the foothills of the western Canadian Rocky Mountains: The Leading Edge, 20,74-79. O'Connell, S. C , 1994, Geological history of the Peace River Arch in Mossop, G. D., and Shetsen, I., Eds., Geological Atlas of the Western Canada Sedimentary Basin:: Canadian Society of Petroleum Geologists and Alberta Research Council, 431-438. Ord, A., and Hobbs, B. E., 1989, The strength of the continental crust, detachment zones and the development of plastic instabilities: Tectonophysics, 158, 269-289. Papasikas, N., and Juhlin, C , 1997, Interpretation of reflections from the central part of the Siljan Ring impact structure based on results from the Stenberg-1 borehole: Tectonophysics, 269, 237-245.  References  210  Parsons, T., Sleep, N. H., and Thompson, G. A., 1992, Host rock rheology controls on the emplacement of tabular intrusions: Implications for underplating of extended crust: Tectonics, 11, 1348-1356. Philpotts, A. R., 1990, Principles of igneous and metamorphic petrology: Prentice Hall. Pilkington, M . , Miles, W. R, Ross, G. M . , and Roest, W. R., 2000, Potential-field signatures of buried Precambrian basement in the Western Canada Sedimentary Basin: Canadian Journal of Earth Sciences, 37, 1453-1471. Pollard, D. D., Muller, O. H., and Dockstader, D. R., 1975, The form and growth of fingered sheet intrusions: Bulletin of the Geological Society of America, 86, 351-363. Pollard, D. D., 1973, Derivation and evaluation of a mechanical model for sheet intrusions: Tectonophysics, 19, 233-269. Porsani, M . J., 1999, Seismic trace interpolation using half-step prediction filters: Geophysics, 64, 1461-1467. Porter, J. W , Price, R. A., and McCrossan, R. G., 1982, The Western Canada Sedimentary Basin: Philosophical Transactions for the Royal Society of London. Series A , Mathematical and Physical Sciences, 305, 169-192. Pratt, T. L., Mondary, J. R, Brown, L. D., Christensen, N. I., and Danbom, S. H., 1993, Crustal structure and deep reflector properties: wide angle shear and compressional wave studies of the midcrustal Surrency bright spot beneath southeastern Georgia: Journal of Geophysical Research, 98, 17723-17735. Price, R. A., and Sears, J. W , 2000, A preliminary palinspastic map of the Mesoproterozoic Belt-Purcell Supergroup, Canada and USA: implications for the tectonic setting and structural evolution of the Purcell Anticlinorium and the Sullivan Deposit in Lydon, J. W , Hoy,  References  211  T., Slack, J. R, and Knapp, M . E., Eds., The geological environment of the Sullivan Deposit, British Columbia:: Geological Association of Canada, Mineral Deposits Division, Special Publication, 61-81. Price, R. A., 1973, Large-scale gravitational flow of supra-crustal rocks, southern Canadian Rockies in Dejong, K. A., and Scholten, R. A., Eds., Gravity and Tectonics:: Wiley, 491502. Pysklywec, R. N., and Mitrovica, J. X., 2000, Mantle flow mechanisms of epeirogeny and their possible role in the evolution of the Western Canada Sedimentary Basin: Canadian Journal of Earth Sciences, 37, 1535-1548. Rabbel, W , and Mooney, W. D., 1996, Seismic anisotropy of the crystalline crust: what does it tell us?: Terra Nova, 8, 16-21. Robinson, E . A., and Treitel, S., 1980, Geophysical signal analysis: Prentice-Hall Inc., first edition. Ronen, J., and Claerbout, J. E , 1985, Surface-consistent residual statics estimation by stackpower maximization: Geophysics, 50, 2759-2767. Ronen, J., 1987, Wave-equation trace interpolation: Geophysics, 52, 973-984. Ross, G. M . , and Eaton, D. W., 1997, Winagami reflection sequence: Seismic evidence forpostcollisional magmatism in the Proterozoic of western Canada: Geology, 25, 199-202. Ross, G. M . , and Eaton, D. W , 2002, Proterozoic tectonic accretion and growth of western Laurentia: results from Lithoprobe studies in northern Alberta: Canadian Journal of Earth Sciences, 39, 313-329.  References  212  Ross, G . M . , Parrish, R., Villeneuve, M . , and Bowring, S., 1991, Geophysics and geochronology of the crystalline basement of the Alberta Basin, western Canada: Canadian Journal of Earth Sciences, 28, 512-522. Ross, G. M . , Eaton, D. W., Boerner, D. E., and Miles, W., 2000, Tectonic entrapment and its role in the evolution of the continental lithosphere: An example from the Precambrian of western Canada: Tectonics, 19, 116-134. Ross, G. M . , 2002, Evolution of Precambrian continental lithosphere in western Canada: results from Lithoprobe studies in Alberta and beyond: Canadian Journal of Earth Sciences, 39, 413^137. Saatcilar, R., and Canitez, N., 1988, A method of ground-roll elimination: Geophysics, 53,894902. Saatcilar, R., and Canitez, N., 1994, The lattice filter in ground-roll suppression: Geophysics, 59,623-631. Sacchi, M . D., and Ulrych, T. J., 1995, High-resolution velocity gathers and offset space reconstruction: Geophysics, 60,1169-1177. Sheriff, R. E., and Geldart, L. P., 1999, Exploration seismology: Cambridge Univ. Press. Sibson, R. H., 1982, Fault zone models, heat flow, and depth distribution of earthquakes in the continental crust of the United States: Bulletin of the Seismological Society of America, 72, 151-163. Sleep, N. H., and Snell, N. S., 1976, Thermal contraction andflexureof mid-continent and Atlantic marginal basins: Geophysical Journal of the Royal Astronomical Society, 45, 125154.  References  213  Spencer, T. W., Sonnad, J. R., and Butler, T. M . , January 1982, Seismic Q - Stratigraphy or dissipation: Geophysics, 47, 16-24. Spitz, S., 1991, Seismic trace interpolation in the f-x domain: Geophysics, 56, 785-794. Spratt, D. A., and Lawton, D. C., 1996, Variations in detachment levels, ramp angles and wedge geometries along the Alberta thrust front: Bulletin of Canadian Petroleum Geology, 44,313— 323. Stiller, M . , 1991, 3-D vertical incidence seismic reflection survey at the K T B location, Oberpfalz in Meissner, R., Brown, L . , Durbaum, H.-J., Franke, W., Fuchs, K., and Seifert, F., Eds., Continental lithosphere: deep seismic reflections:: American Geophysical Union Geodynamics Series 22, 101-114. Stolt, R. H., and Benson, A. K., 1986, Seismic migration - theory and practice: Geophysical Press, London-Amsterdam. Theriault, R. J., and Ross, G. M . , 1991, Nd isotope evidence for crustal recycling in the ca. 2.0 ga subsurface of western Canada: Canadian Journal of Earth Sciences, 28, 1140-1147. Trad, D. O., and Travassos, J. M . , 2000, Wavelet filtering of magnetotelluric data: Geophysics, 65,482-491. Trad, D. O., Sacchi, M . D., and Ulrych, T. J., 2001, A hybrid linear-hyperbolic Radon transform: Journal of Seismic Exploration, 9, 303-318. Trad, D. O., Ulrych, T. J., and Sacchi, M . D., 2002, Accurate interpolation with high-resolution time-variant Radon transforms: Geophysics, 67, 644-656. van der Pluijm, B. A., and Marshak, S., 1997, Earth structure: an introduction to structural geology and tectonics: WCB/McGraw-Hill, first edition.  References  214  Villeneuve, M . E . , Ross, G. M . , Theriault, R. J., Miles, W., Parrish, R. R., and Broome, J., 1993, Tectonic subdivision and U-Pb geochronology of the crystalline basement of the A l berta Basin, western Canada: Geological Survey of Canada Bulletin, 447, 86. Wang, Y., 2002, Seismic trace interpolation in the f-x-y domain: Geophysics, 67,1232-1239. Wang, Y , 2003, Sparseness-constrained least-squares inversion: Application to seismic wave reconstruction: Geophysics, 68, 1633-1638. Warner, M . , 1987, Migration - why doesn't it work for deep continental data?: Geophysical Journal of the Royal Astronomical Society, 89, 21-25. Widess, M . B., 1973, How thin is a thin bed?: Geophysics, 38, 1176-1180. Wright, G. N., McMechan, M . E., and Potter, D. E . G., 1994, Structure and architecture of the Western Canada Sedimentary Basin in Mossop, G. D., and Shetsen, I., Eds., Geological Atlas of the Western Canada Sedimentary Basin:: Canadian Society of Petroleum Geologists and Alberta Research Council, 25-40. Yilmaz, O., 1987, Seismic data processing: Society of Exploration Geophysicists, Tulsa. Yilmaz, O., 2001, Seismic data analysis: Society of Exploration Geophysicists, Tulsa. Yu, Z. Y , McMechan, G. A., Ferguson, J. E , and Anno, P. D., 2002, Adaptive wavelet filtering of seismic data in the wavelet transform domain: Journal of Seismic Exploration, 11, 223246. Zelt, C. A., and Ellis, R. M . , 1989, Seismic structure of the crust and upper mantle in the Peace River Arch region, Canada: Journal of Geophysical Research, 94, 5729-5744. Zhang, R., and Trad, D. O., 2002, Noise attenuation: a hybrid approach based on wavelet transform: Annual Meeting of the Canadian Society of Exploration Geophysicists.  References  215  Zhang, R., and Ulrych, T. J., 2003, Physical wavelet frame denoising: Geophysics, 68, 225231. Zoback, M . L., Zoback, M . D., Adams, J., Assumpcao, M . , Bell, S., Bergman, E. A., Blumling, P., Brereton, N. R., Denham, D., Ding, J., Fuchs, K., Gay, N., Gregersen, S., Gupta, H . K., Gvishiani, A., Jacob, K., Klein, R., Knoll, P., Magee, M . , Mercier, J. L., Muller, B. C., Paquin, C., Rajendran, K., Stephansson, O., Suarez, G., Suter, M . , Udias, A., Xu, Z. PL, and Zhizhin, M . , 1989, Global patterns of tectonic stress: Nature, 341, 291-298.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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.831.1-0052372/manifest

Comment

Related Items