Department of Geophysics and Astronomy, University of British Columbia, Vancouver, Canada
Department of Geophysics, Stanford University, Stanford, California
The accretion of terranes to continental margins is an important process by which large amounts of new continental crust are created. The Canadian Cordillera, with a recognizable history of terrane accretion dating back to the Middle Jurassic (180-185 Ma), is a natural laboratory for the study of this and other processes associated with lithospheric evolution and growth of the continent. Recent multidisciplinary studies within the Lithoprobe Southern Cordillera transect have significantly increased our understanding of the nature, timing, and dynamic evolution of the Cordillera [Clowes et al., 1992].
Seismic refraction/wide-angle reflection data, hereinafter referred to as "refraction data", were acquired by Lithoprobe in two stages: the 1989 Southern Cordillera Refraction Experiment (SCoRE '89), located in southwestern British Columbia, and the 1990 recording program (SCoRE '90), located in southeastern British Columbia with one profile in the Coast Belt (Figure 1, inset). The general objectives of the surveys were (1) to provide quantitative values for velocity variation with depth; (2) to map laterally varying velocity structure; and (3) to map the topography of prominent velocity discontinuities such as the crust-mantle boundary (Moho). An important aspect of the 1989 survey was the triangular recording geometry and abundance of broadside (fan) shots, providing, in addition to two-dimensional (2-D) coverage along three long profiles, three-dimensional (3-D) coverage of a region centered on the Coast-Intermontane Belt (CB-IMB) boundary. This paper presents a new development for the inversion of wide-angle seismic data to determine 3-D velocity structure and depth to reflecting interfaces, and an application of the resulting algorithm to the SCoRE '89 data set.
The Canadian Cordillera is divisible into five physiographically distinct belts [Monger and Price, 1979] (Figure 1, inset). The Foreland, Intermontane and Insular Belts comprise unmetamorphosed and low-grade metamorphic rock and form the Cordilleran superstructure. The Omineca and Coast Belts are major regional, possibly tectonic, welts in which the metamorphic and plutonic infrastructure of the Cordillera is exposed. They comprise high-grade metamorphic and granitic rocks. The Intermontane and Insular Belts, and flanking parts of the Omineca and Coast Belts, contain late Paleozoic and early Mesozoic strata and coeval granitic rocks that form the Intermontane and Insular superterranes, respectively (Figure 1, inset). Rocks of the two superterranes probably formed in ocean basins as oceanic crust and island arcs and later were added (accreted) to the ancient continental margin. In the south, several small terranes of the Coast Belt comprising oceanic and island arc rocks are juxtaposed between the two superterranes.
The times at which the two superterranes were accreted correspond to the times at which the initiation of widespread Mesozoic metamorphism and crustal shortening occurred in the two tectonic welts. This observation, together with the rough correspondence between the locations of the two welts and the eastern boundary regions of the two superterranes (Figure 1, inset), forms the basis of the popular collision model for the origin of the Coast and Omineca Belts [Monger et al., 1982], although other models have been proposed [Gehrels and Saleeby, 1985; van der Heyden, 1992; Monger et al., 1994]. Detailed descriptions of the tectonic evolution of the Canadian Cordillera are given by Monger and Price [1979], Coney et al. [1980], Gabrielse and Yorath [1992], Oldow et al. [1989], Clowes et al. [1992], and Monger et al. [1995].
Sweeney et al. [1991] present a comprehensive listing of geophysical studies in the Cordillera and the most significant results to the early 1980s. However, our understanding of the geophysics and geology of the southern Canadian Cordillera has improved enormously as a result of recent multidisciplinary studies throughout Lithoprobe's Southern Cordillera Transect [Clowes et al., 1992]. Seismic refraction studies in the southwestern Cordillera have provided well- constrained velocity models for the perimeter of the SCoRE '89 triangle [Zelt et al., 1992, 1993; McLean and Spence, 1994] and along-strike through the Coast Belt [O'Leary et al., 1993]. Clowes et al. [1995] provide a synthesis of these and related studies. In terms of general characteristics such as crustal thickness and correspondence with large-scale geological features, the refraction results correlate well with those interpreted from high-quality Lithoprobe multichannel reflection data [Cook et al., 1991, 1992; Varsek et al., 1993]. The crust is highly reflective at all levels, structurally complex and not divisible into transparent upper and reflective lower layers.
Throughout the southern Canadian Cordillera, the Bouguer gravity field is generally low, with values rising significantly around the Coast-Insular belt boundary; the magnetic field is variable. Cook et al. [1995] discuss the relationship of gravity and magnetic variations across the region to specific geological characteristics. Electrical conductivity in the lower crust and/or upper mantle beneath much of the southern Cordillera is high in comparison with the craton [Gough, 1986]. Jones et al. [1992b] present a study of regional variations in crustal conductivity across the southern Cordillera from magnetotelluric studies; Jones and Gough [1995] provide a comprehensive review of the results of recent magnetotelluric experiments in the region. Heat flow throughout most of the southern Cordillera east of the Garibaldi volcanic belt (Figure 1) is relatively high (in comparison with the craton) [Davis and Lewis, 1984; Lewis et al., 1988, 1992]. In the region of the volcanic belt, it is high and variable, but a dramatic decrease in heat flow values occurs immediately west of the belt.
The Lithoprobe 1989 Southern Cordillera Refraction Experiment was designed to image the crustal and upper mantle velocity structure of the southwestern Canadian Cordillera in both two and three dimensions. To meet this objective required both in-line recording and an extensive fan-shot program. To this end, a roughly triangular arrangement of shots and receivers centered over the Fraser River Fault system was employed (Figure 1). This geometry allowed for in-line recording along three strategically located and approximately linear profiles (lines 1, 2, and 3), and a regular distribution of fan shots into these profiles, providing a broad range of offsets and azimuths as required for travel time tomography. In addition to the three main profiles, three short profiles (4, 5, and 6) were recorded in the interior of the triangle.
Data were recorded during three deployments of recording instruments. Average receiver spacing was 1.2 km along lines 1 and 3 and 1.6 km along line 2. Shots were detonated at seventeen locations around the triangle, one site interior to the triangle, and two sites on Vancouver Island. Figure 1 shows the distribution of large (1800 kg), medium (800 kg), and small (200 kg) shots. Average shot spacing was 55 km. In general, all of the 800- and 1800-kg shots were recorded at all receiver sites around and interior to the triangle; that is, these were used for both in-line recording and as fan shots. The exceptions are shot point (SP) 5 which was recorded along line 1 only, and SP 8 which was recorded along lines 2 and 3 only. The 200-kg shots were used for in-line recording only. See Zelt et al. [1990] for further details.
Figure 2 shows the location of shot points and receiver sites used in the analysis. Data from the following shot- receiver combinations are used: (1) in-line data recorded along lines 1, 2, and 3; the line 1 and 3 data sets are described by Zelt et al. [1992, 1993]; McLean and Spence [1994] present the line 2 data set; the data recorded by receivers to the west of SP 5 and data from SPs 8 and 9 (Figure 1) are not used in the analysis; (2) fan data from the seven large shots (SPs 1, 3, 5, 4, 7, 12, and 17) recorded along lines 1, 2 and 3; (3) data recorded along three short interior lines (4, 5, and 6); and (4) in-line data recorded along the southeastern end of line 10 (Figure 1, inset) of SCoRE '90 [O'Leary et al., 1993].
The arrival times of three phases are inverted to determine 3-D crustal and upper mantle structure: (1) first arrival refractions through the crust, Pg; (2) reflections from the Moho, PmP; and (3) first arrival refractions through the upper mantle, Pn. Not all phases are present or recognizable on all sections because of low signal-to-noise ratio or inappropriate offsets. Other arrivals not confidently identified as Pg, PmP, or Pn were not used. Fan data picks were made from various display formats including record sections plotted as a function of shot-receiver azimuth and shot-receiver offset. The latter, which can be used for oblique shot-receiver geometries for far- offset traces, is more useful for trace-to-trace correlation. Azimuthal plots of the data with normal moveout corrections applied were also used. This format is useful for trace-to-trace correlation of the PmP reflection phase. Table 1 displays the number of data recorded along each line.
Elevation corrections were applied to the observed travel times to strip off the transit time through a thin near-surface layer. This layer, as imaged on the 2-D in-line models for lines 1-3 and 10 is typically 0.5-2 km thick, comprising highly variable and relatively low velocities compared to underlying structure [Zelt et al., 1992, 1993; McLean and Spence, 1994; O'Leary et al., 1993]. Corrections were made by calculating the transit time through the near-surface layer assuming a fixed ray angle at all shot and receiver sites. The corrections are approximately equivalent to lowering the shots and receivers to the base of this layer, which forms the datum for the model.
Data quality is, in general, very good with abundant first arrivals and Moho reflections present on most sections. Figures 3-5 show representative fan record sections from shot points located at (1) a corner of the triangle (SP 1 recorded along line 3; Figure 3); (2) a line midpoint (SP 12 recorded along line 1; Figure 4); and (3) the center of the triangle (SP 7 recorded along line 2; Figure 5). Small-scale lateral variations in the arrival times of the picked phases are static shifts due primarily to topographic or near-surface velocity variations mostly accounted for by elevation corrections described above.
Corner shots. First arrivals across the three corner shot record sections are Pn. The arrivals are emergent across most of the line 1 section making trace-to-trace correlation of arrivals difficult. Pn is more impulsive on line 3 (Figure 3) but is obscured by a high noise level across much of the section. PmP is clearly evident on lines 1 and 2, but is not observed on line 3 probably because of the large minimum offset (310 km) and high noise level.
Line midpoint shots. First arrivals across the six sections are a mixture of Pg and Pn, both of which are present on all sections (e.g., Figure 4). The distinction between the two is not definitive in the vicinity of the Pg - Pn crossover point, which on the in-line sections typically occurs at 170- 180 km. In general, first arrivals greater than about 7 s at offsets greater than 170 km are considered to be Pn. PmP arrivals are observed on all sections.
Center shot. This shot produced the best signal-to-noise ratio of all the fan shots. First arrivals are clearly visible across all three sections, comprising Pg and Pn (Figure 5). PmP arrivals are present on all the sections.
A number of methods for imaging 3-D velocity variations exist; see Thurber and Aki [1987] for a brief review. Most of these have been applied to the study of earthquake arrival times from either teleseismic or local events. Popular methods of solution include direct damped least squares [Aki et al., 1977] and algebraic reconstruction [Humphreys and Clayton, 1988]. Phillips and Fehler [1991] compare the popular methods. An important distinction between different techniques is the method used to represent the Earth's velocity structure. The type of parameterization used, along with the amount, type and quality of data, governs the type of structure, and its resolution, that a tomographic inversion procedure will produce [Zhao et al., 1992]. The type of 3-D ray tracing scheme employed is also important, as this too can affect the accuracy of the results.
Most tomographic methods offer advantages and disadvantages which can make the choice of a scheme difficult. Drawbacks can include the use of (1) only first arrival travel times; rarely are both refracted and reflected arrivals used; (2) inaccurate forward modeling (ray tracing) techniques, or accurate methods which are computationally very slow; (3) model parameterizations that limit the resolving power of the data; usually this is governed by the choice of inversion scheme; and (4) inversion schemes which can be computationally very slow for large problems. Our development obviates these disadvantages.
The tomographic inversion procedure of Hole [1992], with modifications to allow the calculation and inversion of reflection travel times [Hole and Zelt, 1995] and layer stripping, is used to determine 3-D velocity structure. This method is computationally efficient and allows for a densely sampled velocity model. Key elements of the procedure include (1) forward modeling of travel times using a 3-D finite difference algorithm; and (2) a choice of inversion parameterization that greatly simplifies and speeds up the linearized inversion by taking advantage of the linear independence of ideal rays.
The usual approach of linearizing the nonlinear relation between travel time t and slowness u is followed, leading to the approximation
Parameterizing the model in terms of functions defined along ray paths for the inversion is not a convenient choice for the forward modeling of travel times. Hole [1992] uses a modified version of the 3-D finite difference algorithm of Vidale [1990] to compute first arrival travel times. The modifications eliminate inaccuracy problems associated with large, sharp velocity contrasts [Hole and Zelt, 1995]. The velocity model is defined at a set of uniformly spaced grid nodes in 3-D space. First arrival travel times are calculated at each node of the grid using finite difference operators based on the eikonal equation of ray tracing. The operators use the average slowness across a grid cell and are equivalent to a plane wave approximation across a grid cell [Hole, 1992]. The travel time at an arbitrarily located receiver is determined using tri-linear interpolation of the travel times at the eight surrounding grid nodes. First arrival ray paths are found by following the steepest gradient in travel time from the receiver to the shot point [Vidale, 1988]. The slowness perturbations in (2) are defined along the ray paths and must be transformed into the same parameterization (i.e., a 3-D grid) used in the forward modeling. Hole [1992] calculates the slowness perturbation at each grid node as the average of all the du terms in some volume surrounding the grid node. This gridding procedure also serves as a smoothing operator which stabilizes the inversion. Additional smoothing of the gridded du terms is performed by applying a 3-D moving average (MA) filter. Synthetic examples demonstrating the capabilities of the procedure and methods for estimating the resolution of resultant models are given by Hole [1992]. An application of the procedure to data from the Queen Charlotte Basin is given by Hole et al. [1993].
Upper crustal velocities are most strongly constrained by shallow phases (e.g., Pg), and it is therefore reasonable to apply the slowness perturbation in (2) to the entire Pg ray path. For deeper phases (e.g., Pn), distributing the slowness perturbation along the entire ray path is generally not desirable because the deeper phases have larger picking errors and additionally are sensitive to a larger degree of uncertainty in deeper velocities due to poorer ray coverage in the lower crust. Including deeper phases in the inversion for shallow structure can actually degrade the solution.
We modify the inversion procedure to include layer stripping, thereby allowing the velocity perturbations of deep phases to be applied only along the part of the ray path that most strongly affects the arrival time. Rather than applying the travel time residual to the entire path length, lj in (2) is taken to be the length of the jth ray path beneath a specified 3- D surface.
The inversion procedure for first arrivals (including layer stripping) is also applied to reflected phases. Thus, for example, PmP travel times can be inverted for lower crustal velocities and/or depth to Moho.
The 3-D finite difference algorithm of Vidale [1990] calculates first arrival travel times to each node of a 3-D grid. Hole and Zelt [1995] present a method for calculating reflection travel times that uses this algorithm and allows the reflector model to vary smoothly with depth, increasing the accuracy compared to methods which require the reflector to lie at a grid node in depth [e.g., Moser, 1991; Hole et al., 1993; Lecomte and Hamran, 1993]. The method relies on an extension to the finite difference algorithm that allows the computations to begin from an arbitrary wall of a previously computed 3-D travel time field. Only the first reflected arrival times are found. Reflection times are computed simultaneously for all receivers, requiring only two applications of the finite difference algorithm for each shot. In comparison, other techniques [e.g., Podvin and Lecomte, 1991; Matsuoka and Ezaka, 1992] based on the reciprocity principal require one pass of the finite difference algorithm for each shot and receiver. The reflector depth model is parameterized by a 2-D grid of depth nodes with the same uniform node spacing as the velocity model.
Hole and Zelt [1995] use the downgoing wave field to analytically calculate the reflected time at one or more grid nodes immediately above, and below, the reflector. The analytically calculated times are used to propagate the reflected time field upward. This method is based on the assumption that the incident and reflected wave fields are locally planar and the reflector surface is locally planar. The algorithm is accurate for reflectors with dips of less than about 35 deg (in 3-D). However, vertical faults are detected and handled properly [Hole and Zelt, 1995].
Here, a simple inversion scheme, which, like the inversion procedure for velocity eliminates the need for matrix inversion, is used to determine reflector depth from reflection times. The procedure is analogous to that described by Hole et al. [1992] for the inversion of refracted first arrivals to determine interface structure. In both procedures a set of depth perturbations are gridded to obtain a reflector (interface) perturbation model, and both assume the velocity model is accurate. Riahi and Juhlin [1994] present a similar inversion scheme; however, their forward modeling procedure is based on the reciprocity principal [Matsuoka and Ezaka, 1992] and their reflector model is required to lie at grid nodes in depth.
For reflections, the partial derivative relating a small change in reflector depth dz to the change in travel time dt is
A set of M data will give M depth perturbations Dzk(x,y), k = 1 to M, which are gridded to obtain a reflector perturbation model Dzij. The choice of gridding scheme is important, especially when there are large gaps in reflection point coverage, or no coverage around the edge of the reflector model. The LaPlace interpolation method, which fits the depth perturbations to a function Dz(x,y) satisfying the equation Ñ2[Dz(x,y)] = 0, is a good choice because it introduces minimal structure between the depth perturbations. For real data with errors, closely spaced reflection sampling points associated with significantly different depth perturbations will produce spikes in the gridded surface. Spikes are removed and the gridded surface is smoothed by applying a 2-D MA filter. Since the relationship between travel time and depth is nonlinear, iterations are required, correctly accounting for the positions of the ray paths and reflection points. Each step of the iterative procedure comprises forward modeling of reflection times, inversion (calculation of the depth perturbation surface), and update of the reflector model.
The resolving power of this method is governed primarily by the density of reflection points, errors in the data, and errors in the assumed velocities [Hole et al., 1992]. These factors control the size of the operator used to smooth the depth perturbation surface. Tests using synthetic data as presented in the following section are the most effective means of evaluating the abilities of this procedure.
Tests using synthetic data. Two examples are presented to illustrate the accuracy with which known structure can be recovered using the inversion procedure outlined above. Both examples use the SCoRE '89 shot and receiver recording geometry and the same model dimensions (nx X ny X nz = 221 X 301 X 56 ~ 3.7 million nodes, or 264 X 360 X 66 km) and grid node spacing (1.2 km) employed in the analysis of the real data presented below (Figure 6). Noise-free synthetic data were generated for the same shot-receiver pairs that recorded real PmP arrivals. Velocities were fixed during the inversions. The velocity model used for both tests is based on the 2-D velocity models for lines 1, 2, and 3 [Zelt et al., 1992, 1993; McLean and Spence, 1994]. Velocities interior to the triangle were determined by linearly interpolating in the strike direction (roughly NW-SE); velocities exterior to the triangle were extrapolated outward from the 2-D models in directions perpendicular to the vertical plane comprising each 2-D model. The 3-D velocity model constructed in this way was smoothed with a 7 X 7 X 3-point MA filter to give the model used in the tests. This model is identical to the starting 3-D velocity model used in the real data analysis.
Sine wave reflector model. The reflector model used to generate the noise-free synthetic data was produced by adding two sinusoidal surfaces (Figure 7a). The primary features of the model within the study area are the shallow dome-like feature centered in the southern part of the triangle and deep bowl-like feature centered on the west side of the northern part of the triangle. The minimum and maximum depths are 29 and 37 km. A total of 1903 data were generated; shot-receiver offsets were limited to a maximum of 250 km. The starting reflector model was a horizontal plane at 35 km depth. Figure 7d shows the travel time residuals for the starting model plotted against the Y coordinate of the reflection point. The rms travel time residual for the starting model is 381 ms. The depth perturbations used in the first iteration are shown in Figure 7f. The roughly sinusoidal pattern in Figures 7d and 7f is an expression of the true reflector geometry.
At each iteration the gridded reflector perturbation model was smoothed with one pass of a 21 X 21-point MA filter. Figure 7c shows the final reflector depth model determined after three iterations (rms travel time residual is 9 ms) and the reflection point locations. Further iterations did not decrease the residuals. Structure between and exterior to the reflection points arises from the gridding and subsequent smoothing of the depth perturbations and is unconstrained.
The final model is in good agreement with the true model. The shape and depth of both primary features are accurately recovered. Other features around the edges of the model outside the region sampled by the data are not recovered. The final model does correctly include an increase in depth at the lower right corner of the model based on reflection points which sample the edge of this feature of the true model. The difference between the final and true model (Figure 7b) shows that errors in the region of coverage are within +/-1 km, i.e., less than the grid node spacing of 1.2 km. The model errors arise in part because of the underdetermined nature of the formulation. This introduces a degree of nonuniqueness in the solution which will be affected by nonlinearity, the starting model, and smoothness constraints. Errors in the data and/or velocities above the reflector will, of course, increase the errors in the final model. Figures 7e and 7g show the travel time residuals and depth perturbations after three iterations. Both are tightly clustered around zero.
Vertically faulted sine wave model. The reflector model comprises a vertically faulted sinusoidal surface similar to the previous example except with a shorter spatial wavelength (Figure 8a). Minimum and maximum depths to the north of the fault are 27 and 35 km and to the south are 32 and 40 km. The fault is located approximately along the surface trace of the CB-IMB boundary and has a vertical throw of 5 km. A number of the 1903 noise-free first arrival reflected data actually represent diffracted arrivals from the corner of the upthrown northeast side of the fault. These nongeometric events will introduce errors in the recovered model because equation (3), which assumes a reflected ray path, will be violated. The starting reflector model was a horizontal plane at 34 km. Smoothing operator dimensions were 11 X 11 grid nodes. The starting travel time residuals and depth perturbations used in the first iteration (Figures 8d and 8f) are significantly greater than the previous example because of the fault.
The final model (after eight iterations; Figure 8c) contains the main features of the true model. The sinusoidal surface north of the fault is recovered accurately with errors of less than 1 km (less than the grid node spacing; Figure 8b). To the south, the recovered model is distorted, but does show the primary features of the true model. The large error at (x,y) = (125,100) (Figure 8b) occurs away from the constraint points. At x = 200 km the vertical fault is imaged as a dipping fault zone of width ~15 km, approximately the width of the operator used to smooth the reflector depth perturbations. The top of the dipping zone corresponds approximately to the true fault location. To the northwest, at x = 150 km, the fault is poorly recovered and the dipping zone occurs ~20 km southwest of the true fault trace. The final travel time residuals and depth perturbations are shown in Figures 8e and 8g. The relatively large scatter on these plots and greater final rms travel time residual (31 ms) in comparison to the previous example arise because it is not possible to exactly reconstruct the true model: the assumption of reflected ray paths is violated for some of the data resulting in distortion in the recovered image. In addition, smoothing the updated reflector model limits the maximum dip with which a vertical fault can be represented. This more difficult synthetic test also serves to stress the importance of having uniform reflection point coverage to avoid the introduction of artifacts that arise solely from the gridding and smoothing operations.
Hole et al. [1992] describe a procedure for inverting first arrival travel time data to find the structure on an interface separating two media of different, and known, velocities. The inversion scheme is similar to that outlined above for the inversion of reflection times for depth, with equation (3) replaced by
The procedure described by Hole et al. [1992] was designed specifically for geometries in which each ray path crosses the interface only once, in which case the depth perturbation for each travel time residual, given by equation (4), is applied at the ray-interface intersection point. For the more usual case where the ray crosses the interface twice, for example, Pn at the Moho, this procedure can be generalized in a straightforward manner.
Rather than applying the entire travel time residual to one crossing, the residual is split between the two crossings. The manner in which the residual is split is arbitrary, but the least biased choice is to apply half of the residual at each crossing, so that
A problem which can affect the performance of the inversion arises because of inaccuracies in the calculated angle of refraction (q 2 in equation (5)) at the first, i.e., closest to the shot point, ray-interface crossing. In the vicinity of the critical point a narrow cone of rays incident on the interface fans out upon refraction through a wide range of angles (typically, 75 deg < theta2 < 90 deg). As represented in the finite difference model with a grid cell size of ~1.2 km, this occurs over a few grid cells. Since only a single ray direction can be represented within each cell (given by the direction of the average travel time gradient), the resolution of the theta2 values is very poor. For a typical model this will lead to errors in the calculated value of dt/dz of ~10-40%. The error is in magnitude only; the direction of the model perturbation will be correct. The accuracy can be increased by regridding the travel time field in the vicinity of the critical point i.e., interpolating the travel times onto a finer grid to more accurately calculate theta2. Note, however, that the errors quoted above are potentially much smaller than that due to the equal split of the residual between the two crossings. Also, because the inversion procedure is iterative, for small perturbations these errors should not seriously affect the accuracy of the inversion assuming isotropic ray coverage. Instead, it is the rate of convergence that is affected.
The complete data set comprises the arrival times of three phases: Pg, Pn (first arrivals) and PmP (Moho reflections). In the southwestern Canadian Cordillera, Pg constrains only the uppermost (<10 km) crustal velocity structure [Zelt et al., 1992, 1993]. There are no first arrival data to constrain deeper crustal velocities. Reasonable options for dealing with the lower crust are: (1) use a fixed, simple lower crustal velocity model, e.g., a 1-D average model based on the 2-D in-line models, or a 3-D model constructed from the 2-D models by interpolation between the models; or (2) use PmP (with layer stripping) to solve for lower crustal velocities (as well as for Moho topography). In the first option the PmP data misfits will only be mapped into Moho depth perturbations, and it is accepted that the lower crust may include gross errors in velocity. In the second option the PmP data misfits are mapped into both depth and velocity perturbations; however, it is not known how the misfits should be partitioned into the two types of perturbations. Regardless of the option chosen it is clear that the classic trade-off between velocity and depth introduces a degree of nonuniqueness that is difficult to assess. We prefer option 2 because it allows for some inference on possible lower crustal velocity variations. The same reasoning can be applied to the use of Pn to constrain both upper mantle velocities and depth to Moho since again there is an ambiguity regarding the source of Pn misfits. Thus both PmP and Pn are used to constrain velocities and depth to Moho.
Steps in the 3-D modeling algorithm are illustrated in the flowchart of Figure 9. Initially, an iterative procedure is used to invert Pg travel times for upper crustal velocities. Upper crustal velocities are then held fixed for all subsequent inversions of deeper phases (layer stripping). To do this, it is necessary to construct a 3-D surface (upper-lower crust boundary) separating the region of the velocity model constrained by Pg from deeper regions. The surface is constructed by gridding a set of points corresponding to the maximum depth of Pg penetration in each vertical column of grid cells comprising the velocity model. After gridding, the surface is smoothed with a 2-D MA filter to produce the final upper-lower crust boundary.
The deeper phases (PmP and Pn) are then inverted quasi- simultaneously for lower crustal velocities, Moho depth, and upper mantle velocities. A true, simultaneous, inversion is not performed. Instead, each step of the quasi-simultaneous inversion comprises: (1) inverting PmP travel times for lower crustal velocities; (2) inverting PmP and Pn travel times simultaneously for Moho depth using the procedures outlined above; and (3) inverting Pn travel times for upper mantle velocities (layer stripping is used to keep crustal velocities fixed during the third step). These three steps are repeated until a stopping criterion is satisfied, e.g., the rms travel time residuals are approximately equal to the estimated travel time uncertainties, or until the rms travel time residual stops decreasing.
In synthetic tests this procedure was sufficient to adequately fit the data. Tests with real data, however, showed that after a few iterations for deeper structure the rms travel time residual stopped decreasing without the data (particularly PmP) being adequately fit. This is apparently caused by convergence to a "final" (although poorly fit) Moho depth model after only a few iterations. The convergence results from the relatively large amount of smoothing applied to the Moho model at each step, which limits the degree to which the rms misfit can be reduced.To further reduce the misfits, the Moho is fixed and PmP and Pn are inverted in two stages to solve for velocities in the lower crust and upper mantle, respectively. In the first stage, PmP travel times are iteratively inverted to solve for lower crustal velocities until a stopping criterion is satisfied. This serves to map (most of) the remaining PmP travel time residuals into velocity variations. In the final stage, Pn travel times are iteratively inverted to solve for upper mantle velocities until a stopping criterion is satisfied. This serves to map (most of) the remaining Pn travel time residuals into velocity variations.
The starting velocity model was constructed from the three 2-D in-line profiles for lines 1, 2, and 3 as described above for the synthetic tests. The starting Moho depth model (Figure 10a) was also constructed from the three in-line profiles. Prior to the inversion of PmP travel times for both lower crustal velocities and Moho depth, velocities below the Moho were stripped off and replaced with velocities immediately above the Moho extended downward with a negative gradient. This ensures that the first arrival travel times at the reflecting surface are from the downgoing wave and not from waves that have turned beneath the reflector. The finite difference model dimensions are identical to the dimensions used in the synthetic tests (Figure 6).
The slowness perturbations were smoothed at two different stages of each iteration. In the gridding operation the slowness perturbation at each grid node was calculated as the average of all the du terms in a cubic grid cell volume centered at the grid node. In the subsequent smoothing operation the gridded du terms were smoothed with a MA filter. The dimensions of the smoothing operators (Table 2) were chosen so that excessive smearing of the velocity anomalies along the ray paths was avoided, while still allowing a good fit to the data (i.e., approximately to within the estimated average travel time uncertainties) by minimizing the operator size. The average estimated uncertainty of the travel time picks is 75 ms (Table 1). This estimate includes errors associated with site mislocations, nonlinear instrument clock drift, and elevation corrections that were applied to strip off the transit times through the 1 to 2-km-thick low-velocity near-surface layer. A large range of smoothing operator dimensions were tested, including fixing the size for all iterations and reducing the size for higher iteration numbers.
A total of 3963 Pg travel times were used in the inversion. Table 2 lists the starting and final rms travel time residuals. After the first few iterations the primary features of the final model were present and the rms travel time residual leveled off. With subsequent iterations the velocity model was only slightly refined and the residual decreased very slowly. After 10 iterations the misfit was 88 ms, and the procedure was stopped. The inability to reduce the misfit to the estimated uncertainty of the data likely indicates the presence of small- scale velocity variations unresolvable by the procedure due to the dimensions of the smoothing operators employed. Figure 11 displays the starting and final misfits for the seven large shot points. Figure 10b shows the upper-lower crust boundary constructed using the method described above. Velocities above this boundary were fixed in subsequent inversions.
The first step of the quasi-simultaneous inversion was the inversion of PmP travel times for lower crustal velocities. The smoothing operator dimensions were increased for the lower crust to reflect the sparser ray coverage and poorer vertical resolution supplied by the PmP travel time constraints (Table 2). A total of 1903 PmP travel times were used in the inversion. Figure 12a displays the starting misfits for the seven large shot points. In the second step, a combined total of 3071 PmP and Pn travel times were inverted for Moho depth. The rms travel time residual for the starting model shown in Table 2 includes the updated lower crustal velocities from the first stage of the quasi-simultaneous inversion procedure. The third step of the quasi-simultaneous inversion was the inversion of Pn travel times for upper mantle velocities. A total of 1168 Pn data were available for the inversion; however, because some of the calculated ray paths did not cross the Moho, not quite all of the data were used in each iteration. The rms travel time residual for the starting model shown in Table 2 includes the updated lower crustal velocities and Moho depth model from the first two steps of the quasi-simultaneous inversion procedure. Figure 12b shows the starting misfits for the seven large shot points.
After three iterations of the quasi-simultaneous inversion procedure, the rms travel time residuals stopped decreasing. Since the residual for PmP (135 ms) was significantly larger than the estimated average uncertainty (110 ms), the model after the third iteration was used as a starting model for the inversion of PmP for lower crustal velocities keeping the Moho fixed. Smoothing operator dimensions were not changed. The final PmP misfits for the seven large shot points are shown in Figure 12c. The "fine tuning" of the lower crustal velocities degraded the rms travel time residual of the Pn data to 125 ms, requiring further modeling of Pn travel times. The final Pn misfits for the seven large shot points are shown in Figure 12d. Table 2 shows the final rms travel time residuals for the inversion of PmP and Pn for lower crustal and upper mantle velocities.
Other inversion techniques, e.g., least squares, provide information on model resolution and parameter uncertainties through a model resolution matrix and covariance matrix. The procedure used here does not provide a quantitative measure of resolution since these matrices are not available. Nevertheless, it is possible to compute resolution kernels (for the linearized problem) which illustrate the degree to which the model at a grid node is derived from the true model. Hole [1992] describes the calculation of resolution kernels and presents an example for a synthetic 3-D survey showing the streaking of the resolution kernels along rays that penetrate a small box (two cubic grid cells) centered at the grid node. Unfortunately, we have found that the resolution kernels are time consuming to calculate and difficult to display and do not offer substantial insight into the overall resolution of a model.
A useful, though qualitative, estimate of resolution and uncertainty can be gauged from synthetic tests, using different starting models, varying the spatial smoothing operator sizes and looking at ray coverage [Hole et al., 1993]. Estimates for absolute uncertainty must also take systematic error sources into account. As pointed out by Hole et al. [1993], relative differences in velocity between neighboring regions are likely to be better resolved because they are less sensitive to systematic errors.
The resolving power of the data can be estimated from synthetic tests using a checkerboard velocity model [Hearn and Ni, 1994]. In these tests, noise-free synthetic data are generated using a velocity model comprising (in plan view) alternating high- and low-velocity squares and a linear vertical velocity gradient. Data are calculated for the same source- receiver combinations as used in the real analysis and are inverted using a 1-D starting velocity model. Results for a test using Pg arrival times and two different true models are shown in Figure 13. In both cases the squares were assigned alternating velocities of 6.0 and 6.3 km/s at the top of the model, beginning with 6.0 km/s in the southwestern corner. Velocities were extended downward with a gradient of 0.02 s-1. In the first model, squares were 25 X 25 km, and in the second, squares were 50 X 50 km. The starting 1-D model in both cases began at 6.15 km/s at the top, extended down with a gradient of 0.02 s-1. At the depth of the slices shown in Figure 13 the velocities of the true models are 6.16 and 6.46 km/s.
Figure 13a shows the results for the model with 25 X 25 km squares obtained by using smoothing operator dimensions that gave the best result. The data are not capable of resolving 25 X 25 km features throughout the study area.
Results for the model with 50 X 50 km squares are shown in Figure 13b. In general, the alternating pattern of low-to-high velocities is well resolved throughout the middle of the triangle. Although velocities in each cell vary by up to 0.3 km/s, the peak values are generally in error by only 0.05 km/s. The smoothing operator dimensions used in this test are identical to those used in the real data analysis and thus it is likely that upper crustal velocities imaged by the real data represent an average of the true Pg velocities over a region roughly 50 X 50 km. Table 3 lists the estimated spatial resolution and absolute model uncertainty for regions of the model that are constrained by ray coverage. The estimates are based on the results of this example and other synthetic tests.
Horizontal slices through the final 3-D velocity model are shown in Plates 1 and 2. Only those portions of the velocity model considered to be well constrained are shown. Uncolored regions represent grid node locations that are further than one half the lateral smoothing operator size from a ray. Although not shown, the velocity model does extend smoothly into these unconstrained regions as indicated by the velocity contours.
The slices show a significant degree of lateral inhomogeneity. The most prominent feature, in both size and magnitude, is the zone of relatively high crustal velocity in the south central region of the model. This feature is present, though small, at 2.8 and 7.6 km depth centered at (x,y) = (125,30) km, with velocities of 0.3-0.4 km/s greater than the mean at these depths. At greater depths the dimensions of this feature increase. Deep in the lower crust it comprises a rectangular zone in the southwest corner of the model extending eastward to the vicinity of the Fraser River Fault. To the west it extends to at least x=50 km; the data do not constrain the deep structure west of this.
The principal feature of the slice at 37.6 km (Plate 2d) is the relatively high upper mantle velocities (8.0 km/s) beneath the southern half of line 2 and at the southern end of line 1 and relatively low velocities (<7.7 km/s) beneath most of line 3 and near the north end of line 1. The lower crustal velocities (<7.0 km/s) in the Coast Belt are present in Plate 2d because this slice cuts through the Moho. Pn velocity throughout most of the interior of the triangle is between 7.7 and 8.0 km/s. The average Pn velocity at this depth is 7.85 km/s.
An indication of how the ray coverage varies with depth is shown in Figure 14a. Percentage-wise, the best coverage occurs near the top of the upper crust, followed by the top of the upper mantle. Coverage is poorest at the base of the upper crust and in the lower crust. Although significantly more PmP rays propagate through the lower crust in comparison to Pn rays through the upper mantle, the PmP ray paths are much closer to vertical and thus penetrate only a few cells at each level in the model. The Pn ray paths are more nearly horizontal, and each ray penetrates many cells at the same depth level.
Plate 3a shows the final Moho depth model. Plate 3c shows the location of the PmP reflection points and Pn Moho intersection points constraining depth to Moho. The points are most dense east of the CB-IMB boundary and are relatively sparse west of the Harrison Fault and along line 3. Overall, however, the coverage is good and provides a well-defined Moho image. The principal feature of the Moho model is the greater depth beneath much of the Coast Belt. Depth to Moho throughout most of the Intermontane Belt is 33-36 km. The minimum strongly constrained depth is ~32 km near the north end of line 1. In contrast, depths are greater than 36 km beneath most of the Coast Belt. Exceptions are the northernmost part of the Coast Belt in the study area, the region in the south that is roughly east of the Fraser River Fault, and the southwestern corner of the triangle where the data indicate a thinning of the crust to ~33 km. The maximum strongly constrained depth is ~39 km in the south central Coast Belt.
Moho depth is converted to two-way vertical travel time (TWTT) in Plate 3b. Throughout most of the Intermontane Belt the TWTT varies between 10.6 and 11.4 s. TWTT in most of the Coast Belt is 11.6-12.2 s. Four Lithoprobe reflection lines (88-11, 88-12, 88-15, and 88-18) overlie regions where the TWTT is well constrained by reflection points. Four other lines (88-13, 88-14, 88-16, and 88-17) are near (<50 km) the well-constrained area. The TWTTs from Plate 3b were compared with the interpreted times on the reflection sections of Cook et al. [1992] and Varsek et al. [1993]. Overall, there is good agreement with differences typically less than 400 ms, although the detailed structure (e.g., dip direction) is not always consistent. Note that differences in the values may be associated with differences in the reflection Moho, which is placed at the base of a zone of reflectivity, and the refraction Moho.
Plate 3d shows the standard deviation of the dz terms used in the final iteration. The values were obtained by dividing the region into a grid of 6 X 6 km cells and taking the standard deviation of the dz terms lying within each cell, provided the cell contained more than one term. The standard deviation is a rough measure of the consistency of the data, small values indicating that the data within a cell require a similar depth perturbation. The values are also a measure of the uncertainty in the reflector depth from errors in the data. The true uncertainty will be much larger because of errors in the velocities above and below the Moho. The majority of the standard deviations are less than 1 km, and there is no obvious trend in the distribution of the values.
Plate 4a shows a vertical slice through the starting velocity model along a north-south profile through the middle of the study area (Figure 2). The final model is shown in Plate 4b, with unconstrained regions uncolored. Note the large unconstrained zones in the upper crust near the middle of the model and the middle crust at the south end. The main differences between the starting and final model include significantly higher lower crustal velocities in the middle of the final model and lower velocities in the lower crust north of 200 km. The final model shows a zone of relatively low lower crustal velocity beneath the CB-IMB boundary. Otherwise, velocities in the lower crust below 20-25 km are significantly higher south of 200 km. Depth to Moho increases from 33 km in the north to 38 km at the south end. Plate 4c shows the difference between the final and starting model. The range of velocity perturbations seen on this slice (+/-0.2 km/s) is typical of the other slices presented below. Plate 4d is a representation of the ray density along a swath 24 km wide, i.e., the lateral dimensions of the smoothing operator used for modeling the upper crust, centered on the profile (see Figure 2). Coverage in the upper crust is incomplete, although where there is coverage the density is high. Coverage in the lower crust is better, but, in general, the ray density is lower. Upper mantle ray coverage is good.
To test the reliability of the results, the modeling procedure was repeated using 1-D starting models for the upper and lower crust and upper mantle, and a horizontal starting Moho. The 1- D velocities were taken from the average 1-D model for the study area (Figure 14b). Starting Moho depth was 34.5 km. The same number of iterations was used to model Pg and Pn. One additional iteration (i.e., four iterations) was used in the quasi-simultaneous inversion step. The final rms travel time residuals were nearly identical to those listed in Table 2.
Representative results for the final model are shown in Plate 5 and can be compared to results for the 3-D starting model in Plates 1d, 2b, 2d, and 3a. Differences in the final models are negligible near the surface (<12.4 km) with close correspondence between local high- and low-velocity anomalies. The differences become slightly more pronounced with increasing depth (Plates 5a-5c). Overall, there is still good agreement in the location of local anomalies but the magnitude of the velocities tend to differ by a greater amount. Maps of depth to Moho (Plates 3a and 5d) show similar features, including a significantly greater depth throughout most of the Coast Belt. One major difference is the much shallower depth within the Coast Belt beneath line 2 (33-34 km for the 1-D starting model case versus 36-37 km). This example of nonuniqueness is a result of the poor ray coverage in the overlying lower crust (Plates 2b, 2c, and 5b). A comparison of TWTTs to Moho with those interpreted by Cook et al. [1992] and Varsek et al. [1993] indicates significantly poorer agreement for several of the reflection profiles in comparison with the results for the 3-D starting model. For this reason, the model determined using the 3-D starting model is preferred.
A striking feature of the 3-D velocity structure shown in Plates 1 and 2 is the lateral heterogeneity indicated by the velocity slices. This heterogeneity is found at all depths, across and along the belts, and does not show any consistent correlation with the Intermontane-Coast belt boundary, which is defined on the basis of morphogeological features. Differences in average velocities between the belts (Figure 14b) are much less than the variations within the belts. Changes in crustal thickness (Plate 3a) also indicate a true 3-D structure is applicable to the region. Only over one small segment of the boundary, south of the shot point interior to the triangle, does a change in crustal thickness coincide with the boundary of the belts. These results indicate that large- scale crustal structure is not related to the morphogeologically defined belts. This conclusion also was reached by Jones and Gough [1995] in their summary of crustal conductivity derived from an extensive array of magnetotelluric stations across the southern Cordillera.
One of the principal features of the 3-D velocity model is the region of high relative velocities in the south central part of the study area, which extends from the upper crust, through the middle crust and into the lower crust (Plates 1 and 2). A significant decrease in these velocities occurs to the east in the vicinity of the Harrison Fault, consistent with the 2-D interpretation of line 3 data [Zelt et al., 1993]. The high velocities correlate with those determined for the upper crust of Vancouver Island [Drew and Clowes, 1990], which comprises Wrangellia, the accreted terrane forming the Insular Belt on the island. The same region includes the strongest Bouguer gravity anomaly in the study area (Plate 6a), a feature which extends westward to Vancouver Island. The relatively high velocity feature also coincides with some known surface exposures of Wrangellia on the mainland, the northeastward limit of which is identified by "uWr" in Plate 6b, and extends eastward beneath Jura-Cretaceous plutonic rocks and part of the Harrison terrane, terminating near the Harrison Lake Fault. Based on the high crustal velocities and high Bouguer anomalies, we infer that Wrangellian crust extends no farther east in the subsurface than the locus of the Harrison Fault and only about 30 km north of the base of the triangle (stippling in Plate 6b).
This interpretation helps to distinguish between two differing interpretations of reflection data along lines 88-14 and 88-17 (see Plate 3b for location) [Monger and Journeay, 1992; Varsek et al., 1993]. The two interpreted models for the collision zone between the Insular superterrane and North America differ significantly. Monger and Journeay [1992] proposed crustal delamination, with the Insular superterrane (Wrangellia) displaced along east vergent faults over rocks of the southeastern Coast Belt. Varsek et al. [1993] proposed crustal wedging, with the Insular superterrane interfingering the Intermontane superterrane. The primary difference is that the interpretation of Varsek et al. [1993] requires the entire crust of Wrangellia to extend at least as far east as line 10 (location in inset, Figure 1), approximately coincident with the Owl Lake/Harrison Faults, and northward to the positions of reflection lines 88-14 and 88-17. In the model of Monger and Journeay [1992], Wrangellia beneath line 10, which bisects the reflection lines, exists only as a narrow wedge between overlying Coast Belt terranes and plutonic assemblages, and underlying rocks of the southwestern Coast Belt. Our results, which indicate limited northward extension of Wrangellian crust, support the conclusion expressed by OÕLeary et al. [1993], who show that velocities along line 10 are similar to those along line 3 east of the Owl Lake/Harrison Faults, that the model of Monger and Journeay [1992] is more consistent with the seismic refraction results.
On the basis of recent geological mapping and geochronological studies, Journeay and Friedman [1993] identified the Coast Belt thrust system, a west vergent contractional belt that formed along the eastern margin of the Insular superterrane in early Late Cretaceous time. It extends between the Thomas Lake fault (TLF) on the southwest and the Bralorne-Kwoiek Creek Fault system on the northeast (Plate 6b). The depth to Moho determined in this study (Plate 3a) indicates the existence of a 3 to 4-km-thick crustal root beneath the Coast Belt thrust system and confirms the limited extent of the feature. Notably, the crustal root does not correspond to a Bouguer gravity low (Plate 6a). Indeed, part of the region with thicker crust has a corresponding gravity high, which is interpreted as the effect of high-velocity Wrangellian crust, as noted above. In general, neither the results from this study nor those from a general summary of refraction results across the southern Canadian Cordillera [Clowes et al., 1995] show a correspondence between variations in crustal thickness and changes in Bouguer gravity values. We infer that strong lateral changes in density within the crust, associated with the lateral changes in velocity, account for the discrepancy.
The Fraser River Fault is a major strike-slip fault with 90- 140 km of documented right-lateral motion [Coleman and Parrish, 1991] in the study area. According to Price and Carmichael [1986], it may be one component of a 2000-km- long intracontinental transform that extends from northern Washington state, where the Fraser Fault is known as the Straight Creek Fault, to the northern Rocky Mountain Trench/Tintina Fault system. However, in their preferred interpretation of the seismic reflection line across the Fraser Fault (line 88-18, location in Plate 3b), Varsek et al. [1993] suggested a listric geometry with the fault soling westward into a midcrustal zone of ductile shear. On the basis of a series of magnetotelluric profiles across the fault zone, Jones et al. [1992a] refuted this suggestion and stated that the Fraser Fault must penetrate the entire crust. Reprocessing and reinterpretation of the reflection data by Perz [1993] provided support for the latter conclusion. Further support is indicated by two features of the seismic velocity structure. In the north- south vertical slice from the velocity model (Plate 4b), a narrow zone of low velocities in the lower crust coincides with the position of the Fraser Fault system, which at that location is coincident with the Coast-Intermontane Belt boundary (Figure 2). Along other profiles in the southern Cordillera, known or inferred faults are sometimes associated with relatively low velocities [e.g., Kanasewich et al., 1994]. From east to west across much of the Fraser Fault, crustal thickness increases by about 3 km (Plates 3a and 6c). The lateral resolution of this study cannot isolate the change to be associated specifically with the fault, but the general location of the change indicates that the fault may be a controlling factor and thus is crustal penetrating.
Plate 6c shows a simple schematic interpretation of the present crustal structure along a southwest-northeast oriented (cross-strike) profile (location in Plate 6b) across the Coast and western Intermontane Belts. Geological boundaries are based on results presented in this study and other geological and geophysical information [Cook et al., 1992; Monger and Journeay, 1992; Perz, 1993; Varsek et al., 1993]. Velocities shown are the average velocities in a +/-25-km-wide swath surrounding the profile (blue rectangle in Plate 6b). As discussed below, a couple of different configurations of the basic units are represented in the presentation.
Beneath the southwestern Coast Belt, the accreted terrane Wrangellia, which forms the Insular Belt to the west, comprises the crust below the pervasive granitic plutons mapped in the region. The easternmost surface exposures of Wrangellia occur about 50 km from the western end of the profile ("uWr" in Plate 6c). Our results suggest that along this profile, deep Wrangellian crust may not extend much further northeast than this position. This preferred interpretation is indicated by the solid black line angling up and to the right from 50 km distance at the Moho. As discussed above, this configuration is consistent with the crustal delamination model of Monger and Journeay [1992]. However, if the association of Wrangellia with relatively high velocities in the deep crust were incorrect (e.g., the high velocities might be due to the more recent Cascade arc magmatism), then Wrangellian crust could extend much further northeast and be interfingered with Intermontane crust as proposed by Varsek et al. [1993]. This configuration is indicated by the jagged, black and white line between distances of 120 and 160 km and depths of 10 and 38 km. However, such an interpretation then requires somewhat lower velocities for Wrangellia in that region compared with the Wrangellian terrane on Vancouver Island [Drew and Clowes, 1990]. Alternatively, if the high velocities in our model were due to recent arc magmatism and Wrangellia also has high velocities, the terrane as a major component of the crust could be restricted to a region west of our study area. This is less consistent with the geological information than our preferred interpretation. Hence we consider the suggestion of the high velocities being due to Cascade arc magmatism a less attractive alternative.
To the east of the "uWr" line, the shallow crust comprises voluminous Jura-Cretaceous plutons that probably intrude Wrangellia. The same plutonic rocks intrude the Harrison terrane farther to the east, linking it to Wrangellia by at least 165 Ma [Monger, 1991]. The eastern limit of the Harrison terrane is the Owl Lake-Harrison Fault system (OLF), which in the south is the locus of lateral changes in the velocity structure, as discussed previously. The velocity model does not distinguish the group of terranes separated by faults of the Coast Belt thrust system but does show, for the first time, a small root below this system. Journeay and Friedman [1993] infer that the Coast Belt terranes east of the Thomas Lake Fault (TLF) are probably relatively shallow features, detached from their basements and stacked along thin-skinned thrust faults. This configuration is shown by the horizontal solid black line at a depth of 10 km below the Coast Belt thrust system. However, reprocessing of reflection profile 88-18 (location in Plate 3b) across the fault by Perz [1993] (location in Plate 3b) indicates that thrust faults of the Coast Belt thrust system may sole out substantially deeper, at depths of 20-25 km (J.M. Journeay (personal communication) as cited by Perz [1993]). Thus the terranes may extend down into the region denoted "CB?" in Plate 6c. In this scenario, the zone marked "Wr?" (above the horizontal dashed white line) may comprise Wrangellia. In the midcrust, the velocity model shows some support for this configuration.
East of the dextral Fraser River Fault, the crust consists primarily of Quesnellia, one of four terranes of the Intermontane superterrane. A zone of low velocities at the top of the crust immediately east of the Fraser River Fault outlines the Middle Jurassic Mount Lytton plutonic complex and Early Jurassic Guichon Creek batholith. Based on reflection and refraction interpretations [Cook et al., 1992; Zelt et al., 1992], the lower lithosphere of Quesnellia is not present in the vicinity of line 1. Presumably the effects of subduction erosion (removal of material from the base of an overlying unit as a result of the process of subduction; i.e., the opposite of subduction underplating) and/or crustal delamination account for this result. Instead, the relatively high velocity deep crust in the Intermontane Belt comprises parautochthonous North American rocks, and possibly a very thin layer of autochthonous North American crust. West of the Fraser River Fault, below the Coast Belt terranes, and east of Wrangellia, the deep crust may comprise highly metamorphosed rocks of the Intermontane superterrane, and possibly deeper elements of the overlying shallow terranes.
The crustal structure derived in this study and depicted in the cross section of Plate 6c is the product of processes that probably span a time interval from the early Mesozoic to the present. In the Intermontane Belt, removal of Quesnellia's lower lithosphere by erosion may have begun in the Late Triassic (prior to the accretion of the Intermontane superterrane to North America) as the oceanic Cache Creek terrane was being subducted beneath it. The disrooting continued in the Middle Jurassic by crustal delamination as Intermontane rocks were thrust eastwards over parautochthonous North American material (J.W.H. Monger, personal communication, 1994). The lower lithosphere of Quesnellia was presumably recycled in the mantle. A similar process may have delaminated the lithosphere of Wrangellia as the Insular superterrane collided with the Intermontane superterrane in the Early Cretaceous. Alternatively, the two superterranes may have wedged together, with Insular material interfingering the Intermontane superterrane at the collision zone. Contraction associated with continuing convergence between the two superterranes into the early Tertiary led to crustal thickening. The Coast Belt thrust system formed in response to convergence in the early Late Cretaceous. Terranes within the Coast Belt thrust system were detached from their basements and stacked along thin-skinned thrust faults, also resulting in a large amount of crustal thickening, the effects of which have been mostly removed by erosion during subsequent uplift. The small crustal root present in the Coast Belt may, in part, be a remnant signature of this contractional activity.
Plutonism, another mechanism of crustal thickening, occurred throughout the entire period of convergence associated with the accretion of the two superterranes. Terranes in the Western Coast Belt are intruded by mainly late Middle Jurassic to mid-Cretaceous plutons. In the Eastern Coast Belt, the intrusives are of mid-Cretaceous to early Tertiary age. Their effects are observed as the extensive surface plutons in the region. Their depth extent is not strongly constrained, but appears to be less than 10 km.
A change from a contractional regime to one of extension, transtension, and strike-slip faulting beginning in the early Eocene is represented in the cross section of Plate 6c by the Fraser River Fault system. The Fraser River and associated faults offset structures by 80-190 km of dextral motion in the middle and late Eocene [Misch, 1977; Monger, 1985; Coleman and Parrish, 1991; Friedman and van der Heyden, 1992].
The most recent phase of crustal evolution, within the last 40 Myr, has been dominated by activity associated with the Cascade magmatic arc. This includes uplift of the Coast Mountains, the development of a northeast trending fault system between 25 to 14 Ma which possibly controlled the emplacement of older volcanic and plutonic rocks, the effusion of young (4-0 Ma) volcanic rocks in the Garibaldi volcanic belt, and arc magmatism [Monger and Journeay, 1994]. The manifestation of effects associated with these younger processes is not readily apparent.
A new procedure for the inversion of wide-angle seismic data to determine 3-D velocity structure and depth to reflecting interfaces was developed to analyze the SCoRE '89 fan shot data set. The procedure is based on the 3-D finite difference algorithms described by Hole [1992] and Hole and Zelt [1995], with further developments for the inversion of reflection times to solve for depth to a reflecting interface and/or velocity structure, inversion of first arrival travel times to solve for depth to a refracting interface, and layer stripping.
Application of the procedure provides a 3-D velocity model for the southwestern Canadian Cordillera, a model that illustrates the pervasive lateral heterogeneity of the Intermontane and Coast Belts in the study area. The results are consistent with the 2-D models of Zelt et al. [1992, 1993], OÕLeary et al. [1993], and McLean and Spence [1994], given differences in modeling techniques and the poorer resolution of the 3-D model. A zone of high-velocity middle and lower crust in the southwest correlates with a strong relative gravity high and outlines the extent of the Insular superterrane (Wrangellia) in the southern Coast Belt, i.e., to the east as far as the Owl Lake-Harrison Fault system and north to ~49.5 deg N latitude. This interpretation favors the crustal delamination model of Monger and Journeay [1992] for the collision zone of the Insular and Intermontane superterranes, rather than the crustal wedging model suggested by Varsek et al. [1993].
The cross section in Plate 6c depicts a simplified interpretation of the present crustal structure in the Coast and Intermontane Belts. Behind the simplicity, however, lies a complex series of events and processes that have been active for an extended period of time. The 3-D velocity model presented in this study provides constraints on the crust and upper mantle velocity structure of the southwestern Canadian Cordillera and therefore helps to better understand the processes that led to its development. The model provides part of the database on which realistic geodynamical models of tectonic processes [Braun, 1993; Beaumont et al., 1994] will be founded.
Acknowledgments. We thank J. W. H. Monger for helpful discussions on Cordilleran geology. We also thank the ~35 field participants for making the recording program a success. Principal financial support for the refraction program was from Lithoprobe, which is funded by a Natural Sciences and Engineering Research Council of Canada (NSERC) Collaborative Project and Programs grant. Supplementary funding was derived from NSERC Research Grants to R.M.C., R.M.E. and E. R. Kanasewich and an NSERC Infrastructure grant to R.M.C. and R.M.E. Analysis was funded by the Energy Mines and Resources Canada Research Agreements Program, an NSERC research grant to R.M.E., and a National Science Foundation grant to J.A.H. B.C.Z. was, in part, financially supported by an NSERC Postgraduate Scholarship and a University of British Columbia Graduate Fellowship. Lithoprobe publication 688.
Aki, K., A. Christoffersson, and E.S. Husebye, Determination of the three-dimensional seismic structure of the lithosphere, J. Geophys. Res., 82, 277-296, 1977.
Beaumont, C., P. Fullsack, and J. Hamilton, Styles of crustal deformation in compressional orogens caused by subduction of the underlying lithosphere, Tectonophysics, 232, 119-132, 1994.
Braun, J., Three-dimensional modeling of compressional orogenies: Thrust geometry and oblique convergence, Geology, 21, 153-156, 1993.
Clowes, R.M., F.A. Cook, A.G. Green, C.E. Keen, J.N. Ludden, J.A. Percival, G.M. Quinlan, and G.F. West, Lithoprobe: New perspectives on crustal evolution, Can. J. Earth Sci., 29, 1813-1864, 1992.
Clowes, R.M., C.A. Zelt, J.R. Amor, and R.M. Ellis, Lithospheric structure in the southern Canadian Cordillera from a network of seismic refraction lines, Can. J. Earth Sci., in press, 1995.
Coleman, M.E., and R.R. Parrish, Eocene dextral strike-slip and extensional faulting in the Bridge River terrane, southwest British Columbia, Tectonics, 10, 1222-1238, 1991.
Coney, P.J., D.L. Jones, and J.W.H. Monger, Cordilleran suspect terranes, Nature, 288, 329-333, 1980.
Cook, F.A., J.L. Varsek, and R.M. Clowes, Lithoprobe reflection transect of southwestern Canada: Mesozoic thrust and fold belt to mid-ocean ridge, in Continental Lithosphere: Deep Seismic Reflections, Geodyn. Ser., vol. 22, edited by R. Meissner, L.D. Brown, H. Durbaum, K. Fuchs, and F. Seifert, pp. 247-255, AGU, Washington, D.C., 1991.
Cook, F.A., J.L. Varsek, R.M. Clowes, E.R. Kanasewich, C.S. Spencer, R.R. Parrish, R.L. Brown, S.D. Carr, B.J. Johnson, and R.A. Price, Lithoprobe crustal reflection cross section of the southern Canadian Cordillera, 1, Foreland Thrust and Fold Belt to Fraser River Fault, Tectonics, 11, 12-35, 1992.
Cook, F.A., J.L. Varsek, and J.B. Thurston, Significance of gravity and magnetic variations along the Lithoprobe Southern Cordillera Transect, Can. J. Earth Sci., in press, 1995. Davis, E.E., and T.J. Lewis, Heat flow in a back-arc environment: Intermontane and Omineca Crystalline Belts, southeastern Canadian Cordillera, Can. J. Earth Sci., 21, 715-726, 1984.
Drew, J.J., and R.M. Clowes, A re-interpretation of the seismic structure across the active subduction zone of western Canada-CCSS Workshop Topic I, onshore-offshore data set, in Studies of Laterally Heterogeneous Structures Using Seismic Refraction and Reflection Data (Proceedings of the 1987 Commission on Controlled Source Seismology Workshop), edited by A.G. Green, Pap. Geol. Surv. Can., 89-13, 115-132, 1990.
Friedman, R.M., and P. van der Heyden, Late Permian U-Pb dates for the Farwell and northern Mount Lytton plutonic bodies; Intermontane Belt, British Columbia, in Current Research, Part A, Pap. Geol. Surv. Can. 92-1a, 137-144, 1992.
Gabrielse, H., and C.J. Yorath (Eds.), Vol. G-2, Geology of the Cordilleran Orogen in Canada, The Geology of North America, 844 pp., Geol. Soc. of Am., Boulder, Colo., 1992.
Gehrels, G.E., and J.B. Saleeby, Constraints and speculations on the displacement and accretionary history of the Alexander-Wrangellia- Peninsular superterrane, Geol. Soc. Am. Abstr. Programs, 17, 356, 1985.
Gough, D.I., Mantle upflow tectonics in the Canadian Cordillera. J. Geophys. Res., 91, 1909-1919, 1986.
Hearn, T.M., and F.F. Ni, Pn velocities beneath continental collision zones: The Turkish-Iranian Plateau, Geophys. J. Int., 117, 273-283, 1994.
Hole, J.A., Nonlinear high-resolution three-dimensional seismic travel time tomography, J. Geophys. Res., 97, 6553-6562, 1992.
Hole, J.A., and B.C. Zelt, Three-dimensional finite-difference reflection travel times, Geophys. J. Int., 121, 427-434, 1995.
Hole, J.A., R.M. Clowes, and R.M. Ellis, Interface inversion using broadside seismic refraction data and three-dimensional travel time calculations, J. Geophys. Res., 97, 3417-3429, 1992.
Hole, J.A., R.M. Clowes, and R.M. Ellis, Interpretation of three- dimensional seismic refraction data from western Hecate Strait, British Columbia: Structure of the crust, Can. J. Earth Sci., 30, 1440- 1452, 1993.
Humphreys, E., and R.W. Clayton, Adaptation of back projection tomography to seismic travel time problems, J. Geophys. Res., 93, 1073-1085, 1988.
Jones, A.G., and D.I. Gough, Electromagnetic images of crustal structures in southern and central Canadian Cordillera, Can. J. Earth Sci., in press, 1995.
Jones, A.G., R.D. Kurtz, D.E. Boerner, J.A. Craven, G.W. McNeice, D.I. Gough, J.M. DeLaurier, and R.G. Ellis, Electromagnetic constraints on strike-slip fault geometry-The Fraser River Fault system, Geology, 20, 561-564, 1992a.
Jones, A.G., D.I. Gough, R.D. Kurtz, J.M. DeLaurier, D.E. Boerner, J.A. Craven, R.G. Ellis, and G.W. McNeice, Electromagnetic images of regional structure in the southern Canadian Cordillera, Geophys. Res. Lett., 19, 2373-2376, 1992b.
Journeay, J.M., and R.M. Friedman, The Coast Belt Thrust System: Evidence of Late Cretaceous shortening in southwest British Columbia, Tectonics, 12, 756-775, 1993.
Kanasewich, E.R., M.J.A. Burianyk, R.M. Ellis, R.M. Clowes, D.J. White, T. C™tŽ, D.A. Forsyth, J.H. Luetgert, and G.D. Spence, Crustal velocity strcuture of the Omineca Belt, southeastern Canadian Cordillera, J. Geophys. Res., 99, 2653-2670, 1994.
Lecomte, I., and S. Hamran, Local plane-wavenumber diffraction tomography in heterogeneous backgrounds, II, Green's functions and finite-difference travel times, J. Seismic Explor., 2, 287-299, 1993.
Lewis, T.J., W.H. Bentkowski, E.E. Davis, and R.D. Hyndman, Subduction of the Juan de Fuca plate: Thermal consequences, J. Geophys. Res., 93, 15,207-15,225, 1988.
Lewis, T.J., W.H. Bentkowski, and R.D. Hyndman, Crustal temperatures near the Lithoprobe southern Canadian Cordilleran transect. Can. J. Earth Sci., 29, 1197-1214, 1992.
Matsuoka, T., and T. Ezaka, Ray tracing using reciprocity, Geophysics, 57, 326-333, 1992.
McLean N.A., and G.D. Spence, Seismic velocity structure from the eastern Insular Belt to the Intermontane Belt, paper presented at Canadian Cordillera Tectonics Workshop, Univ. of Victoria, B.C. Geol. Surv., Victoria, B.C., 1994.
Misch, P., Dextral displacements at some major strike faults in the north Cascades, Geol. Assoc. Can. Program Abstr., 2, 37, 1977.
Monger, J.W.H., Structural evolution of the southwestern Intermontane Belt, Ashcroft and Hope map-areas; British Columbia, in Current Research, Part A, Pap. Geol. Surv. Can., 85-1A, 349-358, 1985.
Monger, J.W.H., Correlation of Settler schist with Darrington phyllite and Shuksan greenschist and its tectonic implications, Coast and Cascade mountains, British Columbia and Washington, Can. J. Earth Sci., 28, 447-458, 1991.
Monger, J.W.H., and J.M. Journeay, Guide to the geology and tectonic evolution of the southern Coast Belt, Field Guide, Penrose Conference on the Tectonic Evolution of the Coast Mountains Orogen, Geol. Surv. Can., Vancouver, 1992.
Monger, J.W.H., and J.M. Journeay, Basement geology and tectonic evolution of the Vancouver region, southwestern British Columbia, in Geology and Geological Hazards of the Vancouver Region, Southwestern British Columbia, Bull. Geol. Surv. Can. 481, 3-25, 1994.
Monger, J.W.H., and R.A. Price, Geodynamic evolution of the Canadian Cordillera-Progress and problems, Can. J. Earth Sci., 16, 770-791, 1979.
Monger, J.W.H., R.A. Price, and D.J. Tempelman-Kluit, Tectonic accretion and the origin of the two major metamorphic and plutonic welts in the Canadian Cordillera, Geology, 10, 70-75, 1982.
Monger, J.W.H., P. van der Heyden, J.M. Journeay, C.A. Evenchick, and J.B. Mahoney, Jurassic-Cretaceous basins along the Canadian Coast Belt: Their bearing on pre-mid-Cretaceous sinistral displacements, Geology, 22, 175-178, 1994.
Monger, J.W.H., R.M. Clowes, D.S. Cowan, C.J. Potter, R.A. Price, and C.J. Yorath, Continent-ocean transitions in western North America between latitudes 46 and 56 degrees, in Phanerozoic Evolution of North American Continent-Ocean Transitions, Centennial Transect Volume CTV-1, edited by R.C. Speed, pp. 357-397,Boulder, Colo., Geol. Soc. of Am., in press, 1995.
Moser, T. J., Shortest path calculation of seismic rays, Geophysics, 56, 59-67, 1991.
Oldow, J.S., A.W. Bally, G.G. AvŽ Lallemant, and W.P. Leeman, Phanerozoic evolution of the North American Cordillera; United States and Canada, in The Geology of North America, Vol. A, The Geology of North America: An Overview, edited by W. Bally and A.R. Palmer, pp. 139-232, Geol. Soc. of Am., Boulder, Colo., 1989.
O'Leary, D.M., R.M. Clowes, and R.M. Ellis, Crustal velocity structure in the southern Coast Belt, British Columbia, Can. J. Earth Sci., 30, 2389-2403, 1993.
Perz, M.J., Characterization of the Fraser River Fault, southwestern British Columbia, and surrounding geology through reprocessing of seismic reflection data, thesis, 191 pp., Univ. of B.C., Vancouver, 1993.
Phillips, W.S., and M.C. Fehler, Travel time tomography: A comparison of popular methods, Geophysics, 56, 1639-1649, 1991.
Podvin, P., and I. Lecomte, Finite difference computation of travel times in very contrasted velocity models: A massively parallel approach and its associated tools, Geophys. J. Int., 105, 271-284, 1991.
Price, R.A., and D.M. Carmichael, Geometric test for late Cretaceous- Paleogene intracontinental transform faulting in the Canadian Cordillera, Geology, 14, 468-471, 1986.
Riahi, M.A., and C. Juhlin, 3-D interpretation of reflected arrival times by finite-difference techniques, Geophysics, 59, 844-849, 1994.
Sweeny, J.F., R.A. Stephenson, R.G. Currie, and J.M. DeLaurier, Part C. Crustal geophysics, in The Geology of North America, vol. G-2, Geology of the Cordilleran Orogen in Canada, edited by H. Gabrielse and C.J. Yorath, Geol. Soc. of Am. pp. 39-59, Boulder, Colo., 1991.
Thurber, C.H., and K. Aki, Three-dimensional seismic imaging, Annu. Rev. Earth Planet. Sci., 15, 115-139, 1987. van der Heyden, P., A middle Jurassic to Early Tertiary Andean-Sierran arc model for the Coast Belt of British Columbia, Tectonics, 11, 82- 97, 1992.
Varsek, J.L., F.A. Cook, R.M. Clowes, J.M. Journeay, J.W.H. Monger, R.R. Parrish, E.R. Kanasewich, and C.S. Spencer, Lithoprobe crustal reflection structure of the southern Canadian Cordillera, 2, Coast Mountains transect, Tectonics, 12, 334-360, 1993.
Vidale, J.E., Finite-difference calculation of travel times, Bull. Seismol. Soc. Am., 78, 2062-2076, 1988.
Vidale, J.E., Finite-difference calculation of travel times in three dimensions, Geophysics, 55, 521-526, 1990.
Zelt, B.C., R.M. Ellis, and R.M. Clowes, SCoRE '89: the Southern Cordillera refraction experiment: Description of data set, Lithoprobe Secr., Lithoprobe Rep., 20, 38 pp., Univ. of B.C., Vancouver, 1990.
Zelt, B.C., R.M. Ellis, R.M. Clowes, E.R. Kanasewich, I. Asudeh, J.H. Luetgert, Z. Hajnal, A. Ikami, G.D. Spence, and R.D. Hyndman, Crust and upper mantle velocity structure of the Intermontane Belt, southern Canadian Cordillera, Can. J. Earth Sci., 29, 1530-1548, 1992.
Zelt, B.C., R.M. Ellis, and R.M. Clowes, Crustal velocity structure in the eastern Insular and southernmost Coast Belts, Canadian Cordillera, Can. J. Earth Sci., 30, 1014-1027, 1993.
Zelt, C.A., and R.B. Smith, Seismic travel time inversion for 2-D crustal velocity structure, Geophys. J. Int., 108, 16-34, 1992.
Zhao, D., A. Hasegawa, and S. Horiuchi, Tomographic imaging of P and S wave velocity structure beneath northeastern Japan, J. Geophys. Res., 97, 19,909-19,928, 1992.
J. A. Hole, Department of Geophysics, Stanford University, Stanford, CA 94305-2215. (e-mail: hole@riel.stanford.edu)
B. C. Zelt, Geotechnology Research Institute, Houston Advanced Research Center, The Woodlands, TX 77381. (e-mail: s22baze@gtri.harc.edu)
(Recieved December 19, 1994; revised August 14, 1995; accepted September 6, 1995.)
1Now at Geotechnology Research Institute, Houston Advanced Research Center, The Woodlands, Texas.
Copyright 1995 by the American Geophysical Union.
Paper number 95JB02807. 0148-0227/95/95JB-02807$05.00
Figure 1. Study area for the 1989 Southern Cordillera Refraction Experiment. Numbered solid triangles are shot points. Black and white dashed lines represent approximate location of receiver sites. Profile numbers are in circles. FRF, Fraser River Fault (dashed line); PF, Pasayten Fault; YF, Yalokom Fault; GVB, Garibaldi volcanic belt; JFP, Juan de Fuca Plate; H, Hope; HM, 100 Mile House; K, Kamloops; P, Princeton; PA, Port Alberni; V, Vancouver; VI, Vancouver Island. Left inset map shows morphogeological belts of the Canadian Cordillera. Dashed gray line within the Coast Belt denotes SCoRE '90 line 10 refraction profile. Right inset map shows the distribution of Cordilleran superterranes.
Figure 2. Shot points and receivers used in 3-D analysis of data. All map view results are displayed using the rotated (7.6 deg) Mercator coordinate system shown here. Dashed line between SPs 1 and 21 denotes the location of the vertical slice shown in Plate 4. IMB, Intermontane Belt; CB, Coast Belt; FRF, Fraser River Fault; HF, Harrison Fault; K, Kamloops; V, Vancouver.
Figure 3. Record section for SP 1 into line 3. Data are trace normalized, band-pass filtered from 3 to 15 Hz and plotted with a reducing velocity of 8 km/s. Horizontal scale is shot-receiver azimuth measured clockwise from north. The approximate locations of the travel time picks are represented by the heavy line. At the offsets for this section, only Pn is observed. Nomogram shows location of shot point (large dot) and recording line (heavy line).
Figure 4. Record section for SP 12 into line 1. See Figure 3 caption for other information.
Figure 5. Record section for SP 7 into line 2. See Figure 3 caption for other information.
Figure 6. Schematic representation of the 3-D finite difference velocity model used in synthetic tests and real data analyses. Grid node spacing is 1.2 km, and model dimensions are 221 X 301 X 56.
Figure 7. (a) Contour plot of sine wave reflector model. White lines represent approximate location of recording lines; triangles represent shot points used in synthetic tests. (b) Difference between final and true models. (c) Final model after three iterations. Dots are PmP reflection points and represent point constraints on depth to reflector. (d) The rms travel time residuals for the starting model plotted as a function of the Y coordinate of the reflection point. (e) The rms travel time residuals after three iterations. (f) Depth perturbations used in first iteration. (g) Depth perturbations used in third iteration.
Figure 8. (a) Contour plot of vertically faulted sine wave reflector model. (b) Difference between final and true models. (c) Final model after three iterations. Dots are PmP reflection points and represent point constraints on depth to reflector. (d) The rms travel time residuals for starting model plotted as a function of the Y coordinate of the reflection point. (e) The rms travel time residuals after eight iterations. (f) Depth perturbations used in first iteration. (g) Depth perturbations used in eighth iteration.
Figure 9. Flowchart of the 3-D modeling algorithm used in the real data analysis. Inset of schematic velocity model illustrates the five basic elements comprising the complete model: three 3-D velocity models and two 3-D surfaces.
Figure 10. (a) Starting Moho depth model for 3-D inversions. The 3-D surface was constructed from the three 2-D Moho depth models for lines 1, 2, and 3 [Zelt et al., 1992, 1993; McLean and Spence, 1994]. Heavy dashed lines represent location of constraint points along the 2-D profiles that were gridded to derive starting model. (b) Upper-lower crust boundary defined by maximum depth penetration of upper crustal rays. IMB, Intermontane Belt; CB, Coast Belt; FRF, Fraser River Fault; HF, Harrison Fault; V, Vancouver. Triangles represent shot points; black and white lines approximate location of receiver sites (see Figure 2 for exact receiver site locations).
Figure 11. Example misfits from the inversion of Pg travel times for upper crustal velocities. (a) Starting and (b) final (after 10 iterations) travel time misfits for the seven large shot points. Spikes along the track for each shot point represent the misfit at each recording site along lines 1, 3, and 2. "Stations" 1-277 are sites along line 1 from north to south; 278-488 are sites along line 3 from east to west; 489-716 are sites along line 2 from southwest to northeast. The rms travel time residual for each shot are given along the upper left axis (see also Table 2).
Figure 12. (a) Starting and (c) final (after three iterations of the quasi-simultaneous inversion procedure and 10 iterations in which PmP were inverted for velocities only) misfits from the inversion of PmP travel times for lower crustal velocities (and depth to Moho). (b) Starting and (d) final (after three iterations of the quasi-simultaneous inversion procedure and six iterations in which Pn were inverted for velocities only) misfits from the inversion of Pn travel times for upper mantle velocities (and depth to Moho). See Figure 11 for other plotting information.
Figure 13. Resolution tests for Pg using checkerboard velocity models with squares of dimension (a) 25 X 25 km and (b) 50 X 50 km (indicated by the white grid lines). At the depth of the slices shown the true velocities alternate between 6.16 and 6.46 km/s, beginning at 6.16 km/s in the lower left square.
Figure 14. (a) Percentage of the study area sampled and number of cells hit by rays plotted against depth. Percentage is calculated as the number of cells penetrated by at least one ray divided by the total number of cells. Study area is defined as triangle formed by SPs 1, 3, and 5. Shaded regions show the range in variation of the upper-lower crust boundary and Moho. Horizontal dashed lines denote average depth to the upper-lower crust boundary (ulcb) and Moho (M). (b) Average 1-D velocity profiles for the entire study area and parts of the Coast and Intermontane Belts within the study area.
Plate 1. Horizontal slices through the final 3-D velocity model at depths of (a) 2.8 km, (b) 7.6 km, (c) 12.4 km (through the upper-lower crust boundary), and (d) 18.4 km. Unconstrained regions are shown as uncolored. CB, Coast Belt; IMB, Intermontane Belt; FRF, Fraser River Fault; HF, Harrison Fault; V, Vancouver.
Plate 2. Horizontal slices through the final 3-D velocity model at depths of (a) 23.2 km, (b) 28.0 km, (c) 32.8 km and (d) 37.6 km. Slices in plates 2c and 2d cut through the Moho. See Plate 1 for other information.
Plate 3. (a) Depth to Moho. (b) Two-way vertical travel time to Moho from a horizontal plane at 1 km above sea level. Purple lines show location of Lithoprobe reflection lines (88- 11 to 88-18). (c) Location of PmP reflection points (blue) and Pn Moho intersection points (yellow) constraining solution in Plate 3a. (d) Standard deviation of dz terms in 6 X 6 km cells.
Plate 4. (a) Vertical slice through starting model along profile through SPs 1 and 21 (see Figure 2 for location). (b) Vertical slice through final model. See Plate 1 for color velocity scale. (c) Difference between final and starting models. (d) Ray density near the slice. Plot shows count of rays in swath (shaded region) shown in Figure 2. CB-IMB, Coast Belt-Intermontane Belt boundary (Fraser River Fault).
Plate 5. Representative results for a 1-D starting model. Horizontal slices through the final 3-D velocity model at depths of (a) 18.4 km (compare with Plate 1d), (b) 28.0 km (compare with Plate 2b), and (c) 37.6 km (compare with Plate 2d. Heavy solid lines are closely spaced contours. (d) Depth to Moho (compare with Plate 3a). See Plate 1 for color velocity scale and Plate 3a for color depth scale.
Plate 6. (a) Bouguer anomaly map for the southwestern Canadian Cordillera. (b) Location map for cross section (heavy solid line) in Plate 6c. Velocities in the cross section represent average values within the blue rectangle. Shaded region represents possible extent of lower crust of Wrangellia based on seismic refraction, reflection and gravity data. uWr, eastern limit of upper Wrangellian surface exposures; HF, Harrison Fault. (c) Schematic cross section across the Coast and western Intermontane Belts. WCB, Western Coast Belt; ECB, Eastern Coast Belt; IMB, Intermontane Belt. Faults depicted are: CBTS, Coast Belt thrust system; TLF, Thomas Lake Fault; OLF, Owl Lake Fault; CCBD, Central Coast Belt detachment; BKFS, Bralorne-Kwoiek Creek Fault system; FRF, Fraser River Fault. Terrane abbreviations: Wr, Wrangellia ("+g" represents the dominant granitic rocks emplaced in Wrangellia); HA, Harrison; CD, Cadwallader; SH, Shuksan; BR, Bridge River; Qn, Quesnellia. Other abbreviations: CB? represents possible depth extent of overlying Coast Belt terranes; Intermontane s.t.? represents possible extent of (deep and highly metamorphosed) Intermontane superterrane rocks; PNA, parautochthonous North America; crosses represent Jura- Cretaceous plutons; M, Moho (white dashed line); Wr? represents possible eastern extent of Wrangellia. See text for other information.
