UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Seismic structure of the crust and upper mantle in the Peace River Arch region Zelt, Colin Andrew 1989

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

Item Metadata

Download

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

Full Text

SEISMIC S T R U C T U R E OF T H E CRUST A N D U P P E R M A N T L E IN THE PEACE RIVER A R C H REGION  By Colin Andrew Zelt B. Sc. Hons., The University of Victoria, 1984  A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E REQUIREMENTS  FOR T H E DEGREE OF  D O C T O R OF P H I L O S O P H Y  in T H E F A C U L T Y OF G R A D U A T E STUDIES (DEPARTMENT OF GEOPHYSICS A N D ASTRONOMY)  We accept this thesis as conforming to the required standard  T H E UNIVERSITY OF BRITISH C O L U M B I A  February 1989 © Colin Andrew Zelt, 1989  In  presenting this  degree at the  thesis  in  University of  partial  fulfilment  of  of  department  this or  thesis for by  his  or  requirements  British Columbia, I agree that the  freely available for reference and study. I further copying  the  representatives.  an advanced  Library shall make it  agree that permission for extensive  scholarly purposes may be her  for  It  is  granted  by the  understood  that  head of copying  my or  publication of this thesis for financial gain shall not be allowed without my written permission.  Department The University of British Columbia Vancouver, Canada  Date  DE-6 (2/88)  Ffl-  / f f  y  Abstract  The Peace River Arch (PRA) is a regional ~ E - W trending geological structure within the Western Canada Sedimentary Basin whose Phanerozoic history of vertical movements is anomalous with respect to the basin as a whole. Four intersecting ~300-km-long reversed refraction lines within the PRA region in northwestern Alberta and northeastern British Columbia have been interpreted for crustal and upper mantle P-wave velocity structure. The data have been analyzed using a new two-dimensional ray-trace forward modelling algorithm to match travel times and amplitudes of first and coherent later arrivals. An inversion of first arrival travel times along a fan shot profile has been performed to constrain crustal thickness northwest of the arch in a region not sampled by the in-line profiles. 5-waves and the observed spectra of the refraction data have been analyzed to infer a regional Poisson's ratio and Q structure, respectively. The consistency of the seismic models with the observed Bouguer gravity data was studied. The new algorithm for tracing rays and calculating amplitudes in two-dimensional media is based on a simple, layered, large-block velocity model parameterization in which velocity is an analytic function of position. This allows for computationally efficient ray tracing. The user's ability to specify kinematically-similar ray families permits practical and rapid forward modelling of refraction data. In addition, the routine allows for 5-wave propagation, converted phases, multiple and surface reflections, approximate attenuation, head waves, a simulation of smooth layer boundaries, and a reverse ray-direction amplitude calculation. Amplitude calculations are based on zero- and first-order asymptotic ray theory. The main attributes of the routine are illustrated with several examples. The major features of the interpreted structural model of the PRA region are (1) ii  weak to moderate lateral variations in crustal structure with no evidence of significant layering or thick low-velocity zones within the crust, (2) an average sub-basement RMS crustal velocity of 6.6 km/s, average upper mantle velocity of 8.25 km/s and average crustal thickness of 40 km, (3) a high-velocity (> 7.0 km/s) lower crust of 5 to 10 km thickness, (4) westward crustal thinning north of the arch, (5) regional variations in structure that appear related to the N-S Precambrian trends as revealed by aeromagnetic data, including crustal thickness, upper crustal and upper mantle velocities and P P m  character, and (6) subtle variations in structure that may be associated with the E - W trending Devonian axis of the PRA, including a shallowing of high lower-crustal velocities, thickening of the crust, and an anisotropic P P character beneath the arch and along-axis m  low-amplitude P arrivals. A high-velocity lower crust and localized shallowing of high n  lower-crustal velocities are commonly observed in continental rift zones. These features and the ~ E - W trend of the arch perpendicular to the ancient western margin suggest that the PRA originated as a Paleozoic failed-rift. The results of supplementary studies show (1) an average crustal Poisson's ratio of 0.25, (2) Q increases with depth from ~500 to ~1000 in the crust and is ~1000 in the upper mantle, and (3) a seismic-gravity relationship that suggests that localized velocity anomalies of the refraction models are not associated with density anomalies. Also, extended-listen-time processing of a 10-km-long industry Vibroseis reflection line coincident with one of the refraction lines shows prominent dipping events that correlate with the zero-offset two-way travel time of a strong intracrustal reflector and the crustmantle boundary of the refraction model. A series of reflections over 1.5 s terminating at the refraction Moho indicates a complex, possibly layered crust-mantle transition zone of 5 km thickness.  in  Table of Contents  Abstract  ii  List of Tables  ix  List of Figures  x  Acknowledgements  1  2  xiv  INTRODUCTION  1  1.1  Thesis Outline . .  2  1.2  Regional Tectonics  3  1.3  Vertical Movements of Continental Crust  9  1.3.1  Passive Models  11  1.3.2  Active Models  13  1.3.3  Complications and Considerations: Peace River Arch  16  1.4  Previous Geophysical Studies  18  1.5  Seismic Refraction and Reflection Methods  : . . . .  21  ACQUISITION, PROCESSING A N D ANALYSIS OF R E F R A C T I O N DATA  22  2.1  Field Experiment: PRASE  22  2.2  Data Presentation  24  2.3  Data Processing  24  2.3.1  24  Site Locations and Timing iv  2.4  3  2.3.2  Field Instrumentation  37  2.3.3  Sedimentary Velocities  37  2.3.4  Topographic Corrections  38  2.3.5  Timing Inconsistencies  42  2.3.6  Amplitude Inconsistencies .  2.3.7  Arrival Picks  43  2.3.8  Filtering  46  :  43  Method of Interpretation  49  2.4.1  Resolution  49  2.4.2  Modelling Procedure  50  2.4.3  Parameter Uncertainties  51  R A Y T R A C I N G IN T W O - D I M E N S I O N A L M E D I A  53  3.1  Introduction  53  3.2  Velocity Model Parameterization  54  3.3  Ray Tracing  57  3.3.1  Ray Tracing Equations and Their Solution  57  3.3.2  Automatic Determination of Ray Families  58  3.3.3  Multiple and Converted Phases and Source and Receiver Locations  62  3.3.4  Shooting Method Versus Two-point Ray Tracing  63  3.3.5  Smooth Layer Boundary Simulation  63  3.4  Amplitude Calculations  64  3.4.1  Turning and Reflected Rays  64  3.4.2  Reverse Ray-direction Amplitude  68  3.4.3  Head Waves  69  3.4.4  Approximate Attenuation  72  v  3.4.5 3.5  3.6 4  Synthetic Seismograms  74  Examples  74  3.5.1  Efficiency and Accuracy  75  3.5.2  PRASE Line C—shot C2  80  3.5.3  Complex Crustal Model  83  3.5.4  Near-surface Velocity Anomalies and the Statics Problem  83  Discussion  86  INTERPRETATION OF R E F R A C T I O N DATA  89  4.1  Introduction  89  4.2  General Results of the In-line Interpretations  90  4.2.1  Simple, One-dimensional Crustal Structure  90  4.2.2  Near-shot Amplitude Inconsistencies  91  4.2.3  Intracrustal Reflections  4.2.4  Prominent Phases of the In-line Data  93  4.2.5  P P Character  94  . .  m  91  4.3  Line A Interpretation  96  4.4  Line B Interpretation  104  4.5  Line C Interpretation  113  4.6  Line D Interpretation  4.7  Interpretation of Fan Shots  125  4.7.1  Line B Fan Shot  126  4.7.2  Line A Fan Shot  128  ..'  119  4.8  S'-wave Analysis  130  4.9  Q Structure of the Crust and Upper Mantle  136  4.9.1  136  Spectral Ratio Method  vi  4.9.2  Linear Programming Inversion  142  4.9.3  Results  146  4.10 Seismic-gravity Consistency  150  4.10.1 Peace River Region Bouguer Gravity Data  152  4.10.2 Calculation of Theoretical Gravity Response  152  4.10.3 Results  154  4.11 Summary: Regional Velocity Structure 5  P R O C E S S I N G A N D I N T E R P R E T A T I O N O F R E F L E C T I O N D A T A 170  5.1  Introduction  170  5.2  Description of Data Set  171  o.o  rrocessmg  174  5.3.1  Extended-listen Correlation  175  5.3.2  Post-correlation Processing  180  5.4 6  159  Results  182  DISCUSSION  189  6.1  Paleozoic Evolution of the Western Canada Sedimentary Basin  189  6.2  Peace River Region Structural Model  190  6.2.1  Precambrian Features  190  6.2.2  Peace River Arch Features  191  6.3  6.4  Tectonic Models  194  6.3.1  Rift-like Features  197  6.3.2  Passive Versus Active Response  201  Supplementary Studies  202  6.4.1  Gravity Modelling  202  6.4.2  Other Studies  206 vii  6.5  Conclusions and Recommendations  References  List of Tables  3.1  CPU time comparison with Spence et al. routine  4.1  Q model of the crust and upper mantle  146  5.1  Parameters of the reflection data set  174  5.2  Stacking velocity function  181  ix  79  List of Figures  1.1  Study area  2  1.2  Precambrian and Paleozoic trends in the Peace River Arch region  4  1.3  Cross-section of the Peace River Arch basement high  5  1.4  Summary of the history of the Peace River Arch.  .7  1.5  NW-SE sedimentary cross-section  1.6  Passive models of vertical crustal movements  12  1.7  Active models of vertical crustal movements  15  1.8  Previous seismic studies near the Peace River region  19  .  8  2.1 Refraction experiment geometry  23  2.2  Line A - shot A l data  25  2.3 Line A - shot A2 data  26  2.4 Line A - shot A3 data  27  2.5 Line A - shot A4 data  28  2.6  Line B - shot B l data. .  29  2.7  Line B - shot B2 data  •'.  30  2.8 Line B - shot B3 data  31  2.9  32  Line B - shot B4 data. . .  2.10 Line C - shot C l data  33  2.11 Line C - shot C2 data  34  2.12 Line D - shot D l data  35  2.13 Line D - shot D2 data  36 x  2.14 Sediment velocity models from sonic logs  39  2.15 Two-dimensional sediment velocity models  40  2.16 Trace energy versus distance  44  2.17 Ray paths of observed P-wave phases  45  2.18 Filtering of the shot C l profile  47  2.19 Filtering of the shot D2 profile  48  3.1  Velocity model parameterization  55  3.2  Three simplest kinematic ray families  60  3.3  Rays traced in search mode  62  3.4  Smooth layer boundary simulation  65  3.5  Velocity-density relationship  67  3.6  Reverse ray-direction amplitude comparison  70  3.7  Comparison of exact and approximate attenuation.  73  3.8  Example to test efficiency and accuracy  77  3.9  Trade-off between accuracy and CPU time.  78  3.10 Velocity model for application to PRASE data  , .  81  3.11 Application of routine to PRASE data.  82  3.12 Complex crustal model example  84  3.13 Exploration statics problem  '  85  4.1  Expanded view of intracrustal reflections  92  4.2  Comparison of sharp and diffuse P P character.  95  4.3  Calculated and observed comparison for Line A - shot A l  97  4.4  Calculated and observed comparison for Line A - shot A2  99  4.5  Calculated and observed comparison for Line A - shot A3  101  4.6  Line A velocity model  103  m  xi  4.7  Calculated and observed comparison for Line B - shot BI  106  4.8  Calculated and observed comparison for Line B - shot B2  108  4.9  Calculated and observed comparison for Line B - shot B3  110  4.10 Line B velocity model  112  4.11 Calculated and observed comparison for Line C - shot C l  114  4.12 Calculated and observed comparison for Line C - shot C2  116  4.13 Line C velocity model  118  4.14 Calculated and observed comparison for Line D - shot D l  120  4.15 Calculated and observed comparison for Line D - shot D2  122  4.16 Line D velocity model  124  4.17 Line B fan shot interpretation  127  4.18 Line A fan shot interpretation  129  4.19 5-phases in the Line A - shot A3 data  131  4.20 5-phases in the Line C - shot C2 data  132  4.21 P S phase in the Line A - shot A2 data  135  4.22 Example of source spectra  140  4.23 Arrivals used in the Q study  141  4.24 Examples of spectral ratio least-squares fit  144  4.25 Observed and predicted  Q  values for  147  4.26 Observed and predicted  Q  values for  4.27 Observed and predicted  Q  values for  m  Q =100 S  Q =250 s  Q =500 s  148 149  4.28 Q model  151  4.29 Bouguer gravity anomaly map  153  4.30 Gravity modelling of Line A  155  4.31 Gravity modelling of Line B  156  4.32 Gravity modelling of Line C  157 xii  4.33 Gravity modelling of Line D  158  4.34 One-dimensional velocity model  160  4.35 Basement velocity  162  4.36 RMS crustal velocity  163  4.37 Depth to 6.8 km/s  165  4.38 Thickness of greater than 7.2 km/s crust  166  4.39 Velocity at the base of the crust  167  4.40 Crustal thickness and P P character  168  4.41 Upper mantle velocity  169  5.1  Location of reflection line within Line B model  171  5.2  Velocity-depth profile at position of reflection line  172  5.3  Geometry, stack fold and elevation profile  173  5.4  Normal, fixed and self-truncating sweep correlation.  176  5.5  Sweep length and high frequency as a function of time  178  5.6  Representative common-shot gathers. . .  179  5.7  Stack of the reflection section  183  5.8  Stack of the sedimentary reflection section  184  5.9  Expanded view comparison of deep near-vertical and wide-angle events. . 186  6.1  Isostatic gravity anomalies of western Canada  193  6.2  Proterozoic and Paleozoic rifts of North America  196  6.3  Generalized velocity structure of continental rifts  198  6.4  Model of transverse anisotropy required by gravity modelling  204  m  xm  Acknowledgements  I wish to thank my thesis supervisor, Bob Ellis, for his guidance and confidence in me throughout my years at UBC. I have benefited from the experienced advice of Ron Clowes on numerous matters and Ron, along with Bob, encouraged me to enter a Ph.D. program. I am sincerely grateful to Randell Stephenson for it is perhaps his efforts which are most responsible for this thesis. Randell proposed and played the key role in acquisition of funding for the refraction experiment and provided leadership in the field operations. Subsequently, Randell's continued interest and participation extended far beyond his formal role as scientific authority on the analysis and interpretation contracts. He was keenly interested in my work and was always willing to discuss my research at each stage. I learned much from these discussions and Randell must be partly credited with some of the ideas presented in this thesis, particularly with regard to matters geological and tectonic. Randell's approach was always straightforward and honest and for this I am most thankful. I acknowledge the significant contributions of Zoli Hajnal, Patrick Morel-a-l'Huissier, Bob Mereu, Doug Northey, Gordon West and Ernie Kanasewich to the experimental design and/or field program. In addition, this thesis would not have been what it is without the dedicated efforts of all those who assisted in the field experiment ensuring the high quality of the refraction data: Andy Boland, Doug Brown, John Cassidy, Ben Ciammaichella, Connie Cudrak, Jeff Drew, Michael Ehling, Bob Meldrum, Geoff Coleman, Brian Reilkoff, Doug Scott, Dave Wilkonson, Tina Hewak, Karin Michel, Karen Morrissey, Andrew Gunstensen, Trevor Shortt, David Anholt, Sean Lavery and Blake Wright.  xiv  I thank Bob Krider for his help in the unpleasant task of preparing the final PRASE data tape, Bob Meldrum for helping me to understand the instrumentation, and John Amor for helping me to process the Vibroseis data. I also thank Lydia Wong, Susan Clegg and Carole Rennie for their help and accurate secretarial skills, usually on very short notice. I wish to thank all of the people in the Department of Geophysics and Astronomy I have met over the years for their friendship. With the many wine and cheeses, cake parties, golf tournaments, baseball and hockey games, and the occasional attempt at darts and badminton, I think the last four years have been the most enjoyable of my life. A number of people deserve special mention for the help or advice they have given me during my studies: Dave Aldridge, Andy Boland, Andy Calvert, John Cassidy, David Dalton, Sonya Dehler, Stan Dosso, Dave Lumley, Tim Scheuer, Don White, Ken Whittall and Matt Yedlin. Some of you may have forgotten your "contribution", but I have not. An extra special thanks to Jeff Drew and my office mate John Cassidy, and his wife Lynn, for their friendship and the chance to meet my fiancee: Jeff accidentally injured my field partner, John, during the field experiment so Bob had to bring in Connie as my new field assistant. I thank Connie Cudrak for her willingness to work with an "intense" field person and her constant friendship, support and confidence in me since then. I also thank my brother, Barry, for his friendship and help while we have been university students together and, more importantly, for his humour in the special way that he sees things. I thank Lorraine and Ron MacLeod for making my stay in Vancouver an even more enjoyable (nutritious) experience. Finally, I thank my parents and brothers for their constant love and support. Principal funding for the data acquisition and analysis has been through contracts from the Institute of Sedimentary and Petroleum Geology, Geological Survey of Canada.  xv  Supplementary support for the refraction data acquisition was provided by Canterra Energy Limited and Chevron Canada Resources Limited. Bill Davitt suggested the use of uncorrected Vibroseis data and its existence within the industry. Shell Canada Limited provided the Vibroseis reflection data set and I thank Clive Watson and Keith Wilkinson for allowing me to bother them with so many questions about the data. Personal financial support was provided principally by a Postgraduate Scholarship from the Natural Sciences and Engineering Research Council of Canada.  xvi  Chapter 1  INTRODUCTION  In 1985, a large-scale crustal refraction survey was conducted in the Peace River region of northwestern Alberta and northeastern British Columbia (Figure 1.1) to provide a regional P-wave seismic velocity model of the crust and upper mantle and thereby permit a better understanding of the tectonic processes that generated the anomalous vertical movements of the Peace River Arch (PRA). In-line shots were recorded along four intersecting lines, each approximately 300 km in length. Fan shots were recorded along two of these lines. The resulting data set is of exceptional quality—it has high signal-to-noise ratio, near-uniform receiver spacing, reliable amplitude control and numerous coherent later arrivals. This thesis presents the results of a two-dimensional travel time and amplitude analysis of the complete in-line data set through ray-trace forward modelling. In addition, an interpretation of the first arrival times of the two fan shots is included to constrain crustal thickness northwest of the arch in a region not sampled by the in-line profiles and determine upper mantle velocity variations in the centre of the survey region. The travel time of the S S phase is analyzed to infer an average crustal Poisson's m  ratio and a spectral ratio method is applied to the refraction data to obtain a regional one-dimensional crustal and upper mantle Q model. Finally, extended-listen-time processing of a 10-km-long industry Vibroseis reflection line has been performed and results compared with a coincident refraction model. Regional gravity and aeromagnetic data are considered as both a check on the refraction models and as an aid in interpreting their structural trends.  1  Figure 1.1: Location of the study area within the Western Canada Sedimentary Basin. Thickness of the Phanerozoic sediments east of the Cordilleran foreland thrust belt is indicated by contours in km. Adapted from Porter et al. (1982). 1.1  Thesis Outline  The remainder of chapter 1 provides the background for the PRA seismic experiment, discussing both the regional tectonics and models for vertical crustal movements as well as previous geophysical studies and the advantages and disadvantages of the seismic refraction and reflection methods. Chapter 2 describes the refraction experiment, data processing and the method of interpretation. Chapter 3 describes the new two-dimensional ray tracing algorithm developed and applied in the travel time and amplitude forward modelling of the refraction data. Chapter 4 presents the results of the refraction modelling including the qualitative analysis of two fan shot profiles and the estimation of Poisson's ratio and Q structure as well as a check of the seismic-gravity consistency. Chapter 5 presents the results of processing industry Vibroseis reflection data coincident with one of the refraction lines. Chapter 6 relates the results of the thesis to the tectonic  Chapter 1. INTRODUCTION  3  history of the PRA.  1.2  Regional Tectonics  The Precambrian continental craton that forms the basement of the Western Canada Sedimentary Basin consists mainly of rocks of the Churchill province of the Canadian Shield (Porter et ai., 1982). The Precambrian tectonics of the PRA region have been interpreted from regional aeromagnetic trends extrapolated from the exposed Slave Province of the Canadian Shield to the northeast (Hoffman, 1987; 1988) and from radiometric studies of basement cores (Ross et al., 1988). Figure 1.2 shows the three dominant, roughly N-S trending basement terranes of the survey region, referred to, from east to west, as the Taltson Magmatic Arc (TMA), the Chinchaga Low (CL) and the Wilson Island Terrane (WIT). These are zones of predominantly positive, negative and positive residual aeromagnetic anomaly, respectively. The Taltson Magmatic Arc consists of plutons associated with an Early Proterozoic collision and subduction beneath the Archean Rae Province, lying to the east of the PRA region, and the north-south trending Chinchaga Low may represent a suture zone associated with this eastward-dipping subduction (Ross et al., 1988). The Wilson Island Terrane is the same age as, and may be a part of, the Taltson Magmatic Arc (Ross et al., 1988). The Western Canada Sedimentary Basin is a northeasterly tapering wedge of sedimentary rocks up to 6 km thick extending from the Cordilleran foreland thrust belt to the Canadian Shield. In the tectonic development of the Western Canada Sedimentary Basin, two subsidence mechanisms are relevant: thermal contraction associated with the establishment of a rifted continental margin in the late Precambrian, and flexural downwarp of the lithosphere associated with the loading of the adjacent fold and thrust belt beginning in the Middle Jurassic (e.g., Beaumont et ai., 1982). Thus, the Western  Chapter 1. INTRODUCTION  4  120°  119°  118"  7  H °  LONGITUDE (W) Figure 1.2: The study area showing the location of the Peace River Arch as defined by an anomalous sedimentary cover. The Paleozoic trends (arches and embayments) are from McCrossan and Glaister (1964) as compiled by J.E. Barclay (pers. comm., 1986). The Devonian trend continues ~100 km further to the northeast and the Carboniferous trend continues ~50 km westward. The negative (dotted) and positive (elsewhere) residual aeromagnetic domains define the Precambrian basement terranes (TMA - Taltson Magmatic Arc; CL - Chinchaga Low; WIT - Wilson Island Terrane). Aeromagnetic data are only available north of the dashed line (Geological Survey of Canada, 1986). Precambrian basement topography with reference to sea level where sampled by drilling is indicated by contours in km (U. S. Geological Survey, 1972). Thickening of the sediments southwest of the 3 km contour is assumed based on extrapolation of observed basement trends. See Figure 1.5 for NW-SE cross-section.  Chapter  1.  5  INTRODUCTION  rianVJJ;::/^s^Reef'T^I^i^dic ^.-~..y.. ^ . . . ~ „  ^.  Devonian  . \-f^*eq^gLower D e v o n i a n  Granite Wash  Figure 1.3: Diagrammatic cross-section of the Devonian stratigraphic units onlapping the Precambrian basement high. The Devonian section missing over the arch from the Lower to Upper Devonian is estimated to be 800-1000 m thickness. The uplift averages about 140 km in width measured at the level of the middle Devonian. Adapted from Cant (1988). Canada Sedimentary Basin represents a foreland basin superimposed on a rift margin basin. The PRA is situated within the Western Canada Sedimentary Basin and has a Phanerozoic history of anomalous vertical motions with respect to the basin as a whole (Porter et al., 1982). The ENE-WSW (hereafter referred to as E-W) trending PRA (~350 km long, 140 km wide) is perpendicular to the craton margin and depositional trend of the Western Canada Sedimentary Basin (Figure 1.2). The arch is superimposed on the pre-existing cratonic structure and is identified on the basis of an anomalous Phanerozoic sedimentary cover within the region: no Cambro-Ordovician sediments are preserved in the Peace River region, hence its stratigraphic designation as an arch. The granitic basement along the arch axis stands above the regional trend (Figure 1.2); Cant (1988) estimates 800 to 1000 m of uplift. The basement is onlapped by Mid to Upper Devonian sediments, with only the uppermost Devonian strata extending across and covering the arch (Figure 1.3). It is thought that the arch was uplifted and eroded after early Paleozoic sedimentation (deMille, 1958), although it is possible that the arch was a positive feature throughout the Early Paleozoic (Stelck et al., 1978).  Many of the Devonian, Mississippian and  Chapter 1. INTRODUCTION  6  Pennsylvanian units above the arch are cut by major normal faults approximately parallel and perpendicular to the arch axis, forming horsts and grabens (Cant, 1988). The Peace River region contains the only major sequence of Permian strata east of the present deformed belt. Therefore, near the end of the Devonian Period, through normal faulting, the arch became a depocentre with greater subsidence and a thicker accumulation of sediment than elsewhere (Lavoie, 1958; Procter and Macauley, 1968) and/or became a site of favored preservation of these sediments (Sikabonyi and Rodgers, 1959). The Paleozoic to Mesozoic history of the PRA is summarized in Figure 1.4. Until the Early Cretaceous, the arch's fault system continued to control sedimentation patterns (Williams, 1958; Stott, 1984). During the foreland basin phase of the Western Canada Sedimentary Basin, the PRA region may have been a favored site of Early Cretaceous subsidence (Rudkin, 1964), but the arch may have again inverted to become a positive feature by the close of the Cretaceous Period (Zolnai, 1982). A NW-SE sediment cross section across the PRA is shown in Figure 1.5. The P R A was discovered by exploratory drilling in 1948 after the Precambrian Churchillian basement was reached below an abbreviated Devonian section (Stelck et al., 1978). In the 1950's, numerous oil and gas discoveries were made and the PRA area continues to be an important producing region. Faulting associated with the arch plays a major role in influencing the petroleum geology of the area (Cant, 1988). Hydrocarbons occur in Paleozoic reef and structural traps and structurally influenced stratigraphic traps in the Mesozoic (Cant, 1988). No metallic minerals have as yet been discovered in the subsurface of the PRA region (Podruski, 1988).  Chapter 1. INTRODUCTION  7  Cambrian to Early Devonian  Mid to Late Devonian  A  ^  ^  Devonian *2Z/ \ " ~ _  ~  Mississippian  fMississippians rr.:: Pc  d  Pennsylvanian  Pennsylvanian  s  sMississippiang •Devonian',^, -rvT: :  :  -  \  1  =—zr -  Mesozoic • Cretaceous;  Figure 1.4: Diagrammatic summary of the history of the Peace River Arch. Arrows indicate direction of motion and Pc is the Pecambrian basement, (a) Cambrian to Early Devonian—tectonic uplift of 800-1000 m; (b) Mid to Late Devonian—normal faulting, graben formation, onlap by Devonian units; (c) Mississippian—block subsidence; (d) Pennsylvanian—grabens reactivated; (e) Mesozoic—general sagging of the arch. Figure from Cant (1988).  Chapter 1. INTRODUCTION  8  VERTICAL : HORIZONTAL = 1 0 0 : 1  Figure 1.5: Northwest-southeast sedimentary cross-section across the Peace River Arch (see Figure 1.2). Only the major depositional units are shown. Pc—Precambrian; C—Cambrian; D—Devonian; M—Mississippian; P—Permian; T—Triassic; J—Jurassic; iK—Lower Cretaceous; uK—Upper Cretaceous. Figure from Podruski (1988).  Chapter 1. INTRODUCTION  1.3  9  Vertical Movements of Continental Crust  The PRA is of particular interest because of its anomalous vertical movements during the Phanerozoic development of the Western Canada Sedimentary Basin. Therefore, a review of the tectonic processes thought to be responsible for such movements within continents is presented. Movements within the PRA region are classified as epeirogenic— broadscale, mainly vertical movements of the crust that produce plateaus, basins and arches— as opposed to orogenic—mountain building associated with compressional tectonics. Orogeny is relatively well understood within the plate tectonic framework; however, within continents, epeirogeny is generally not well understood (Press and Siever, 1986). The PRA is an example of relatively small-scale epeirogeny since it can be contrasted with larger-scale features such as, for example, the Colorado Plateau or the Western Canada Sedimentary Basin, within which the PRA represents one of many intrabasinal features (Porter et al, 1982). The present day location of the PRA is within the Western Canada Sedimentary Basin of the Interior Platform of the North American continent; however, during its development, it has been situated within or adjacent to a number of different tectonic environments (Porter et ai., 1982). Its early history is least certain, but the PRA may have been an intracratonic feature prior to rifting and the establishment of the passive margin off the western margin of the craton in Late Precambrian (Bond and Kominz, 1984; Bond et ai., 1985). Prior to commencement of the foreland basin phase of the Western Canada Sedimentary Basin and the development of the Rocky Mountains in the Jurassic, the PRA was situated on the western edge of the North American craton. During the formation of the Rocky Mountains, the western North American margin, and hence the PRA region, was subjected to predominantly compressional forces (Monger and Price, 1979). Thus, to understand the history of the PRA requires an understanding and  Chapter 1. INTRODUCTION  10  consideration of the large-scale processes acting since at least Late Precambrian along the western margin of the North American craton. The processes of classical plate tectonic theory are capable of explaining the development of the observed features of the ocean basins and continental margins. Within continents, the age and structural complexity of the crust requires an extension or modification of plate tectonic theory to explain features such as basins, plateaus, rifts and arches. The Wilson cycle (Wilson, 1966) provides a link between continental evolution and plate tectonics (e.g., Burke, 1980). Meissner (1986, p. 241) summarizes the Wilson cycle as (1) continental rifting, (2) development of passive margin and oceanic spreading ridge, (3) subduction of oceanic lithosphere at former passive margin, (4) closing of ocean and continent-continent collision, (5) continental accretion and mountain building, (6) erosion of mountain belt and disappearance of mountain root, and (7) remaining zone of weakness. The Wilson cycle and its associated driving forces—mantle convection, plumes and contraction by cooling—will be important when considering continental epeirogeny. Two broad classes of models, designated passive and active, are proposed for the origin and development of intracratonic features and are separated on the basis of whether the associated tectonic processes are acting laterally from a distance (passive) or directly beneath the feature within the lithosphere or asthenosphere (active). A combination of the two models can be envisaged. For example, lateral tensional forces may lead to asthenospheric upwelling or injection of mantle material into the lower crust, a passive rift (Morgan and Baker, 1983). The two basic models will now be considered in further detail. This discussion will serve to define terms and introduce ideas that will be used when the tectonic processes associated with the evolution of the PRA are considered in chapter 6. Included in this review are the geophysical signatures associated with passive and active processes which will allow an assessment of the details of these processes that might be resolved by the geophysical data available in the PRA region.  Chapter 1. INTRODUCTION  1.3.1  11  Passive M o d e l s  Passive models of small-scale epeirogeny suggest that anomalous horizontal or vertical forces acting locally are due to distant sources within or upon the lithosphere. Consequently, a local anomalous crustal, upper mantle, or deeper structure or process is not necessarily implied, although these may initially or subsequently modify the epeirogenic movements. Passive models are perhaps more easily understood and accepted since the relevant driving mechanisms are confined to forces acting within the lithosphere, where processes and properties are reasonably well understood. There are two subclasses of passive models: those associated with lithospheric loading and unloading, and those associated with plate boundary interactions. The lithospheric loading model is based upon the concept of flexure (Walcott, 1970). The distant load may either be deposited on the surface or intruded at depth from below the lithosphere. If the lithosphere is loaded or unloaded by deposition or erosion, the lithosphere will respond flexurally and a peripheral "bulge" or "sag" is formed in the vicinity of the loading or unloading (Figure 1.6a). The peripheral response is due to the flexural rigidity of the lithosphere and its position, elevation and timing is associated with the flexural properties of the lithosphere, such as its flexural rigidity and relaxation time for elastic and viscoelastic models. The usual source of loading is sediment, but the load may also be ice or water. A complication of this model involves the consideration of two or more points of loading so that a flexural interaction between the associated peripheral features is possible (Figure 1.6b). These qualitative ideas form the basis of quantitative flexural models that account for the broad-scale features of stratigraphy within sedimentary basins and include a network of interrelated basins, arches and domes (e.g., Quinlan and Beaumont, 1984). The second subclass of passive models considers the horizontal stress field within the  Loading  Flexural moat  -3-  I Peripheral bulge  J  I  Unloading  Peripheral sag  a  Figure 1.6: Passive models of vertical crustal movements (qualitative and not to scale), (a) Flexural response of lithosphere to loading and unloading. Figure from Beaumont et al. (1988). (b) Flexural interaction between two lithospheric loads; cratonic basin on left, foreland basin on right. Upper panel is a plan view prior to coupling and corresponds to first cross-section. Subsequent cross-sections show the interaction as the basins move closer together. Figure from Quinlan and Beaumont (1984). (c) Deformation due to compressional forces. Figure from Press and Siever (1986).  Chapter 1. INTRODUCTION  13  lithosphere due to interactions with other plates around its margin. These interactions are those of conventional plate tectonics: compression at collisional boundaries, extension at rift margins, and shear at transform boundaries. As a result of these stresses, the lithosphere may deform through warping, buckling, faulting or flow resulting in regions of relative uplift and downwarp (Bott, 1982, pp. 314-315) (Figure 1.6c). An important model of this type is one in which extensional forces initiate continental rifting through stretching, thinning or faulting of the lithosphere (Morgan and Baker, 1983). A combination of the plate interaction model with the loading model is the development of a foreland basin at a collisional margin: the load of the over-riding adjacent fold-thrust belt causes a flexural downwarp of the lithosphere at its margin (Beaumont, 1981). The application of a passive model requires knowledge of the distant loading history which may be relevant to the study area and/or an understanding of the historical plate tectonic environment around the lithospheric plate of interest. This information is often difficult or impossible to obtain in detail and generally only concerns large-scale processes. By the nature of the passive model, no inherent structural or compositional anomaly is required within a region suspected of past epeirogenic movement. Thus, the present day response to all geophysical testing may yield nothing anomalous within the crust or upper mantle with respect to large-scale regional averages. On the other hand, the presence of anomalous structure does not negate a proposed passive model since prior or subsequent tectonic processes may be responsible. Instead, it is necessary to consider the surrounding loading and plate tectonic environment to validate a passive model. 1.3.2  Active Models  Active models of epeirogeny suggest that the driving mechanism is situated within the underlying lithosphere or asthenosphere. Active models are perhaps viewed with more  Chapter 1. INTRODUCTION  14  skepticism since the driving mechanism is often within the more "mysterious" asthenosphere and there may at present be no evidence for the mechanism's existence beneath or adjacent to the study area if the epeirogeny occurred sufficiently long ago. Numerous active models have been proposed to explain continental epeirogeny—many are interrelated but they can be divided into two subclasses: thermal expansion and contraction due to igneous intrusion in the lower lithosphere or positioning over a mantle plume, and isostatic warping due to igneous intrusion or a phase change in the lower lithosphere. Thermally induced forces are generally recognized as the principal origin of stress within the lithosphere; these forces are due to thermal expansion and contraction (Hinze et al., 1980). Many different models have been proposed (e.g., McGetchin et al., 1980); for example, thermal anomalies due to igneous intrusions in the lower lithosphere (e.g. Artyushkov et al., 1980) or mantle penetrative convection resulting in regional heating and tensional and compressional patterns within the lower lithosphere (Hinze et al., 1980). In the latter model, stationary continental lithosphere overlies a mantle plume and rising material spreads out in the upper asthenosphere producing stresses in the overlying plate.  Regional expansion and contraction occurs over rising (hotter) and  sinking (cooler) mantle material, respectively (Figure 1.7a). A profound mantle heat source may serve to initiate continental rifting (Figure 1.7b) (e.g., Burke and Dewey, 1973). In the isostatic warping model, density variations due to intrusions or phase changes in the lower lithosphere will be counter-balanced by isostatic adjustment causing surface uplift or depression (Figure 1.7c). The intrusion model generally proceeds with a lower density intrusion of basic magma into the lower crust or upper mantle due to gravitational instability, causing surface uplift (e.g. McKenzie, 1984). Subsequent cooling forms a high density body and results in subsidence. The phase change model can be related to the thermal model whereby a hot mantle diapiric plume causes a conversion from gabbro to  Chapter 1. INTRODUCTION  15  CRUSTAL RIFTING THERMAL EXPANSION AND CONTRACTION  .  A  Figure 1.7: Active models of vertical crustal movements (qualitative and not to scale), (a) Thermal expansion and contraction due to sub-crustal heat source, (b) A profound heat source initiating continental rifting, (c) Isostatic warping due to lower crustal intrusion or phase change. Figure from Hinze et al. (1980). more dense eclogite in the crust (Haxby et ai., 1976). If the mantle source is removed or cools by conduction, subsidence will occur under the load of the eclogite. This discussion has shown that active models of epeirogeny result in a local anomalous lithospheric structure. In particular, a high density lower crust and/or variations in crustal thickness or the crust-mantle transition might be expected. By analogy with the passive model, the presence of an anomalous structure does not in itself validate an active model because of the possibility of pre-existing structure or the effects of contemporaneous plate dynamics. The geophysical signature associated with active models is complicated by the inferred timing of the epeirogenic activity. For example, heat flow anomalies associated with Pre-Cenozoic continental rifting events have now decayed to immeasurable levels •(Morgan, 1982). Thus, for a thermal driving mechanism, surface heat flow measurements may not serve as a strong indicator. A similar argument can be made for local seismicity  Chapter 1. INTRODUCTION  16  levels. Present day seismicity is controlled by the current local stress regime, although the distribution and focal mechanism pattern may delineate faults that were once associated with the epeirogenic event (Hinze et al., 1980). Aeromagnetic data is capable of evaluating the magnetic properties of the crust down to depths at which the Curie point is reached, 15-50 km for continental crust (Zietz, 1980). Since crystalline rocks usually have a magnetic signature associated with them, the magnetic method can be used in geological mapping for both subsurface lithology and structure, including faults. Thus, for relatively shallow anomalous structure due to an active model, magnetic data may be relevant. If an active process has altered the lower crust or upper mantle, gravity and seismic data are prime indicators. Density anomalies are most likely within the lower lithosphere or due to crustal thickness variations. Density anomalies may correspond to seismic velocity anomalies which can be revealed best through wide-angle seismic reflection/refraction data. Crustal thickness variations could be detected by both seismic refraction and reflection techniques. 1.3.3  Complications and Considerations: Peace River A r c h  In the study of any epeirogenic feature, two complications arise: (1) the influence of preexisting lithospheric structure and overprinting by subsequent plate tectonic processes, and (2) a combination of passive and active processes acting simultaneously and/or at different times. The former complication can be lessened only through knowledge of both the prior tectonic processes which might have formed and shaped the region of interest and the plate tectonic environment surrounding the region during and after the epeirogeny. Of special importance is the presence of pre-existing faults or zones of weakness within the crust at which epeirogenic movements might be concentrated (Hinze et al., 1980; Morgan and Baker, 1983). The latter complication is partially due to the  Chapter 1. INTRODUCTION  17  fact that the models presented in this discussion are oversimplified in that it is often impossible to apply one model to an observed epeirogenic feature without another model being of equal or secondary importance by necessity. For example, extensional forces (passive model) may initiate rifting at a point in the lithosphere where it is weakest over an existing thermal anomally (active model). Thus, it may be important to consider more than one driving mechanism and then establish, if possible, which was dominant. The PRA is a likely candidate for multiple driving mechanisms since its anomalous history is traced from as early as Late Precambrian to Late Mesozoic (over 500 Myr period). During this time, the P R A region may have been an intracratonic feature at first, then adjacent to a rift margin, and finally subjected to the effects of an adjacent collisional boundary. Throughout these different tectonic environments, the PRA region continued to be an anomalous feature within the Western Canada Sedimentary Basin. The relevant geophysical data for the PRA region includes Bouguer gravity and residual aeromagnetic data.  As discussed in section 1.4, both data sets contain no  obvious trends that might be associated with the E-W trending axis of the PRA. Instead, the gravity data suggest a near-uniform platformal crust depressed by an adjacent isostatically-compensated fold-thrust belt as interpreted further south by Berry and Forsyth (1975). The N-S aeromagnetic trends are interpreted to represent pre-existing Precambrian basement features (Hoffman, 1987; 1988). This thesis presents the interpretation of crustal refraction data within the PRA region. The areal extent of the PRA (350 X 140 km) suggests that the refraction data will be able to distinguish between a uniform or normal crust and upper mantle structure or anomalous structure.  This may permit the dominant epeirogenic model(s) to be  established. In particular, an anomalous crustal and/or upper mantle structure related to the PRA would support an active model and the details of the anomaly may allow a specific tectonic process to be established.  Chapter 1. INTRODUCTION  1.4  18  Previous Geophysical Studies  Within the PRA region, gravity and aeromagnetic data are available. The Bouguer gravity field within the region varies smoothly from about -70 mGal in the northeast to -90 mGal in the southwest with a rapid decrease to less than -140 mGal in the extreme southwest associated with the adjacent Rocky Mountains (Figure 4.29). The predominantly N-S trending residual aeromagnetic anomalies were discussed in section 1.2. Much of Alberta, including the PRA region, appears to be aseismic (Wetmiller, 1986), although the ability to detect earthquakes within the PRA region of magnitude greater than 3 has only been possible since 1965 (Milne et al., 1978). A NE-SW compressive stress field exists today thoughout the Interior Platform of western Canada (Fordjor et al., 1983). Heat flow within the Western Canada Sedimentary Basin is high in comparison with other Precambrian platform areas. Majorowicz and Jessop (1981) suggest that the heat flow pattern is controlled by basin-wide groundwater flow patterns, while Bachu (1988) suggests that the dominant controlling mechanism is crustal thickness and/or radiogenic basement heat generation. The PRA region is a site of anomalously high basement rock heat generation (Majorowicz and Jessop, 1981). A relatively low geothermal gradient is observed in the British Columbia PRA region. Lam and Jones (1984) note that magnetic highs within the PRA region correspond to zones of high basement heat generation. No previous seismic studies have been performed in the PRA region except that of industry reflection profiling. The nearest deep crustal seismic studies include: Forsyth et al. (1974), Bennett et al. (1975) and Berry and Forsyth (1975), to the west and south in British Columbia; Ganley and Cumming (1974) and Sommerville and Ellis (1972), approximately 400 km southeast; and Mereu et al. (1977), approximately 200 km south of the PRA region (Figure 1.8).  Bennett et al. (1975) estimate a crustal thickness  of between 50 and 60 km and a low-velocity layer about 8 km below the Mohorovicic  Chapter 1. INTRODUCTION  19  U . S . A . Figure 1.8: Previous seismic studies near the PRA region. The study area of this thesis is indicated. The long straight lines are the approximate refraction profiles of Mereu et al. (1977)—1, Bennett et al. (1975)—2, and Chandra and Cumming (1972)—3. The short lines are the reflection profiles of Ganley and Cumming (1974)—4, Clowes et al. (1968), Kanasewich et al. (1969) and Clowes and Kanasewich (1970)—5, and Cumming and Chandra (1975)—6. The triangles are the seismograph stations used in the teleseismic studies of Sommerville and Ellis (1972)—7, and EUis and Basham (1968)—8.  Chapter 1. INTRODUCTION  20  discontinuity (Moho) from a refraction survey 400 km south along the southern Rocky Mountain Trench. Ganley and Cumming (1974) report an intracrustal boundary at 20 km depth dipping 15° to 20° to the southeast, a tentative "Riel" discontinuity at 32 km depth, and the base of the crust at about 35 km from reflection profiles near Edmonton, Alberta.  This intracrustal dip is consistent with the results of Ellis and  Basham (1968) on the basis of anomalies in the direction of approach of teleseismic P arrivals. Sommerville and Ellis (1972) present evidence from a study of teleseismic Pcoda spectral ratios near Leduc, Alberta, for a low-velocity zone within the crust at 12 km depth. Mereu et ai. (1977) conclude the following from a refraction survey across the Rocky Mountains near Jasper, Alberta: (1) roughly uniform crustal thickness of 50 km, (2) Bouguer gravity anomalies associated with lateral density variations within the crust, (3) crustal velocities increasing from 6.0 to 7.3 km/s from the basement to the base of the crust with any possible intracrustal discontinuities being irregular and not continuous, (4) Q values within the upper crust, lower crust, vicinity of the crust-mantle boundary, and upper mantle of less than 500, greater than 500, less than 200, and greater than 500, respectively, and (5) a high velocity (8.5-8.8 km/s) layer within the upper mantle at 60 km depth. In southern Alberta, deep seismic studies include those of Clowes et al. (1968), Kanasewich et ai. (1969), Clowes and Kanasewich (1970), Cumming and Chandra (1975) and Chandra and Cumming (1972). Kanasewich et ai. (1969) analyzed deep reflection, gravity and magnetic data and proposed the existence of an ~ E - W trending Precambrian rift. Depths to the intracrustal Riel discontinuity and crust-mantle boundary were 2637 km and 38-47 km, respectively. Clowes and Kanasewich (1970) estimated the average Q value in the sediments and crust to be 300 and 1500, respectively, from near-vertical incidence deep reflection data.  Cumming and Chandra (1975) determined depths of  35 and 46 km to the Riel and crust-mantle boundary, respectively, from near-vertical  Chapter 1. INTRODUCTION  21  reflections. Chandra and Cumming (1972) estimated crustal thickness to be 45-50 km from a refraction survey across and to the east of the Rocky Mountains. To summarize the results of the seismic studies, crustal thickness south of the PRA region is ~50 km beneath and immediately east of the Rocky Mountains and is ~45 km in southern Alberta. The depth to the Riel discontinuity is ~35 km in central and southern Alberta.  1.5  Seismic Refraction and Reflection Methods  Recent studies have shown the advantage of performing both seismic refraction and reflection experiments in crust and upper mantle studies (e.g., Mooney and Brocher, 1987) because both methods have inherent limitations in their ability to provide a complete subsurface structural image (e.g., Braile and Chiang, 1986). The refraction method suffers from poor lateral resolution due to relatively large station spacing and the averaging of structure between source and receiver along predominantly horizontal ray paths. The reflection method requires good subsurface velocity control and is dependent on the presence of sharp impedance contrasts. The advantage of the refraction method is high resolution of vertical velocity structure and high signal-to-noise ratio due to wideangle reflection/transmission coefficients and the use of first arrivals. The advantage of the reflection method is high lateral resolution and the processing of redundant data to enhance the signal-to-noise ratio. One additional reason for the application of both methods is the comparatively low and high frequency content of refraction and reflection data, respectively, so that each images crustal features of different scale.  Chapter 2  ACQUISITION, PROCESSING A N D ANALYSIS O F REFRACTION  2.1  Field Experiment:  DATA  PRASE  The 1985 Peace River Arch Seismic Experiment (PRASE) was a cooperative effort involving personnel and instrumentation from the Geological Survey of Canada and five universities: Alberta, British Columbia, Saskatchewan, Toronto and Western Ontario. A total of 56 vertical-component portable seismograph units were used to record 19 shots along 4 lines, each line approximately 300 km in length. The complete refraction data set comprises over 1000 seismograms. Lines A and D are roughly parallel about 100 km apart trending N-S across the axis of the PRA, Line B trends roughly E - W over the axis of the arch and Line C trends E-W about 100 km south of the arch axis (Figure 2.1). Lines C and D were each recorded with the 56 recording units distributed along the lines to record end shots ( C l and C2 and D l and D2, respectively). The resulting average receiver spacing is 4.8 km and 5.3 km for lines C and D, respectively. Line A was recorded in two stages with the 56 recorders first distributed along the southern half of the line to record end shots ( A l and A3), a centre shot (A2) and a fan shot (A4). The recorders were then moved to the northern half of the line and the same shot sequence repeated, resulting in an average receiver spacing of 3.1 km along Line A. Recording of Line B was similar to that of Line A except that for the fan shot (B4), the 56 recorders were distributed along the entire line. The average receiver spacing for the in-line shots of line B is 2.9 km and for the fan shot is 5.7 km.  22  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  L O N G I T U D E (W) Figure 2.1: Shot points and receiver sites along the four refraction lines. The shot points labelled A l , A2 and A3 identify Line A and similarly for Lines B, C and D. Fan shot profiles were recorded along Lines A and B from shots points A4 and B4, respectively. The location of wells from which sonic logs were used to obtain the sediment velocity model are indicated as well as the Devonian axis of the Peace River Arch.  23  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  24  Charge sizes of 600 and 2000 kg were used as sources; the smaller size was reserved for lines A and B in which all receivers were within half the length of the line. The explosives were arranged in shot-hole arrays, each hole 30 m deep and a minimum of 25 m apart containing 200 kg of Geogel 60%. Further details of the field experiment are described by Ellis et al. (1986).  2.2  Data Presentation  The complete data set comprising 10 in-line and 2 fan shot profiles is presented in Figures 2.2 to 2.13 in trace normalized, reduced travel time format. The data in these figures are "raw"; no filtering or topographic corrections have been applied, only the removal of unusable traces. For the in-line profiles, distance is measured along each line with respect to shots A l , BI, C l and D l for lines A, B, C and D, respectively. Since the prominent F-phases are best observed in trace normalized format, examples of each have been labelled and are described in section 2.3.7. 2.3 2.3.1  Data Processing Site Locations and Timing  Shot and receiver sites were located on 1:50,000 scale topographic maps providing a nominal accuracy of 50 m in surface location and 10 m in elevation. An accuracy of 10 ms in timing was achieved through the use of WWVB and WWV radio signals for both shots and receivers, either recorded simultaneously with the seismic signal or used to rate internal clocks at instrument deployment and retrieval. A linear clock drift between deployment and retrieval was assumed.  Chapter 2. DATA ACQUISITION,  PROCESSING  AND  ANALYSIS  o  o Z\  6  (S)  9  8/a-i  C  0  Figure 2.3: Line A - shot A2 data. Each trace is scaled to a common maximum amplitude. See section 2.3.7 for a description of the phases labelled.  to Ob  Figure 2.4: Line A - shot A3 data. Each trace is scaled to a common maximum amplitude. See section 2.3.7 for a description of the phases labelled.  Figure 2.5: Line A - shot A4 data. Each trace is scaled to a common maximum amplitude. See section 2,3.7 for a description of the phases labelled.  to CO  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  Figure 2.7: Line B - shot B2 data. Each trace is scaled to a common maximum amplitude. See section 2.3.7 for a description of the phases labelled.  CO  o  Figure 2.8: Line B - shot B3 data. Each trace is scaled to a common maximum amplitude. See section 2.3.7 for a description of the phases labelled.  155  165  175  185 A Z I M U T H  195  205  215  ( D E G )  225  % £ Co  Figure 2.9: Line B - shot B4 data. Each trace is scaled to a common maximum amplitude, See section 2.3.7 for a description of the phases labelled.  CO  to  Chapter 2.  DATA  ACQUISITION,  PROCESSING  AND  ANALYSIS  LO  Z\  6  9  (S)  8/a-i  £  0  33  Figure 2.11: Line C - shot C2 data. Each trace is scaled to a common maximum amplitude. See section 2.3.7 for a description of the phases labelled.  CO  Figure 2.12: Line D - shot D l data. Each trace is scaled to a common maximum amplitude. See section 2.3.7 for a description of the phases labelled.  CO  Chapter 2.  DATA ACQUISITION, PROCESSING AND ANALYSIS O  Z\  6  9  (S)  8/Q-l  £  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  2.3.2  37  Field Instrumentation  The 56 portable seismographs were of six different types, four digital and two analog, each with its own characteristic gain and frequency response. All recording units were calibrated to determine the full instrument frequency response to allow the conversion of each seismogram to velocity units of ground motion. The average frequency range of uniform response was 2.5-14 Hz, but varied from 1-3 Hz at the low end to 10-20 Hz at the high end. All data have been reduced to a common 60 Hz sample rate. Variations in analog tape speed were taken into account. Further details are described by Ellis et al. (1986).  2.3.3  Sedimentary Velocities  The interpretation of the refraction data is sensitive to the relatively low velocity and variable thickness of the 2 to 5 km of sediments covering the survey region (Figure 1.2). However, only within about 10 to 20 km of shot points can velocity information for the sediments be extracted from the refraction data. Therefore, a set of 30 industry sonic logs from wells close to the lines, 26 reaching the Precambrian basement, were obtained and analyzed to provide velocity control within the sedimentary section (Figure 2.1). Although sonic velocities are relatively high-frequency measurements (kHz), other studies have observed good agreement with seismically derived velocities (e.g., Fuis et al., 1984). The sonic logs were smoothed and fit with straight lines because of the requirements of the ray tracing algorithm. This is justified since the modelling of sub-basement arrivals is mainly sensitive to the travel time of near-vertical propagation through the sediments which was preserved in the sediment velocity models. The sonic logs were digitized at 100 m depth intervals and smoothed using a three-point filter.  The two best-fitting  Chapter  2. DATA  ACQUISITION,  PROCESSING AND  ANALYSIS  38  least-squares straight lines were then determined using a least-squares criterion so that a two-layer, linear velocity-gradient velocity model was obtained for each log (Figure 2.14). The sonic logs were smoothed prior to fitting straight lines to reduce the effects of the rapidly fluctuating sonic velocities. The sonic log models were interpolated to form continuous two-dimensional sediment profiles for each of the four refraction lines using the parameterization described in chapter 3 (Figure 2.15).  The four models in  Figure 2.15 represent the upper two layers of the four velocity models through which the ray-trace modelling of the refraction data was performed to obtain crust and upper mantle structure. The velocity structure in the near surface could not be obtained from the sonic logs (Figure 2.14) and so was extrapolated using the straight-line fit to the data below. The four sediment models contain strong lateral variations in both thickness and -1  velocity. High vertical velocity gradients (> 1 s ) in the first layer are common and in some regions the second layer contains a negative velocity gradient. Small adjustments to these models was necessary near a few shots to match the observed first arrival times of near-offset traces. 2.3.4  Topographic Corrections  Topographic corrections were made to avoid the inconvenience of incorporating the surface topography into the ray tracing models. A datum of 0.65 km above sea level was used, the mean elevation of the survey region. A near-surface velocity of 2 km/s and a ray-angle of 20° from the vertical were assumed in the correction, resulting in a time and offset adjustment of each trace. The near-surface velocity is a lower bound obtained by extrapolating the sonic velocities to the surface and a 20° incidence angle is appropriate for sub-basement turning and reflecting rays assuming a horizontally-layered model and applying Snell's law. Topographic corrections were not applied to the near-offset traces in which the first arrival energy propagated solely within the sedimentary section since  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  W  -20  39  E  i  30  i  80  i  130  DISTANCE  i  180  i  230  i  280  i  330  (KM)  Figure 2.14: Sediment velocity models from sonic logs along Line B. Distance is measured from shot BI at 0 km and sea level corresponds to 0 km depth. The solid lines directly above and beneath the sonic logs are the surface and sediment-basement boundary, respectively. Portions of two logs were interpolated from adjacent logs to provide uniform coverage from the near-surface to the basement, (a) The seven digitized sonic logs are shown as dashed and solid lines before and after smoothing, respectively. The location of each well is indicated along the surface, (b) The smoothed logs (solid) have been fit with the two best-fitting least-squares straight lines (dashed). Note that for two of the logs the best fit was achieved by a single line.  Chapter 2. DATA ACQUISITION,  -20  30  O-B-GO-IOW.'.  J  W -20  30  G-21-05-10W5  i  130  U  0-17-71-19W5.  80  L_l  180  230  *l T . 1i  1(1-32-70-20 W5  9-13-80-21W5  ,1-1G-85-21WS  DISTANCE (KM)  •-10-A 93 P/10 H-36-7B-15W6 16-28-79-11W6  o  AND ANALYSIS  DISTANCE (KM)  80  .L _ . - _ _ y ,  PROCESSING  13012-24-82-6W6 180 , 2-31-83-1W6 T  I  280  4.0  J T  230  ,  11-5-01-21W5  N 330380  l T  , 12-5-05-22W5  280  E 330  T  >—-C\J  n: h— n cn LU CD  b  X  LINE B  Figure 2.15: Isovelocity contour representation of the sediment velocity models; contours in km/s. Models were smoothed before contouring. The projected position onto the line and location (legal subdivision-section-township-range W meridian) of the wells from which sonic logs were used to obtain the models are indicated along the surface of each model. The surface of the models corresponds to a datum level of 0.65 km above sea level. The dashed line represents the sediment-basement boundary, (a) Line A. (b) Line B.  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  S 6-4-65-5W6 -20  d  | 30 6-17-71-5W6 80  11-22-67-5W6  DISTANCE (KM) 7-30-78-4W6 . 6-24-85-5W6 10-20-90-5WC 13012-24-82-0W6180 230  10-23-05-5W6 280  41  N 330  LINE D  Figure 2.15: Isovelocity contour representation of the sediment velocity models; contours in km/s. Models were smoothed before contouring. The projected position onto the line and location (legal subdivision-section-township-range W meridian) of the wells from which sonic logs were used to obtain the models are indicated along the surface of each model. The surface of the models corresponds to a datum level of 0.65 km above sea level. The dashed line represents the sediment-basement boundary, (c) Line C. (d) Line D.  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  42  the assumptions of the topographic correction are not valid at close offset—horizontal propagation independent of the datum level was assumed. Most corrections were small (~0.1 s and 0.1 km) since 85% of all traces were recorded at sites within 150 m of the datum level. 2.3.5  Timing Inconsistencies  Two major timing inconsistencies were discovered in the data set. First, a 0.2 s difference in first arrival time in the shot D l and D2 profiles at 295 km offset, the shot separation distance. The probable cause of a discrepancy this large is an error in either shot location (~1.6 km) or the timing of either shot. Ray tracing through the sedimentary velocity model obtained from a sonic log near shot D2 showed a difference between theoretical and observed first arrival times of roughly 0.2 s for the nearest offset trace. Thus, a static shift of +0.2 s was applied to all traces of the shot D2 profile. The second timing discrepancy concerns the Line A fan shot (shot A4) profile. Shots A4 and B l are approximately 200 m apart. A receiver at the intersection of Lines A and B recorded both shots (Figure 2.1) and the first arrivals of these two traces are of relatively high signal-to-noise ratio. There is a travel time difference of ~0.6 s between the first arrivals. Neither trace appears to have anomalous timing when compared with adjacent traces of each profile. After the modelling of all in-line profiles was complete, an inversion of the first arrival times of the Line A fan shot profile was performed to determine the upper mantle velocity in the centre of the survey region (section 4.7.2). A time advance of ~0.6 s of the fan shot profile was necessary to match the upper mantle velocities as obtained from the in-line modelling. Therefore, a 0.6 s time shift was applied to the Line A fan shot traces.  43  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  2.3.6  Amplitude Inconsistencies  The recording success rate was 97% with ~90% of these data containing reliable amplitude information where reliable is a qualitative assessment made on the basis of examining the record sections in true relative amplitude format and identifying traces with clearly anomalous low or high amplitude. To extract amplitude information from the traces with amplitude inconsistencies, an amplitude scale factor for each of these traces was calculated so that when applied the total trace energy followed the smoothed energy-distance trend of the section (Figure 2.16). To generate composite profiles in true relative amplitude format consisting of two shots of differing charge size (Lines A and B, see section 2.1), a scale factor proportional to the total charge weight raised to the 0.72 power was applied to the trace amplitudes (Michel, 1986). This power law was derived empirically by examining the amplitudes of traces recorded along Line C at the same site but for differing charge sizes from the same shot point; the smaller shots were misfires. 2.3.7  Arrival Picks  Travel time picks were made by hand from all record sections. A semi-automated picking routine (Zelt et al., 1987) was applied to a portion of the Line A-shot A3 record section as an objective check of the hand picks. Numerous record section formats were used for making the arrival picks, including various reducing velocities, small and large scale plots, unfiltered or a range of bandpasses, trace normalized or true relative amplitude, and time varying gain. The arrivals picked include all first arrivals as well as numerous later arrivals that provide additional constraints on the velocity models. Thefirstarrivals correspond to refracted energy turning within the sedimentary column (P ), the upper 3  crust (P ) and the upper mantle (P ) (Figure 2.17). The later arrivals are interpreted g  n  to be intracrustal reflections (R\,R2, etc.), a free-surface multiple generated within the  DATA ACQUISITION, PROCESSING AND ANALYSIS  E  0  44  W  50  100  150  200  250  300  DISTANCE (KM)  Figure 2.16: Total trace energy versus source-receiver distance before (solid) and after (dashed) smoothing. The dots identify the position of traces whose amplitude was adjusted so that its total energy matched the smooth trend of the section, (a) Line B shot B3. (b) Line C - shot C l .  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  0  p  50  .*"  i  DISTANCE (KM)  100  150  200  i  i  45  250  •  300  •  o CD  Figure 2.17: Ray paths of observed P-wave phases. The model is a generalized one-dimensional crust and upper mantle. For purposes of exposition a single reflecting horizon within both the crust and upper mantle is shown. A representative ray path for each of the nine phases is labelled. P„, P , P and P are rays which turn within the sediment, upper crust, mid or lower crust, and upper mantle, respectively. Ri, P P and P' are rays which reflect off an intracrustal boundary, the crust-mantle boundary, and an upper mantle boundary, respectively. The value of i may vary between one and five depending on the number of intracrustal reflectors identified in a particular record section, one corresponding to the shallowest, two the next shallowest, and so on. P propagates at the surface sedimentary velocity and is represented as a high-order surface-reflected multiple of the P phase. Mi is a first-order peg-leg multiple of the P phase with the peg-leg occurring beneath the receiver as shown, or beneath the source. g  c  n  m  w  e  g  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  46  sedimentary section (Mi), mid and lower crustal refracted arrivals (P ), and wide-angle c  reflections off the crust-mantle boundary (P P) and a boundary within the upper mantle m  (P~) (Figure 2.17). Intracrustal reflected phases with the same label on different record sections do not necessarily refer to reflections from the same boundary within the crust. Not all phases have been picked on all record sections due to variable data quality or the absence of a particular phase. A high-amplitude, low apparent velocity branch is present on most record sections and is interpreted to be the whispering gallery phase, P —a w  combination of the direct wave  and high-order free-surface multiples (Hill et al., 1982). This phase propagates directly beneath the free surface and may be used to obtain the near-surface sedimentary velocity and Q structure (Hwang and Mooney, 1986). The P phase will not be considered since w  it cannot be modelled using ray theory and the focus of this thesis is on sub-basement velocity structure.  2.3.8  Filtering  To aid arrival picking, the data were filtered with a zero-phase, eight-pole Butterworth filter (Kanasewich, 1981, pp. 237-280) to remove noise outside the seismic frequency band. For presentation it is convenient to apply the same filter characteristics over an entire section; 2-12 Hz and 1-8 Hz bandpasses, were selected for the in-line and fan shot profiles, respectively. Figure 2.18 shows the shot C l profile before and after filtering as well as the power spectra of each trace. This profile has the lowest signal-to-ratio of the PRASE data set. Figure 2.19 shows the filtering of the shot D2 profile.  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS w  Line C - shot CI  47 E  CO  a  DISTANCE (KM) Figure 2.18: Filtering of the shot CI profile, (a) The unfiltered shot CI profile in trace normalized format. The shot point is at 0 km. (b) The power spectrum of each trace calculated over the 12 s window in (a) as a function of distance with power increasing to the right. The power spectra have been scaled to a common maximum amplitude. The dashed lines show the 2-12 Hz bandpass applied, (c) The filtered shot C l profile, trace normalized after filtering.  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  48  Line D - shot D2  330  T 80  130  180  DISTANCE (KM!  330  Figure 2.19: Filtering of the shot D2 profile, (a) The unfiltered shot D2 profile in trace normalized format. The shot point is at 295 km. (b) The power spectrum of each trace calculated over the 12 s window in (a) as a function of distance with power increasing to the right. The power spectra have been scaled to a common maximum amplitude. The dashed lines show the 2-12 Hz bandpass applied, (c) The filtered shot D2 profile, trace normalized after filtering.  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  2.4  49  Method of Interpretation  All in-line profiles were interpreted using two-dimensional ray-trace forward modelling of both travel times and amplitudes of the first and later arrivals discussed earlier. The ray tracing algorithm is based on asymptotic ray theory (Cerveny et ai., 1977) and is described in chapter 3. Ray theory is inaccurate or invalid in regions about caustics, critical points, and transitions between strongly and weakly illuminated zones. Although the PRASE data set appears to possess reliable amplitude information, it generally exhibits large amplitude fluctuations within phases. Probable causes of these amplitude fluctuations include variations in near-surface geology at receivers, seismometer coupling problems and small-scale crustal heterogeneities not resolvable with the refraction data (e.g., Ojo and Mereu, 1986). Therefore, only the general trends in the amplitude character of major phases and the relative amplitude of phases are matched in the modelling. In these circumstances ray theory is sufficiently accurate. 2.4.1  Resolution  A more fundamental limitation than the use of ray theory is that associated with the refraction experiment geometry. The use of large source-receiver offsets results in lateral averaging of structure between source and receiver. Furthermore, the receiver spacing and number of shots used in this experiment prevents the resolution of small-scale crustal structure.  In addition, the frequency content of the data must be considered when  observations are compared with infinite-frequency ray-theory calculations and in assessing model resolution^ Resolution as referred to here is defined as the length scale over which features of the sampled medium are averaged. The dominant frequency of the data set is about 5 Hz and the receiver spacing varies between 3 and 5 km. Therefore, the minimum vertical resolution is a few hundred metres (one-quarter of the observed wavelength) and  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  50  the lateral resolution is at best 5 to 10 km. To improve upon the often simplistic and poorly resolved models obtained from refraction analysis, the quality and geometry of this data set has permitted the following three aspects of analysis: (1) the consideration of later phases, (2) amplitude modelling, and (3) a tie-in of each line's model at the intersection points with two other models. 2.4.2  Modelling Procedure  The modelling procedure used for each line proceeded in the following steps. First, sonic logs from wells close to the line were analyzed to obtain a sedimentary velocity model as described in section 2.3.3. The sub-basement crustal and upper mantle starting velocity model was obtained by interpolating the best-fitting laterally homogeneous models that satisfied the travel times of the profiles from both end shots and, where available, the centre shot. Through a trial-and-error scheme, this model was adjusted to match the travel times of arrivals from all shots simultaneously from progressively deeper layers within the crust and upper mantle. A comparison of the amplitudes obtained from the best-fit travel time model with the observed record sections showed that, in most cases, only a few minor adjustments were necessary to fit the observed amplitude trends. The velocity gradients of layers were adjusted to better match turning ray amplitudes and velocity contrasts across layer boundaries were adjusted to improve reflection amplitudes. Finally, the tie-in at intersection points with previously modelled lines was checked and the necessary adjustments were made. Most of the steps in the above sequence were repeated several times for further model refinement. Ray diagrams and comparisons between observed and theoretical travel times and amplitudes will be presented in chapter 4. The ray diagrams indicate the regions constrained by the phases modelled. The travel time picks of first arrivals generally have been fit to within ±0.1 s. This tolerance was chosen by taking into account both the  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  51  uncertainties associated with the picks and the rapid fluctuations observed in the arrival times along branches. Later arrivals were fit with less accuracy since the picking of these events, even when of high amplitude, is less certain. The constraint provided by the amplitude modelling of later arrivals is often as important as the matching of the corresponding travel times, particularly for the P and P P phases. For the intracrustal c  m  reflections, the observed amplitudes were generally not matched for the reasons discussed in section 4.2.3 so that only the travel times were used to constrain the positions of the reflectors. The requirement of model tie-in at the four line intersection points (Figure 2.1) was valuable in providing additional constraints on the models, particularly in the central portion of the Line C and D models which lacked centre shots. The agreement in crustal velocities at the four intersection points is generally within 0.15 km/s. The depths to the crust-mantle boundary agree to within 1 km and the upper mantle velocity to within 0.1 km/s. This level of agreement is considered satisfactory in view of the uncertainties discussed in section 2.4.3 and the possible presence of anisotropy. 2.4.3  Parameter Uncertainties  Since the models presented in this thesis are derived through forward modelling, it is not possible to quantify the uncertainties associated with the model velocities and boundary positions. Perhaps the most significant uncertainties are introduced by the interpretive process of phase identification. An estimate of the uncertainty associated with a parT  ticular parameter can be obtained by varying its value until a satisfactory match w ith travel times (±0.1 s) or amplitudes is lost. Using this approach, velocities and layer boundary positions have an uncertainty of roughly 0.05 km/s and 0.5 km, respectively. However, this is probably an underestimate since the simultaneous adjustment of other  Chapter 2. DATA ACQUISITION, PROCESSING AND ANALYSIS  52  model parameters will allow a particular parameter to be varied outside the uncertainties listed above and the model may still match the observed data. More appropriate model uncertainties for velocities and layer boundary positions are 0.1 km/s and 1 km, respectively, in the upper crust and 0.2 km/s and 2 km in the mid and lower crust. The additional constraint offered by intersecting lines may reduce these uncertainties.  Chapter 3  R A Y T R A C I N G IN T W O - D I M E N S I O N A L  3.1  MEDIA  Introduction  The interpretation of crustal seismic refraction data is generally carried out using a trialand-error forward modelling approach based on two-dimensional ray tracing  (e.g.,  Spence  et al, 1985; Mereu et al., 1986; Catchings and Mooney, 1988). The theoretical travel time and amplitude response of a laterally inhomogeneous medium is repeatedly compared with observed record sections until a model is constructed which provides a satisfactory match between calculation and observation. This iterative, forward modelling approach is necessary because no inversion scheme exists which is capable of providing geologically reasonable models in laterally varying media for a typical crustal refraction data set. Thus, an appropriate algorithm for forward modelling must be capable of accurately calculating the travel times and amplitudes associated with a general two-dimensional model in an efficient and practical manner to allow many trials. A number of algorithms based on asymptotic ray theory (Cerveny et al., 1977) have been developed and widely used, each satisfying the above mentioned requirements to varying degrees. These algorithms include those of McMechan and Mooney (1980), Cerveny and Psencfk (1981), Cassell (1982) and Spence et al. (1984). Numerous more accurate, but less efficient techniques for forward modelling in one-, two-, and threedimensional media exist (e.g., reflectivity and finite difference) and are necessary for  53  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  54  detailed studies of particular wave propagation problems and as a final check on models obtained using the ray method; however, routines based on the ray method are at present the only economical and practical choice for the interpreting seismologist. This chapter presents a routine based on the ray method which possesses similar advantages and limitations of other algorithms, but offers greater practicality. The fundamental feature of the new routine is its simple, layered, large-block velocity model parameterization. This allows for efficient ray tracing and since a minimum number of parameters are necessary to completely define the model, changes to it can be made quickly by the user. Rapid forward modelling is achieved by allowing the user to select simple kinematic ray families for exploring a part of the model of interest. In addition, a number of features have been incorporated into the routine to expand its flexibility and range of applicability. These features include S-waves, converted phases, multiples, approximate attenuation, head waves, a simulation of smooth layer boundaries, and a reverse ray-direction amplitude calculation suitable for the interpretation of commonreceiver profiles.  3.2  Velocity M o d e l Parameterization  Ray tracing algorithms are most strongly characterized by their velocity model parameterization. The accuracy, practicality and efficiency of a particular routine is ultimately linked to the parameterization, with some compromise between these three key attributes. The velocity model of the routine presented here can be described as a simple, layered, large-block, two-dimensional parameterization. In addition, the medium is assumed to be isotropic and laterally homogeneous in a direction normal to the plane of the model. The model is composed of a sequence of quasi-horizontal layers, the layers separated by boundaries which must cross the model from left to right without crossing another  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  55  X (DISTANCE) i  i i i  I  v,  1—  Q_  UJ  •v (x .z ) 0  v  _  r\  0  0  2  z=m x+b2 2  M  Figure 3.1: An example velocity model consisting of 5 layers and 27 trapezoidal blocks. The velocity, v , at the point (x ,z ) is given by equation (3.1). 0  0  0  boundary. Each boundary is defined by an arbitrary number of points connected by straight line segments of arbitrary dip. Within each layer, the P-wave velocity structure is defined by specifying a single upper and lower layer velocity for each straight line segment of the upper layer boundary. The upper and lower velocity pair may vary laterally within the layer (e.g., Figure 3.10a). Vertical boundaries separate each layer into large blocks and the emplacement of these boundaries is performed automatically by the routine as a requirement of the model parameterization. Vertical boundaries are necessary at each of the following points within a layer: (1) changes in slope of the line segments of the upper layer boundary, (2) changes in slope of the line segments of the lower layer boundary, and (3) changes in the upper or lower velocity within the layer. As a result of the vertical boundaries, each layer is broken up laterally into a series of large trapezoidal blocks with vertical left and right sides and upper and lower boundaries of arbitrary dip (Figure 3.1).  Chapter  RAY  3.  TRACING  IN  TWO-DIMENSIONAL  MEDIA  56  The velocity at a point within a trapezoid is determined by linear interpolation between the upper and lower velocity along a vertical path. For a trapezoid bounded above and below by line segments whose equations in the x-z plane are +  z = mix  the P-wave velocity, v  0  =  {(vim  -  2  VQ,  v m )xz 2  x  at the point + (v  2  -  and  bx  x  0  0  2  -+-  b  2l  within the trapezoid is given by  (x ,zo)  v )z  z = mx  + (v b x  2  -  v bi)]/[(m 2  2  -  m^xo  + [b  2  -  &i)]  (3.1)  where V\ and v are the upper and lower velocities, respectively, within the trapezoid (Fig2  ure 3.1). The S'-wave velocity is obtained from the P-wave velocity given in equation (3.1) using an assigned Poisson's ratio which may vary throughout the model, but is constant within each trapezoidal block. There are a number of advantages to having a large-block, trapezoidal model parameterization with constant velocities along the upper and lower boundaries: (1) a wedge-shaped structure with a vertical and horizontal velocity gradient can be modelled, (2) velocity discontinuities across layer boundaries can be readily incorporated into the model, (3) layer boundaries can be modelled as horizontal, dipping, or with considerable topographic relief through an appropriate combination of straight line segments, (4) layers whose boundaries have considerable topographic structure can be assigned a single constant velocity along the full extent of the upper and/or lower layer boundary, (5) a trapezoid can be of sufficiently general shape so that a minimum number are required to represent typical earth models, and, therefore, the trapezoids will, on average, be large blocks within the model, and (6) the layered, large-block character of the model results in a non-uniform but consistent grid so that the position of a single boundary point can be adjusted and the velocity model remains consistent. The algorithm allows for layers to shrink to zero thickness so that pinch-outs or isolated bodies can be considered. However, within a trapezoid adjacent to the pinch-out  Chapter  3.  RAY  TRACING  IN  TWO-DIMENSIONAL  57  MEDIA  at the point where the layer thickness shrinks to zero, the velocity will be undefined if the upper and lower velocities within the trapezoid are unequal. To avoid an undefined velocity, an equal upper and lower velocity within the adjacent trapezoid is specified.  3.3  R a y Tracing  3.3.1  R a y Tracing Equations and T h e i r Solution  To trace rays through the velocity model, the ray tracing equations are solved numerically (cf., Cerveny et al., 1977; McMechan and Mooney, 1980). The two-dimensional ray tracing equations solved by the routine are a pair of first-order ordinary differential equations; there are two sets: dz „ — = cotanfl.  dQ (v, — i>,cotan#) — = ^ '-  dx  dx  v  (3.2a)  anc ^  =  t  j  m  j  d8  ?  dz  (v tan6-v ) z  =  dz  x  v  with initial conditions x = x, Q  Z  =z, 0  9 =  8  0  (Cerveny et al., 1977, equations 3.19 and 3.19'). The variable 6 is the angle between the tangent to the ray and the z-axis. v is velocity, and v and v are partial derivatives of x  velocity with respect to  x  and  z,  respectively. The point  z  (x ,z ) 0  0  is the source location  and 6 is the ray take-off angle. To ensure stability, system (3.2a) is solved with x as the 0  integration variable when the ray path is near horizontal, and system (3.2b) is solved with z as the integration variable when the ray path is near vertical. To solve either system, the routine uses a Runge-Kutta method (Sheriff and Geldart, 1983, pp. 157-158) with error control as suggested by Cerveny et al. (1977). To complete the ray tracing algorithm, Snell's law must be satisfied at each point of intersection of a ray with a model boundary.  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  58  The ray step length, A, used in solving system (3.2a) or (3.2b) is an increment in either the x or z direction, respectively, and is given by the relationship  \v \ + \v \ x  z  '  where a is a user specified constant which will be discussed further with reference to the examples. An examination of equation (3.3) reveals the following: if the local derivative of the velocity field is large, the corresponding step length will be small, and vice versa. This approach is based on the realization that strong ray bending will occur if the ray is under the influence of a strong velocity gradient, whereas ray paths are essentially straight lines if the local velocity gradient is near zero. Thus, equation (3.3) is used to adjust the ray step length at each point to maximize efficiency while maintaining accuracy by avoiding unnecessarily small step lengths when the ray bending is small and the potential inaccuracy due to the use of a large step length when the ray bending is large. If the velocity within a trapezoid is constant, no ray bending occurs and the ray path is calculated analytically and the numerical solution of system (3.2a) or (3.2b) is avoided. Since the velocity, given by equation (3.1), and its partial derivatives are analytic functions of position, systems (3.2a) and (3.2b) can, in conjunction with equation (3.3), be solved efficiently. Once a ray has been traced through the model, it is defined by a series of points, the spacing and number of points being dependent on the value of a used in equation (3.3). The travel time at the ray end point is then evaluated by simple numerical integration along the ray path using the trapezoidal rule.  3.3.2  A u t o m a t i c Determination of Ray Families  Rapid forward modelling is accomplished through the routine's ability to search for the take-off angles of specific kinematic ray families. For a specific model layer, three types  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  59  of ray families can be searched for: (1) rays which turn (refract) within the layer, (2) rays which reflect off the bottom of the layer, and (3) the ray which generates head waves along the bottom of the layer (Figure 3.2). For a turning ray family, the search mode seeks the take-off angles of the shallowest and deepest ray to turn within the layer. Once the take-off angles of these two rays have been determined within a pre-specified distance, or a maximum allowable number of rays have been traced in the search mode, the entire ray family consisting of a specified number of rays with take-off angles between the shallowand deep-ray take-off angles is traced.  For a reflected ray family, the search mode  seeks the take-off angle of the ray with the smallest take-off angle (measured from the horizontal) which reflects off the bottom of the specified model layer. Once this take-off angle is appropriately determined, the complete family consisting of a specified number of rays with take-off angles between the smallest reflected take-off angle and a pre-specified maximum take-off angle, typically near 90°, is traced. To generate head waves along a layer boundary, the search mode seeks the ray which intersects the bottom of the layer at the critical angle, within a pre-specified tolerance (cf., Whittall and Clowes, 1979). To trace the complete head wave family, rays emerging at the critical angle are traced upward from the layer boundary and at a spacing specified by the user. In an analogous manner, a general diffracted ray family, in which the rays traced upward from the layer boundary emerge at the angle of incidence of the generating ray, may be traced. This family may be used, for example, to model the kinematics of diffractions which travel along the top of a low-velocity layer and occur as first arrivals within the geometrical shadow. The generating ray in this case is the first ray of the reflected ray family in the layer above the low-velocity zone. The manner in which the search mode operates in its search for any one of the four basic ray types is similar, so only one will be considered in detail. Consider the case of determining the take-off angle of the shallowest ray to turn within a specified layer.  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL  60  MEDIA  D I S T A N C E (KM) 0  60  120  180  240  300  0  60  120  180  240  300  0  60  120  180  240  300  LT)  LD  Figure 3.2: Example of the three simplest kinematic ray families for a single layer (layer 4) of a 5-layer model, (a) Turning rays (ray code 4.1); (b) reflected rays (ray code 4.2); (c) head waves (ray code 4.3).  61  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  If the medium were laterally homogeneous, this take-off angle, <£ , measured from the 0  horizontal in degrees, would be given by <t> = 90° - sin-Vo/u)  (3.4)  0  where VQ is the velocity at the source and v is the velocity at the top of the specified layer. This take-off angle is used as a starting value in the search mode where v is the velocity at the top of the layer directly beneath the source. If the value of v is not the maximum velocity between the source and the top of the layer along a vertical path directly beneath the source, the value of v in equation (3.4) is replaced by this maximum velocity. If all velocities between the source and the top of the layer are less than or equal to the velocity at the source, a pre-specified minimum take-off angle is used as the initial take-off angle in the search mode. An initial ray with a take-off angle cp is traced to its turning or reflection point. A 0  second ray is traced with a take-off angle of layer specified, or  <f> — 8(p 0  fraction of the quantity  cj> + 8<f> 0  if the initial ray did not enter the  if the ray entered the layer. The quantity  \<f>' — cp \ Q  0  where  <f>'  0  8cp  is a pre-specified  is the initial take-off angle used in the search  for the deepest turning ray in the layer. A third ray will be traced with a take-off angle of  cp ± 28<p 0  if both of the first two rays either did not enter the layer (-f) or did enter  the layer ( —), or (p ± \8(j) if the first ray did not enter the layer, but the second did 0  (+) or vice versa ( —). Thus, when two ray paths bracket the upper layer boundary, a bisection of angles begins, but until this point the take-off angles are incremented by 8</> (Figure 3.3). Within the search mode, rays are generally traced only to their turning or reflection point. However, for models with strong lateral variations, it may be necessary to trace rays to their completion to correctly determine the appropriate take-off angles; a switch allows the user to select this option.  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  62  Figure 3.3: The first seven rays traced in a hypothetical search for the shallowest ray to turn in the lower layer. The bold horizontal line is the model boundary separating the upper and lower layer. The maximum number of rays traced in the search mode is specified by the user, but will be terminated before this number if the distance between the end points of two successive rays is less than a pre-specified distance. Typically, 10 to 20 rays are sufficient to accurately define a specific take-off angle; however, for models with strong lateral variations, more rays may be necessary. The search routine has been found through experimentation to work successfully and efficiently for models with strong lateral variations. Take-off angles can also be supplied manually by the user.  3.3.3  Multiple and Converted Phases and Source and Receiver Locations  The three basic kinematic ray families obtained by the search mode can each be modified through simple numerical codes and switches to include any number or combination of multiple reflections, free surface reflections, and P-S or S-P conversions at any layer boundary. There is no restriction on the location of shot points within the model. Rays  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  63  from multiple shots may be traced in a single run so that the corresponding fit of the theoretical to observed data can be monitored simultaneously for all shot points along a particular line. For shot points located below the model surface, rays may also be traced upwards. Receiver locations are assumed to be at the top of the model. 3.3.4  Shooting M e t h o d Versus Two-point R a y Tracing  The search routine is based on the determination of ray families composed of kinematically similar rays, as opposed to two-point ray tracing in which rays are traced to specific receiver locations (cf., Cassell, 1982). A ray family search is a superior approach when forward modelling since a particular region of interest within the model can be examined through specification of the appropriate ray families. Therefore, when constructing a velocity model, the interpreter can begin model construction at the earth's surface and proceed downwards, matching arrivals from progressively deeper layers within the earth. The travel time (and amplitude) associated with an arbitrary point on the earth's surface is determined by linearly interpolating across the endpoints of the two closest rays which bracket the point of interest. For models with strong lateral variations, it may be necessary to trace relatively more rays to obtain accurate interpolated travel times and amplitudes. 3.3.5  Smooth Layer Boundary  Simulation  A simulation of smooth layer boundaries is possible whereby each boundary, consisting of a series of straight line segments, is uniformly sampled at many points (typically a few hundred) and smoothed through the application of a running-mean filter a specified number of times. The ray tracing inside each trapezoidal block remains unaffected; however, the angle of incidence and reflection or transmission at layer boundaries is calculated assuming the local slope of the smooth boundary. This simulation of smooth  Chapter 3.  RAY  TRACING  IN  TWO-DIMENSIONAL  MEDIA  64  or curved boundaries is simple and efficient to apply and can be important for rays paths which cross layer boundaries near a point at which the boundary consists of two line segments of differing slope. The accuracy of ray tracing in this situation will depend upon how the position of the original and smooth boundary differ—this can be checked visually by the user and differences reduced by increasing the number of points defining the layer boundary. Travel time curves will be smoothed in accordance with the amount of model smoothing and geometrical shadows caused by model blockiness will be reduced. An example of this effect is shown in Figure 3.4. The additional CPU time required for the smooth-layer ray tracing was ~10% for this example.  3.4  Amplitude  3.4.1  Calculations  Turning and Reflected Rays  The complex amplitude of turning and reflected rays, possibly multiply reflected and/or converted, is calculated by zero-order asymptotic ray theory (Cerveny et al. 1977, equation 2.56) as given by the expression A = A q(-iy/L  (3.5)  0  where A is the amplitude of the ray at its endpoint, A is the initial ray amplitude, 0  q is a factor which accounts for the energy partitioning at model boundaries, and L is the geometrical spreading. The factor e is equal to the number of times the direction of positive displacement changes along SV segments of the ray path with respect to the ray coordinate system (Cerveny and Ravindra, 1971, pp. 72-73). A point source with uniform directional characteristics is assumed. The initial ray amplitude is set equal to unity so that the amplitude, A , is a relative amplitude; however, the ratio of SV- to P-wave energy generated by the source can be  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  65  Figure 3.4: The effect on ray tracing and travel time curves of the smooth layer boundary simulation, (a) Rays traced through velocity model without smooth layer simulation. The velocity in km/s is indicated in each layer. The upper dashed layer boundary is horizontal and the lower dashed boundary is comprised of 10 straight-line segments. The maximum change in slope of adjacent line segments is 30° at 3.25 km model distance, (b) Travel time curves corresponding to rays traced in (a), (c) Rays traced through velocity model with smooth layer simulation. For the purposes of ray tracing, the lower dashed boundary was uniformly sampled at 250 points and smoothed with 100 iterations of a 3-point running-mean filter. The maximum deviation between the dashed boundary shown and the smooth boundary is 30 m at 3.25 km. (d) Travel times curves corresponding to rays traced in (c).  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  66  specified. The quantity q is equal to (Cerveny et ai. 1977, equation 2.58)  (  ) 2  VoPo  9 =  (3.6)  Zi  VrPr  where n — number of model boundaries encountered by ray = velocity and density at the source whose product is impedance  Vo,Po  vp T  T  ViPi  —  impedance at the receiver  = impedance at the point of incidence at the z'th boundary  v'^ = impedance at the point of emergence at the ith boundary Zi — Zoeppritz amplitude displacement coefficient at the zth boundary. 3  The relationship between density, p (g/cm ), and P-wave velocity, v (km/s), is specp  ified as p = -0.6997 + 2.230v - 0.5980^ + 0.07036^ - 0.002831i£ p  This represents a 4th-order polynomial fit to a collection of velocity-density values presented by Ludwig et ai. (1970) and provides a reasonable fit for 1 < v < 9 (Figure 3.5). p  The coefficients of the polynomial may also be specified by the user. For a water layer. Poisson's ratio is specified as 0.5 and the density is automatically set to one. The complex Zoeppritz displacement amplitude coefficients for incident P- and SVwaves are calculated using a routine based on that described by Young and Braile (1976) which allows for the calculation of surface reflection coefficients, but has been modified here to allow for the calculation of surface conversion coefficients for incident P and S type rays of vertical and horizontal component type (Cerveny and Ravindra. 1971, equations 2.89 and 2.90).  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  1  P-WRVE V E L O C I T Y  67  (KM/S)  Figure 3.5: The default velocity-density relationship is a 4th-order polynomial fit to data presented by Ludwig et ai. (1970). The geometrical spreading is given by ( cos Qi ,=i V c o s ^ where J± is an out-of-plane spreading term, J|| is an in-plane spreading term, and the product term represents a compensation for the discontinuity in the spreading functions across model boundaries (Cerveny et ai. 1977, p. 38). The angles Qi and Q\ are taken at the point of incidence and emergence at the ray's intersection with the ith model boundary. The out-of-plane spreading term can be expressed as J  1  1J  =  v  0  vdl  (3.7)  J  where / is the ray path length (Cerveny, 1981, equation 31; Cerveny and Hron, 1980). The integral in equation (3.7) is evaluated using the trapezoidal rule. The in-plane spreading term is given by 8r  Jw = cos Q,  (3.8)  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  68  where 9 is the angle between the vertical and the ray path at the receiver, r is the distance T  between the source and receiver (range), and 0 is the angle between the vertical and O  the ray path at the source (take-off angle) (Cerveny et al. 1977, equation 3.69). The evaluation of the derivative in equation (3.8) can be made through an examination of the kinematics of a set of similar rays. The approach used is most like that used by Marks (1980). To evaluate the derivative, a separate cubic spline is fit to each prograde and retrograde segment of the curve defined by range versus take-off angle for all rays of the same kinematic family and the derivative for each ray in the family can be obtained directly from the splines. The advantage of this approach is that additional rays, beyond those already traced for the family, do not have to be traced, as is necessary with other routines. A limitation is that for models with strong lateral variations, the slope of the cubic spline may be sensitive to the number of rays which constitute the family. The greater the lateral inhomogeneity of the model, the greater the number of rays which must be traced to fully sample all aspects of the lateral variation of the model. For most models which have been studied, 10 to 20 rays per family were sufficient. 3.4.2  Reverse Ray-direction  Amplitude  For the interpretation of common-receiver profiles, often encountered in marine ocean bottom seismometer (OBS) studies, it is necessary to calculate amplitudes in the reverse direction to which rays are traced. A reverse ray-direction amplitude calculation is accomplished by switching the source and receiver quantities and incident and emergent quantities in equation (3.6) and adjusting the geometrical spreading by a factor dependent on the ratio of the source and receiver velocities (Richards, 1971, equation 18). These adjustments are made since the reciprocity of amplitudes is not necessarily preserved by the approximations of ray theory (Razavy and Lenoach, 1986). In addition, the surface conversion coefficient must be calculated at the ray source.  Chapter  3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  69  Figure 3.6 shows a comparison of amplitudes calculated in both directions for a onedimensional model of shallow crustal structure as determined by an OBS survey (Cudrak, 1988). The maximum difference in amplitudes is a factor of ~2 for the first arrival at 2 km offset after the reverse ray-direction profile was scaled up by a factor 4 to account for the absolute amplitude difference between sections. Arrivals beyond 6 km offset agree to with one percent. This example shows that the effect of the reverse ray-direction calculation is not significant unless an accurate near-offset amplitude analysis is performed.  3.4.3  H e a d Waves  In most practical situations, the modelling of head waves is not required as turning rays in the lower medium dominate. However, in some instances topographic structure on a layer boundary can result in a shadow zone with respect to turning rays in the lower medium. In this case, head waves will represent the first arrivals at certain receiver locations. The routine allows for the calculation of P and SV head wave amplitudes through first-order asymptotic ray theory assuming the layer along the top of which the wave propagates has zero velocity gradient. The routine requires that the ray paths toward and away from the head wave boundary do not contain reflection or conversion points. The calculation of head wave amplitudes is made in the same manner as that of Spence et ai. (1984); however, no explicit reparameterization of the velocity model into a series of thin homogeneous layers is required since the region between each point defining the ray path toward and away from the head wave boundary can be considered a homogeneous layer if the average velocity between points is used as the layer velocity. The head wave H  amplitude, A , is given by (Cerveny and Ravindra, 1971, equations 5.22 and 5.29) H  A =  vVk  tan 6'  k  (3.9)  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  70  DISTANCE (KM] 10  15  to •_ 00 \ Q  O  20  25  20  25  o  CO  DISTANCE  15 (KM:  Figure 3.6: Reverse ray-direction amplitude comparison, (a) Turning and reflected rays traced through one-dimensional velocity model of shallow oceanic crustal structure (Cudrak, 1988). The ray source at the bottom of the water layer is the location of the OBS which recorded airgun shots from just below the water surface. The upper and lower velocity in km/s is indicated in each layer, (b) Synthetic section for rays traced in (a) from left-to-right, (c) Synthetic section for rays traced in (a) from right-to-left (reverse direction). A vertical component seismometer and surface conversion coefRient was applied at the water surface for (b) and at the OBS position for (c). Amplitudes in both sections are scaled by a factor proportional to distance. The major discrepancies between (b) and (c) occur in the strong low-apparent-velocity branch at near offset which corresponds to turning rays in the first crustal layer.  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  where  and  (  COS  71  $i  cos dl  JV = number of line segments comprising the complete head wave ray path angle of incidence and emergence at the zth boundray Vi = velocity along the ith ray segment li = length of ith ray segment k = head wave boundary number v = velocity below head wave boundary I = length of ray path along head wave boundary UJ = dominant frequency of head waves Tfc = head wave coefficient (Cerveny and Ravindra, 1971, pp. 108-109). The Zoeppritz routine has been modified to include the computation of the head wave coefficient. If the velocity below the head wave boundary varies laterally, the value of v is set equal to the appropriate average value. At the critical point, the value of I equals zero so that the corresponding ray amplitude is undefined. To avoid this problem, an ad hoc approach used by Spence (1983) is incorporated which assures that the value of I used in equation (3.9) remains greater than a pre-specified minimum value. Therefore, the head wave amplitude within the interference zone with the reflected ray is not valid. The routine does not allow for the calculation of amplitudes for the general diffracted ray family described earlier; however, this would be possible if an angle-dependent diffraction coefficient were incorporated into a formula similar to equation (3.9).  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  3.4.4  72  Approximate Attenuation  The effect on amplitudes due to attenuation can be considered in an efficient manner by the routine through a major simplifying assumption: the attenuation is independent of frequency. The attenuation is calculated at a specified dominant frequency and this attenuation is applied equally at all frequencies (e.g., Chiang and Braile, 1984). This results in the incorporation of a single multiplicative scale factor into equations (3.5) and (3.9) so that the coupled phenomenon of dispersion is neglected. The scale factor,  AQ,  is given by (Aki and Richards, 1980, pp. 168-169) N  A  Q  = IJ exp[-u;*,/2Q,]  (3.10)  i= l  where N — number of line segments comprising the complete ray path UJ = dominant frequency ti = travel time along the ith ray segment Qi = Q-attenuation value along ith ray segment. P- and 5-wave Q values can be assigned to each model block. This approximate method is a good approximation to a more accurate frequencydependent calculation if the source energy is concentrated within a narrow frequency band or the attenuation is not large. Figure 3.7 is a comparison of an exact dispersive calculation with the approximate method for propagation through a homogeneous constant-Q medium. The exact calculation was made using equation (3.10) for each value of UJ (N = 1) and the dispersion relationship presented by Kanamori and Anderson (1977, equation (15)). Figure 3.7 shows that for an initial 10 Hz dominant-frequency waveform, and, for example, Q = 1000 and a propagation time of the reference frequency  Chapter 3. RAY TRACING  IN TWO-DIMENSIONAL  MEDIA  73  in  0.0  0.02  0.04  TIME/Q  0.06  0.08  0.1  (S)  TIME/0 (S) Figure 3.7: (a) Comparison of exact (solid) and approximate (dashed) waveforms as a function of the time of propagation through a homogeneous constant- Q medium divided by Q. The initial waveform was a low-passed Pucker wavelet whose power spectrum is shown in (b). The initial central frequency is 10 Hz. The exact waveform shows dispersion and was calculated using a frequency-dependent attenuation and a dispersion relationship (see text for description). The approximate waveform was calculated using the non-dispersive method described in the text, (c) The ratio of the exact to approximate peak amplitudes, (d) The time difference between the central peaks of the exact and approximate waveforms. The "choppiness" of the curves in (c) and (d) is an artifact of the sampling interval.  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  74  (1 Hz) of 40 s (TIME/Q = 0.04), the ratio of the approximate to exact amplitude is 0.92 and the time delay of the central peak of the exact waveform with respect to the approximate waveform is 0.055 s. 3.4.5  Synthetic Seismograms  Once the amplitude calculations for all rays which reached the model surface are complete, synthetic seismograms for a set of specified surface locations can be generated.  The  approach used is similar to that used by McMechan and Mooney (1980) and Spence et al. (1984). The time and complex amplitude of all arrivals on a particular seismogram are obtained by linearly interpolating across the end points of the two closest rays which bracket the seismogram location and are of the same kinematic ray type. A superposition of all arrivals on the seismogram results in a trace consisting of a series of impulses, each with appropriate amplitude and phase rotation. The 90° phase rotation associated with passage through a caustic can be applied automatically to retrograde turning ray branches or supplied manually by the user. The final seismogram is obtained through a convolution with an apparent source function. 3.5  Examples  The routine presented in this chapter is based on the ray method and as such the accuracy and limitations of its dynamic calculations are similar to all other ray tracing routines. For a thorough discussion of the limitations of the ray method and its comparison to other more accurate methods, the reader is referred to Cerveny (1985). The following examples are not intended to show the weaknesses of the ray method, but rather the positive aspects of the new routine as a practical forward modelling algorithm. All computations were performed on a 12 Mips Amdahl 5860 computer.  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  3.5.1  75  Efficiency and Accuracy  The first example illustrates the routine's efficiency and accuracy as compared to the algorithm of Spence et al. (1984) in which the model is parameterized in terms of large polygonal blocks within which the velocity gradient vector is constant. Since the Spence et ai. algorithm uses analytic expressions to trace rays and calculate travel times, these quantities may be considered to be exact with respect to ray theory. Therefore, a comparison of these values with the new routine's will reveal the inaccuracy associated with solving the ray tracing and travel time equations numerically. As a check on the routine's amplitude calculations, numerous models were tested and the results showed close agreement with those obtained from the Spence et al. algorithm which has been compared with the McMechan and Mooney (1980) routine and the reflectivity method of Fuchs and Muller (1971) by Spence et al. (1984). Due to the manner in'which the new routine and that of Spence et al. parameterize velocity models, only models in which the boundaries are horizontal or vertical can be compared. The model chosen for the comparison study is a generalized earth model consisting of three layers (Figure 3.8a). The three layers represent sediment, crust and upper mantle and the velocity structure within each layer has been broken laterally into six blocks. Two different numerical tests were made with this model: (1) a study of the tradeoff between CPU time and RMS travel time error, measured against the Spence et al. algorithm, as a function of the ray step length parameter (a), and (2) a comparison of the CPU time required for a number of different types of runs of the two routines. In all cases, 100 rays were traced through the model—25 turning rays in each layer and 25 rays reflected off the crust-mantle boundary. For some runs of the new routine, additional rays were traced in the search mode to demonstrate the extra CPU time required. Figure 3.8b shows the 100 rays traced through the model by the new routine and the corresponding  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  76  travel time curves and synthetic section are shown in Figures 3.8c and 3.8d. Figure 3.9 shows the tradeoff between accuracy and CPU time. The value of the ray step length parameter, a, used in equation (3.3) was varied between 0.4 and 0.0075. In all runs, amplitudes were not calculated and no plots were generated. Also, a maximum of 20 additional rays were traced in the search mode to define the take-off angles of each ray family. One might expect the CPU time to increase monotonically and the RMS error to decrease monotonically as the value of a decreases. This is the general trend observed; however, a minimum CPU time occurs at a = 0.15 and a minimum RMS error occurs at a = 0.015. The C P U time minimum is due to two factors. First, regardless of the value of a used, a constant single-step tolerance is used in the Runge-Kutta routine. Therefore, for relatively large step lengths, the number of sub-intervals required by the Runge-Kutta routine may become large. Second, the search mode does not operate effectively when too large a value of a is used. This results in relatively more rays traced in the search mode. The minimum in RMS travel time error is also due to two factors. The fixed single-step tolerance used in the Runge-Kutta routine does not allow for increased accuracy when relatively small step lengths are used. Also, the accumulated round-off error of this single precision routine prevents increased accuracy with decreased step length. Figure 3.9 shows that the CPU time increases by a factor of about five for a corresponding decrease in the value of a of about fifty. The variations in the RMS travel time error are perhaps more important and should be considered when selecting an appropriate value of a. Typically, crustal refraction data will have a timing uncertainty of between 10 and-20 ms. Therefore, a value of a less than about 0.15 should be sufficient in most cases (Figure 3.9). Also, values of a less than about 0.05 increase the C P U time required without significantly improving accuracy. Thus, the CPU time and RMS error minimums suggest that for this model, and probably all models that are normally  Figure 3.8: Example in which the efficiency and accuracy of the new routine is studied, (a) Velocity model consisting of three layers with six blocks in each layer. The upper and lower velocity in km/s is indicated within each model block, (b) One reflected and three turning ray families traced through the model, (c) Reduced travel time curves corresponding to the rays traced in (b). (d) Synthetic section corresponding to the rays traced in (b). The source wavelet is a single-cycle sine wave. The trace amplitudes have been scaled by a factor proportional to distance. —i  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  78  RMS T R A V E L TIME ERROR  50  CPU TIME  ti o ti ti w w 10 H  2  a  a  $ ti  O •O O o«  .01  .04  .1  .4  RAY STEP LENGTH PARAMETER (a)  Figure 3.9: Trade-off between accuracy and CPU time as a function of the ray step-length parameter, a , used in equation (3.3). considered, the value of a should be between 0.015 and 0.15. The number of points defining rays increases as the value of a decreases. For this example, the ray defined by the greatest number of points was, in each run, the deepest turning ray in the second layer. For a = 0.4, 0.1, 0.03 and 0.01, the number of points defining this ray was 13, 33, 94 and 267, respectively. A comparison of the CPU time required for the new routine with that of Spence et ai. was the second test performed with the model shown in Figure 3.8a. Table 3.1 summarizes the results of eleven runs. For all runs of the new routine, a was equal to 0.1. For the first five runs listed in Table 3.1, the take-off angles for the new routine were supplied manually. For the last six runs the take-off angles for the new routine were determined through the search mode with the maximum number of rays traced to define each take-off angle equal to 10, 25 or 50. The additional number of rays traced in each case was 55, 77 and 102, respectively. For all runs, the take-off angles for the Spence et ai. routine were supplied manually.  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  79  CPU time (s) Description CPU time (s) of runt (this routine) (Spence et ai.) 1 0.70 0.65 0.74 1.24 1,2 1,3 0.81 1.07 1,2,3 1.50 0.95 1,2,3,4 2.04 6.08 0.97 1.02 1.08 1,3,5 1.07 1,2,3,5 1.50 1.05 1.10 1.14 Tl - calculate ray paths and travel times 2 - calculate amplitudes 3 - generate plot of rays and travel time curves 4 - calculate synthetic section and generate plot 5 - search mode for new routine with maximum number of rays traced to define each take-off angle equal to 10, 25 and 50 Table 3.1: CPU time comparison with Spence et ai. routine for model and rays traced in Figure 3.8 on 12 MIPS Amdahl 5860 computer. Table 3.1 shows that the new routine required less CPU time for all runs except one in which the search mode was used and amplitudes not calculated. There are three reasons for the new routine's favourable comparison. First, the new routine operates in single precision whereas the Spence et ai. routine requires double precision because of the way in which travel times and the in-plane geometrical spreading are calculated. Second, the manner in which the next block is determined when a ray intersects a model boundary is generally more efficiently performed by the new routine due to its layered, grid-like model parameterization. Finally, to calculate the in-plane spreading, the routine of Spence et ai. traces a second ray for every ray traced in order to estimate the derivative in equation (3.8). The effect of this additional ray tracing is evident in the CPU times for the runs in which amplitudes are calculated. The relative efficiency of the new routine for the run in which the synthetic section was generated is due in part to the incorporation  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  80  of a time/cost saving feature in which zero amplitude portions of traces are plotted as two points at either end connected by a straight line. Even with the additional ray tracing in the search mode, the new routine generally required less CPU time. With the maximum number of rays traced to define each take-off angle equal to 10, the routine was able to determine the take-off angles of each of the ray families. Setting this maximum to 25 and 50 resulted in additional refining of the take-off angles. 3.5.2  P R A S E Line C—shot C 2  The second example is an application of the routine to the modelling of the PRASE data set. Figure 3.10a shows the final block velocity model for Line C in which profiles from shots at 0 and 270 km model distance were interpreted through a comparison of observed and theoretical travel times and amplitudes. Figure 3.10b is a smoothed isovelocity contour representation of the model and shows more clearly that the model possesses significant lateral variations. The model consists of 9 layers—2 sedimentary, 5 crustal and 2 upper mantle—and 51 blocks. Figure 3.11a shows the complete set of all rays traced through the model from the shot point at 270 km. Fourteen ray families consisting of 150 rays were traced. An additional 252 rays were traced in the search mode. The ray families include turning rays in each layer and reflections off the four intracrustal boundaries, the crust-mantle boundary, and an upper mantle boundary at a depth of about 56 km. The corresponding travel time curves are shown in Figure 3.11b compared with the observed travel time picks. The synthetic record section is shown in Figure 3.11c. With a value of a = 0.1, the C P U time required to calculate and plot all rays, travel time curves, and the synthetic section was 4.26 s. The corresponding observed record section is shown in Figure 3.lid. The agreement between the synthetic and observed data is good; however, the relatively high  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  81  Figure 3.10: Velocity model interpreted for PRASE Line C. (a) Block model representation in which the upper and lower velocity in km/s is indicated within each trapezoid. The surface velocity in the first sedimentary layer varies from 4.25 to 1.75 km/s from left to right and the lower velocity in the second sedimentary layer varies from 6.0 to 4.0 km/s form left to right, (b) A smoothed, isovelocity contour representation of the velocity model in (a). The velocities of contours are indicated in km/s.  Figure 3.11: Application of the new routine to the interpretation of PRASE Line C. (a) Rays traced through velocity model, (b) Comparison of the observed (symbols) and calculated (lines) travel times, (c) Synthetic section corresponding to the rays traced in (a). The source wavelet is a low-passed Ricker wavelet, (d) Observed record section bandpass filtered from 2-12 Hz with- the calculated travel time curves superimposed. The amplitudes in (c) and (d) have been scaled by a factor proportional to the source-receiver distance.  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  83  amplitudes of the intracrustal reflected phases have not been reproduced by the model. Strong reflections of this type could not be explained by large step increases in velocity suggesting that they may be due to a zone of thin, alternating high- and low-velocity layers (e.g., Sandmeier and Wenzel, 1986) caused by variations in composition, texture and/or anisotropy (Christensen and Szymanski, 1988). 3.5.3  C o m p l e x Crustal M o d e l  For the third example, a generalized crustal model of an oceanic-continental subduction zone is presented (Figure 3.12a). This model is a simplified version of that developed by Spence et ai. (1985) for the subduction zone of western Canada. The purpose of this example is to illustrate the amount of lateral inhomogeneity possible with the model parameterization of the new routine. The layer pinchout feature discussed earlier has allowed the obliquely dipping oceanic plate to be incorporated into the model. Five continental layers on the right side of the model have been pinched out above the subducting plate and four oceanic layers on the left side of the model have been pinched out at the bottom right side of the model. Figure 3.12b shows turning and reflected rays traced through the subducting oceanic upper mantle and the corresponding synthetic section is shown in Figure 3.12c.  3.5.4  Near-surface Velocity Anomalies and the Statics Problem  The final example is an application of the routine to the study of the effects of nearsurface velocity anomalies on common-midpoint (CMP) data, i.e., the statics problem. The routine was modified so that a collection of CMP gathers could be obtained in a single run. The model is presented in Figure 3.13a and consists of a horizontal reflector underlying a low-velocity layer of varying thickness. A common-shot gather and CMP gather for the same point in the model are shown in Figures 3.13b and 3.13c, respectively.  150  190  230  270  310  350  D I S T A N C E (KM) Figure 3.12: Application of the new routine to a model with strong lateral inhomogeneity. (a) Velocity model with only layer boundaries shown. The velocities in km/s are indicated within each layer, (b) Some turning and reflected rays traced through the model, (c) Synthetic section corresponding to the rays traced in (b). The source wavelet is a single-cycle sine wave. The two turning ray families are identified in (b) and (c) by the letters A and C and the two reflected families by B and D.  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  85  Figure 3.13: Application of the routine to the study of the effects of near-surface velocity anomalies on CMP data, (a) Reflected rays in second layer traced through model. Upper and lower velocities in km/s are indicated with each layer. Vertical block boundaries are not shown, (b) Common-shot gather at 3 km model distance for turning and reflected rays traced through first two layers of model. Each trace is scaled to a common maximum amplitude, (c) CMP gather at 3 km model distance. Only the reflected phase from the second layer is included, (d) CMP gather in (c) after NMO correction using the rms velocity to the second layer reflector (2.38 km/s). Distance in (c) and (d) is half offset and the corrected relative amplitude of each trace is presented.  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  86  The CMP gather after NMO correction is presented in Figure 3.13d. The misalignment of arrivals in Figure 3.13d is due primarily to the lateral variations in thickness of the upper layer, but is also due to the increase in velocity with depth within the model. 3.6  Discussion  This routine offers greater practicality than existing forward modelling ray tracing algorithms. In addition, the efficiency and range of model types which can be considered is comparable to current, widely-used routines. Fundamental to the new routine is its simple, layered, large-block velocity model parameterization.  The layered nature of the  model is valuable in that it allows the user to select specific layers within which, or layer boundaries at which, turning, reflection and conversion are to occur. The large-block parameterization is vital for efficiency since the velocity and its derivatives necessary for the solution of the ray tracing equations are simple analytic functions of position and the number of model boundaries crossed by each ray is minimized. Another important advantage of large blocks is that the number of zero- and first-order discontinuities in velocity within the model is minimized. This reduces the difficulties which can arise in the application of the ray method to "blocky" models, in particular, with regard to amplitude calculations (Cormier and Spudich, 1984, equations 4.1 to 4.3). The parameterization requires only a single upper and lower velocity to define the velocity field within each trapezoid and layer boundaries are a series of straight line segments. This allows the user to alter the model easily and quickly so that the new velocity structure is always well-defined and predictable. A number of relatively simple features incorporated into the algorithm, and which are often necessary in refraction interpretation, give the routine flexibility and a wide range of applicability. These features include layer pinchouts, S-waves, converted phases,  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  87  multiple and surface reflections, approximate attenuation, head waves, no restriction on source location, multiple sources in a single run, a simulation of smooth layer boundaries, and a reverse ray-direction amplitude calculation. The latter capability is necessary for the interpretation of marine data in which common-receiver profiles are involved. As a practical algorithm for forward modelling interpretation of crustal refraction data, the only significant limitation of the routine is that of the ray method itself. In particular, a blocky velocity structure can introduce false highs and lows into amplitude calculations with respect to finite frequency wave propagation due to discontinuities in velocity gradients at model boundaries (Cerveny, 1985, Figure 11). Velocity discontinuities also cause lows near the critical point where the reflected, refracted and head wave arrivals interfere. The lows are reduced somewhat if the head wave amplitude is included; however, this improvement is arbitrary since the head wave amplitude is invalid near the critical point. This problem is in addition to the fundamental inapplicability of the ray method in regions about caustics, critical points and geometrical shadows. The incorporation of Gaussian beam amplitudes into the routine would be a straightforward procedure and would remove the above mentioned difficulties associated with the dynamic properties of ray theory. However, this modification does not appear to be necessary or practical at this point given the amplitude control of most refraction data sets and the somewhat unpredictable nature of Gaussian beam seismograms (Miiller, 1984; Cerveny, 1985). To reduce the artificial highs and lows associated with ray amplitudes, allowance within the routine is made for a specified amount of smoothing of the amplitude curves for each ray family. In addition, layer boundaries may also be smoothed as described in section 3.3.5. The routine was developed for use in the interpretation of crustal refraction data; however, a number of other applications are possible—for instance, the CMP example presented earlier. The routine would also be suitable as a forward modelling algorithm  Chapter 3. RAY TRACING IN TWO-DIMENSIONAL MEDIA  88  for iterative inversion of refraction/wide-angle reflection data or migration of reflection data. The appropriate conversion of coding would be straightforward to allow for the routine's implementation on a microcomputer, as was done for the Spence et al. (1984) algorithm (Crossley, 1987). A possible alternate model parameterization could be implemented in which the velocity field within each trapezoid is defined by velocity values specified at each of the four corners of the block. This would allow the velocity structure within a layer to be free of velocity discontinuities at vertical block boundaries and the parameterization would require only a single additional upper and lower velocity pair to be specified for each layer. It may be possible to improve the efficiency of the routine by automatically breaking each trapezoid into two or more triangular blocks within which the velocity gradient is a constant (e.g., Chapman and Drummond, 1982; Miiller, 1984). However, this would at least double the number of blocks within the model and result in each ray crossing about twice as many model boundaries. This may reduce efficiency and degrade the amplitude calculations due to the increased number of discontinuities in the velocity gradient.  Chapter 4  INTERPRETATION OF REFRACTION  4.1  DATA  Introduction  This chapter presents (1) observations and aspects of the in-line interpretations common to the four lines, (2) features unique to the interpretation of each line, (3) the two fan shot interpretations, (4) an examination of the S-wave data, (5) an estimation of crustal and upper mantle Q structure, (6) a check of the seismic-gravity consistency, and (7) a summary of the regional velocity structure. For each of the ten in-line shots, the ray diagram, observed and theoretical travel times, and the true relative amplitude observed and synthetic section will be presented. The synthetic sections were calculated assuming the velocity-density relationship shown in Figure 3.5 and a Poisson's ratio of 0.25; the calculations are relatively insensitive to these parameters given the degree to which the observed and theoretical amplitudes are matched. To compensate for anelastic attenuation, amplitudes were calculated assuming crustal and upper mantle Q values determined from a study of the observed spectra of the PRASE data set as presented in section 4.9. The incorporation of attenuation in the amplitude calculations is important to obtain more reliable estimates of vertical velocity gradients. The source wavelet used in the synthetic sections is a low-passed, time-shifted (causal) Ricker wavelet with 4.7 Hz dominant frequency. This wavelet resembles the P  g  waveform on a number of high signal-to-noise sections. Finally, layer boundaries and amplitude-distance curves of each ray family were smoothed as described in chapter 3,  89  Chapter 4. INTERPRETATION OF REFRACTION DATA  90  although the effect of this smoothing was minimal since the refraction models do not contain strong lateral variations. The four profile velocity models are presented in cross-sectional P-wave isovelocity contour format. The models were smoothed prior to contouring to remove the " block" nature of the model parameterization required by the ray tracing routine. This presents the models in a way that is more in keeping with the resolution of the refraction data. The positions of all reflecting boundaries constrained by travel time picks are indicated as well as the intersection points with other lines. In addition, velocity-depth profiles taken from within the models at points 50 km from either end shot are included on either side of the velocity model diagrams. The letters N, S, E and W at model ends indicate approximate line orientation. The top of each model at 0 km depth corresponds to the datum level of 0.65 km above sea level. Unless otherwise stated, "distance" refers to model distance.  4.2  4.2.1  General Results of the In-line Interpretations  Simple, One-dimensional Crustal Structure  A preliminary inspection of the ten observed record sections (Figures 2.2 to 2.13) reveals an important result of the complete seismic interpretation: the PRA region has a simple near-laterally-homogeneous crustal structure. The structure is simple in the sense that there are no strong breakovers or offsets in first arrivals with associated reflected branches suggesting significant layering or thick low-velocity zones within the crust. Lateral homogeneity is expected due to the similarity in first arrival apparent velocities and intercept times of the P and P phases. It is the high quality of the data set in terms of reliable g  n  amplitude information and coherent later phases that has allowed the detection of the generally small perturbations from the near-laterally-homogeneous structure.  Chapter 4. INTERPRETATION OF REFRACTION DATA  4.2.2  91  Near-shot Amplitude Inconsistencies  The first arrival amplitude of traces near shot points have not been reproduced well by the synthetics for a number of profiles; for some the calculated amplitudes are too high, for others, too low. These discrepancies are expected for a number of reasons. First, the near-shot first arrival energy has propagated solely within the sedimentary section and is sensitive to the variable velocity gradients and low Q structure of the sediments. As mentioned earlier, the sonic logs, from which the sediment velocity models were obtained, were smoothed and fit with straight lines and a constant Q value was assigned to the sedimentary layers. Second, the rapid change in material properties within the sediments suggests that the application of asymptotic ray theory may not be valid for the modelling of these arrivals. Finally, the near-shot first arrival energy is partially composed of the direct wave, which is not modelled.  4.2.3  Intracrustal Reflections  Between two and five intracrustal reflected phases have been picked on all profiles, except for Line B-shot B2, on the basis of the characteristic moveout of these branches. These phases have generally weak lateral coherency and are best observed in trace-normalized format (Figure 4.1). The uncertainty in the corresponding time picks is increased due to their location in the high amplitude wavetrain following the first arrivals. Often, a reflected phase will appear absent from a number of traces within the branch, perhaps due to destructive interference of many later phases. Through forward modelling it was found that the amplitude of these phases, often stronger than the first arrivals, could not be explained using simple step increases in velocity. The discontinuities in velocity would have to be so large (>0.5 km/s) that there would be an inconsistency with the observed travel time data. Two other structural zones capable of generating high reflectivity within  92  Chapter 4. INTERPRETATION OF REFRACTION DATA  Figure 4.1: Example of intracrustal reflections in the Line B - shot B3 section. Distance is from source to receiver and the traces are bandpass filtered from 2-12 Hz and scaled to a common maximum amplitude. The arrowheads identify the four intracrustal phases, Ri through i? . Reflections from the crust-mantle boundary, P P , are also labelled. 4  m  Chapter 4. INTERPRETATION OF REFRACTION DATA  93  the crust are: (1) a low-velocity zone with a significant velocity contrast and thickness on the order of or greater than one observed wavelength (approximately one kilometre)  (e.g.,  Morel-a-l'Huissier et al., 1987), or (2) a zone of thin (less than one wavelength thickness) alternating high- and low-velocity layers  (e.g.,  Berry and Mair, 1977; Sandmeier and  Wenzel, 1986; Wenzel and Sandmeier, 1988). The latter mechanism is favored here since neither the travel time or amplitude character of any profile suggests the presence of thick low-velocity zones. Another possibility for strong crustal reflectivity is random velocity variations (Gibson and Levander, 1988). A layered transition zone may be due to variations in composition, texture and/or anisotropy (Christensen and Szymanski, 1988). The theoretical response of a finely layered zone is best modelled using the reflectivity method of Fuchs and Miiller (1971). Therefore, although the positions of the intracrustal reflectors are constrained by the travel time picks, the precise structure of the reflectors cannot be determined through the present analysis and the calculated amplitudes are not expected to show good agreement with the observed amplitudes. The use of small step increases in velocity to represent reflectors has resulted in a monotonic increase in velocity with depth in the sub-basement of all models. Only the most obvious reflected phases have been picked; other weaker phases are present.  4.2.4  P r o m i n e n t P h a s e s of t h e I n - l i n e  Data  The P crossover point occurs between 160 and 180 km source-receiver distance on all n  end-shot profiles (Figures 2.2 to 2.13). The first arrivals at smaller offsets correspond to refracted arrivals in the sedimentary column, P , and the upper and mid crust, P . a  g  The P phase is a later arrival beyond the P crossover point and is composed of mid c  n  and lower crustal refracted arrivals; however, due to the relatively high vertical velocity gradients in the lower crust in all models, the lower crustal P phase is focused within c  ~50 km beyond the P crossover and is of relatively high amplitude. Beyond this point, n  Chapter-4. INTERPRETATION OF REFRACTION DATA  94  P corresponds to mid crustal arrivals of lower amplitude. The generally strong P P c  m  phase appears as a later arrival behind P before the P crossover point and after P up g  n  c  to ~50 km beyond the P crossover point. The critical point, where a P P ray path is n  m  coincident with the shallowest-penetrating P ray path, occurs between 110 and 140 km n  offset. 4.2.5  P P Character m  The strength and/or coherency of P P in all profiles is exceptional for a crustal refraction m  data set. As a result, the character of P P near the critical point has been qualitatively m  analyzed for each profile to better understand the nature of the crust-mantle boundary in the PRA region. The characterizations have been divided into three classes. First, the phase is described as sharp if it is represented by a short coda of ~0.25 s duration and is laterally coherent (Figure 4.2a). Second, the phase is described as diffuse if it is contained within a wavetrain of greater than 1 s duration that may represent an extended P P coda or the P P phase may comprise only a portion of the wavetrain and m  m  it is not readily identifiable therein (Figure 4.2b). In either case, the phase lacks lateral coherency. The third class is referred to as intermediate and describes a phase that has characteristics between or a combination of the sharp and diffuse classes. Analysis of the frequency content of all profiles suggests that these P P characters are not associated m  with source characteristics alone. Instead, a sharp P P may be associated with a laterally m  continuous, uniform crust-mantle boundary that may be a transition zone no more than about two kilometres thick. A diffuse P P character may indicate a disturbed or irregular m  boundary or a transition zone of several kilometres thickness. It is not possible to match the diffuse character of P P in the ray method synthetics, only the general amplitude m  trend. The tectonic implications of these boundary types, both of which are present in the PRA region, may be fundamental to the understanding of the history of the area.  Chapter 4. INTERPRETATION OF REFRACTION DATA P  m  80  95  P Character  140  170  DISTANCE  KM)  200  230  s  N  cn CO  cn 'r  '  OO  o  CD  i  in  80  140  170  DISTANCE (KM)  200  230  Figure 4.2: Comparison of sharp and diffuse P P character. Distance is from source to receiver and the traces have been bandpass filtered from 2-12 Hz and are presented in true relative amplitude format scaled by a factor proportional to distance. The arrowheads identify the near-critical P P phase, (a) Sharp P P character from Line B - shot B3. (b) Diffuse P P character from Line A - shot A3. Note that these P P phases sample roughly the same region of the crust-mantle boundary in approximately perpendicular directions (Figure 2.1). m  m  m  m  m  Chapter 4. INTERPRETATION OF REFRACTION DATA  96  The maximum amplitude of the observed P P is expected to occur at an offset ~20 km m  greater than the location of the ray theoretical critical point given the 5 Hz dominant frequency  Line A Interpretation  4.3  Ray diagrams, observed and theoretical travel times, and observed and synthetic sections for shots at 0, 177 and 342 km ( A l , A2 and A3) along Line A are shown in Figures 4.3, 4.4 and 4.5 and the Line A model in Figure 4.6.  Line A is the eastern N-S line that  crosses the axis of the PRA at 240 km (Figure 2.1). The observed data are of excellent quality except for shot A3 at distances less than 80 km where P arrivals are difficult n  to pick. This segment of the northern shot profile has the poorest quality of the data set and hence the upper mantle structure along Line A is less certain than for the other lines. The sedimentary section decreases in thickness from 3.75 to 1.9 km toward the north. The velocity of the underlying basement is roughly constant along the line, ranging from 6 to 6.1 km/s. A weak multiple generated within the sedimentary column (M ) is present x  on the middle-shot profile to the south where the velocity contrast at the sedimentbasement boundary is 4.5/6.1 km/s. This phase is most prominent on Line B. The upper crustal velocity anomaly centred at 140 km is required to match a slight decrease in apparent velocity of the first arrival along the A2 profile between 40 and 90 km. The most important feature of the Line A data is the difference in the apparent velocity of P and the character of P P between the southern and northern shot profiles. c  m  The higher apparent velocity of P from shot A3 requires a shallowing of high-velocity c  lower crust between 200 and 280 km, roughly beneath the axis of the PRA (Figure 4.6). For shots A l and A3, the strong later energy beyond the P crossover point is modelled n  Chapter 4. INTERPRETATION OF REFRACTION DATA  S  97  Line A - shot A l  N  DISTANCE (KM)  -20  30  80  130  180  230  280  330  380  -20  30  80  130  180  230  280  330  380  DISTANCE (KM) Figure 4.3: Ray diagram and comparison of calculated and observed travel times and amplitudes for Line A - shot A l . (a) Rays traced through velocity model. Dashed lines are model boundaries, (b) Comparison of calculated (lines) and observed (symbols) travel times. Symbol height is ±0.1 s.  o oC,O„ H  -20  30  80  130  180  230  280  330  DISTANCE (KM) Figure 4.4: Ray diagram and comparison of calculated and observed travel times and amplitudes for Line A - shot A2. (a) Rays traced through velocity model. Dashed lines are model boundaries, (b) Comparison of calculated (lines) and observed (symbols) travel times. Symbol height is ±0.1 s.  380  Chapter  -20  4.  INTERPRETATION  30  80  OF REFRACTION  130  180  DATA  230  100  280  330  380  DISTANCE (KM) Figure 4.4: (c) Synthetic section, (d) Observed section bandpass filtered from 2-12 Hz. Traces in (c) and (d) are presented in true relative amplitude format scaled by a factor proportional to the source-receiver distance.  30 8 0  bit*  sue!  U  ^ -I s.  of  330  ^ 6 ,  Chapter 4. INTERPRETATION OF REFRACTION DATA  -20  30  80  130  180  230  102  280  330  380  DI STANCE (KM) Figure 4.5: (c) Synthetic section, (d) Observed section bandpass filtered from 2-12 Hz. Traces in (c) and (d) are presented in true relative amplitude format scaled by a factor proportional to the source-receiver distance.  N  D I S T A N C E (KM) VELOCITY ( K M / S ) - 2 0 . . 13 5 7 9 Al  30  80  1301  £  180  230  280  A2|  330 „ A  A3  380 VELOCITY IKM/SJ 13 5 7 9  8 fa "a  bq  Pa S3  2 o  % o  S3  Chinchaga Low | Taltson Magmatic Arc A Peace River Arch  Figure 4.6: Smoothed isovelocity contour representation of the Line A model. Contours are P-wave velocity in km/s with 0.2 km/s contour interval. The shaded region is not constrained by the refraction data. Velocity-depth profiles at 50 and 290 km are shown at the left and right of the model, respectively. The upper dashed line at 2-5 km depth is the sediment-basement boundary. The short dashed line segments are the position of reflectors constrained by travel time picks. The crust-mantle boundary near 40 km depth appears as a gradient zone due to model smoothing before contouring. At the top of the model the location of the shot points (Al, A2 and A3) are indicated as well as the intersection with Lines B and C. The intersection of Line A with the Devonian axis of the PRA and the boundary between the Chinchaga Low and the Talston Magmatic Arc is indicated below the model.  O  % tl  2  o CO  Chapter 4. INTERPRETATION OF REFRACTION DATA  104  as lower crustal arrivals, P , refracted in the relatively strong gradients within the lower c  crust, resulting in the focusing of rays at caustics. The P P character of data from m  shot A l near the critical point is sharp, whereas for shot A3 it is diffuse. This suggests that the crust-mantle boundary is a sharp, uniform discontinuity south of the arch, but is more disturbed and non-uniform beneath the arch. The P P character of data from m  shot A2 toward the south is also sharp, but again diffuse toward the north. The absence of a post-critical P P on the shot A l profile has been modelled well by a crust-mantle m  transition zone 1.5 km thick from 7.85 to 8.2 km/s south of 100 km. Depths to the crust-mantle boundary are constrained by P P and P . At the southern m  n  end of the line the crustal thickness is 40.5 km and in the north it is 42.5 km with a downwarp at 260 km to 44 km depth; this location is roughly coincident with the axis of the PRA (Figure 4.6). The upper mantle velocity for the southern half of the line is 8.2 km/s, but the high apparent velocity of the P phase for the southern shot requires n  an upper mantle velocity of 8.45 km/s for the northern half of the line. The P phase n  on the southern shot profile between 180 and 230 km is relatively strong and has been modelled with a high velocity gradient in the upper 4 km of the mantle between 80 and _1  160 km. The gradient here, 0.02 s , is about twice that of the region's average upper mantle gradient. The P" phase in the shot A l profile (Figure 2.2) can be modelled as reflections off a horizontal boundary within the upper mantle at a depth of 74 km and associated velocity increase of 8.5 to 8.7 km/s. 4.4  Line B Interpretation  Ray diagrams, observed and theoretical travel times, and observed and synthetic sections for shots at 0, 156 and 312 km (BI, B2 and B3) along Line B are shown in Figures 4.7,  Chapter 4. INTERPRETATION OF REFRACTION DATA  4.8 and 4.9 and the Line B model in Figure 4.10.  105  Line B trends E-W and overlies the  axis of the PRA (Figure 2.1). The observed data on all profiles are of excellent quality. The thickness of the sediments along Line B decreases from 4.90 km to 1.95 km toward the east. The basement velocity (5.9 km/s) is roughly constant along the line except for a high of 6.2 km/s between 50 and 90 km. This upper crustal velocity anomaly is required to match a slight decrease in apparent velocity along the middle shot profile to the west between 0 and 40 km. It was not possible to pick intracrustal reflections from data on the shot B2 profile due to the predominance of strong multiples generated within the sedimentary section. (Figure 4.8d). The presence of these multiple phases suggests a reflective zone at the sediment-basement boundary, perhaps a sharp velocity discontinuity or a series of thin layers capable of generating constructive interference. A velocity discontinuity varying between 0.1 and 0.5 km/s at the sediment-basement boundary was sufficient to match observed amplitudes west of shot B2. The reason for the predominance of the multiples on the B2 profile overlying the PRA axis but not on other profiles of the data set, except shot A2 to the south, is unclear. One possibility is the unconformity at the sedimentbasement boundary defining the Devonian axis of the PRA. There is a significant difference in the P P character on both the end shot profiles and m  the middle shot profile of Line B. P P from shot B l is diffuse, whereas the shot B3 profile m  displays a sharp PmP. From shot B2, P P appears intermediate in character toward m  the west, but is sharp toward the east. This suggests that the crust-mantle boundary is disturbed or transitional in the west, but becomes uniform and sharp toward the east. To model the corresponding P P amplitude behavior, the lower crust in the west has a m  strong velocity gradient in a region extending 10 km above the crust-mantle boundary and a relatively small velocity discontinuity of 7.6 to 8.05 km/s at the boundary. At the east end of Line B the crust-mantle discontinuity is large with a contrast of 7.4 to  Chapter 4. INTERPRETATION OF REFRACTION DATA  107  Figure 4.7: (c) Synthetic section, (d) Observed section bandpass filtered from 2-12 Hz. Traces in (c) and (d) are presented in true relative amplitude format scaled by a factor proportional to the source-receiver distance.  Chapter 4. INTERPRETATION OF REFRACTION DATA  108  Line B - shot B 2  W  -20  30  80  -20  30  80  DISTANCE (KM)  E  130  180  230  280  330  130  180  230  280  330  DISTANCE (KM) Figure 4.8: Ray diagram and comparison of calculated and observed travel times and amplitudes for Line B - shot B2. (a) Rays traced through velocity model. Dashed lines are model boundaries, (b) Comparison of calculated (lines) and observed (symbols) travel times. Symbol height is ±0.1 s.  Chapter 4. INTERPRETATION OF REFRACTION DATA  w  109  Line B - shot B2  -20  30  80  130  180  230  280  -20  30  80  130  180  230  280  c  330  330  DISTANCE (KM) Figure 4.8: (c) Synthetic section, (d) Observed section bandpass filtered from 2-12 Hz. Traces in (c) and (d) are presented in true relative amplitude format scaled by a factor proportional to the source-receiver distance.  Chapter 4. INTERPRETATION OF REFRACTION DATA  110  Line B - shot B 3  W  E  DISTANCE (KM)  -20  30  80  130  180  230  280  330  -20  30  80  130  180  230  280  330  DISTANCE (KM) Figure 4.9: Ray diagram and comparison of calculated and observed travel times and amplitudes for Line B - shot B3. (a) Rays traced through velocity model. Dashed lines are model boundaries, (b) Comparison of calculated (lines) and observed (symbols) travel times. Symbol height is ±0.1 s.  Chapter 4. INTERPRETATION  OF REFRACTION  DATA  111  Line B - shot B3  w  E  cn CO  a  i cn  CD  -20  30  80  130  180  -20  30  80  130  180  230  280  330  230  280  330  C\J  cn cn CO  CO  CD  DISTANCE (KM) Figure 4.9: (c) Synthetic section, (d) Observed section bandpass filtered from 2-12 Hz. Traces in (c) and (d) are presented in true relative amplitude format scaled by a factor proportional to the source-receiver distance. '  DISTANCE VELOCITY (KM/S) "20 3  5  7  9  Bl  tg to TJ to  s O o  to to Wilson Island Terrane |Chinchaga| Taltson Magmatic A r c Low  Figure 4.10: Smoothed isovelocity contour representation of the Line B model. Contours are P-wave velocity in km/s with 0.2 km/s contour interval. The shaded region is not constrained by the refraction data. Velocity-depth profiles at 50 and 260 km are shown at the left and right of the model, respectively. The upper dashed line at 2-5 km depth is the sediment-basement boundary. The short dashed line segments are the position of reflectors constrained by travel time picks. The crust-mantle boundary near 40 km depth appears as a gradient zone due to model smoothing before contouring. At the top of the model the location of the shot points (Bl, B2 and B3) are indicated as well as the intersection with Lines A and D. The intersection of Line B with the boundaries between the Wilson Island Terrane, Chinchaga Low and the Taltson Magmatic Arc is indicated below the model.  Q  o  >»  B  to  Chapter 4. INTERPRETATION OF REFRACTION DATA  113  8.35 km/s. As on Line A, the strong later energy near and beyond the P crossover point is n  modelled as refracted P arrivals that are focused by strong velocity gradients in the c  lower crust. The depth to the crust-mantle boundary varies slightly around 40 km in the west and centre of Line B, but drops to about 43.5 km at the extreme eastern end of the line. The P* phase in the shot B3 profile (Figure 2.8) can be modelled as reflections off a horizontal boundary within the upper mantle at a depth of 69 km and associated velocity increase of 8.3 to 8.5 km/s.  4.5  Line C  Interpretation  Ray diagrams, observed and theoretical travel times, and observed and synthetic sections for shots at 0 and 271 km ( C l and C2) along Line C are shown in Figures 4.11 and 4.12 and the Line C model in Figure 4.13.  Line C trends E - W and roughly parallels the  Devonian axis of the PRA about 100 km to the south (Figure 2.1). The quality of the shot C2 data is excellent while that of C l is the poorest overall of any profile. The low signal-to-noise ratio of this profile is likely due to the poor transfer of energy between the shot and surrounding media. Reliable picks for all phases were possible. The thickness of the sediments along Line C decreases from 4.5 to 2.3 km toward the east. The basement velocity is ~6 km/s along the entire line except for an increase to 6.1 km/s near 50 km. This velocity anomaly accounts for a slight decrease in first arrival apparent velocity along the shot C l profile at about 100 km. The strong P refraction c  in the shot C2 profile is modelled by ray focusing in strong velocity gradients within the lower crust. A comparison of P P character near the critical point between the C l and C2 profiles m  is difficult due to the poor quality of the C l profile. P P character on the C2 profile m  -15  35  85  135  185  235  285  DISTANCE (KM) Figure 4.11: Ray diagram and comparison of calculated and observed travel times and amplitudes for Line C - shot C l . (a) Rays traced through velocity model. Dashed lines are model boundaries, (b) Comparison of calculated (lines) and observed (symbols) travel times. Symbol height is ±0.1 s.  Figure 4.11: (c) Synthetic section, (d) Observed section bandpass filtered from 2-12 Hz. Traces in (c) and (d) are presented in true relative amplitude format scaled by a factor proportional to the source-receiver distance.  Chapter 4. INTERPRETATION  OF REFRACTION  DATA  116  Line C - shot C 2  J o  n  -15  i  35  i  85  i  \b i  135  i  185  235  DISTANCE (KM) Figure 4.12: Ray diagram and comparison of calculated and observed travel times and amplitudes for Line C - shot C2. (a) Rays traced through velocity model. Dashed lines are model boundaries, (b) Comparison of calculated (lines) and observed (symbols) travel times. Symbol height is ±0.1 s.  2  Chapter 4. INTERPRETATION OF REFRACTION DATA  w  -15  117  Line C - shot C2  35  85  135  DISTANCE  185  235  285  (KM)  Figure 4.12: (c) Synthetic section, (d) Observed section bandpass filtered from 2-12 Hz. Traces in (c) and (d) are presented in true relative amplitude format scaled by a factor proportional to the source-receiver distance.  w VELOCITY (KM/S) 1  3  5  7  9  _  15 , r  85  D I S T A N C E (KM) 135  185.  285vEL0CITY  (KM/S)  tn to to to tg  2 H •—I  O to tg ^g  to >>  Wilson Island Terrane | Chinchaga | Taltson Magmatic Arc Low  Figure 4.13: Smoothed isovelocity contour representation of the Line C model. Contours are P-wave velocity in km/s with 0.2 km/s contour interval. The shaded region is riot constrained by the refraction data. Velocity-depth profiles at 50 and 220 km are shown at the left and right of the model, respectively. The upper dashed line at 2-5 km depth is the sediment-basement boundary. The short dashed line segments are the position of reflectors constrained by travel time picks. The crust-mantle boundary near 40 km depth appears as a gradient zone due to model smoothing before contouring. At the top of the model the location of the shot points (Cl and C2) are indicated as well as the intersection with Lines A and D. The intersection of Line C with the boundaries between the Wilson Island Terrane, Chinchaga Low and the Taltson Magmatic Arc is indicated below the model.  O o  GO  Chapter 4. INTERPRETATION OF REFRACTION DATA  119  is sharp, whereas the phase on the C l profile appears intermediate. The crust-mantle velocity contrast is 7.15 to 8.15 km/s in the west and 7.35 to 8.2 km/s in the east. The depth to this boundary is 41 km at the extreme west end of Line C and 41.5 km at the east end with an upwarp to 37 km in the centre of the line. An important feature of the shot C2 profile is the strong, high-apparent-velocity P~ phase ~0.5 s after P (Figure 2.11). The time and amplitude of this phase has been n  modelled as reflections off an upper mantle boundary at a depth of 56 km and associated velocity increase of 8.3 to 8.6 km/s. A relatively high velocity gradient of 0.025 s  _1  has  been assigned to the layer below this boundary and turning rays have been traced through it; however, it is the reflected phase that contributes most of the observed energy. The P" phase of the shot C2 profile is significantly higher in amplitude than the shot A l and B3 P" phases and corresponds to a shallower boundary (56 km as compared to 74 and 69 km depth). This suggests that these boundaries are relatively localized and not of related origin.  4.6  Line D Interpretation  Ray diagrams, observed and theoretical travel times, and observed and synthetic sections for shots at 0 and 295 km (Dl and D2) along Line D are shown in Figures 4.14 and 4.15 and the Line D model in Figure 4.16.  Line D trends N-S approximately 100 km west  of Line A and crosses the PRA axis at about 150 km (Figure 2.1). Both profiles are of excellent quality. The sedimentary section along Line D thins from 4.05 to 2.4 km towards the north. The basement velocity is 6.05 km/s for the southern half of the line and 5.65 km/s for the northern half. A local decrease of upper crustal velocity is required near 200 km to produce a slight decrease in the apparent velocity of the first arrivals along the shot D2  Chapter 4. INTERPRETATION  OF REFRACTION  DATA  120  Line D - shot Dl S  -20  30  80  -20  30  80  DISTANCE (KM)  N  130  180  230  280  330  130  180  230  280  330  DISTANCE (KM) Figure 4.14: Ray diagram and comparison of calculated and observed travel times and amplitudes for Line D - shot D l . (a) Rays traced through velocity model. Dashed lines are model boundaries, (b) Comparison of calculated (lines) and observed (symbols) travel times. Symbol height is ±0.1 s.  Figure 4.14: (c) Synthetic section, (d) Observed section bandpass filtered from 2-12 Hz. Traces in (c) and (d) are presented in true relative amplitude format scaled by a factor proportional to the source-receiver distance.  Chapter 4. INTERPRETATION  s  -20  OF REFRACTION  122  DATA  Line D - shot D2 30  80  DISTANCE (KM) 130  180  230  N  280  330  CM  cn CO  Q  I CO  o  \ •20  30  80  130  180  230  b  280  DISTANCE (KM) Figure 4.15: Ray diagram and comparison of calculated and observed travel times and amplitudes for Line D - shot D2. (a) Rays traced through velocity model. Dashed lines are model boundaries, (b) Comparison of calculated (lines) and observed (symbols) travel times. Symbol height is ±0.1 s.  330  Figure 4.15: (c) Synthetic section, (d) Observed section bandpass filtered from 2-12 Hz. Traces in (c) and (d) are presented in true relative amplitude format scaled by a factor proportional to the source-receiver distance.  S  DISTANCE  Wilson Island Terrane  A  |  (KM)  N  Chinchaga Low  Peace River A r c h  Figure 4.16: Smoothed isovelocity contour representation of the Line D model. Contours are P-wave velocity in km/s with 0.2 km/s contour interval. The shaded region is riot constrained by the refraction data. Velocity-depth profiles at 50 and 250 km are shown at the left and right of the model, respectively. The upper dashed line at 2-5 km depth is the sediment-basement boundary. The short dashed line segments are the position of reflectors constrained by travel time picks. The crust-mantle boundary near 40 km depth appears as a gradient zone due to model smoothing before contouring. At the top of the model the location of the shot points (Dl and D2) are indicated as well as the intersection with Lines B and C. The intersection of Line D with the Devonian axis of the PRA and the boundary between the Wilson Island Terrane and the Chinchaga Low is indicated below the model.  Chapter 4. INTERPRETATION OF REFRACTION DATA  125  profile between 150 and 200 km. Relatively strong lower crustal gradients have been used to generate high amplitude P on both profiles. The P P character on both profiles c  m  is intermediate. The weak amplitude of the shot D l P P phase near the critical point m  as compared to the high amplitude at post-critical distances between 170 and 220 km has not been reproduced well in the synthetics. This suggests that the crust-mantle boundary in the south is a velocity-gradient transition zone of several kilometres thickness according to the reflectivity modelling of Braile and Smith (1975). This feature has not been incorporated into the velocity model. The crust-mantle velocity contrast is 7.1 to 8.25 km/s in the south and 7.25 to 8.35 km/s in the north at depths of 37.5 km and 38.5 km, respectively.  4.7  Interpretation of Fan Shots  Due to the thinner than expected crust based on previous nearby refraction studies, strong _1  lower-crustal velocity gradients (> 0.03 s ) and large source-receiver offsets (> 250 km), the later arrivals observed in the Line A and B fan shot profiles are not P P, which could m  be used to extract crustal thicknesses  (e.g.,  Mereu  et al.,  1986), but are most likely mid  or lower crustal refracted arrivals, P . These arrivals are prominent (Figures 2.5 and 2.9), c  but have not been modelled. Therefore, the interpretation of data from the fan shots consists of a simple conversion of first arrival travel times of P into apparent crustal n  thickness variations (Line B) or upper mantle velocity variations (Line A) along profiles approximately 50 km north of Line B or 50 km west of Line A where the upcoming P  n  ray path samples the crust-mantle boundary (Figure 2.1). A weakness of this approach is that the travel time of P is sensitive to the upper mantle velocities (Line B) or crustal n  thicknesses (Line A) assumed, quantities which may vary laterally across the fan shot profile in regions not constrained by the in-line modelling. For this reason, the analysis is  Chapter 4. INTERPRETATION OF REFRACTION DATA  126  more qualitative than quantitative and the obtained parameter variations are "apparent". 4.7.1  Line B Fan Shot  The observed record section for the Line B fan shot (shot B4) is shown in Figure 4.17a. Using a reducing velocity of 8.25 km/s, an average obtained from the in-line modelling, the variation in first arrival time is about 1 s from west to east. The source-receiver offset varies from 340 km in the west to 280 km in the east with a minimum of 260 km near the centre of the profile. The ray-trace model used in the inversion of each travel time measurement along the profile is shown in Figure 4.17b. For all ray paths, a constant RMS crustal velocity (6.6 km/s), upper mantle velocity (8.25 km/s) and crustal thickness beneath the shot (40 km) was assumed, each value being an average obtained from the in-line models further south. In addition, the known and variable sediment velocity structure beneath each receiver along Line B was incorporated into the time-to-depth conversion. Crustal thickness at the point where the ray path leaves the crust-mantle boundary along an E-W profile approximately 50 km north of the PRA axis was then solved for to match each travel time measurement. The smoothed result of the inversion is shown in Figure 4.17c. An apparent thinning of the crust toward the west by as much as 10 km is observed. Errors in the assumed structure beneath shot B4 would not affect this thickness trend. Thinning of the crust toward the west is in agreement with the E - W Line B in-line model and the northern ends of the Line A and D models. However, the amount of thinning suggested by the fan shot interpretation is much larger than the 4 km obtained from the in-line models. The absolute amount of thinning cannot be determined from the fan shot interpretation due to the presence of possible E-W lateral variations in RMS crustal and upper mantle velocity as well as out-of-plane dip not accounted for by the model in Figure 4.17b. To better quantify the amount of thinning and an associated  Chapter 4. INTERPRETATION OF REFRACTION DATA  127  Figure 4.17: Line B fan shot (shot B4) interpretation, (a) Observed data bandpass filtered from 1-8 Hz and in true relative amplitude format scaled by a factor proportional to the source-receiver distance, (b) Simple ray-trace model used in the inversion of first arrival times. The crustal thickness at the point where the ray path leaves the crust-mantle boundary is determined separately for each travel time measurement, (c) The smoothed apparent crustal thickness variation along an E - W profile approximately 50 km north of Line B assuming a constant crustal and uppermantle velocity, (d) Smoothed apparent crustal thickness variations using the variable Line B RMS crustal velocity for the upward travelling ray path and the upper mantle velocities indicated. The intersection points of this profile with the north end of Lines A and D is labelled along with the upper mantle velocities required to match the crustal thicknesses at the north end of the Line A and D models. The upper mantle velocity at any point along the profile is obtained by linearly interpolating between these values and linearly extrapolating outside these values to those indicated at either end.  Chapter 4. INTERPRETATION OF REFRACTION DATA  128  uncertainty, a similar but more complicated inversion was performed in which the known RMS crustal velocity variations along Line B were incorporated into the upward travelling ray path and the upper mantle velocity was varied from east to west to provide agreement with the crustal thickness at the north ends of Lines A and D at azimuth 210° and 188°, respectively. The average RMS crustal velocity and crustal thickness from the Line B model was used for the ray path beneath the shot. Using a range of plausible upper mantle velocities for the ray paths west of 188° (8.1-8.5 km/s), a crustal thickness in the west of 36±4 km was obtained (Figure 4.17d). 4.7.2  Line A Fan Shot  In contrast to the Line B fan shot, the P ray paths of the Line A fan shot (shot A4) n  travel beneath the centre of the survey region that is sampled by the in-line modelling of the four lines. Therefore, this interpretation serves as a check on the in-line models as opposed to providing new structural information. The observed record section is shown in Figure 4.18a. Using a reducing velocity of 8.25 km/s, a variation in first arrival times from south to north of about 1 s is observed. The source-receiver offset varies from 320 km in the south to 290 km in the north with a minimum of 250 km near the centre of the profile. A drop in amplitude north of azimuth 270° is apparent. Figure 4.18b is the smoothed trace energy versus azimuth and shows a minimum at ~295°. This azimuth corresponds roughly to the Devonian axis of the PRA and may suggest that the energy along the P ray paths beneath the arch have been attenuated or scattered more strongly. n  The variation in source-receiver offset across the profile cannot account for this amplitude reduction. Increased scattering within the more-faulted sediments above the PRA is a possibility, although the dominant frequency of the P arrivals is 3 to 4 Hz. n  A similar inversion of first arrivals of the Line A fan shot was performed; however, the known velocity structure and crustal thickness below both the shot and receivers was  Chapter 4. INTERPRETATION OF REFRACTION DATA  129  Figure 4.18: Line A fan shot (shot A4) interpretation, (a) Observed data bandpass filtered from 1-8 Hz and in true relative amplitude format scaled by a factor proportional to the source-receiver distance, (b) The smoothed trace energy versus azimuth. The azimuth of the Devonian axis of the PRA is indicated, (c) Simple ray-trace model used in the inversion of first arrival times. The upper mantle velocity, V , is determined separately for each travel time measurement assuming the velocity structure and crustal thickness from the Line B and A in-line models for the downward and upperward travelling ray path, respectively, (d) The smoothed apparent upper mantle velocity variation along a N-S profile through the centre of the survey region. m  Chapter 4. INTERPRETATION OF REFRACTION DATA  130  obtained from the Line B and Line A models, respectively (Figure 4.18c). Therefore, the upper mantle velocity was determined for each arrival to match the observed first arrival travel time. The smoothed result of the inversion are shown in Figure 4.18d. An increase in velocity from about 8.1 to 8.4 km/s toward the north is observed and this is in agreement with the trend in upper mantle velocities obtained from the in-hne modelling in the centre of the survey region.  4.8  5-wave Analysis  Two features of the PRASE data set decrease the likelihood of observing shear phases. First, explosive sources may not generate significant shear energy. Second, vertical component seismometers have a sensitivity to compressional motion of roughly three times that of shear motion for the near-vertical approach (~20° from the vertical) of subbasement arrivals. Nevertheless, 5-phases have been observed in wide-angle, vertical component, explosive source data and used to infer crustal and upper mantle shear wave velocity structure (e.g., Tarkov et a l , 1981; Hall and Ali, 1985; Kullinger and Lund, 1986; Grad and Luosto, 1987). A search of the data set for pure and converted 5-phases was made by plotting all profiles at the equivalent S-wave velocity for a P-wave velocity of 8 km/s and a Poisson's ratio (cr) of 0.25 (4.62 km/s). The equivalent S'-phases corresponding to the P-phases used in the in-line modelling (S , S , S , S S , S and S") as e  g  c  m  n  well as the converted phases P S and S P were traced through the four velocity models m  m  using values of a near 0.25. A comparison of the theoretical travel times of these phases was then made with the observed data for all profiles. The S S phase was the only pure 5-phase identified in the observed data. Figures m  4.19a and 4.20a show the shot A3 and C2 data using a reducing velocity of 4.62 km/s. The S S phase occurs at ~13 s between 130 and 250 km on the A3 profile and between m  Chapter 4. INTERPRETATION  OF REFRACTION  131  DATA  Line A - shot A3  130  130  180  180  230  230  N  280  280  330  330  380  380  DISTANCE (KM) Figure 4.19: S-phases in the Line A - shot A3 record section, (a) Observed data bandpass filtered from 2-12 Hz and in common maximum amplitude format, (b) Calculated travel times of the S-phases for rays traced through the P-wave velocity model of Line A using a constant Poisson's ratio of 0.25 (dashed) and 0.251 (solid). A Poisson's ratio of 0.251 provides a best fit of the observed S S phase. m  Chapter 4. INTERPRETATION  OF REFRACTION  DATA  132  Figure 4.20: S-phases in the Line C - shot C2 record section, (a) Observed data bandpass filtered from 2-12 Hz and in common maximum amplitude format, (b) Calculated travel times of the S-phases for rays traced through the P-wave velocity model of Line C using a constant Poisson's ratio of 0.25 (dashed) and 0.24 (solid). A Poisson's ratio of 0.24 provides, a best fit of the observed S S phase. m  Chapter 4. INTERPRETATION OF REFRACTION DATA  133  60 and 170 km for shot C2. The ratio of P P to S S amplitude on all profiles varies m  m  from ~10 at wide-angle to between 2 and 5 near the critical point. The S„, S and S g  n  phases are not observed on any profiles. Given the presence of S S , the absence of S m  g  is puzzling since the corresponding take-off angles are similar. A contributing factor to this observation may be that at the predicted arrival time of S , the amplitudes of later g  P and converted phases are substantially above the noise level. The absence of S is n  commonly observed in refraction data. Two possible explanations given by Gajewski et al. (1987) are an increasing Poisson's ratio with depth or the effects of anisotropy in the upper mantle. The S phase may contribute part of the observed energy beyond 200 km c  source-receiver offset on some profiles (e.g., Figure 4.20a). The S S phase on all profiles was picked, except for the poor quality C l profile TO  in which S S is not observed. These picks have uncertainties of ~ ±0.3 s due to the m  low amplitude of S S relative to the later P and converted phases. A constant cr for m  the crust was determined for each of the eleven S S phases by matching observed and m  theoretical travel times for rays traced through the corresponding velocity model derived by the P-wave refraction modelling. The values of <r obtained range from 0.240 to 0.257 with an average of 0.249 and standard deviation of 0.005. The sensitivity of S S travel m  time at 140 km offset to variations in a (At/ACT) is ~55 s. Given the estimated arrival pick uncertainty of ±0.3 s, this corresponds to a change in a of ~0.005, the same as the standard deviation of the calculated a values. Figures 4.19b and 4.20b show the theoretical travel time curves of the S-phases for shots A3 and C2 for a = 0.25 and for the value of cr determined for that particular shot (0.251 and 0.240, respectively). Values of cr are greater than 0.25 in the northeast of the survey region where the north end of Line A intersects the east end of Line B (Figure 2.1), and less than 0.25 elsewhere. Given the estimated ±0.005 uncertainty associated with cr, this regionalization of cr is tentative, but may suggest a less rigid and/or more mafic crust in the northeast. Note  Chapter 4. INTERPRETATION OF REFRACTION DATA  134  that to obtain the crustal Poisson's ratios, only the travel times of the relatively high amplitude S S near the critical point was matched. The S phase in Figure 4.20 suggests m  c  a higher cr is appropriate as a mid crustal Poisson's ratio, implying a depth-dependent cr. It is possible that the onset of the phase picked as S S corresponds to a P to S m  converted phase at a discontinuity within the sediments or at the sediment-basement boundary.  Ray tracing tests indicate that a greater than 1.5 km/s velocity contrast  would be required to generate significant converted energy and the P-wave velocity above the discontinuity could be no more than ~4.5 km/s or the converted phase would not be observed beyond 120 km offset. Sedimentary velocity structure similar to this is observed in the sonic logs near some shot points; however, reflectivity modelling is required to test whether these zones of rapid velocity increase would appear as first-order discontinuities or smooth transition zones to the seismic wavefield. Depending on the depth to and the velocity above such a discontinuity, the converted phase would arrive 0.3 to 0.6 s before This corresponds to an underestimate of a by 0.005 to 0.01 if the converted phase  SmS.  is assumed for S S. m  A search of the data set for the crust-mantle converted phase, P S (or S P), unm  m  covered a weak series of arrivals on several profiles with the identification based on ray tracing through the P-wave models using the value of cr determined for the corresponding S S phase of the same shot. Figure 4.21 is the best observed example of this phase m  and is for shot A2 to the south. The P S phase has been observed in explosion crustal m  studies (e.g., Jacob and Booth, 1977; Luetgert et al, 1987; Holbrook et ai., 1988). Fuchs (1975) has shown that a strong P 5 requires a sharp first-order Moho discontinuity or m  a gradient transition zone of less than one seismic wavelength thickness. The strongest P S phase observed in the PRASE data set is shown in Figure 4.21; the corresponding m  P P character for the same shot is sharp. The absent or generally weak P S observed m  m  on all profiles suggests that the PRA region has a crust-mantle transition zone of greater  135  Chapter 4. INTERPRETATION OF REFRACTION DATA  L i n e A (south) — shot  60  80  100  120  DISTANCE  140  A2  160  180  (KM)  Figure 4.21: P S phase in the southern half of the Line A - shot A2 record section. A 2-12 Hz bandpass has been applied and each trace is scaled to a common maximum amplitude. Distance is measured from source to receiver. Both the PmS and P P phases are identified. m  m  Chapter 4. INTERPRETATION  OF REFRACTION  DATA  136  than a few hundred metres thickness. 4.9  Q Structure of the Crust and Upper Mantle  The quality factor, Q, is a dimensionless quantity that is a measure of the internal friction or anelasticity experienced by a propagating wave. Given an initial spectrum A (UJ), 0  a wave which has propagated through a homogeneous anelastic medium will have  a spectrum, A (uj), after a travel time t given by (Aki and Richards, 1980, pp. 168-169) t  A (w) = AoHexp[-u;i/2Q].; t  (4.1)  Two methods are commonly used for estimating Q from seismic data: modelling of amplitude variations with distance or time, or analysis of spectral ratios. Each technique has been applied to crustal refraction data (e.g., Braile, 1977; Hashizume, 1979). A requirement for the application of the amplitude method is a data set having precise amplitude control or numerous coherent phases on each seismogram.  A weakness of  this approach is the sensitivity of amplitudes to both velocity gradients and attenuation, which operate on the data in a similar manner (Hill, 1971). To allow application of the spectral ratio technique, the frequency response of receivers must be known and the data should have a broad frequency band. Amplitude modelling the PRASE data set suggested that it would be impossible to decouple the effects of velocity gradients and attenuation. Therefore, the spectral ratio technique was applied. 4.9.1  Spectral Ratio Method  The assumptions made in the application of the spectral ratio method to be presented here are described by Jacobson et ai. (1984). To summarize, it is assumed that the frequency-dependent effects upon the spectra are due to attenuation and that attenuation can be described by Q as stated in equation (4.1) in which Q is independent of frequency.  Chapter 4. INTERPRETATION OF REFRACTION DATA  137  Attenuation, as determined in this study, refers to effective attenuation, a combination of intrinsic and apparent attenuation. Intrinsic attenuation accounts for anelasticity as given in equation (4.1). Apparent attenuation includes inhomogeneous and structural effects such as scattering and tuning. For the PRASE data set, it would not be possible to separate the effects of intrinsic and apparent attenuation (Richards and Menke, 1983). Intrinsic attenuation can be assumed to be independent of frequency over the seismic frequency range (e.g., Knopoff, 1964). Apparent attenuation appears to be frequency dependent, but Menke (1983) shows that for narrow band data the dependence is weak. A 2.5-12 Hz frequency band was chosen taking into account both the observed spectra and instrument responses. The assumption of frequency-independent Q implies dispersion and according to a model presented by Kanamori and Anderson (1977), the amount of dispersion for a 2.5-12 Hz band is less than 0.05% and hence is not resolvable with the refraction data (Jacobson et al, 1981). The spectrum of a seismic wave that has propagated through an attenuating medium a distance r is given by (e.g., Hashizume, 1979)  'A (u>) = A {w)I(u)R(r)exp[-wt/2Q]. r  0  This is an extension of equation (4.1) to include the recording instrument response,  (4.2) I(UJ),  and the geometrical spreading losses R(r), assumed to be independent of frequency. For this experiment, the instrument response is essentially uniform for all receivers over the 2.5-12 Hz frequency band considered (Ellis et al, 1986) and can be set equal to one. Taking the logarithm of equation (4.2) and the ratio of the measured spectrum to the source spectrum to obtain units of dB gives 20log[A (u)/A (uj)} = 20log R(r) + 20log e[-ut/2Q]. T  o  (4.3)  Chapter 4. INTERPRETATION OF REFRACTION DATA  138  The slope, s, from equation (4.3) of the spectral ratio versus frequency in Hz is s ~ -27.3-f-  (4.4)  which allows the calculation of the raypath-averaged Q, i.e., Q. Several practical considerations concerning the application of equations (4.3) and (4.4) are necessary. First, the manner in which the spectra Ao(u>) and A (u) are calculated r  must be considered. Since equations (4.3) and (4.4) are applied to a single arrival with known ray path, the spectrum A,(a;) must be of that arrival only. An examination of the observed data to estimate the nominal wavelet length suggests that a 0.5 s time windowis appropriate for all calculations of  A (UJ). t  The Fourier spectra for these short time  series were calculated using a method described by Dziewonski et al. (1969) in which the maximum amplitude in the time window is picked after applying a Gaussian filter in the frequency domain. More specifically, the seismogram s(t) is Fourier transformed to a Gaussian weighting function W(u>,o> )=exp[ — p(uj 0  inverse Fourier transform of  W(ui, UJ )S(UJ) q  — u> ) /u) \ 0  2  2  is imposed on  S(UJ),  S(u>),  an  is taken to obtain s(t,aj ) and the spectrum 0  at the frequency u> is given by the maximum value of s(t,u) ) in the time window. The 0  0  value of a>o is varied over the bandwidth considered and p=15 was used as the Gaussian width parameter (Hashizume, 1979). The estimation of the source spectrum AQ(U) is required for each shot. One possibility is to use an analytic functional form chosen on the basis of its performance in the subsequent determination of the Q structure (Hashizume, 1979). A more desirable approach is to extract the source spectra directly form the data. The method used here was to consider the nearest few traces to each shot in which the first arrival energy propagates solely within the sedimentary section; the offsets of these traces were generally less than 10 km. From the measured spectra of the first arrivals and the known travel times of each, an estimate of A (UJ) 0  was made using equation (4.1) assuming sedimentary  Q  values,  Q, T  Chapter 4. INTERPRETATION OF REFRACTION DATA  139  of 100, 250 and 500. The average of these spectra for each value of Q, was determined for each shot (Figure 4.22). Three source spectra were obtained for each shot to assess the sensitivity of the inversion for crustal and upper mantle Q to the assumed value of Q. s  The 100 to 500 range in Q is assumed to bracket the "true" average Q„ value in a  the PRA region based on estimates of Q in southern Alberta (Clowes and Kanasewich, a  1970) and the Great Valley, California (Hwang and Mooney, 1986). The final consideration is the selection of data to be used which in turn leads to the consideration of the appropriate Q model parameterization. Only those phases within the data that are positively identified through the ray trace forward modelling can be considered since ray path information is required. This includes P , P , P m P , P and P". g  c  n  Where these phases intersect (i.e., the P crossover point), a mixture of arrivals is present n  and must be avoided. Therefore, only that portion of a phase isolated from all others by about 1 s was considered. The source spectra of the two shots used for each of the Line A and B profiles were assumed to be the same and the shot B l source spectrum was assumed for the Line A fan shot P arrivals. This assumption was necessary to increase n  the number of lower crustal and upper mantle observations. The Line B fan shot data was not used since no reliable estimate of the source spectrum was possible and these P  n  ray paths sample the crust and upper mantle outside of the main survey region.  Consideration of the ray paths shows P can provide upper crustal Q estimates, P and g  PmP  c  provide mid and lower crustal Q estimates, and P and P* provide upper mantle n  Q estimates (Figure 4.23). The wide scatter of Q values obtained for each phase from all shots suggested that a 3-layer (sub-basement), constant-Q model would represent the limit of resolution possible with this data set. Nevertheless, a comparison of the average Q values in the upper crust (Q ), lower crust (Q ) and upper mantle (Q ) is possible. g  c  m  In presenting the regional Q model, the upper 3 km represents the sedimentary layer in which Q was assumed, the boundary between the upper crust and the lower crust is 3  Chapter  4.  INTERPRETATION  OF  REFRACTION  140  DATA  o .  CM  Source Spectrum  Line C - shot C2  CO' CD  cr^ I—  CJ  CL.2 CO i  Qs  o  =  250  a  i  8  o  10  12  14  C\J  COCD Q  a  Q  3  = 100 = 250  <9, = 500  o  CM .  i  i  6  8  :  i  10  12  14  FREQUENCY (HZ) Figure 4.22: Example of source spectra for shot C2. (a) Source spectra of the three nearest-offset traces (solid) and the corresponding average spectrum (dashed) for C?,=250. (b) The average source spectrum for shot C2 for values of Q, of 100, 250 and 500.  Chapter 4. INTERPRETATION  OF REFRACTION  DATA  DISTANCE (KM) 90  120  1.41  150  180  210 J  240  270  150  180  210  240  270  Q model and Ray Paths  60  90  120  DISTANCE (KM)  Figure 4.23: Arrivals used in the Q study, (a) Representative ray paths of the six phases. The Q model is parameterized into four layers. For purposes of presentation the boundary between the upper and lower crust is placed at 15 km and that between the crust and mantle at 40 km based on the average depth of penetration of the P ray paths and the average crustal thickness, respectively, (b) Example of the windows of data corresponding to the six phases used in the Q study for the shot C2 section. A 0.5 s window about each arrival was used to estimate the observed spectra. Arrivals near phase intersection points were avoided. The windows as illustrated are only approximate since the actual pick of each arrival was used as the start of the 0.5 s window. For this presentation a 2-12 Hz bandpass filter was applied and each trace is scaled to a common maximum amplitude and distance is measured from source to receiver. G  Chapter 4. INTERPRETATION OF REFRACTION DATA  142  positioned at 15 km depth, and the crust-mantle boundary at 40 km. The crust below 15 km cannot be resolved into more than one layer since the arrivals from turning rays within the lower 10 km of the crust are often coincident with mid crustal arrivals. 4.9.2  Linear Programming Inversion  For a ray that has passed through a medium parameterized in terms of constant- Q cells (Dorman, 1969),  where t- and Qi are the travel time and Q value in the z'th cell, t is the total travel time, t  and m is the number of cells traversed by the ray path. Equations (4.4) and (4.5) show that the spectral ratio is a linear combination of the attenuation values in each model cell. Thus, it is possible to use linear inverse theory to obtain the Qi's assuming there is sufficient ray coverage in each model cell. From equation (4.5) the following system of linear equations is obtained for an m-layer Q model and the Q values for n arrivals (ray paths):  tuQi'  + t Q? 12  + ... + t Q lm  i  m  = t<?r  1  1  (4.6)  t^Qi  1  + t Q n2  l  2  + ... + t^Q-  1  x  = t Qn  1  where m~4 for the application to be presented here and the attenuation Q^ corresponds x  to Q~ , Q  X  2  1  corresponds to Qg , and so on; tij represents the travel time of the xth arrival 1  in the jth layer; t; is the total travel time of the z'th arrival; and Qi" is the raypathaveraged attenuation of the ith arrival. Equation (4.6) can be expressed more compactly in matrix form as Ta = s  (4.7)  143  Chapter 4. INTERPRETATION OF REFRACTION DATA  where T is an n x m matrix of travel times, a is the unknown attenuation vector, and s is the vector of measured spectral ratios adjusted by a scale factor given in equation (4.4). The solution of equation (4.7) through linear inverse methods yields the attenuation structure of the model. A total of 1205 arrivals was considered, representing 507 from the upper crust (P ), g  399 from the lower crust (P and P P) and 299 from the upper mantle (P and P"). Of c  m  n  these 1205 arrivals, 400 are from Line A, 420 from Line B, 148 from Line C, 154 from Line D and 83 from the Line A fan shot. The spectral ratio for each arrival was obtained from a least-squares-fit straight line within the 2.5-12 Hz frequency band. The quality of each least-square fit was measured and used later as a weighting applied to the spectral ratio measurements in the inversion. The fit of the z'th arrival, Wi, is defined as  Wi  where  (4.8)  = exp  \  p  is the observed spectra, r - is the predicted spectra from the straight Une fit, and  N is the number of frequency samples within the 2.5-12 Hz range. The value of w will be one for a perfect fit. The highest observed w was 0.679 while w=0.1 was the lowest allowable fit for a spectral ratio to be included in the subsequent inversion (Figure 4.24). The weights of the arrivals were scaled so that the total weight of all arrivals from each layer was equal. The solution of equation (4.7) was made through linear programming (Cooper and Steinberg, 1970, pp. 183-285) in which an L\ norm of the weighted misfit errors was minimized using the objective function  u  =  H  W  i  UJi\U \ ° i a  {  — ~  a?  i\  where a° and af are the observed and predicted average attenuation of the z'th arrival, respectively, and Wi is the weight of the zth arrival given by equation (4.8). The advantage  Chapter 4. INTERPRETATION OF REFRACTION DATA  144  Figure 4.24: Six examples of spectral ratios and their corresponding least-squares fit for P arrivals recorded north of shot A2 and C7„=250. The upper and lower solid curve in each case is the source and arrival spectrum, respectively. The dashed curve is the spectral ratio and the straight line is the corresponding least-squares fit over 2.5-12 Hz. The raypath-averaged Q value and fit parameter given by equation (4.8) is indicated in the lower right for each. g  Chapter 4. INTERPRETATION  OF REFRACTION  145  DATA  of the Li norm is that "outliers" will not strongly effect the solution. As mentioned earlier, the observed Qi's were widely scattered within phases and for all shots. An Li norm is also appropriate here since the Q model is parameterized in terms of constantQ layers that are suitable for solution by linear programming. An advantage of linear programming is that upper and lower bounds can be imposed on model parameters. The manner in which the source spectra were obtained required an assumed average sedimentary Q value. When solving equation (4.7), the value of Q^  1  1  (=Q7 )  was fixed at  1  the value used to obtain the source spectra and the Q^ values. If the true sedimentary Q structure in the PRA region is a constant value of Q , the solution of equation (4.7) s  would be independent of.the Q„ value assumed. The true Q structure is more complex and so equation (4.7) was solved three times assuming different Q, values to assess the sensitivity of the solution to the assumption of a constant Q,. The variation in observed Q^s for each phase was large. Of the 1205 spectral ratios measured, 200 to 300 Q°'s were less than zero depending on the assumed value of Q . s  Of the remaining positive Q°'s, a greater than one order-of-magnitude variation was observed, most values falling between 200 and 2000. Likely causes of this scatter as well as negative Q values are (1) the simple Q model assumed, (2) the spectral estimation of numerous arrivals within the short time window or an incorrectly picked arrival, and (3) deviations of the actual wave propagation from the assumed anelastic propagation described by equation (4.1)— in particular, the effects of scattering, tuning and noise. To minimize the effect on the inversion of the extreme low and high <5°'s, two steps were taken. First, only arrivals whose Q° values were between 250 and 2000 were considered. This restriction, along with the requirement that w > 0.1 was satisfied, reduced the number of acceptable arrivals to between 380 and 425. Second, to remove outliers, equation (4.7) was solved iteratively so that arrivals whose value of  8  =  io  \Q\ s[Q°/^i^  was more than 2.5 were removed from the inversion on the subsequent iteration. This  Chapter 4. INTERPRETATION OF REFRACTION DATA  Q. Layer upper crust lower crust upper mantle  100 650 786 1204  250 480 1008 897  146  Average 500 Q model 467 532 1017 937 858 986  Table 4.1: Q model of the crust and upper mantle elimination process required 3 iterations for each value of Q resulting in an average 8 of 3  between 1.43 and 1.47 for all 355 to 390 arrivals used on the final iteration. 4.9.3  Results  The results of inverting all spectral ratios for a 3-layer crustal and upper mantle Q model assuming three different values of Q are presented in Table 4.1. The first observation is B  the relative sensitivity of the solution for an assumed Q —100 versus <5«=250 and 500. s  This is probably due to the scatter of Q° values since the number of arrivals used in the inversion increases as Q decreases. The variation of the sub-basement Q values over the s  three models is a factor of 1.4, 1.3 and 1.4 for the upper crust, lower crust and upper mantle, respectively. For the average model, Q increases from ~500 to ~900 in the crust and is ~1000 in the upper mantle. To assess the significance of this Q model it is necessary to examine the fit of the observed Q value for each arrival to that predicted by the model. Figures 4.25a, 4.26a and 4.27a show the comparison of the observed and predicted Q values for all arrivals used and for Q, equal to 100, 250 and 500, respectively. For a perfect fit of the model, all points would he along the diagonal line defined by  P  Q =Q°.  Note that the  0,  Qs  fall  between 250 and 2000 since values outside this range were not used in the solution. Figures 4.25b, 4.26b and 4.27b show for each arrival, the average Q value along the  147  Chapter 4. INTERPRETATION OF REFRACTION DATA  Q  = 100  s  Upper Crust  Upper Mantle  O o CL<£)-  -\—i 2 3  0.0  0.2  i  a  r  15 6  10'  2  0 OBSERVED  0.4  0.6  3  156  156  10'  2  3 4 5 6  3  0 OBSERVED  0.8  SPECTRRL RATIO LST 50 FIT  1.0  0.0  0.2  0.4  0.6  156  10'  2  3  4 56  10*  0 OBSERVED'  0.8  SPECTRAL RATIO LST SQ FIT  1.0  0.0  0.2  0.4  0.6  0.6  SPECTRAL RATIO LST SO FIT  Figure 4.25: Comparison of observed and predicted Q values for Q = 100. (a) Logarithm of predicted Q values versus the logarithm of observed Q values for the upper crust, lower crust and upper mantle (each point is represented by a number from one to ten that corresponds to the fit given in equation (4.8) multiplied by ten), (b) Logarithm of the raypath-corrected Q values versus the fit given by equation (4.8) for the upper crust, lower crust and upper mantle. The dashed line represents the predicted Q for each layer. 3  1.0  Chapter 4. INTERPRETATION OF REFRACTION DATA  Q  s  Upper Crust  148  = 250 Upper Mantle  Lower Crust  Q  Q t— <_>  OoLUCE Q_co-  2  3 4 5 6  10'  2  3  ce Q-UJ-  -1 2  a 5 6  i i i i • •1 3 •» 5 6 10' 1  Q OBSERVED  0.0  ~i 1 1 r 0.2 0.4 0.6 0.8 1.0 SPECTRRL RATIO LST SQ FIT  r~i 2 3  i i i ' 4 56  a 10'  2  "i 0.2  1 0.4  1 0.6  10'  3156  0 OBSERVED  0.0  /  /  2  3 4 5 6  0 OBSERVED  r 0.6  SPECTRAL ROT 10 LST SO FIT  0.0  I  :  1  r  0.2  0.4  0.6  0.8  SPECTRAL RATIO LST SO FIT  Figure 4.26: Comparison of observed and predicted Q values for Q,=250. (a) Logarithm of predicted Q values versus the logarithm of observed Q values for the upper crust, lower crust and upper mantle (each point is represented by a number from one to ten that corresponds to the fit given in equation (4.8) multiplied by ten), (b) Logarithm of the raypath-corrected Q values versus the fit given by equation (4.8) for the upper crust, lower crust and upper mantle. The dashed line represents the predicted Q for each layer.  149  Chapter 4. INTERPRETATION OF REFRACTION DATA  = 500  Q  s  Upper Mantle  1—i—I  - i — i — i 111 " i 2  0 OBSERVED  3 4 5 6  10'  2  i—i M ' " i — i — i — M Ia'  IM"  3 4 5 6  3  1  ~ i — i i 11 " i 3 4 5 6  1—i—n~r  10'  2  0 OBSERVED  0.0  0.2  0.4  0.6  0.8  SPECTRAL RATIO LST 50 FIT  1.0  0.0  0.2  0.4  0.6  0 OBSERVED  4 56  10'  2  3  4 5 6  3 4 5 6  0.8  SPECTRAL RATIO LST SO FIT  1.0  0.0  0.2  0.4  0.6  o.  SPECTRAL RATIO LST SOFIT  Figure 4.27: Comparison of observed and predicted Q values for <5«=500. (a) Logarithm of predicted Q values versus the logarithm of observed Q values for the upper crust, lower crust and upper mantle (each point is represented by a number from one to ten that corresponds to the fit given in equation (4.8) multiplied by ten), (b) Logarithm of the raypath-corrected Q values versus the fit given by equation (4.8) for the upper crust, lower crust and upper mantle. The dashed line represents the predicted Q for each layer.  10'  Chapter 4. INTERPRETATION OF REFRACTION DATA  150  portion of the raypath in the deepest layer sampled by that particular arrival. The comparison of the observed and predicted Q's indicates the wide scatter of the observations. The raypath-corrected Q values exhibit an even larger spread and in some cases the mean of the points does not appear to correspond to the predicted Q value for that layer. However, the logarithmic dependence of the plots and the effect of the weights given by equation (4.8) must be considered. The reason for the large scatter of points in Figures 4.25b, 4.26b and 4.27b is that the L norm was applied to the raypathx  averaged Q values of each arrival, not the raypath-corrected Q value of the deepest layer sampled. When evaluating the raypath-corrected Q values, a phenomenon similar to that experienced in the inversion of RMS velocities for interval velocities is observed: small errors in shallower model layers generate large errors when evaluating a feature of the model within an interval. In fact, negative raypath-corrected Q values were obtained. Figure 4.28a shows the three Q models for each of the assumed Q, values presented in Table 4.1. Figure 4.28b shows the average Q model along with the range of Q values for each layer from the three calculated models. This figure indicates that the value of Q increases with depth from 550 ± 100 to 950 ± 150 in the crust and is 1000 ± 200 in the upper mantle. These uncertainties are required to bracket the three models in Figure 4.28a and may be underestimates since only the sensitivity to the assumed Q„ value is considered.  4.10  Seismic-gravity Consistency  Regional Bouguer gravity data for the PRA region is available (Geological Survey of Canada, National Data Base) and many studies have shown that there exists an empirical relationship between seismic velocity and density (e.g., Birch, 1964; Gardner et al., 1974), although a wide scatter of observations is common. Thus, it is appropriate  Chapter 4. INTERPRETATION OF REFRACTION DATA Q  Q 500  151  1000  1500  500  0  1000  1500  _J CD  -—  <—-C\J  0  ~ -J Q_cn  Q, = 100 Q. = 250 Q, = 500  : 550 ± 100  CD CM  _  C D .Qc = 950 ± 150 CO  CD o .  Q  m  in  = 1000 ± 200 ]  CD .  in  Figure 4.28: One-dimensional Q model for the PRA region, (a) Q models for values of Q, equal to 100, 250 and 500. (b) Average Q model (solid) and range of Q models for the three values of Q„ in (a) (dashed). to consider the relationship between the interpreted velocity structure of the PRA region and the observed Bouguer gravity data. It is not possible to infer unique density structure within the crust and upper mantle due to the ability of many different density models to reproduce a given observed gravity anomaly (Skeels, 1947). The purpose of this study is to check the seismic-gravity consistency with respect to the typical scatter of points upon which empirical velocity-density relationships are based (e.g., Ludwig et al, 1970). Only the qualitative assessment of plausible velocity-density relationships within the PRA region will be possible. The conversion of velocity models derived by the seismic refraction method to density for the purpose of assessing the seismic-gravity consistency is common practice (e.g., Catchings and Mooney, 1988). Often, to match the observed gravity profiles, anomalous low and/or high density regions are incorporated into the models with little regard as to  Chapter 4. INTERPRETATION OF REFRACTION DATA  152  their geological significance (e.g., Zucca et ai., 1986). Barton (1986) points out that given the range of densities possible for rocks of each seismic velocity and vice versa, the use of a seismic velocity as the only indication of density does not provide a useful constraint when reproducing observed gravity variations. Thus, the merits of such an exercise is perhaps questionable; however, the possibility of discovering major discrepencies with the seismic results and exploring the possible implications is sufficient justification  (e.g..  Spence et ai., 1985).  4.10.1  Peace River Region Bouguer G r a v i t y D a t a  The Bouguer gravity in the survey region dips smoothly from -70 mGal in the northeast to -90 mGal toward the southwest with a rapid decrease to less than -140 mGal in the extreme southwest reflecting the presence of the adjacent Rocky Mountains (Figure 4.29). 3  A reduction density of 2.67 g/cm was applied in the Bouguer correction and the average 2  station spacing within the region is one in every 10 km . No terrain correction has been applied to the data which has an estimated uncertainty of ±2 mGal. For comparison with theoretical gravity calculations, the areal distribution of the gravity data was sampled along each of the four refraction fines. 4.10.2  Calculation of Theoretical Gravity Response  Each of the four velocity models was converted to density using the velocity-density relationship shown in Figure 3.5. The corresponding theoretical two-dimensional gravity response was calculated using the technique of Talwani et ai. (1959) in which a two-dimensional cross-section is represented by n-sided polygons of constant density. A constant density was assigned to each model block (e.g., Figure 3.10a) using the average block velocity. To reduce end effects, the density models were extended a sufficient distance (>300 km) using the one-dimensional density structure at either model end.  Chapter 4. INTERPRETATION OF REFRACTION DATA B O U G U E R  G R A V I T Y  A N O M A L Y  153  (mgal)  LONGITUDE (W) Figure 4.29: Bouguer gravity anomaly map of the PRA region; contours in mGal. The location of the four refraction lines, Devonian axis of the PRA and edge of the Rocky Mountain disturbed belt is indicated. There is insufficient data in the southwest for contouring.  Chapter 4. INTERPRETATION  OF REFRACTION  154  DATA  Two-dimensional gravity calculations in which lateral homogeneity is assumed in the out-of-plane direction are acceptable when the in-plane length of an anomalous body is about five times less than its out-of-plane width (Grant and West, 1965, p. 230). Slight 3  perturbations (< ±0.05 g/cm ) were made to the densities of blocks through a trialand-error approach until the agreement between the theoretical and observed gravity was satisfactory (±5 mGal). The position of model boundaries were not altered in this procedure except for slight adjustments to the dip of the crust-mantle boundary outside of the region constrained by the refraction modelling. Only the relative shape of the observed gravity profiles were fit. 4.10.3  Results  Initial discrepancies between the observed and theoretical gravity responses along each line were generally no greater than ±20 mgal. For each line, Figures 4.30 to 4.33 show the isovelocity contour model, the position of blocks within the model whose densities 3  were significantly altered from their starting values (> ±0.02 g/cm ), and the final fit of the theoretical to observed gravity. The level of agreement between the observed and theoretical is reasonable considering the wide scatter of empirical measurements upon which the velocity-density relationship is based and suggests that the seismic models and observed gravity are consistent. However, the localized velocity anomalies within the seismic models that represent deviations from lateral homogeneity, generally correspond to regions where adjustments to the densities were required to provide a satisfactory fit (compare Figures 4.30a and 4.30b,  etc.).  This implies that the velocity anomalies do  not in general give rise to density anomalies and results in a more laterally homogeneous density structure than velocity structure. Of course, the manner in which the densities were perturbed is not unique, but the need to reduce densities in high-velocity regions and vice versa is likely significant.  155  Chapter 4. INTERPRETATION OF REFRACTION DATA  CD  i  "1  -20  I  30  r  80  i  130  i  180  i  230  DISTANCE (KM)  i  280  i  330  380  Figure 4.30: Gravity modelling of Line A. (a) Isovelocity contour representation of Line A. See Figure caption 4.6 for details, (b) Location of model blocks that required an adjustment from the initial density of greater than ±0.02 g/cm (less than ±0.05 g/cm ) to match the observed Bouguer gravity anomaly (< —0.02 g/cm —dots; > +0.02 g/cm —pluses), (c) Comparison of observed (symbols) and theoretical (Une) Bouguer gravity anomaly. Symbol height is ±2 mGal. 3  3  3  3  156  Chapter 4. INTERPRETATION OF REFRACTION DATA  -20  30  80  130  180  230  280  330  T r + + + + *+ + + + + + + + .... > •,+r+ + + +++  -20  30  80  130  4+  180  230  280  330  DISTANCE (KM) Figure. 4.31: Gravity modelling of Line B. (a) Isovelocity contour representation of Line B. See Figure caption 4.10 for details, (b) Location of model blocks that required an adjustment from the initial density of greater than ±0.02 g/cm (less than ±0.05 g/cm ) to match the observed Bouguer gravity anomaly (< —0.02 g/cm —dots; > +0.02 g/cm —pluses), (c) Comparison of observed (symbols) and theoretical (line) Bouguer gravity anomaly. Symbol height is ±2 mGal. 3  3  3  3  Chapter 4. INTERPRETATION OF REFRACTION DATA  -15  CD i  "1 -15  35  85  135  185  '"I  I  I  35  85  135  DISTRNCE (KM)  I  185  157  235  285  I  235  285  Figure 4.32: Gravity modelling of Line C. (a) Isovelocity contour representation of Line C. See Figure caption 4.13 for details, (b) Location of model blocks that required an adjustment from the initial density of greater than ±0.02 g/cm (less than ±0.05 g/cm ) to match the observed Bouguer gravity anomaly (< —0.02 g/cm —dots; > +0.02 g/cm —pluses), (c) Comparison of observed (symbols) and theoretical (line) Bouguer gravity anomaly. Symbol height is ±2 mGal. 3  3  3  3  Chapter 4. INTERPRETATION OF REFRACTION DATA  -20  30  80  130  180  230  158  280  330  280  330  . + •). +4- + + -H-'+ + 4 +  +++ ++>+• + ++++ +•+++•+++++++•++4 +++++++++•+++++++ ++•++++++•+++++++ •++++++•••++•+++4 * + 4- + 4- + +- + + + 4- + + + 4-4.4  Mo. •  C\J  CL.cn ' LU  Q  o.  130  180  DISTANCE (KM]  230  Figure 4.33: Gravity modelling of Line D. (a) Isovelocity contour representation of Line D. See Figure caption 4.16 for details, (b) Location of model blocks that required an adjustment from the initial density of greater than ±0.02 g/cm (less than ±0.05 g/cm ) to match the observed Bouguer gravity anomaly (< —0.02 g/cm —dots; > ±0.02 g/cm —pluses), (c) Comparison of observed (symbols) and theoretical (line) Bouguer gravity anomaly. Symbol height is ±2 mGal. 3  3  3  3  Chapter 4. INTERPRETATION OF REFRACTION DATA  4.11  159  S u m m a r y : R e g i o n a l Velocity Structure  The crust in the PRA region contains weak to moderate lateral variations in structure: the maximum dip on a sub-basement model boundary is 5°. The refraction data do not suggest that the crust is composed of a few layers separated by significant velocity discontinuities or thick low-velocity zones. Instead, the crust is interpreted to consist of zones of positive velocity gradient extending from the basement to the crust-mantle boundary with reflective horizons scattered throughout. If these reflectors are associated with discontinuities, they can be no more than ~0.3 km/s increases in velocity. The velocity gradient variation with depth in the PRA region is quite consistent and may be summarized as follows: moderate gradients of between 0.02 and 0.03 s 15 km of the crust, weak gradients of between 0.005 and 0.01 s  -1  -1  in the upper  in the mid 15 to  30 km of the crust, and strong gradients of between 0.035 and 0.055 s  _1  in the lower  10 km of the crust. In this sense, the crust may be "divided" into three zones. Average velocities within these zones range from 5.95-6.35 km/s, 6.45-6.75 km/s and 6.90-7.35 km/s, respectively. The average velocity and gradient below the crust-mantle boundary is 8.25 km/s and 0.01 s  _1  (Figure 4.34).  The thickness and range of velocities and gradients in the three crustal zones are very similar to those reported for the East European Platform from DSS data (Pavlenkova, 1988). Pavlenkova suggests that the upper and middle crust are of the same granitegneiss composition and so are separated by mechanical or physiochemical changes under certain P - T conditions. The structure of the upper crust is more laterally inhomogeneous and the mid crust represents a more homogeneous weakened layer. The lower crust is formed as a result of material differentiation separating crustal rocks from more basic mantle rocks. Having presented the average one-dimensional velocity structure of the PRA region,  160  Chapter 4. INTERPRETATION OF REFRACTION DATA  1  VELOCITY 3 5  8.25 km/s  (KM/S) 7 9  1  O.Ol s"  CD LO  Figure 4.34: One-dimensional velocity model of the PRA region. The average upper and lower velocity and velocity gradient within the sediments (0-3 km), upper crust (3-15 km), mid crust (15-30 km) and lower crust (30-40 km) is indicated as well as the average velocity and gradient directly below the crust-mantle boundary. These averages were obtained by laterally averaging the four in-line refraction models. Note that the velocity gradients within each layer as calculated from the average upper and lower layer velocity may not correspond to the gradients indicated since each of the three crustal layers are generally composed of two or more thinner layers separated by small velocity discontinuities for the purpose of ray tracing and amplitude calculations.  Chapter 4. INTERPRETATION OF REFRACTION DATA  161  we now examine the regional variations in velocity structure.  This is best achieved  through a series of contour plots of the variation in a particular parameter (velocity, depth or thickness) over the survey region as sampled along the four refraction lines. However, there are a number of complications with this approach. First, the contours must be a smooth interpolation or extrapolation of the trends established along the refractions lines in off-line regions. Second, at the four intersection points of the lines (Figure 2.1), exact agreement was not achieved for reasons discussed in section 2.4.2. Therefore, the average of the two parameter values at each intersection point was used by adding a Gaussian-shaped adjustment term to the parameter values along each line centred about each intersection point and with a half-width of 20 km and peak amplitude of half the difference between the two parameter values at the intersection point. This correction is necessary to present regional variations in a consistent format. The contour surfaces are not sensitive to the width of the Gaussian smoothing given the scale of regional structure resolved by the refraction data. The effect of this correction must be considered when interpreting the contour structure since artifacts can be introduced near intersection points where the most significant inconsistencies occur. In particular, the lower crustal structure along Line D is the most inconsistent with its two intersecting lines and results in artificial structure in the corresponding contour surfaces. Given the relatively poor constraint in the centre of Line D, it may be possible to improve the consistency at the intersection points and thereby remove these artifacts. The basement velocity varies between 5.6 and 6.2 km/s with an average of 5.95 km/s for the PRA region (Figure 4.35). Basement velocities are generally less than 6 km/s in the northwest of the survey region and greater elsewhere. Localized upper crustal velocity anomalies will be discussed in chapter 6. The average RMS crustal velocity in the region is 6.6 km/s with values ranging from 6.45 to 6.75 km/s (Figure 4.36). The RMS crustal velocity is greater than 6.6 km/s in the northeast and generally less  Chapter 4. INTERPRETATION  OF REFRACTION DATA  162  59°,  58°  S57° LU Q  < 56°  55° 100  km Devonian axis of Peace River Arch  54°  121°  120° 119° 118°  117 116  c  c  L O N G I T U D E (W) Figure 4.35: Basement velocity in the PRA region as determined from the four in-line refraction models; contour interval is 0.1 km/s. The negative (dotted) and positive (elsewhere) aeromagnetic domains and the Devonian axis of the PRA are shown (see Figure 1.2).  Chapter 4. INTERPRETATION OF REFRACTION DATA  L O N G I T U D E (W) Figure 4.36: RMS crustal velocity in the PRA region as determined from the four in-line refraction models; contour interval is 0.05 km/s. The negative (dotted) and positive (elsewhere) aeromagnetic domains and the Devonian axis of the PRA are shown (see Figure 1.2).  Chapter 4. INTERPRETATION OF REFRACTION DATA  164  elsewhere. The average depth at which a velocity of 6.8 km/s is reached is ~30 km with a localized shallowing to 26 km in the east (Figure 4.37). There is a significant thickness of high-velocity (> 7.2 km/s) lower crustal material throughout most of the PRA region, the greatest thickness in the extreme midwest and mideast (Figure 4.38). The average velocity at the base of the crust is 7.35 km/s, but varies between 7.1 and 7.75 km/s with the highest velocities present in the northeast and northwest (Figure 4.39). The average crustal thickness of the P R A region is 40 km with variations between 37 and 44 km (Figure 4.40). The crust is greater than 40 km thick in the eastern half of the survey region and generally less in the west. The thickest and thinnest crust occurs in the northeast and northwest, respectively. P P character is generally sharp in the east, m  but intermediate or diffuse in the west (Figure 4.40). The average upper mantle velocity is 8.25 km/s with local deviations of ±0.2 km/s (Figure 4.41). The upper mantle velocity is greater than 8.2 km/s in the northeast and less elsewhere.  165  LONGITUDE (W) Figure 4.37: Depth to 6.8 km/s in the PRA region as determined from the four in-line refraction models; contour interval is 1 km. The negative (dotted) and positive (elsewhere) aeromagnetic domains and the Devonian axis of the PRA are shown (see Figure 1.2).  Chapter 4. INTERPRETATION OF REFRACTION DATA  L O N G I T U D E (VV) Figure 4.38: Thickness of the lower crust with velocity greater than 7.2 km/s in the PRA region as determined from the four in-line refraction models; contour interval is 1 km. The negative (dotted) and positive (elsewhere) aeromagnetic domains and the Devonian axis of the PRA are shown (see Figure 1.2).  Chapter 4.  INTERPRETATION  OF REFRACTION  DATA  167  L O N G I T U D E (W) Figure 4.39: Velocity at the base of the crust in the PRA region as determined from the four in-line refraction models; contour interval is 0.1 km/s. The negative (dotted) and positive (elsewhere) aeromagnetic domains and the Devonian axis of the PRA are shown (see Figure 1.2).  Chapter 4. INTERPRETATION  OF REFRACTION  DATA  168  L O N G I T U D E (W) Figure 4.40: Crustal thickness in the P R A region as determined from the four in-line refraction models (solid contours) and the Line B fan shot (dashed contours); contour interval is 1 km. The P P character near the critical point is indicated along each line at the point where the crust-mantle boundary is imaged. The negative (dotted) and positive (elsewhere) aeromagnetic domains and the Devonian axis of the PRA are shown (see Figure 1.2). m  Chapter 4. INTERPRETATION OF REFRACTION DATA  169  L O N G I T U D E (W) Figure 4.41: Upper mantle velocity in the PRA region as determined from the four in-line refraction models; contour interval is 0.1 km/s. The negative (dotted) and positive (elsewhere) aeromagnetic domains and the Devonian axis of the PRA are shown (see Figure 1.2).  Chapter 5  PROCESSING A N D INTERPRETATION OF REFLECTION  5.1  DATA  Introduction  Section 1.5 discussed the advantages of combined near-vertical reflection and wide-angle reflection/refraction surveys to provide more detailed structural and velocity information than is possible from either technique alone. Interpretations of the two data sets can present conflicting pictures of the nature of the crust  (e.g.,  Klemperer et al., 1986; Okaya,  1986; Phinney, 1986). Of course, the structure imaged by the two methods is the same, it is simply being "looked at" in two different ways. Braile and Chiang (1986) clearly demonstrate the different response of each technique to a model of the crust-mantle boundary and show features of the model to which the techniques respond. Within the PRA region, no deep near-vertical reflection data exists. However, through a simple processing procedure, it is possible to extend the listen-time of industry-acquired uncorrelated Vibroseis reflection data to image deep sub-basement structure. A single ~10-km-long line coincident with Line B of the refraction data set was obtained for such a purpose. The processing and interpretation of the data set to be presented here is viewed as a test of the potential of extended-listen processing for deep crustal imaging in the PRA region. Although the reflection line coincides with less than 1% of the total length of the refraction lines, if the extended-listen processing is successful, a comparison of the reflection data with the refraction model where coincident will (1) produce a different view of refraction boundaries and layers modelled as step-discontinuities and gradient  170.  Chapter  5.  PROCESSING  AND  DISTANCE  W -20  INTERPRETATION  30  80  130  OF  REFLECTION  DATA  (KM) 180  171  E 230  280  330  reflection line  Figure 5.1: Location of the reflection line within the Line B refraction model is indicated. See Figure 4.10 caption for details. zones, (2) serve as a check on the depth to intracrustal reflectors and the crust-mantle boundary by converting the refraction model to two-way, zero-offset travel time, (3) elucidate detailed structure not resolvable with refraction data, and (4) provide a better understanding of and confidence in the refraction interpretation of all four lines. 5.2  Description of Data Set  A Vibroseis reflection data set shot in 1983 was obtained from Shell Canada Limited. The line is 9.36 km long and is oriented E-W, roughly coincident with the refraction Line B at 100 km model distance (Figure 5.1). The west and east ends of the reflection line are approximately 8 and 12 km south, respectively, of Line B. Figure 5.2 shows the velocity-depth profile at 100 km distance as determined from the Line B refraction model.  Chapter 5.  PROCESSING  AND  INTERPRETATION  VELOCITY  o _l  OF REFLECTION  DATA  172  (KM/S)  1  1  in  Figure 5.2: Velocity-depth profile taken from Line B at position of reflection line (100 km model distance). The depth and zero-offset two-way travel time of the basement, intracrustal reflectors and the crust-mantle boundary are indicated. The depth and two-way travel time of the basement, intracrustal reflectors and the crustmantle boundary are indicated for later comparison with the processed reflection data. The two-way travel times were obtained by performing zero-offset ray tracing through the refraction model taking into account the dip on each boundary. The parameters describing the data set are presented in Table 5.1. Most importantly, the source was a 13 s linear Vibroseis upsweep from 15 to 80 Hz with a total listen time of 16 s. A 96-channel, split-spread configuration was used with a 2 ms sampling interval. The receiver spacing was 30 m and the source interval was variable, but was nominally 90 m, resulting in a nominal stack fold of 16. Figure 5.3a shows the source and spread geometry and Figure 5.3b shows the realized fold along the line due to irregular source spacing. Figure 5.3c shows the relatively smooth elevation variation along the line from approximately 0.7 to 0.8 km above sea level from west to east.  Chapter'5. PROCESSING  AND INTERPRETATION  DATA  DISTANCE (m)  W 0  OF REFLECTION  1000  — i  2000 1  3000  4000  1  1  5000 1  6000 1  173  E 7000 1  8000  9000 i  STACKING D I A G R A M source  Y 4 CMP's  end of spread  o IS  «  Q g<  Figure 5.3: The source-receiver geometry, stack fold and elevation variations of the reflection line, (a) The stacking diagram in which the source location increases upward and the receiver points increase to the right. The source location, extent of the spread and CMP points are defined in the upper left inset. The horizontal and vertical axes are distance along the line, (b) Stack fold, (c) Elevation above sea level.  Chapter 5. PROCESSING AND INTERPRETATION OF REFLECTION DATA 174  Shell Line 5-925 Source Vibroseis Linear 13 s upsweep 15-80 Hz Split-spread Receiver array 96-channel 1.5 km max offset Source interval 90 m (nominal) Receiver interval 30 m Stack fold 16 (nominal) Sample rate 2 ms Total listen time 16 s Table 5.1: Parameters of the reflection data set. 5.3  Processing  The demultiplexed uncorrelated field tapes were acquired so that reprocessing at the correlation step could be performed. It was possible through a special correlation procedure to extend the final correlated traces from 3 s to 15 s two-way travel time. This record length corresponds to subcrustal depths as predicted by the refraction model (Figure 5.2). The processing of the data after correlation was entirely routine and in one aspect is simplified over conventional processing: accurate velocity analysis in the sub-basement is not required since the relatively short maximum offset (1.5 km) results in low sensitivity of the normal moveout of deep events to significant changes in the stacking velocity (Okaya et al., 1988). The normal moveout at 1.5 km offset below 4 s two-way travel time is less than 10 ms. The data and sweep were resampled to 4 ms prior to correlating.  Chapter 5.  5.3.1  PROCESSING  AND INTERPRETATION  OF REFLECTION  DATA 175  Extended-listen C o r r e l a t i o n  Deep crustal reflections may be extracted from shallow industry data provided a Vibroseis source has been used and the original uncorrelated field gathers are available. Extendedlisten processing has been summarized by Okaya (1986). Given a Vibroseis sweep of duration t  sweep  , a reflected event at  t  e v e  will end at i  nt  s w e e p  +t  eve  nt-  After correlation,  the uncorrelated event will be compressed to a simple, non-causal, zero-phase wavelet centred at £  event  (Figure 5.4a). The latest possible reflected event with a fully recorded  sweep is given by ^event  where t  recorc  ^record  ^sweep  i is the total listen time. In other words, the length of the final record section  after correlation is ^section  Events after t  scctlon  ^record  ^swccp •  may be extracted if the tail-end of the sweep is truncated by At  seconds. The resulting correlated profile is lengthened to ^extend  —  ^record  (^sweep  At)  —  ^record  ^swecp ~l~ At  —  ^section H~ At.  For an upsweep from low to high frequency, shortening the sweep by At will remove high frequencies from events occurring after  £  s e c  tion;  however, for deep reflections, this is ac-  ceptable since higher frequencies are not returned from the lower crust due to attenuation effects. Two methods are available for extended-listen processing. First, a fixed or truncated sweep, that has a constant length and frequency bandwidth lower than the original sweep, can be correlated with the data so that the entire section has a bandwidth reduced from  Chapter 5.  PROCESSING  AND  B fixed  A normal  >  >  INTERPRETATION  -4  OF  REFLECTION  DATA  176  self-truncating  V-  J4-  f  > >  4  'r  >  > >  v  s w e e p'  s e c t i o n  >  -  te x t e n d  ^  ^ r e c o r d  Figure 5.4: A comparison of (a) normal, (b) fixed and (c) self-truncating sweep correlation. R is the field record and C is the result of correlating R with the Vibroseis sweep, shown between R and C at various lags. In the normal method (a), the last possible correlated time is t ection- All frequencies contained in the sweep are present in the correlated trace. Wavelet centred at t ection is missing lower half because correlation stopped at t i . In the fixed method (b), more correlation of the shorter sweep with the field recording is possible, but at the expense of a decreased frequency content of the correlated trace. Wavelet centred at t a is missing lower half because correlation stopped at t d- In the self-truncating method (c), normal correlation occurs until Section- Subsequent correlation is performed with the sweep truncated at t cord resulting in a decreased frequency content with increasing time (as seen in the deeper correlated wavelets). Figure from Okaya (1986). S  S  s e c t  o n  exten  exten  re  Chapter 5. PROCESSING AND INTERPRETATION OF REFLECTION DATA 177  the original section (Figure 5.4b). Second, a self-truncating sweep, in which the full sweep is used after the data have been tapered and padded with zeros to the desired additional time, results in a diminishment of high frequencies with time after  t  section  (Figure 5.4c).  Most importantly, the self-truncating method allows for maximum use of the original sweep at each point of the extended data after £  sect  sweep could be used to identify events as late as  t  ; . In theory, the self-truncating on  r e C 0  rd;  however, the data would be  extemely band-limited near this time. The fixed sweep method was first used by Oliver and Kaufman (1977) to add an additional 3.5 s to a 15 s deep reflection profile. Finckh et al. (1984), Finckh et ai. (1986) and Rotstein and Trachtman (1986) applied the fixed sweep method to examine deep crustal structure using data originally obtained for shallow studies. Okaya (1986) applied the self-truncating method to a conventional 4 s data set to extend the record lengths to 12 s. The self-truncating method as described by Okaya (1986) was applied to the data set of this study to extend the correlated record lengths from 3 s to 15 s. Figure 5.5 shows the decrease in sweep length and high frequency with time for this study. This figure shows that reflections from the crust-mantle boundary, predicted to occur at ~13 s from the refraction modelling, will be imaged by 3 s of the sweep with a bandwidth of 15-30 Hz. The only special processing step applied in the self-truncating correlation was a 0.2 s cosine taper to the end of the uncorrelated data before padding with zeros to 28 s. Figure 5.6 shows six common-shot gathers from different points along the fine after extended-listen correlation selected to illustrate different features of the data set. Each of the shot gathers in Figure 5.6 has been bandpass filtered (corner frequencies 5-15-4060 Hz) and a time varying gain control was applied to balance the gathers and enhance the deeper events. The character of each gather is quite different. After the strong refracted first arrivals on all gathers, the reflections from within the sedimentary section are variable with only  Chapter 5. PROCESSING AND INTERPRETATION OF REFLECTION DATA 178  ~i 4  5  6  7  8  9  10  i  i  i  r  11  12  13  14  15  TWO-WAY TRAVEL TIME (S)  Figure 5.5: The sweep length used and drop in high frequency as a function of time for a 13 s linear upsweep and application of the self-truncating method. The full sweep is used from 0 to 3 s. gathers from 7.2 km and 8.46 km having clear events between 0.5-2.0 s (event A). Gathers from 5.94 km and 7.2 km contain a steeply westward-dipping event between 3-5 s and 2-4 s, respectively, having an apparent velocity of ~2.2 km/s (event B). These events give rise to a prominent steeply-dipping artifact in the stack and may represent a direct wave reflected from a near-surface high impedance contrast to the east of these gathers. The eastern side of the 7.2 km gather below 2 s appears banded in character with a uniform hyperbolic moveout centred outside of the spread (event C). This feature is believed to be environmental noise (e.g., due to an oilwell pump). This type of impulsive-like source would cause a replication of the sweep into the gather after sweep correlation. A similar banded noise column is present in the right-centre of the gather from 8.46 km (event C). Gathers from 1.59 km and 2.61 km between 10-14 s appear smeared-like (event D) and the 5.94 km gather between 12-14 s contains events which mimic the first arrivals (event E) (cf., Okaya, 1986). The precise origin of these features is uncertain.  Chapter 5. PROCESSING AND INTERPRETATION OF REFLECTION DATA 180  Karageorgi et al. (1988) suggest that these artifacts are related to resonance of the vibrator sources and present a method for their removal in which filtering is applied in the frequency-uncorrelated time domain. Harmonic distortion is not considered as a cause since a low-to-high frequency upsweep was used (Seriff and Kim, 1970). Eastwarddipping events near 9.4 s and 12.2 s in the 5.94 km and 7.2 km gathers, respectively, give rise to prominent events in the stack (events F and G). The 4.17 km gather corresponds to the central part of the reflection line in which the stack has a relatively transparent appearance within the sediments and the sub-basement.  5.3.2  Post-correlation Processing  A conventional pre-stack processing sequence was followed. First, a thorough examination of the shot gathers was made to assess the need for manual trace editing. While prominent noise features exist throughout the data set (Figure 5.6), the relatively low signal-to-noise ratio below 2 s and the low stack fold suggested that trace editing would not be advantageous and so was not performed.  Elevation corrections and refraction  statics were performed using the method of Coppens (1985) with a replacement velocity of 2 km/s and a datum of 0.75 km above sea level. The refraction models presented in chapter 4 are corrected to a 0.65 km datum so that a time delay of 0.1 s of events on the reflection profile as compared to the same event in the refraction model is expected. Residual statics and velocity analysis within the upper 2 s of the section were performed alternatively for two iterations. Residual statics were applied over the most prominent reflectors within the sediments between 0.5-1.8 s in which the maximum time shift permitted in the iterative stack-power maximization technique of Ronen and Claerbout (1985) was 0.1 s. Velocity analysis was through the method of constant-velocity stacks over the upper 2 s. A single velocity function was found to be sufficient for the sedimentary section and is presented in Table 5.2. The purpose of velocity analysis was  Chapter 5. PROCESSING AND INTERPRETATION OF REFLECTION DATA 181  Time (s) 0.0 0.6 0.9 1.3 1.6  Sediment Velocity (km/s) 2.7 2.7 3.6 4.5 4.8  Sub-basement Time (s) Velocity (Km/s) 2.9 5.3 4.8 5.7 6.3 5.9 9.4 6.1 12.8 6.4 15.0 6.7  Table 5.2: Sediment and sub-basement stacking velocity function. to improve the estimation of residual statics; the enhancement of sedimentary reflectors was of secondary importance. Velocity analysis in the sub-basement was not necessary since RMS velocities were obtained from the Line B refraction model and tests showed the insensitivity of the normal moveout of deep events to the stacking velocity. Table 5.2 gives the stacking velocity function for the sub-basement. The final pre-stack processes involved a CMP-sort, automatic gain control of each trace over the 15 s window (L  x  norm, 1.0 s operator length), a top mute of first arrivals, and an energy balance of each trace over the sub-basement window (2-15 s). Post-stack processing included only bandpass filtering and additional gain control (0.5 and 3.0 s operators in the sediments and sub-basement, respectively). The bandpass filtering was segmented over four windows with a 1 s overlap between each window for subsequent trace merging and decreasing high cut to account for the narrowing bandwidth of the sweep with time. This filtering was not significant since the sweep correlation acts as a time-varying bandpass filter; however, it was performed as a precaution against the possible generation of noise during the processing steps.  Chapter 5. PROCESSING AND INTERPRETATION OF REFLECTION DATA 182  5.4  Results  The stack is presented in Figure 5.7. Traces near 1 and 2 km appear muted above 0.7 s due to source gaps at these points (Figure 5.3a). Horizontal reflectors within the sedimentary section between 0.6 and 1.6 s are present. A close-up of the sedimentary section between 0 and 2 s is presented in Figure 5.8. This figure shows clearly that a number of the sedimentary reflectors appear to fade-out between 2.5-5.0 km. A statics problem was suspected for this gap in reflectivity; however, an examination of the shot and CMP gathers in this part of the line suggests that a source coupling or near-surface transmission problem has resulted in a reduced signal-to-noise ratio. In fact, the vibrators were operated at 30% full drive rate west of 5 km due to nearby pipelines (it is possible that the vibrators were returned to full rate west of 2.4 km where the sedimentary reflectors appear to strengthen again, but this has not been confirmed). This may explain the reduced reflectivity in certain parts of the section west of 5 km in both the sediments and sub-basement (Figure 5.7). The location of the intracrustal reflectors and crust-mantle boundary as determined from the refraction modelling of Line B are indicated in Figure 5.7. Note that from Figure 5.1, the intracrustal reflectors at 13, 18 and 29 km depth are imaged approximately 65, 30 and 15 km to the west of the location of the Vibroseis Une. Extrapolating the existence and depth of these reflectors to the east becomes more uncertain as the distance increases. Nevertheless, there is a weak horizontal event in the stack between the 13 km and 18 km wide-angle boundaries at 5.2 s between 5-7.5 km. This event may correspond to one of the wide-angle boundaries imaged at a sUghtly shallower (~1 km) or deeper (~5 km) depth further to the west. A prominent eastward-dipping event between 9.2-9.4 s and 5-7.5 km is observed in Figure 5.7. It is also prominent in the 5.94 km shot gather between 9.3-9.5 s (event F,  Q  i  n  »i Ox  Pa —1  O o  tq  Co  o1  CO  >  Q  -<  >>  —1 XI  >  rn i—  H  tq •n  <=o —i  m CO  rn  o  Pa s H lb.  !j  O  o S3  5 | g  t-< ts  O H o  to  Figure 5.7: Stack of the reflection section with an average horizontal exaggeration of 8. Distance in metres is indicated along the top. The position of the basement at 3.5 km, intracrustal reflectors at 13, 18 and 29 km and the crust-mantle boundary at 41 km from the refraction model are indicated.  > > »—' OO CO  1200  1800  2400  3000  3600  4200  4800  5400  6000  6600  7200  7800  8400  0.00  9000  J7J •XT  0.00  0.25  to to  o o  fe Co CO  g la to to to  2.00  2.00  o o to  Figure 5.8: Upper 2 s of the stack with an average horizontal exaggeration of one. Distance in metres is indicated along the top. The near-horizontal reflections originate within the 3.5-km-thick sedimentary cover.  to fe  a fe  o  H O  to >  oo  4-  Chapter 5. PROCESSING AND INTERPRETATION OF REFLECTION DATA 185  Figure 5.6). This corresponds well with one of the strongest intracrustal reflected phases present in the refraction data set (Figure 5.9b). Within the Line B model this boundary is at 29 km depth or 9.4 s zero-offset two-way travel time and represents the transition between the low-velocity-gradient mid crust and high-gradient lower crust with an associated velocity jump of 0.2 km/s. This boundary may correspond to the Riel discontinuity observed on near-vertical data at similar depths as a strong multi-cyclic event by Kanasewich et al. (1969) and Cumming and Chandra (1975) in southern Alberta and Ganley and Cumming (1974) ~400 km to the southeast. Within the refraction model this boundary is imaged 15 km to the west with a 4° westward dip, but extends horizontally beneath the Vibroseis line. The corresponding near-vertical reflection event has an apparent eastward dip of 14° using the velocity at 29 km depth from the Line B refraction model (6.7 km/s) and measured over a reflecting segment of 2.4 km horizontal extent (Figure 5.9a). After migration this dip would increase approximately 0.4° according to the constant-velocity depth migration equation s i n a = tana m  where a is the recorded dip angle and a r  m  r  (°-l)  is the migrated dip (Chun and Jacewitz, 1981).  This event extends westward across the reflection section where it reaches ~9.0 s at the extreme western end, although its strength drops significantly west of 5 km; this is the section where the vibrators operated at a reduced rate. A strong, 8° eastward-dipping event at 11.4 s between 5.7 and 7.4 km is observed in the reflection section and has no clear counterpart in the refraction data or Line B model. Between 12.0 and 12.8 s and 5.7 and 8.0 km is a strong series of eastward dipping events, most strong near 12.2 s. The apparent dip of the steepest dipping segment about 1.4 km long at 12.2 s is 11°. This event is visible on the west side of the 7.2 km shot gather at 12.2 s (event G, Figure 5.6). The bottom of this series of refections at  Chapter 5. PROCESSING AND INTERPRETATION OF REFLECTION DATA 186  Figure 5.9: An expanded view comparison of deep near-vertical and wide-angle reflection events, (a) Prominent events between 5-8 km and 8.5-13.5 s from the stack. The apparent dip of some short reflecting segments are indicated. Distance in metres is indicated along the top. (b) An intracrustal reflected phase (R) and the crust-mantle reflection ) from shot B l of the refraction data. Traces are scaled to a common-maximum amplitude and bandpass filtered between 2-12 Hz and source-receiver distance is indicated The near-vertical event between 9.2-9.5 s is interpreted to correspond to the wide-angle event labelled R. The base of the reflective zone in (a) between 12-13 s corresponds well to the zero-offset time of the crust-mantle boundary in the refraction model.  Chapter 5. PROCESSING AND INTERPRETATION OF REFLECTION DATA 187  ~12.8 s corresponds well with the 1.5° eastward-dipping refraction Moho at 41 km depth below the Vibroseis line. The character of the wide-angle P P from shots B l and B2 m  is diffuse and intermediate, respectively. Figure 5.9 is a close-up comparison of the 29km-depth intracrustal and crust-mantle events on the shot B l refraction profile and the reflection line. The diffuse character of the wide-angle P P suggests that the crustm  mantle boundary is a locally complex or disturbed transition zone. The 0.8 s series of near-vertical reflections terminating at 12.8 s suggests that locally the Moho is a complex, possibly layered transition zone of ~3 km thickness assuming a lower crustal velocity of 7 km/s. If the 11.4 s near-vertical event represents the top of the crust-mantle transition zone, a thickness of ~5 km is suggested. Within the Line B refraction model the lower crust immediately above the Moho has been modelled as a strong velocity-gradient layer _1  (0.06 s ) of 10 km thickness (Figure 5.2). The near-vertical reflection data suggest that locally, the lowermost 5 km of this layer is structurally complex. This may explain the diffuse and intermediate character of the wide-angle P P phases. The termination of m  lower crustal near-vertical reflectivity at a depth corresponding to the refraction Moho is commonly observed in continental areas (Mooney and Brocher, 1987). Steeply dipping, multi-cyclic events between 3 and 5 s and 5.0 and 8.4 km are present in the stack. A strong, relatively low frequency event (~20 Hz) dipping 70° westward intersects a weaker 70° eastward dipping event at approximately 3.1 s and 6.6 km. These events may be due to near-surface, near-vertical, out-of-plane structures of high impedance contrast since an apparent dip of 45° and a migrated dip of 90° is obtained using equation (5.1) if a reasonable near-surface velocity of 2 km/s is assumed. The shot gathers from 5.94 and 7.2 km in Figure 5.6 contain the westward dipping event between 2 and 5 s with an apparent velocity of 2.2 km/s (event B). The apparent transparency of the sub-basement within the 3.6-5.0 km column may be a result of the reduced Vibroseis source energy west of 5.0 km, and corresponds to  Chapter 5. PROCESSING AND INTERPRETATION OF REFLECTION DATA 188  a fading of some sedimentary reflectors within the same zone (Figure 5.8). The subbasement reflectivity appears to increase west of 3.6 km, as does the strength of the sedimentary reflectors above. However, the sub-basement reflectivity west of 3.6 km is not as sharp or focussed as it is east of 5 km, but appears randomly scattered. A relatively strong deep event at 13 s is present between 2.5-3.6 km and may represent the base of the crust further to the west of the 12.0-12.8 s series of reflectors. Since no parallel or cross lines have been analyzed, it is possible that the correspondence between the deep near-vertical events and the wide-angle intracrustal reflector and crust-mantle boundary is a coincidence and the near-vertical events may be an artifact of shallower out-of-plane structure. The character and apparent dip of these events suggests that they cannot be multiples of shallower events. At present, the correspondence of these events with the refraction boundaries suggests that it is more likely that their origin is within the lower crust rather than a result of shallower or out-of-plane structure or a processing artifact. Thus, this test of extended-listen processing has shown the potential for deep reflection imaging in the PRA region using industry-acquired, uncorrelated, upsweep Vibroseis data.  Chapter 6  DISCUSSION  6.1  Paleozoic Evolution of the Western C a n a d a Sedimentary Basin  The Phanerozoic development of the Western Canada Sedimentary Basin is dominated by thermal contraction of thinned, heated lithosphere and subsequent flexural downwarp of the lithosphere as a result of transverse (~N-S) loading along the western North American margin. The Paleozoic evolution of the Western Canada Sedimentary Basin, which is relevant to the uplift and inversion of the PRA, is often ascribed to a thermal subsidence driving mechanism related to the formation of a pre-Cordilleran passive continental margin beginning in the late Precambrian (Bond and Kominz, 1984; Bond et ai., 1985). The Mesozoic evolution of the Western Canada Sedimentary Basin is controlled by the Columbian and Laramide syn-orogenic phases of foreland basin development along the western North American margin (e.g., Beaumont, 1981). The PRA was structurally active in the Devonian and into the Mesozoic (Cant, 1988). Since it is normally assumed that thermal subsidence does not take place later than about 200 Myr after rifting (e.g., McKenzie, 1978), the Devonian and later subsidence of the Western Canada Sedimentary Basin is puzzling. In contrast with Bond and Kominz (1984) and Bond et al. (1985), evidence for an extensional event along the western margin at about 370 Ma has been presented by Gordey et al. (1987), Struik (1987) and Thompson et al. (1987). Struik (1987) suggests that this rifting is responsible for the creation of the ancient western North American passive margin and the Carboniferous  189  Chapter 6. DISCUSSION  190  and later subsidence. Beaumont et al. (1989) suggest to the contrary that the middle Paleozoic subsidence of the Western Canada Sedimentary Basin was a flexural response due to loading along the continental margin. The source of loading would be either continued progradation of the existing sedimentary wedge or due to thrusting of loads onto the continental margin during an early compressional phase of mountain building related to the Antler Orogeny further to the south. Thus, the development of the PRA can be viewed as either related to isostatic flexure and relaxation due to the emplacement and removal of loads on the lithosphere (e.g., Beaumont et al, 1989) or thermal uplift and decay due to extensional rifting (e.g., Cant, 1988). In a third class of model the crust beneath the PRA would have anomalous properties or structure at which vertical movements would be concentrated: a basement zone-of-weakness model. The results of the refraction interpretation will now be considered to test these opposing models,  6.2  6.2.1  Peace River Region Structural  Precambrian  Model  Features  Much of the interpreted crustal and upper mantle structure of the PRA region appears to be related to the pre-existing predominantly N-S trending Early Proterozoic cratonic structure (Figure 1.2). These relationships include trends in crustal thickness and P P m  character, upper mantle velocities, and the upper crustal velocity anomalies of Lines B and C and the basement velocity along Line D. Crustal thickness is greater than 40 km within the Taltson Magmatic Arc and generally less than or equal to 40 km in the Chinchaga Low and Wilson Island Terrane (Figure 4.40). P P character is generally m  sharp in the Taltson Magmatic Arc but intermediate or diffuse in the Wilson Island Terrane (Figure 4.40). The upper mantle velocity is generally greater than or equal to 8.2 km/s in the Taltson Magmatic Arc and the Chinchaga Low and less than or equal  Chapter 6. DISCUSSION  191  to 8.2 km/s in the Wilson Island Terrane (Figure 4.41). The upper crustal velocity anomalies at about 70 km and 50 km along Lines B and C, respectively, correlate well with localized aeromagnetic lows within the Wilson Island Terrane (Figure 4.35). The basement velocity along Line D is greater than 6 km/s within the Wilson Island Terrane and less than 5.7 km/s within the Chinchaga Low (Figure 4.35).  6.2.2  Peace River A r c h Features  Superimposed upon the cratonic structure are a number of subtle crustal features that may be associated with the Devonian axis of the PRA. These include the presence of high lower-crustal velocities at shallower depths beneath the arch, crustal thickening along the arch axis, relatively low amplitude P arrivals for propagation along the arch axis, inn  tracrustal reflecting horizons that appear to dip away from the arch axis and centre of the PRA region, and an anisotropic P P character beneath the arch. A shallowing of the m  6.6-7.0 km/s contours is observed in the Line A model (Figure 4.6) roughly beneath the arch axis (Figure 4.37). A similar shallowing of high velocities beneath the arch is not observed in the Line D model (Figure 4.16). However, the lack of a centre shot for Line D reduces the structural constraint in the lower crust here. Beneath the arch axis, a crustal thickening of 1.5 km and 1.0 km is observed in the Line A and D models, respectively (Figure 4.40). Additional evidence supporting crustal thickening beneath the arch is a short wavelength relative isostatic gravity low aligned with the axis of the PRA (Stephenson et ai., 1989). Figures 6.1a and 6.1b show the long wavelength (700-1200 km) and short wavelength (400-700 km) isostatic gravity anomalies of western Canada, respectively. There is a trade-off between the size and depth of a density anomaly when interpreting a feature of a particular wavelength. The long wavelength features are indicative of localized density anomalies at upper mantle depths or shallower anomalies of area! extent  Chapter 6. DISCUSSION  192  comparable to the wavelength of the gravity anomaly. Similarly, short wavelength features are indicative of localized lower crustal anomalies or shallower anomalies of greater areal extent. The long wavelength anomalies suggest that the crust in northern Alberta is thin with respect to regions north and south in the Western Canada Sedimentary Basin, a result consistent with the refraction interpretation when compared with other seismic studies in central and southern Alberta. However, the short wavelength anomalies show that there is also evidence in the gravity data for a slight amount of thickening roughly coincident with the arch axis. The long wavelength anomalies may alternatively suggest a relatively high density (velocity) lower crust, again consistent with the interpreted velocity structure. The short wavelength anomalies may alternatively suggest that the sediments overlying the PRA are of relatively low density, perhaps due to the younger than average age of these sediments as a result of an absence of Early Paleozoic strata and thickened Late Paleozoic sequence. The north-south P amplitude variations across n  the Line A fan shot suggest that energy propagating along the arch axis is strongly attenuated or scattered within the upper mantle or within the crust below Line A. This may suggest a structurally complex or disturbed upper mantle and/or crust-mantle boundary beneath the arch axis, or within the crust where it intersects Line A in the east. For the four models, approximately 60% of the intracrustal reflectors are near horizontal while about 40% dip away from the centre of the PRA axis or the centre of the PRA region. A few of the reflectors modelled as dipping are capable of satisfying the travel time data of the reflected phase equally well if replaced by a horizontal boundary due to the relatively large uncertainty of the arrival picks; the dip on these boundaries is required by the modelling of refracted rays. The weakly constrained trend of these dips may be a remnant expression of crustal uplift localized in the centre of the survey region. The two PP m  phases which sample the crust-mantle boundary at the north end of Line A are  diffuse, whereas in the two profiles along the eastern end of Line B they are sharp. These  Chapter 6. DISCUSSION  193  Figure 6.1: Isostatic gravity anomalies of western Canada. Contours in mGal. (a) Wavelengths of 700-1200 km. (b) Wavelengths of 400-700 km. Figure from Stephenson et al. (1988).  Chapter 6. DISCUSSION  194  phases sample the same region of the crust-mantle boundary in perpendicular directions (Figure 4.40). The anisotropic P P character suggests an E - W structural grain at the m  base of the crust beneath the arch axis in the east perhaps due to irregular topography or E - W faulting of the crust-mantle boundary. The inferred structural grain would act as a sharp reflector for E - W propagating energy but scatter N-S propagating energy. Further evidence for an irregular crust-mantle boundary beneath the PRA is the diffuse P P character of the shot B l profile along the arch axis in the west (Figure 4.40). m  6.3  Tectonic  Models  The presence of a crustal expression of the PRA would suggest that its uplift and subsequent inversion do not have a passive, flexural origin but are related to an anomalous thermal history during the Paleozoic Era. In particular, the Middle Paleozoic rifting event proposed by Struik (1987) suggests an origin of the PRA related to subtle thermoextensional rifting, perhaps as a Paleozoic failed-rift as suggested by Cant (1988). A number of the ideas associated with the models proposed by Cant (1988), Beaumont et al. (1989) and others will now be considered in more detail and the identification of specific epeirogenic processes suitable for the PRA will be addressed. Stelck et al. (1978) suggest that a mantle hot spot may have caused potassium metasomatism or enrichment in the basement, rocks reducing their density and thereby allowing them to maintain their position as a high. The drift of the continent over a hot spot would result in a lineal trace of altered basement rock. Presumably, the removal of the mantle heat source would cause the collapse of the high. Zolnai (1982) suggests that the vertical movements within the PRA region are due to a rejuvenation of Precambrian basement features as a result of tensional stresses set up within the North American continent as it was "pulled down" under the load of the Appalachian  Chapter 6. DISCUSSION  195  orogenic foreland thrust belt. Thus, no thermal domes or plumes are required (Zolnai, 1986). Cant (1988) points out that the uplift of domes and arches can result from rifting, associated with a hot spot or zone of crustal weakness. In addition, the orientation of the PRA into the continent in a manner similar to a failed third arm (e.g., Burke and Dewey, 1973) leads Cant (1988) to suggest that the arch originated as an incipient rift in the Early Paleozoic. This rift may have developed along pre-existing Precambrian basement faults, as described by Zolnai (1982): the rift oriented ~N-S and the failed rift oriented ~ E - W (Figure 6.2). Cant (1988) makes a comparison of the development of the PRA and the New Madrid Rift Complex of the Mississippi Embayment as described by Braile et al. (1986) on the basis of their similar, long and complex histories. The length and complexity of their development suggests the importance of the reactivation of failed rifts along established lines of weakness, perhaps due to changing intraplate stresses (e.g., Zolnai, 1982) as a result of the uplift and subsidence of adjacent arches and basins or the accretion of allocthonous terranes and lithospheric loading at the western end of the arch (Cant, 1988). Beaumont et al. (1989) have modelled the major Paleozoic sedimentary sequences of the Western Interior Basin in terms of flexural models driven by loading in the Williston Basin and on the passive continental margin. However, Beaumont et aJ. (1989) also point out that the Paleozoic evolution of the Western Canada Sedimentary Basin is least well understood and that their modelling of features such as the PRA can test whether certain concepts are viable but cannot distinguish between equally viable mechanisms. The results of this thesis suggest that the crust beneath the PRA is anomalous. In addition, there is little evidence to support the necessary loading for the flexural model in the Early Paleozoic along the western continental margin (Cant, 1988). Therefore, the purely flexural model is considered less likely. The passive model proposed by Zolnai (1982) is a basement zone-of-weakness model with reactivation a result of a regional stress  Chapter 6. DISCUSSION  Figure 6.2: Major Proterozoic and Paleozoic rifts of North America with the Peace River Arch added. Adapted from Burke (1980).  Chapter 6. DISCUSSION  197  field. This model may have played a role in the later development of the PRA, however, there are no disruptions in basement trends marked by distinctive aeromagnetic and/or gravity anomalies as might be expected for a profound basement feature. The anomalous lower crust beneath the arch suggests a thermo-extensional model is appropriate. The Stelck et al. (1978) thermal model is one possibility. However, no E-W trending basement feature as required is observed in the gravity, magnetic or seismic data. The E - W trend of the PRA perpendicular to the Paleozoic western margin supports the failed-rift model as proposed by Cant (1988) and best fits the interpreted seismic structure. The relatively weak features of the velocity models related to the PRA and the absence of a gravity and magnetic signature, which is observed in other rift zones (e.g., southern Alberta and Mississippi Embayment), indicates that the PRA is not a profound crustal feature and may suggest minimal extension or rift failure at an early stage. The only constraint on the timing of the initiating thermo-extensional event from the refraction data is that it is post-Early Proterozoic. It has been suggested that the arch existed as early as the latest Precambrian or as late as the Devonian. A thermal origin for the Middle Devonian PRA implies a thermal subsidence origin for the subsequent collapse of the arch by the Permian. This supports a Devonian initiating event given the thermal diffusion properties of the continental lithosphere (McKenzie, 1978). It now remains to work out the appropriate driving mechanism, i.e., a passive or active response. 6.3.1  Rift-like Features  A number of the interpreted structural features beneath the PRA are similar to those observed in continental rift zones. The relatively high lower-crustal velocities (> 7.2 km/s, Figure 4.38), and more localized shallowing of high-velocity lower-crustal material (Figure 4.37) is reported, for example, by Mooney et al. (1983) for the Mississippi Embayment and by Catchings and Mooney (1988) for the Columbia Plateau (Figure 6.3). Morgan  Chapter 6.  DISCUSSION  S -20  30  l  198  DISTANCE  80  130  l  l 2.5-fi.n  o  180  l  (KM)  280  230  l  l  330  l  N 380  o _  6.0-7.0 DZ o  > 7.0  LU Q o —  o  LINE  A  8.2  8.4  b  PEACE RIVER A R C H Figure 6.3: A comparison of the generalized velocity structure of known continental rifts where seismic refraction measurements have been made with the PRA. Anomalous high-velocity (> 7.0 km/s) lower crustal layer is indicated by dotted pattern, (a) Known continental rifts. Figure from Catchings and Mooney (1988). (b) Line A from the PRA region which crosses the arch at roughly 240 km model distance.  Chapter 6. DISCUSSION  199  and Ramberg (1987) present a detailed summary of geophysical properties commonly observed in "paleorifts" and most are consistent with observations of the PRA; in particular, a high-velocity lower crust. In addition, high upper mantle velocities (8.4 km/s) have been interpreted beneath the Limagnegraben (Him and Perrier, 1974), the Rhinegraben (Zucca, 1984) and the Columbia Plateau (Catchings and Mooney, 1988). The high-velocity lower crustal layer beneath the Mississippi Embayment is interpreted by Mooney et ai. (1983) to represent a crust-mantle mix formed by the injection of mantle material into the lower crust. Furlong and Fountain (1986) point out that a mafic crustal underplate is common in continental rift and extension zones and suggest that, in general, high lower crustal velocities (> 7.1 km/s) may indicate that continental crust thinned by extension has been replaced by mafic intrusions and underplating from the mantle. Strong anisotropy does not appear to be the cause of the observed high upper mantle velocity in the northeast of the PRA region since similar high values of 8.45 and 8.35 km/s are measured along the perpendicular Lines A and B, respectively. A mechanism for generating high upper mantle velocities is the alteration of mineralogy due to the removal of basalt that may leave a higher-velocity upper mantle (Catchings and Mooney, 1988) as part of the thermal differentiation process associated with the intrusions and/or underplating from the mantle. If a balancing of material between the lower crust and uppermost mantle is required by this process, a thickening of the crust is implied: a thickening of no more than 2 km is observed beneath the PRA. Three additional features of the interpreted structural model may provide weak evidence for, or are consistent with, extensional processes. First, the relatively thin average crustal thickness and the observed thinning of the crust towards the west, most pronounced north of the PRA, may be related to the stretching of the lithosphere associated with the rifting event which established the passive continental margin in the late Precambrian or Paleozoic. The amount of thinning is 4 km from the in-line modelling and  Chapter 6. DISCUSSION  200  as much as 10 km from the Line B fan shot north of the arch. This significant crustal thinning toward the west is not supported by the gravity data which may suggest strong E-W lateral variations in velocity and/or density north of the arch in addition to possible crustal thinning. The presence of numerous sub-horizontal intracrustal reflectors may be a remnant of ductile crustal extension, such as proposed by Valasek et ai. (1987) in the Basin and Range, and associated with a rifting event. However, intracrustal reflectors are not unusual within the continental crust and so their correlation with extension is tentative. Finally, upper mantle boundaries at depths of 74, 69 and 56 km are observed along Lines A, B and C, respectively, and all within the Taltson Magmatic Arc. The extent of these boundaries imaged in each case suggests that they are laterally continuous over at least 35 to 50 km. Refraction results indicating boundaries within the upper mantle are not uncommon (e.g., Fuchs, 1986; Mueller and Ansorge, 1986) and indeed Mereu et ai. (1977) describe a high velocity (8.5-8.8 km/s) upper mantle layer at a depth of ~60 km along a refraction profile approximately 200 km south of the PRA region. Fuchs (1986) suggests that strong horizontal shear deformation in the subcrustal lithosphere could produce nearly horizontal layers within the upper mantle. Another possibility is that intrusions from the lower lithosphere or asthenosphere of the appropriate density could form sills within the upper mantle (e.g., McKenzie, 1984). The upper mantle boundaries have been modelled as increases in velocity of 0.2 to 0.3 km/s. An alternative model is a decrease in velocity implying a low-velocity zone within the upper mantle similar to that proposed by Bennett et ai. (1975) in the southern Rocky Mountain Trench at a depth of 63 km or ~8 km below the Moho. This model is less likely since a decrease in velocity to less than 7 km/s cannot generate more than half the amplitude of the increasing-velocity model. A low-velocity model would suggest a layer of partial melt within the upper mantle (e.g., Catchings and Mooney, 1988).  Chapter 6. DISCUSSION  6.3.2  201  Passive Versus A c t i v e Response  Given that the failed-rift model accounts best for the interpreted crustal structure, can the distinction between a passive or active response be determined? A passive rift would proceed with lithospheric extension first, followed by lower crustal and upper mantle alteration through magmatic intrusions and/or underplating. An active response would most likely imply an underlying mantle plume or thermal anomaly that would heat, alter and finally extend the lithosphere. The anisotropic P P character in the northeast of m  the PRA region (Figure 4.40) may suggest that the Moho is disrupted by a series of approximately E - W trending faults and thereby imply that the lower crust and upper mantle responded to possible extension in a brittle manner. This supports the passive rift model in which there was heating and modification of the lower crust subsequent to extension. The dominant extensional orientation may have been approximately E - W along the western continental margin during the Paleozoic (e.g., Struik, 1987) with only mild extension in the P R A region as a failed third arm (e.g., Cant, 1988). The plume-generated rift model (Burke and Dewey, 1973) can also be supported by the observations. A number of structural patterns appear focused about the northeast of the PRA region, including a shallowing of high-velocity lower crust (Figures 4.37 and 4.38), thickening of the crust (Figure 4.40) and high upper mantle velocities (Figure 4.41). In addition, the intracrustal reflectors may dip away from the centre of the PRA region, not exactly coincident with the anomalous structure in the northeast, as opposed to the axis of the arch. This apparent focusing of anomalous structure is tentative given the limited areal extent of the refraction data further to the north, northeast and east, but may support a profound mantle heat source as a rift-initiating mechanism. Similar anomalous structure may also be focused at the west end of Line B, perhaps due to the motion of the plate at an irregular rate over the mantle plume, similar to that proposed  Chapter 6. DISCUSSION  202  in the Stelck et al. (1978) model. If E-W compression was dominant along the western margin in the Paleozoic, it is possible that the plumes associated with the active model may correspond to the rising magmas of a continental arc at an ancient subduction zone. This model is quite distinct from the extensional model, but again involves a thermal process that would be expected to alter the lower crust and upper mantle structure and/or composition. Thus, although the PRA may be associated with a failed rift, other thermal models are possible and a dominant passive or active response is not resolvable with the geophysical data alone, but requires further understanding of the Paleozoic development of the Western Canada Sedimentary Basin and western margin as a whole.  6.4 6.4.1  Supplementary Studies Gravity Modelling  The results of the gravity modelling suggest that most of the lateral velocity anomalies in the refraction models do not correspond to density anomalies since no gravity signatures are observed. As a result, regions of the initial density models obtained from the seismic modelling were adjusted primarily within zones of anomalous velocity structure (i.e., upwarps or downwarps). The magnitude of these velocity anomalies is ±0.1 to 0.3 km/s when compared with lateral averages. Densities were lowered by as much as 0.05 g/cm  3  within velocity upwarps and vice versa resulting in all models having essentially laterally homogeneous density structure.  One possible explanation for this observed velocity-  density behaviour is lateral variations in composition; in particular, a compositionally heterogeneous crust in isostatic equilibrium. Zones of relative low and high velocity could correspond to enriched felsic and mafic rock, respectively. The aeromagnetic data show that the basement high-velocity zones along Line B and C and the low-velocity northern end of Line D correspond to magnetic lows. Thus, a direct correspondence between  Chapter 6. DISCUSSION  203  composition and these velocity anomalies is not supported by the magnetic data since felsic and mafic rocks generally produce magnetic lows and highs, respectively; however, the opposite magnetic-composition relationship is observed in the exposed Canadian Shield (Zietz, 1980). If the heterogeneities which give rise to the velocity upwarps (and downwarps) are due to intrusions, as might be associated with a thermo-extensional process, their emplacement into the crust could be in a manner similar to that described by McKenzie (1984): intruding material will rise upward until it reaches a zone of equidensity and then moves laterally to form sills. If the intrusions are thermally driven, subsequent cooling may alter their density. Thus, the compositional model is probably overly-simplistic since compositional variations are likely to give rise to density variations. Another possible model is one incorporating transverse anisotropy. Anomalous zones of high and low velocity within the crust could represent regions of relatively fast or slow horizontal velocity, respectively. Therefore, these anomalies may not be compositional variations but rather structural variations whereby high-velocity propagation is aligned with the layering that is often observed in deep reflection data, particularly in the lower crust.  Zones of relative low velocity may suggest an absence of layering or may be  due to thermal or structural alteration or disruption. If a layer within the lower crust is produced by underplating or intrusions, the manner in which material is added or injected into the crust will depend upon the densities involved and whether the process is a passive upwelling associated with extension or forced intrusion driven by a density instability. A compositionally homogeneous layer may result but may be disrupted by or consist entirely of directional dyke- and sill-like features giving rise to lateral variations in horizontal velocity. Of course, a combination of the composition and anisotropy models may be involved; one may be prevalent in the upper crust, the other in the lower crust. Figure 6.4 illustrates a possible model for the lower crustal upwarp and Moho downwarp below the PRA along Line A. Figure 6.4a is the simplified interpreted velocity structure and Figure  Chapter 6.  204  DISCUSSION  s  N  DISTANCE (KM)  -20  30  80  130  180  230  280  o  o  CL.cn LU O CD .  330  380  LINE A _  \-V  p  =  6.8-  -MOHO a  o . in  -20  30  80  130  180  230  280  330  380  •—'CM  rr CL.cn LU  p = 2.8  a a. XT  o in  Figure 6.4: Qualitative model of transverse anisotropy required by gravity modelling of Line A. (a) Simplified lower crustal and Moho velocity structure, (b) Laterally homogeneous density structure and transverse anisotropic body required by gravity modelling based on qualitative arguments. The densities indicated are averages for the mid and lower crust and upper mantle.  Chapter 6. DISCUSSION  205  6.4b is the density structure required by the gravity modelling: laterally homogeneous density structure with a zone of transverse anisotropy superimposed. Note that the Moho would be "flattened" in this interpretation since the P P m  a n  d P ray paths travel more n  vertically than the P ray paths through the anisotropic body. The short wavelength c  isostatic gravity anomaly (Figure 6.1b) suggests that the Moho is slightly depressed below the arch, but could also imply a low-density anomaly within the sediments or lower crust. Layering in the lower crust often observed in deep reflection data may be produced by sill intrusions or underplating (e.g., McKenzie, 1984; Furlong and Fountain, 1986) and suggests that transverse anisotropy may be important when considering the consistency between wide-angle and near-vertical seismic data. The agreement between the twoway travel time of deep refraction and reflection events presented in chapter 5 suggests that strong anisotropy is not present in this location. However, the amount of tranverse anisotropy required to produce the prominent lower crustal "upwarp" beneath the PRA along Line A (Figure 6.4) would result in a discrepancy between the wide-angle and nearvertical Moho of about 0.35 s, an amount within the uncertainties of the Moho depth estimates (~1 km). A compositional anomaly may be substituted for the anisotropic body in Figure 6.4b, but the Moho would remain depressed as shown in Figure 6.4a. The anisotropy model best explains the observed velocity-density behaviour in the PRA region since the regional Bouguer gravity data is relatively flat and compositional variations would most likely involve velocity and density variations. The effect of anisotropy may be particularly important for continental rift models that often have an upwarped highvelocity lower crust (Figure 6.3). This idea could be quantitatively tested by anisotropic ray tracing, although a data set containing a range of ray paths from horizontal to nearvertical in the region of interest would be required to achieve an unambiguous isotropicversus-anisotropic interpretation.  Chapter 6. DISCUSSION 6.4.2  206  Other Studies  The result of the 5-wave analysis suggests that the average crustal Poisson's ratio is 0.25. There is weak evidence to suggest that a is greater than 0.25 in the northeast of the PRA region, as determined from the four S S phases at the intersection of Lines A and B, and m  less than 0.25 elsewhere. Assuming the relatively high a in the northeast is significant, two explanations are offered. First, a localized compositional anomaly would imply a more mafic crust in the northeast. Second, assuming the same gross crustal composition throughout the PRA region, a slightly less rigid crust in the northeast is implied. The former model is consistent with the lower crustal up warp along Line A beneath the PRA representing a mafic intrusion or underplate. The latter model suggests that a slight thermal anomaly may exist beneath the eastern end of the PRA. This model is less likely given the high lower crustal and upper mantle velocities in the northeast (Black and Braile, 1982). In either case, a localized cr anomaly in the northeast supports the active plume-generated rift model discussed earlier. An alternative explanation is that the relatively high Poisson's ratio in the northeast is an expression of the Proterozoic basement terrane of the eastern PRA region, the Taltson Magmatic Arc. The one-dimensional Q model of the crust and upper mantle in the PRA region suggests a Q within the crust increasing with depth from about 500 to 1000 and an upper mantle Q of 1000. This model is generally consistent with a number of crustal Q studies. From an amplitude analysis of refraction data, Mereu et al. (1977) estimated a Q of less than 500 in the upper 9 km of the crust, greater than 500 in the mid 35 km of the crust, and greater than 500 13 km below the Moho. They also claim to resolve a zone of Q of less than 200 in a 20 km band about the Moho at 51 km depth on the basis of the relatively low observed amplitude of reflections from an upper mantle boundary (/?*). This result is not consistent with observations of P" in the PRASE data set.  Chapter 6. DISCUSSION  207  Clowes and Kanasewich (1970) obtained an average cmstal Q of 1500 from a spectral analysis of near-vertical reflection data in southern Alberta. Singh and Herrmann (1983) regionalized crustal coda Q in the continental United States and obtained values of 800 to 1300 in the mid-continent. An analysis of spectral ratios in a manner similar to that presented in section 4.9 for refraction data from southern Japan yielded a crustal Q of between 1000 and 2000 (Hashizume, 1979). The main conclusion to be drawn from the Q model of the PRA region is that there is no evidence for a regional, thick, low-<5 layer within the crust or uppermost mantle that would suggest a heated or partial melt zone. The main results of the analysis of the reflection data presented in chapter 5 are (1) a consistency between the zero-offset two-way travel time of a prominent intracrustal reflector and the crust-mantle boundary of the Line B refraction model with deep events in the extended-listen Vibroseis data set, and (2) a band of reflectors 1 to 1.5 s in duration terminating at the refraction Moho. The intracrustal reflector of the refraction model represents the transition from the mid crust to the high-velocity lower crust. The near vertical reflections are multi-cyclic over approximately 0.35 s suggesting a layered transition zone of about 1 km thickness. This zone may represent alternating thin layers of granitic mid crustal material and more mafic lower crustal material and is perhaps a result of an intrusion or underplating process described by McKenzie (1984) and Furlong and Fountain (1986). Alternatively, this boundary may have formed as a result of fusion and differentiation processes that separated the weakened mid crust from the more basic lower crust (Pavlenkova, 1988). The crust-mantle boundary of the refraction model is represented by a step discontinuity of 0.6 km/s at the base of a high vertical-velocitygradient lower crust while the reflection data suggest a complex layered transition of 3 to 5 km thickness. Similarly, Edel et ai. (1975) and Fuchs et ai. (1981) report a high gradient lower crust over a transitional crust-mantle boundary of 4 to 5 km thickness localized  Chapter 6. DISCUSSION  208  beneath the Rhinegraben. They propose that the transitional Moho is a zone of crustmantle interaction formed by intrusions from the upper mantle, differentiation, or phase transformation. The banded appearance of lower crustal reflections can best be explained by alternating layers of high- and low-velocity material, each layer typically 100-200 m thick (Fuchs, 1969; Clowes and Kanasewich, 1970). Common geological interpretations of such zones include a series of acidic lower crustal and more basic mantle material, metasediments and metavolcanics, partially molten and solid matter, sill-like mantle intrusions, metamorphic layering, and anisotropy combined with the action of shear stresses (Meissner, 1967; 1973; Clowes and Kanasewich, 1970; Berry and Mair, 1977; Smithson and Brown, 1977). Given that wide-angle refraction and near-vertical reflection data are sensitive to different velocity and structural features, a consistency between the Line B refraction model and the short reflection Une is observed.  6.5  Conclusions and Recommendations  It is proposed that the PRA originated as a Paleozoic failed-rift, although a dominant passive or active response cannot be differentiated.  This conclusion is based upon an  assumed Middle Paleozoic extensional event along the western North American margin and the foUowing features of the P R A structural model: 1. The approximate east-west orientation of the arch perpendicular to the Paleozoic western margin. 2. The regional high-velocity lower crust and shaUowing of high-velocity lower crust beneath the eastern end of the P R A axis. 3. A sUght crustal thickening beneath the arch axis. 4. An anisotropic P P character beneath the eastern end of the arch axis. m  Chapter 6. DISCUSSION  209  5. The along-axis low-amplitude P ,arrivals. n  These features are subtle suggesting that the arch is a less profound crustal feature than other known continental rifts, perhaps an indication of minimal extension or early rift failure. The weak seismic expression of the PRA is perhaps not unexpected and is consistent with the absence of a gravity and magnetic signature. An explanation of the weak crustal expression is the inferred Paleozoic extension as compared to the generally more recent extension of the documented rifts shown in Figure 6.3a. The evolution of these young stronger crustal features with time is unknown but may involve a "smearing" or smoothing of lower crustal structure through ductile flow or other adjustment mechanisms due to intraplate stresses. The later Paleozoic and Mesozoic evolution of the PRA may have been a result of, or modified by, changing intraplate stresses or flexure due to lithospheric loading. Further crustal research in the PRA region should address the following issues: • The weak crust and upper mantle features associated with the PRA need to be better resolved to confidently propose a particular epeirogenic process. • Localized structural anomalies in the northeast PRA region need to be better resolved by data further to the north, northeast and east. Similarly, data in the northwest is required to study the apparent westward crustal thinning as revealed by the Line B fan shot. • A better understanding of the relationship between the observed P P character m  and the structure of the crust-mantle boundary is needed. • S is not observed when the corresponding S S is recorded. g  m  • A puzzling seismic-gravity relationship exists in which the most prominent structural features of the refraction models do not appear as gravity anomalies.  Chapter 6. DISCUSSION  210  To address these issues, high spatial density refraction and deep reflection data both within the PRA region and the surrounding area is needed. A data set containing a dense coverage of rays with paths from horizontal to near-vertical could be used in the assessment of possible transverse anisotropy and to obtain the velocity structure using an inverse or tomographic technique. A portion of these data should be concentrated in the northeast to better resolve the three-dimensional velocity structure, particularly the shallowing of high-velocity lower crust beneath the arch. Increased'shot density along refraction lines would better resolve the size and magnitude of the upper crustal velocity anomalies inferred from the current refraction data and thereby resolve or confirm the apparent seismic-gravity conflict. Three component recording may reveal or explain the "missing" 5-phases of the PRASE data set and shear wave splitting is an indicator of anisotropy. The near-vertical reflection character of the crust-mantle boundary could be analyzed to assess a possible correlation with the observed sharp and diffuse wideangle P P phases andfinite-differencemodelling of two- and three-dimensional structures m  may aid in relating the near-vertical and wide-angle reflection characters. The results of chapter 5 suggest that the recorrelation of existing industry Vibroseis data could be used for deep reflection imaging both within the PRA region and throughout the Western Canada Sedimentary Basin where available. The potential data volume is considerable and it could provide a wealth of information at relatively little expense. Geodynamic modelling of the Paleozoic evolution of the PRA region and surrounding area involving thermo-extensional processes should be performed to check the validity of the tectonic models proposed in this thesis. This may be especially important since a thermal evolutionary process could have significant implications for hydrocarbon maturation. To fully understand the development of the PRA, it will be necessary to learn more about the Paleozoic evolution of the Western Canada Sedimentary Basin and the western North American margin as a whole.  References  Aki, K. and P.G. Richards, 1980, Quantitative Seismology, W.H. Freeman Co., San Francisco. Artyushkov, E.V., A.E. Shlesinger, and A.L. Yanshin, 1980, The origin of vertical crustal movements within lithospheric plates, in Dynamics of Plate Interiors, Geodyn. Ser., vol. 1, A.W. Bally, P.L. Bender, T.R. McGetchin, and R.I. Walcott, ed., AGU, Washington, D.C., 37-51. Bachu, S., 1988, Analysis of heat transfer processes and the geothermal pattern in the Alberta Basin, Canada, J. Geophys. Res., 93, 7767-7781. Barton, P.J., 1986, The relationship between seismic velocity and density in the continental crust—a useful constraint?, Geophys. J. R. Astron. Soc, 87, 195-208. Beaumont, C , 1981, Foreland basins, Geophys. J. R. Astron. Soc, 65, 291-329. Beaumont, C , C E . Keen, and R. Boutilier, 1982, A comparison of foreland and rift margin sedimentary basins, Philos. Trans. R. Soc. London, Ser. A, 305, 295-317. Beaumont, C , G. Quinlan, and G.S. Stockmal, 1989, The evolution of the western interior basin: causes, consequences and unsolved problems, in Geodynamic Evolution of the Western Interior Basin, Geol. Assoc. Can., Spec. Paper, in press. Bennett, G.T., R.M. Clowes, and R.M. Ellis, 1975, A seismic refraction survey along the southern Rocky Mountain Trench, Canada, Bull. Seism. Soc. Am., 65, 37-54. Berry, M.J., and D.A. Forsyth, 1975, Structure of the Canadian Cordillera from seismic refraction and other data, Can. J. Earth Sci., 12, 182-208. Berry, M.J., and J.A. Mair, 1977, The nature of the earth's crust in Canada, in The Earth's Crust, Geophys. Mono. 20, J.G. Heacock, ed., AGU, Washington, D.C., 319-348. Birch, F., 1964, Density and compostion of mantle and core, J. Geophys. Res., 69, 4377-4388.  211  References  212  Black, P.R., and L.W. Braile, 1982, P velocity and cooling of the continental lithosphere, J. Geophys. Res., 87, 10557-10568. n  Bond, G.C., and M.A. Kominz, 1984, Construction of tectonic subsidence curves for the early Paleozoic miogeocline, southern Canadian Rocky Mountains: implications for subsidence mechanisms, age of breakup, and crustal thinning, Geol. Soc. Am. Bull., 95, 155-173. Bond, G.C., N. Christie-Blick, M.A. Kominz, and W.J. Devlin, 1985, An early Cambrian rift to post-rift transition in the Cordillera of western North America, Nature, 315, 742-746. Bott, M.H.P., 1982, The Interior of the Earth, its Structure, Constitution and Evolution, E. Arnold, London. Braile, L.W., 1977, Interpretation of crustal velocity gradients and Q structure using amplitude-corrected seismic refraction profiles, in The Earth's Crust, Geophys. Mono. 20, J.G. Heacock, ed., AGU, Washington, D.C., 427-439. Braile, L.W., and C.S. Chiang, 1986, The continental Mohorovicic discontinuity: results from near-vertical and wide-angle seismic reflection studies, in Reflection Seismology: A Global Perspective, Geodyn. Ser., vol. 13, M. Barazangi and L. Brown, ed., AGU, Washington, D.C., 257-272. Braile, L.W., W.J. Hinze, G.R. Keller, E.G. Lidiak, and J.L. Sexton, 1986, Tectonic development of the New Madrid Rift Complex, Mississippi Embayment, North America, Tectonophysics, 131, 1-21. Braile, L.W., and R.B. Smith, 1975, Guide to the interpretation of crustal refraction profiles, Geophys. J. R. Astron. Soc, 40, 145-176. Burke, K., 1980, Intracontinental rifts and aulacogens, in Continental Tectonics, Studies in Geophysics, National Academy of Sciences, Washington, D.C., 42-49. Burke, K., and J.F. Dewey, 1973, Plume-generated triple junctions: key indicators in applying plate tectonics to old rocks, J. Geol, 81, 406-433. Cant, D.J., 1988, Regional structure and development of the Peace River Arch, Alberta: a Paleozoic failed-rift system?, JBuli. Can. Petrol. Geol., 36, 284-295. Cassell, B.R., 1982, A method for calculating synthetic seismograms in laterally varying media, Geophys. J. R. Astron. Soc, 69, 339-354.  References  213  Catchings, R.D., and W.D. Mooney, 1988, Crustal structure of the Columbia Plateau: evidence for continental rifting, J. Geophys. Res., 93, 459-474. Cerveny, V., 1966, On dynamic properties of reflected and head waves in the n-layered Earth's crust, Geophys. J. R. Astron. Soc, 11, 139-147. Cerveny, V., 1981, Computation of geometrical spreading by dynamic ray tracing, Stanford Exploration Project, Rept. No. 28, J. Claerbout, ed., Stanford Univ., Dept. of Geophysics, 49-60. Cerveny, V., 1985, Gaussian beam synthetic seismograms, J. Geophys., 58, 44-72. Cerveny, V. and F. Hron, 1980, The ray series method and dynamic ray tracing system for three-dimensional inhomogeneous media, Bull. Seism. Soc. Am., 70, 47-77. Cerveny, V., I. Molotkov, and I. Psencik, 1977, Ray Method in Seismology, Charles Univ. Press, Prague. Cerveny, V. and I. Psencik, 1981, 2-D seismic ray package, Res. Rept., Institute of Geophysics, Charles Univ., Prague. Cerveny, V. and R. Ravindra, 1971, Theory of Seismic Head waves, Univ. of Toronto Press, Toronto. Chandra, N.N., and G.L. Cumming, 1972, Seismic refraction studies in western Canada, Can. J. Earth Sci., 9, 1099-1109. Chapman, C H . and R. Drummond, 1982, Body-wave seismograms in inhomogeneous media using Maslov asymptotic theory, Bull. Seism. Soc. Am., 72, S277-S317. Chiang, C.S., and L.W. Braile, 1984, An example of two-dimensional synthetic seismogram modeling, Bull. Seism. Soc. Am., 74, 509-519. Christensen, N.I., and D.L. Szymanski, 1988, Origin of reflections from the Brevard fault zone, J. Geophys. Res., 93, 1087-1102. Chun, J.H., and C A . Jacewitz, 1981, Fundamentals of frequency domain migration, Geophysics, 46, 717-733. Clowes, R.M., E.R. Kanasewich, and G.L. Cumming, 1968, Deep crustal seismic reflections at near-vertical incidence, Geophysics, 33, 441-451. Clowes, R.M., and E.R. Kanasewich, 1970, Seismic attenuation and the nature of reflecting horizons within the crust, J. Geophys. Res., 75, 6693-6705.  References  214  Cooper, L., and D. Steinberg, 1970, Introduction to Methods of Optimization, W.B. Saunders Co., Toronto. Coppens, F., 1985, First arrival picking on common-offset trace collections for automatic estimation of static corrections, Geophys. Prosp., 33, 1212-1231. Cormier, V.F. and P. Spudich, 1984, Amplification of ground motion and waveform complexity in fault zones: examples from the San Andreas and Calaveras Faults, Geophys. J. R. Astron. Soc, 79, 135-152. Crossley, D.J., 1987, The development of geophysics software for PC's (abstract), paper SW1-2, IUGG, Abstracts, vol. 1, Vancouver, 378. Cumming, G.L., and N.N. Chandra, 1975, Further studies of reflections from the deep crust in southern Alberta, Can. J. Earth. Sci., 12, 539-557. Cudrak, C.F., 1988, Shallow crustal structure of the Endeavour Ridge segment, Juan de Fuca Ridge, from a detailed seismic refraction survey, M.Sc. thesis, Univ. of British Columbia, Vancouver. deMille, G., 1958, Pre-Mississippian history of the Peace River Arch, J. Alta. Soc. Petrol. Geol., 6, 61-68. Dorman, L.M., 1969, Anelasticity and the spectra of body waves, J. Geophys. Res., 73, 3877-3307. Dziewonski, A., S. Bloch, and M. Landisman, 1969, A technique for the analysis of transient seismic signals, Bull. Seism. Soc. Am., 59, 427-444. Edel, J.B., K. Fuchs, C. Gelbke, and C. Prodehl, 1975, Deep structure of the southern Rhinegraben area from seismic refraction investigations, J. Geophys., 41, 333-356. Ellis, R.M., and P.W. Basham, 1968, Crustal characteristics from short-period P waves, Bull. Seism. Soc. Am., 58, 1681-1700. Ellis, R.M., Z. Hajnal, and R.A. Stephenson, 1986, PRASE 85: Crustal seismic refraction profiles in the Peace River Arch region, northwestern Alberta and northeastern British Columbia, Open File Rept. 1317, Geol. Surv. Can. Finckh, P., J. Ansorge, St. Mueller, and Chr. Sprecher, 1984, Deep crustal reflections from a Vibroseis survey in northern Switzerland, Tectonophysics, 109, 1-14. Finckh, P., W. Frei, B. Fuller, R. Johnson, St. Mueller, S. Smithson, and Chr. Sprecher,  1986, in Reflection Seismology: A Global Perspective, Geodyn. Ser., vol. 14, M.  References  215  Barazangi and L. Brown, ed., AGU, Washington, D.C., 43-54. Fordjor, C.K., J.S. Bell, and D.I. Gough, 1983, Breakouts in Alberta and stress in the North American plate, Can. J. Earth Sci., 20, 1445-1455. Forsyth, D.A., M.J. Berry, and R.M. Ellis, 1974, A refraction survey across the Canadian Cordillera at 54° N, Can. J. Earth Sci., 11, 533-548. Fuchs, K., 1969, On the properties of deep crustal reflectors, Z. Geophys., 35, 133-149. Fuchs, K., 1975, Synthetic seismograms of PS-reflections from transition zones computed with the reflectivity method, J. Geophys., 41, 445-462. Fuchs, K., 1986, Reflections from the subcrustal lithosphere, in Reflection Seismology: The Continental Crust, Geodyn. Ser., vol. 14, M . Barazangi and L. Brown, ed., AGU, Washington, D.C., 67-76. Fuchs, K., K.P. Bonjer, and C. Prodehl, 1981, The continental rift system of the Rhinegraben— structure, physical properties and dynamical processes, Tectonophysics, 73, 79-90. Fuchs, K., and G. Midler, 1971, Computation of synthetic seismograms with the reflectivity method and comparison with observations, Geophys. J. R. Astron. Soc, 23, 417-433. Fuis, G.S., W.D. Mooney, J.H. Healy, G.A. McMechan, and W.J. Lutter, 1984, A seismic refraction survey of the Imperial Valley region, California, J. Geophys. Res., 89, 1165-1189. Furlong, K.P., and D.M. Fountain, 1986, Continental crustal underplating: thermal considerations and seismic-petrologic consequences, J. Geophys. Res., 91, 8285-8294. Gajewski, D., K. Fuchs, C. Prodehl, and R. Stangl, 1987, Anomalous disappearance of 5 -phases near the continental Moho (abstract), paper U7-13, IUGG, Abstracts, vol. 1, Vancouver, 58. n  Ganley, D.C., and G.L. Cumming, 1974, A seismic reflection model of the crust near Edmonton, Alberta, Can. J. Earth Sci., 11, 101-109. Gardner, G.H.F., L.W. Gardner, and A.R. Gregory, 1974, Formation velocity and density— the diagnostic basics for stratigraphic traps, Geophysics, 39, 770-780. Geological Survey of Canada, 1986, Magnetic anomaly map, NO-10-M (Fort Nelson), N O - l l - M (Peace River).  References  216  Gibson, B.S., and A.R. Levander, 1988, Lower crustal reflectivity patterns in wide-angle seismic recordings, Geophys. Res. Lett., 15, 617-620. Gordey, S.P., J.G. Abbott, D.J. Tempelman-Kluit, and H. Gabrielse, 1987, "Antler" elastics in the Canadian Cordillera, Geology, 15, 103-107. Grant, F.S., and G.F. West, 1965, Interpretation Theory in AppHed Geophysics, McGrawHill Book Co., Toronto. Grad, M., and U. Luosto, 1987, Seismic models of the crust of the Baltic shield along the SVEKA profile in Finland, Ann. Geophys., 5B, 639-650. Hall, J., and M . Ali, 1985, Shear waves in a seismic survey of Lewisian basement: an extra control on lithological variation and porosity, J. Geol. Soc. London, 142, 677-688. Hashizume, M., 1979, Q of the crust beneath southwestern Honshu, Japan, derived from explosion seismic waves, Phys. Earth Planet. Inter., 20, 25-32. Haxby, W.F., D.L. Turcotte, and J.M. Bird, 1976, Thermal and mechanical evolution of the Michigan Basin, Tectonophysics, 36, 57-75. Hill, D.P., 1971, Velocity gradients and anelasticity from crustal body wave amplitudes, J. Geophys. Res., 76, 3309-3325. Hill, D.P., W.D. Mooney, G. Fuis, and A. Walter, 1982, The whispering-gallery phase in deep sedimentary basins (abstract), Earthquake Notes, 53, 26. Hinze, W.J., L.W. Braile, G.R. Keller, and E.G. Lidiak, 1980, Models for midcontinent tectonism, in Continental Tectonics, Studies in Geophysics, National Academy of Sciences, Washington, D.C., 73-83. Hirn, A., and G. Perrier, 1974, Deep seismic soundings in the Limagnegraben, in Approaches to Taphrogenesis, H. Lilies and K. Fuchs, ed., Schweizerbart, Stuttgart, 563-567. Hoffman, P.F., 1987, Continental transform tectonics: Great Slave Lake shear zone (ca. 1.9 Ga), northwest Canada, Geology, 15, 785-788. Hoffman, P.F., 1988, United plates of America, the birth of a craton: Early Proterozoic assembly and growth of Laurentia, Ann. Rev. Earth Planet. Sci., 16, 543-603. Holbrook, W.S., D. Gajewski, A. Krammer, and C. Prodehl, 1988, An interpretation of wide-angle compressional and shear wave data in southeast Germany: Poisson's  References  217  ratio and penological implications, J. Geophys. Res., 93, 12081-12106. Hwang, L.J., and W.D. Mooney, 1986, Velocity and Q structure of the Great Valley, California, based on synthetic seismogram modeling of seismic refraction data, Buli. Seism. Soc. Am., 76, 1053-1067. Jacob, A.W.B., and D.C. Booth, 1977, Observations of PS reflections from the Moho, J. Geophys., 43, 687-692. Jacobson, R.S., G.G. Shor, and L.M. Dorman, 1981, Linear inversion of body wave data— Part II: attenuation versus depth using spectral ratios, Geophysics, 46, 152-162. Jacobson, R.S., G.G. Shor, and M. Bee, 1984, A comparison of velocity and attenuation between the Nicobar and Bengal Deep Sea Fans, J. Geophys. Res., 89, 6181-6196. Kanamori, H., and D.L. Anderson, 1977, Importance of physical dispersion in surface wave and free oscillation problems: review, Rev. Geophys. Space Phys., 15, 105112. Kanasewich, E.R., 1981, Time Sequence Analysis in Geophysics, Univ. of Alberta Press, Edmonton. Kanasewich, E.R., R.M. Clowes, and C H . McCloughan, 1969, A buried Precambrian rift in western Canada, Tectonophysics, 8, 513-527. Karageorgi, E., D. Okaya, T. McEvilly, and P. Malin, 1988, Vibroseis uncorrelated seismogram decomposition and filtering in the frequency-uncorrelated time domain to suppress narrow-band Vibroseis correlation artifacts (abstract), Eos, 69, 1327. Klemperer, S.L., T.A. Hauge, E . C Hauser, J.E. Oliver, and C J . Potter, 1986, The Moho in the northern Basin and Range province, Nevada, along the COCORP 40°N seismic-reflection transect, Geol. Soc. Am. Bull., 97, 603-618. Knopoff, L., 1964, Q, Rev. Geophys., 2, 625-660. Kulhnger, B., and C E . Lund, 1986, A preliminary interpretation of S-wave traveltimes from Fennolora data, Tectonophysics, 126, 375-388. Lam, H.L., and F.W. Jones, 1984, Geothermal gradients of Alberta in western Canada, Geothermics, 13, 181-192. Lavoie, D.H., 1958, The Peace River Arch during Mississippian and Permo-Pennsylvanian time, J. Alta. Soc. Petrol. Geol, 6, 69-74.  References  218  Ludwig, W.J., J.E. Nafe, and C.L. Drake, 1970, Seismic refraction, in The Sea, vol. 4, pt. 1, A.E. Maxwell, ed., New York, 53-84. Luetgert, J.H., C.E. Mann, and S.L. Klemperer, 1987, Wide-angle deep crustal reflections in the northern Appalachians, Geophys. J. R. Astron. Soc, 89, 183-188. Majorowicz, J.A., and A.M. Jessop, 1981, Regional heat flow patterns in the Western Canadian Sedimentary Basin, Tectonophysics, 74, 209-238. Marks, L.W., 1980, Computational topics in ray seismology, Ph.D. thesis, Univ. of Alberta, Edmonton. McCrossan, R.G., and R.P. Glaister (Eds.), 1964, Geological History of Western Canada, Alia. Soc. Petrol. Geol, Calgary. McGetchin, T.R., K.C. Burke, G.A. Thompson, and R.A. Young, 1980, Mode and mechanisms of plateau uplifts, in Dynamics of Plate Interiors, Geodyn. Ser., vol. 1, A.W. Bally, P.L. Bender, T.R. McGetchin, and R.L Walcott, ed., AGU, Washington, D.C., 99-110. McKenzie, D., 1978, Some remarks on the development of sedimentary basins, Earth Planet. Sci. Lett., 40, 25-32. McKenzie, D., 1984, A possible mechanism for epeirogenic uplit, Nature, 307, 616-618. McMechan, G.A. and W.D. Mooney, 1980, Asymptotic ray theory and synthetic seismograms for laterally varying structures: theory and application to the Imperial Valley, California, Bull. Seism. Soc. Am., 70, 2021-2035. Meissner, R., 1967, Exploring deep interfaces by seismic wide angle measurements, Geophys. Prosp., 15, 598-617. Meissner, R., 1973, The 'Moho' as a transition zone, Geophys. Surv., 1, 195-216. Meissner, R., 1986, The Continental Crust: a Geophysical Approach, Academic Press, Inc., Toronto. Menke, W., 1983, A formula for the apparent attenuation of acoustic waves in randomly layered media, Geophys. J. R. Astron. Soc, 75, 541-544. Mereu, R.F., S.C. Majumdar, and R.E. White, 1977, The structure of the crust and upper mantle under the highest ranges of the Canadian Rockies from a seismic refraction survey, Can. J. Earth Sci., 14, 196-208.  References  219  Mereu, R.F., D. Wang, 0. Kuhn, D.A. Forsyth, A.G. Green, P. Morel, G.G.R. Buchbinder, D. Crossley, E. Schwarz, R. duBerger, C. Brooks, and R. Clowes, 1986, The 1982 COCRUST seismic experiment across the Ottawa-Bonnechere graben and Grenville Front in Ontario and Quebec, Geophys. J. R. Astron. Soc, 84, 491-514. Michel, H.K., 1986, Prehminary results of the 1985 Peace River Arch crustal refraction survey, B.Sc. thesis, Univ. of Western Ontario, London. Milne, W.G., G.C. Rogers, R.P. Riddihough, G.A. McMechan, and R.D. Hyndman, 1978, Seismicity of western Canada, Can. J. Earth Sci., 15, 1170-1193. Monger, J.W.H., and R.A. Price, 1979, Geodynamic evolution of the Canadian Cordillera— progress and problems, Can. J. Earth. Sci., 16, 770-791. Mooney, W.D., M.C. Andrews, A. Ginzburg, D.A. Peters, and R.M. Hamilton, 1983, Crustal structure of the northern Mississippi Embayment and a comparison with other continental rift zones, Tectonophysics, 94, 327-348. Mooney, W.D., and T.M. Brocher, 1987, Coincident seismic reflection/refraction studies of the continental lithosphere: a global review, Rev. Geophys., 25, 723-742. Morel-a-l'Huissier, P., A.G. Green, and C.J. Pike, 1987, Crustal refraction surveys across the Trans-Hudson orogen/Williston Basin of south central Canada, J. Geophys. Res., 92, 6403-6420. Morgan, P., 1982, Heat flow in rift zones, in Continental and Oceanic Rifts, Geodyn. Ser., vol. 8, G. Palmason, ed., AGU, Washington, D.C., 107-122. Morgan, P., and B.H. Baker, 1983, Introduction—Processes of continental rifting, Tectono-  physics, 94, 1-10. Morgan, P., and LB. Ramberg, 1987, Physical changes in the lithosphere associated with thermal relaxation after rifting, Tectonophysics, 143, 1-11. Mueller, St., and J. Ansorge, 1986, Long-range seismic refraction profiles in Europe, in  Reflection Seismology: A Global Perspective, Geodyn. Ser., vol. 13, M. Barazan and L. Brown, ed., AGU, Washington, D.C., 167-182. Miiller, G., 1984, Efficient calculation of Gaussian beam seismograms for two-dimensional inhomogeneous media, Geophys. J. R. Astron. Soc, 79, 153-166. Ojo, S.B., and R.F. Mereu, 1986, The effect of random velocity functions on the travel times and amplitudes of seismic waves, Geophys. J. R. Astron. Soc, 84, 607-618.  References  220  Okaya, D.A., 1986, Seismic profiling of the lower crust: Dixie Valley, Nevada, in Reflection Seismology: The Continental Crust, Geodyn. Ser., vol. 14, M . Barazangi and L. Brown, ed., AGU, Washington, D.C., 269-279. Okaya, D.A., T.L. Henyey, and E.G. Frost, 1988, Deep crustal imaging in southwest U.S. using industry seismic profiles (abstract), in Seismic Probing of Continents and their Margins, J.C. Dooley, ed., Canberra, Australia, 110-111. Oliver, J., and S. Kaufman, 1977, Complexities of the deep basement from seismic reflection profiling, in Tie Earth's Crust, Geophys. Mono. 20, J.G. Heacock, ed., AGU, Washington, D.C., 243-253. Pavlenkova, N.I., 1988, The nature of seismic boundaries in the continental lithosphere, Tectonophysics, 154, 211-225. Phinney, R.A., 1986, A seismic cross section of the New England Appalachians: the orogen exposed, in Reflection Seismology: The Continental Crust, Geodyn. Ser., vol 14, M. Barazangi and L. Brown, ed., AGU, Washington, D.C., 157-172. Podruski, J.A., 1988, Contrasting character of the Peace River and Sweetgrass Arches, Western Canada Sedimentary Basin, Geosci. Can., 15, 94-97. Porter, J.W., R.A. Price, and R.G. McCrossan, 1982, The western Canada sedimentary basin, Philos. Trans. R. Soc. London, Ser. A., 305, 169-192. Press. F., and R. Siever, 1986, Earth, W.H. Freeman Co., New York. Procter, R.M., and G. Macauley, 1968, Mississippian of western Canada and Williston Basin, Am. Assoc. Petrol. Geol. Bull., 52, 1956-1968. Quinlan, G.M., and C. Beaumont, 1984, Appalachian thrusting, lithospheric flexure, and the Paleozoic stratigraphy of the eastern interior of North America, Can. J. Earth Sci., 21, 973-996. Razavy, M . and Lenoach, B., 1986, Reciprocity principle and the approximate solution of the wave equation, Bull. Seism. Soc. Am., 76, 1776-1789. Richards, P.G., 1971, An elasticity theorem for heterogeneous media, with an example of body wave dispersion in the earth, Geophys. J. R. Astron. Soc, 22, 453-472. Richards, P.G., and W. Menke, 1983, The apparent attenuation of a scattering medium, Bull. Seism. Soc. Am., 73, 1005-1021. Ronen, J., and J.F. Claerbout, 1985, Surface-consistent  residual statics estimation by  References  221  stack-power maximization, Geophysics, 50, 2759-2767. Ross, G.M., R.R. Parrish, S.A. Bowring, and A.J. Tankard, 1988, Tectonics of the Canadian Shield in the Alberta subsurface (abstract), Geol. Assoc. Can., Program with abstracts, 13, A106. Rotstein, Y., and P. Trachtman, 1986, On the use of seismic reflection surveys from oil exploration in deep crustal studies, Phys. Earth Planet. Inter., 43, 236-243. Rudkin, R.A., 1964, Lower Cretaceous, in Geological History of Western Canada, Alta. Soc. Petrol. Geol, R.G. McCrossan and R.P. Glaister, ed., Calgary, 156-168. Sandmeier, K.-J., and F. Wenzel, 1986, Synthetic seismograms for a complex crustal model, Geophys. Res. Lett., 13, 22-25. Seriff, A.J., and W.H. Kim, 1970, The effect of harmonic distortion in the use of vibratory surface sources, Geophysics, 35, 234-246. Sheriff, R.E., and L.P. Geldart, 1983, Exploration Seismology, vol. 2: Data-processing and interpretation, Cambridge Univ. Press. Sikabonyi, L.A., and W.J. Rodgers, 1959, Paleozoic tectonics and sedimentation in the northern half of the west Canadian Basin, J. Alta. Soc. Petrol. Geol., 7, 193-216. Singh, S., and R.B. Herrmann, 1983, Regionalization of crustal coda Q in the continental United States, J. Geophys. Res., 88, 527-538. Skeels, D.C., 1947, Ambiguity in gravity interpretation, Geophysics, 12, 43-56. Smithson, S.B., and S.K. Brown, 1977, A model for lower continental crust, Earth. Planet. Sci. Lett., 35, 134-144. Somerville, P.G., and R.M. Ellis, 1972, P-coda evidence for a layer of anomalous velocity in the crust beneath Leduc, Alberta, Can. J. Earth Sci., 9, 845-856. Spence, G.D., 1983, RAYAMP: An algorithm for tracing rays and calculating amplitudes in laterally varying media, Program Documentation, Univ. of British Columbia, Vancouver. Spence, G.D., R.M. Clowes, and R.M. Ellis, 1985, Seismic structure across the active subduction zone of western Canada, J. Geophys. Res., 90, 6754-6772. Spence, G.D., K.P. Whittall, and R.M. Clowes, 1984, Practical synthetic seismograms for laterally varying media calculated by asymptotic ray theory, Bull. Seism. Soc.  References  222  Am., 74, 1209-1223. Stelck, C.R., R.A. Burwash, and D.R. Stelck, 1978, The Vreeland High: a Cordilleran expression of the Peace River Arch, Bull. Can. Petrol. Geol., 26, 87-104. Stephenson, R.A., C A . Zelt, R.M. Ellis, Z. Hajnal, P. Morel-a-l'Huissier, R.F. Mereu, D.J. Northey, G.F. West, and E.R. Kanasewich, 1989, Crust and upper mantle structure and origin of the Peace River Arch, BuL". Can. Petrol. Geol., in press. Stott, D.F., 1984, Cretaceous sequences of the Foothills of the Canadian Rocky Moun-  tains, in The Mesozoic of Middle North America, Can. Soc. Petrol. Geol., Mem. 9, D.F. Stott and D.J. Glass, ed., 85-107. Struik, L . C , 1987, The ancient western North American margin: an Alpine rift model for the east-central Canadian Cordillera, Geol. Surv. Can, paper 87-15. Talwani, M., J.L. Worzel, and M. Landisman, 1959, Rapid gravity computations for twodimensional bodies with application to the Mendocino submarine fracture zone, J. Geophys. Res., 64, 49-59. Tarkov, A.P., LP. Basula, V.G. Generalov, A.L Dubyansky, and V.V. Chernykh, Composite travel times of seismic waves and general velocity models of the Voronezh Shield crust and upper mantle, Geophys. J. R. Astron. Soc, 67, 137-143. Thompson, B., E. Mercier, and C. Roots, 1987, Extension and its influence on Canadian Cordilleran passive-margin evolution, in Continental Extensional Tectonics, Geol. Soc Spec. Publ. No. 28, M.P. Coward, J.F. Dewey and P.L. Hancock, ed., 409-417. U. S. Geological Survey, 1972, Tectonic map of North America, compiled by P.B. King, Washington, D . C Valasek, P.A., R.B. Hawman, R.A. Johnson, and S.B. Smithson, 1987, Nature of the lower crust and Moho in eastern Nevada from "wide-angle" reflection measurements, Geophys. Res. Lett, 14, 1111-1114. Walcott, R.L, 1970, Flexural rigidity, thickness, and viscosity of the lithosphere, J. Geophys. Res., 75, 3941-3954. Wenzel, F., and K.-J. Sandmeier, 1988, Reflectivity method for dipping layers, J. Geophys. Res., 93, 15046-15056. Wetmiller, R.J., 1986, Earthquakes near Rocky Mountain House, Alberta, and their  References  223  relationship to gas production facilities, Can. J. Earth Sci., 23, 172-181. Whittall, K.P., and R.M. Clowes, 1979, A simple, efficient method for the calculation of travel times and ray paths in laterally inhomogeneous media, J. Can. Soc. Expl. Geophys., 15, 21-29. Williams, G.K., 1958, Influence of the Peace River Arch on Mesozoic strata, J. Aita. Soc. Petrol. Geol, 6, 74-81. Wilson, J.T., 1966, Did the Atlantic close and then re-open?, Nature, 211, 676-681. Young, G.B., and L.W. Braile, 1976, A computer program for the application of Zoeppritz's amplitude equations and Knott's energy equations, Bull. Seism. Soc. Am., 66, 1881-1885. Zelt, C.A., J.J. Drew, M.J. Yedlin, and R.M. Ellis, 1987, Picking noisy refraction data using semblance supplemented by a Monte Carlo procedure and spectral balancing, Bull. Seism. Soc. Am., 77, 942-957. Zietz, I., 1980, Exploration of the continental crust using aeromagnetic data, in Continental Tectonics, Studies in Geophysics, National Academy of Sciences, Washington, D.C., 127-138. Zolnai, G., 1982, Late Mississippian-Pennsylvanian orogenic movements of the west Canadian platform and adjacent areas—their role in sedimentation and hydrocarbon accumulations (abstract), Am. Assoc. Petrol. Geol., 66, 646-647. Zolnai, G., 1986, Les aulacogenes du continent nord-americain, Bull. Soc. Geol. France, 8, 809-818. Zucca, J.J., 1984, The crustal structure of the southern Rhinegraben from re-interpretation of seismic refraction data, J. Geophys., 55, 13-22. Zucca, J.J.; G.S. Fuis, B. Milkereit, W.D. Mooney, and R.D. Catchings, 1986, Crustal structure of northeastern California, J. Geophys. Res., 91, 7359-7382.  

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-0052504/manifest

Comment

Related Items