Reservoir Characterization    Gashydrate    Tomography    Near Surface    Casadia Subduction Zone

Regional P wave velocity structure of the Northern Cascadia Subduction Zone

Citation: Ramachandran, K., R. D. Hyndman, and T. M. Brocher (2006), Regional P wave velocity structure of the Northern Cascadia Subduction Zone, J. Geophys. Res., 111, B12301, doi:10.1029/2005JB004108.

Abstract

This paper presents the first regional three-dimensional P wave velocity model for the Northern Cascadia Subduction Zone (SW British Columbia and NW Washington State) constructed through tomographic inversion of first-arrival traveltime data from active source experiments together with earthquake traveltime data recorded at permanent stations. The velocity model images the structure of the subducting Juan de Fuca plate, megathrust, and the fore-arc crust and upper mantle. Beneath southern Vancouver Island the megathrust above the Juan de Fuca plate is characterized by a broad zone (25–35 km depth) having relatively low velocities of 6.4–6.6 km/s. This relative low velocity zone coincides with the location of most of the episodic tremors recently mapped beneath Vancouver Island, and its low velocity may also partially reflect the presence of trapped fluids and sheared lower crustal rocks. The rocks of the Olympic Subduction Complex are inferred to deform aseismically as evidenced by the lack of earthquakes within the low-velocity rocks. The fore-arc upper mantle beneath the Strait of Georgia and Puget Sound is characterized by velocities of 7.2–7.6 km/s. Such low velocities represent regional serpentinization of the upper fore-arc mantle and provide evidence for slab dewatering and densification. Tertiary sedimentary basins in the Strait of Georgia and Puget Lowland imaged by the velocity model lie above the inferred region of slab dewatering and densification and may therefore partly result from a higher rate of slab sinking. In contrast, sedimentary basins in the Strait of Juan de Fuca lie in a synclinal depression in the Crescent Terrane. The correlation of in-slab earthquake hypocenters M > 4 with P wave velocities greater than 7.8 km/s at the hypocenters suggests that they originate near the oceanic Moho of the subducting Juan de Fuca plate.

1. Introduction

Figure 1.  Location map showing the SHIPS temporary land-based receiving stations (blue triangles) and air gun shot positions (red lines) of the active source data used in the present study. Bottom left inset shows the earthquakes (blue stars) and permanent recording stations (red triangles) used in this study. Inset to the right top shows the plate tectonic regime of the study area.

Northern Cascadia subduction zone in the Pacific Northwest of North America is a region of high earthquake hazard. Megathrust earthquakes, such as the one that occurred in 1700 [Satake et al., 1996], have an average recurrence of ~500 to 600 years [Atwater and Hemphill-Haley, 1997] and nearly 10 earthquakes of magnitude 6 and above have occurred in the past 125 years [e.g., Clague, 1997]. Cenozoic fore-arc basins such as the Seattle, Tacoma, and the Georgia basins may amplify ground motion due to earthquake energy in high population centers [e.g., Pratt et al., 2003]. An accurate estimate of the subsurface velocity structure is necessary to model the shaking effects of earthquakes and to understand the origin of the earthquakes. The main objective of this study is to construct a regional P wave velocity model for the northern Cascadia subduction zone through tomographic inversion of active source and earthquake data from SW British Columbia, Canada and NW Washington State (Figure 1).

The velocity models from previous tomographic studies [Zelt et al., 2001; Ramachandran et al., 2004; Parsons et al., 1999; Brocher et al., 2001; Van Wagoner et al., 2002; Zhao, 2001; Graindorge et al., 2003; Ramachandran et al., 2005; Lees and Crosson, 1990; Symons and Crosson, 1997; Stanley et al., 1999; Tréhu et al., 2002] varied in regional coverage, some with overlapping borders, and the velocity models were constructed with different traveltime data sets, inversion parameters, and their resolutions were different. In the present study, a joint tomography of active source and earthquake data is employed to image the active subduction zone and the fore-arc structure beneath SW British Columbia and NW Washington State from 46.5°N to 49.0°N, and from 120.0°E to 125.5°E (Figure 1). The tomographic inversion included 147,000 first arrival times recorded at 225 temporary stations widely distributed over SW British Columbia and NW Washington from the Seismic Hazards Investigation in Puget Sound (SHIPS) experiment in 1998 [Fisher et al., 1999; Brocher et al., 1999]. In addition, 72,000 arrival times from approximately 3,000 local earthquakes recorded at 91 permanent recording stations located in SW British Columbia and NW Washington were also employed in the inversion. The P wave velocity model constructed in this study provides the first regional image of the the northern Cascadia subduction zone obtained from a single tomographic inversion. The velocity model helps constrain the position of the subducting plate along a 300 km north-south stretch of the northern Cascadia subduction zone representing nearly one fourth of the entire zone. The relocated in-slab earthquake positions and the associated velocity at the hypocenters provide evidence that the in-slab earthquakes occur near the oceanic Moho of the subducting Juan de Fuca plate.

2. Geology, Tectonics, and Seismicity

Figure 2.  Sedimentary basin and fault map. CFTB, Cowichan Fold and Thrust Belt; CH, Chuckanut subbasin; CLB, Clallam basin; CPC, Coast Plutonic Complex; CRBF, Coast Range Boundary fault; CR, Crescent terrane; DDMF, Darrington-Devils Mountain fault; EB, Everett basin; HCF-Hood Canal fault; HRF-Hurricane Ridge fault; KA, Kingston Arch; LIF, Lummi Island fault; LRF, Leech River fault; B, Muckleshoot Basin; NA, Nanaimo subbasin; OF, Olympia fault; OIF, Outer Islands fault; PB, Possesion Basin; PR, Pacific Rim terrane; PTB, Port Townsend basin; SB, Seattle basin; SF, Seattle fault; SJF, San Juan fault; SMF, Survey Mountain fault; SQB, Sequim basin; SQF, Sequim fault; SU, Seattle uplift; SWIF, southern Whidbey Island fault; TB, Tacoma basin; TF, Tacoma fault; WA, Whatcom subbasin. AB, CD, EF, GH, IJ, KL, MN, NO, PQ, and QR show the location of the vertical cross sections shown in Figures 6, 7, and 8. ALB and LAS are receiver function study locations from Cassidy and Ellis [1993].

The Intermontane superterrane (Figure 1, top right inset), made up mostly of sedimentary and volcanic rocks, collided with the North America plate, about 200 My ago. The last major collisional episode, around mid-Cretaceous time, emplaced the Insular superterrane against the Intermontane superterrane (Figure 2). The two superterranes form the Intermontane and Insular belts, respectively. The mid-Cretaceous to early Tertiary intrusive rocks of the Coast Belt were emplaced during the collisional episode [Monger et al., 1982; Monger, 1990]. In the SW British Columbia margin, Vancouver Island is dominated by the Wrangellia terrane, emplaced during the middle Cretaceous [Smith and Tipper, 1986] and thought to represent a largely Jurassic island arc [Jones et al., 1977; Muller, 1977]. The mainly metasedimentary Mesozoic Pacific Rim terrane and the volcanic Eocene Crescent Terrane lie along the west coast and southern end of Vancouver Island and were the last to accrete to the continent and reached their present locations during late Cretaceous and Tertiary periods [Johnson, 1984]. To the south of Vancouver Island, the Strait of Juan de Fuca lies in a synclinal depression formed in the Crescent Terrane. To the east of Vancouver Island, the Strait of Georgia is a fore-arc basin that straddles the boundary of the Insular and Coast belts.
 

The NW Washington margin is comprised of the Coast Range province to the west and the Cascade Range province to the east (Figure 2). The Coast Range province includes the Olympic subduction complex and Puget Sound. The Olympic Subduction Complex is an exposed accretionary wedge, metamorphosed and uplifted as a result of subduction of the Juan de Fuca Plate [e.g., Brandon and Calderwood, 1998]. To the east of the Olympic Subduction Complex, the Puget Sound is a fore-arc basin. Coast Range basement beneath Puget Sound is comprised of the Eocene Crescent terrane rocks, probably formed in a continental margin rift setting [Wells et al., 1984; Babcock et al., 1992]. The Crescent terrane along the eastern margin of the Olympic Subduction Complex is tilted to the east by the underthrusting and uplift of the accretionary sediments [e.g., Tabor and Cady, 1978; Brandon and Calderwood, 1990]. The unexposed Coast Range Boundary fault (CRBF) lies on the eastern margin of the Coast Range Province [Johnson, 1984, 1985]. In northern Puget Sound the CRBF becomes the northwest trending southern Whidbey Island fault (SWIF) [Johnson et al., 1994, 1996]. The Leech River fault in southern Vancouver Island is postulated to connect to SWIF [Johnson et al., 1996]. The Cascade Range province to the east of Puget Lowland comprises of a variety of igneous, sedimentary and metamorphic rocks that comprise several distinct crustal terranes [e.g., Tabor, 1994].
 

The Juan de Fuca oceanic plate converges with the North America continental plate at a relative rate of ~46 mm/yr directed N56°E [Riddihough and Hyndman, 1991]. Subduction of the Juan de Fuca plate has resulted in a number of damaging earthquakes [e.g., Rogers, 1998]. Earthquakes occur in three distinctive source regions: (1) earthquakes in the North American continental crust, (2) earthquakes in the subducting oceanic plate, and (3) megathrust earthquakes at the interface of the oceanic and the North American plates. The relatively large number of earthquakes in the continental crust are driven by a north-northwest compressive stress parallel to the continental margin [e.g., Rogers, 1998]. The deeper earthquakes, in the depth range of 40–70 km, occur due to a tensional stress regime within the subducted oceanic plate. Paleoseismic evidence shows that megathrust earthquakes have occurred at intervals of 500–600 years [Atwater and Hemphill-Haley, 1997]; the most recent event occurred in 1700 [Satake et al., 1996].


3. Data

[7] Data from active source experiments and earthquake recordings were used to construct a three-dimensional P wave minimum structure velocity model. The active source data constrains the upper crustal velocities that have large lateral variations between the sedimentary basins and the basement rocks made up of accreted terranes. The well-constrained upper crustal velocity structure improves the estimation of deeper velocities and earthquake hypocentral parameters.

3.1. Active Source Data

[8] The 1998 SHIPS experiment [Fisher et al., 1999; Brocher et al., 1999] recorded arrivals from a total of 33,000 air gun shots fired on 11 shot lines in the waterways of the Strait of Georgia, the Strait of Juan de Fuca, and Puget Sound (Figure 1). The shots were recorded widely over southwestern British Columbia and northwestern Washington (Figure 1) at 257 temporary land-based, Reftek stations [Incorporated Research Institutions for Seismology, 1991], and 15 ocean bottom seismometers at offsets from 1 to 370 km [Brocher et al., 1999]. Approximately 147,000 first arrival time picks from 225 temporary recording stations were employed in the inversion. In addition, 1000 first arrival traveltimes from a refraction line acquired in 1991 [Miller et al., 1997] were also included. No wide-angle reflection traveltimes were included in the inversion.

3.2. Earthquake Data

[9] Approximately 70,000 first arrival times from nearly 3000 earthquakes that occurred beneath the SW British Columbia and NW Washington margin in the past 25 years and were recorded at 91 permanent recording stations (Figure 1, bottom left inset) were used in the tomographic inversion. Selection criteria included earthquakes from any depth that were recorded by at least six stations within the study region. Hypocentral data from the earthquake catalogs were used as initial hypocentral parameters for the tomographic inversion and were updated during tomographic inversion.

4. Tomographic Inversion

Figure 3. (a) Traveltime misfit of active source data for initial 1-D velocity model, (b) traveltime misfit of active source data for final velocity model, (c) traveltime misfit of earthquake data for initial 1-D velocity model, and (d) traveltime misfit of earthquake data for final velocity model.

The velocity model was parameterized in the forward and inverse steps by a node and cell spacing of (1.5 × 1.5 × 1.5) km and (4.5 × 4.5 × 1.5) km, respectively. The velocity model dimensions in (x, y, z) directions are (460 × 490 × 96) km. The top of the model is set to 3 km above sea level to allow positioning the receivers at their actual elevations in the velocity model. An initial 1-D velocity model (Table 1) was constructed by forward modeling the first arrival traveltime data. The traveltime misfit for the initial model as a function of offset for the SHIPS data and the earthquake events is shown in Figures 3a and 3c, respectively.

The RMS traveltime residual for this model for approximately 210,000 observations was 762 ms for a normalized c² of 52. During the tomographic inversion, the hypocentral parameters and the velocity model were updated after each iteration. Ray tracing was performed during each iteration to account for the change in the hypocentral parameters between iterations. A second-order smoothing constraint in the horizontal and vertical direction was applied during the inversion. The ratio of the vertical to horizontal smoothing parameter was set to 0.2 (i.e., three times more smoothing in the horizontal direction than vertical direction) and was held fixed throughout the inversion. After 30 iterations, a stable minimum was obtained with an RMS misfit of ~132 ms and normalized c² of ~1.1, which represents a 97% variance reduction. The traveltime misfit with offset for the final model for the SHIPS data and the earthquake events is shown in Figures 3b and 3d, respectively.

4.1. Ray Hit Count and Checkerboard Tests

Figure 4. Depth slices at 3, 9, 15, 21, and 27 km depth of (a) ray hit count, (b) checkerboard test recovered anomaly pattern, and (c) semblance values, for 30 km grid size.

Figure 5. Depth slices at 21, 39, 45, 51, and 57 km depth of (a) ray hit count, (b) checkerboard test recovered anomaly pattern, and (c) semblance values, for 50 km grid size.

The ray hit count method offers a quick and easy way to assess the velocity constraint for each cell. Rays were traced though the final velocity model with the actual receiver geometry, relocated earthquake positions, and active source positions to determine the number of rays that pass through each cell. The ray hit counts for the 3–27 km and 21–57 km depth ranges are shown in Figures 4a and 5a, respectively. The results show significant ray coverage in the upper and middle crust down to 20 km, where most of the rays are from the SHIPS active source data set. The ray coverage in the mid crust and below is primarily from the earthquake data and is sparser. However, since these rays from deeper earthquakes travel through the well constrained upper crustal velocity model, the velocities modeled for the lower crust and upper mantle are inferred to be better constrained than they would be if only earthquake data were available.

Checkerboard tests were carried out with grid sizes of 30, 40 and 50 km to access the ability of the data to resolve model features of these sizes at different depth levels. Semblance values measuring the correlation between the input and recovered patterns, as discussed by Zelt and Barton [1998], were used to classify model volumes having reasonable lateral resolution. Semblance values of 1 and 0.5 indicate total recovery of the velocity perturbation and no recovery of the velocity perturbation, respectively. Semblance values less than 0.5 indicate negative correlation. Semblance values between the initial and recovered checkerboard anomaly patterns were computed with a window size of 16.5 × 16.5 × 3.0 km.

The recovered checkerboard anomaly pattern for a grid size of 30 km between depths of 0–21 km (Figure 4b) is characterized by semblance values of 0.7 and above (Figure 4c), indicating adequate resolution. The sedimentary basins in the Straits of Georgia and Juan de Fuca, and the Puget Lowland show good recovery of the checkerboard anomaly patterns in the depth slices at 3 and 9 km. The low semblance values on the fringe of the model (Figure 4c) are inferred to be due to limited ray directivity and fewer rays at the edges. The recovered checkerboard anomaly pattern and semblance values for a grid size of 40 km indicate adequate resolution down to 45 km depth for features of this size. The recovered checkerboard anomaly pattern for a grid size of 50 km (Figure 5b) suggests that the ray coverage is sufficient to recover features of this size to 57 km depth. Semblance plots (Figure 5c) indicate adequate lateral resolution, in regions with ray coverage, at all depth levels west of 122.25°W.

5. Results and Discussion

Figure 6. Vertical cross sections (a) AB, (b) CD, and (c) EF. The cross sections are approximately in the margin-perpendicular direction. Location of the cross sections is shown on Figure 2. Abbreviations are as in Figure 2.

Figure 7. Vertical cross sections (a) GH, (b) IJ, and (c) KL. The cross sections are approximately in the margin-perpendicular direction. Location of the cross sections is shown on Figure 2. Abbreviations are as in Figure 2.

Figure 8. Vertical cross sections (a) PQ and QR and (b) MN and NO. The cross sections are approximately in the margin-parallel direction. Location of the cross sections is shown on Figure 2. Abbreviations are as in Figure 2.

Figure 9. Horizontal cross sections at (a) 27 km, (b) 33 km, (c) 36 km, (d) 45 km, (e) 51 km, and (f) 57 km depth.

Figure 10. a) Horizontal cross section at 3 km depth showing the outline of sedimentary basins. (b) Gravity anomaly map (data from National Geophysical Data Center [1999]). The thick red dashed line marks the approximate position of the junction of the subducting Juan de Fuca crust and the fore-arc mantle wedge. Abbreviations are as in Figure 2.

Vertical cross sections along selected orientations (Figures 6, 7, and 8) and horizontal (depth) slices (Figure 9) define the structure of the subducting Juan de Fuca plate. A horizontal velocity slice at 3 km depth showing the inferred outline of the younger sedimentary basins (Figure 10a), and a gravity anomaly map of the region (Figure 10b) indicate the setting of the sedimentary basins in relation to the forearc mantle wedge.

5.1. Lower Crustal Low-Velocity Zone

A broad zone of low velocities (6.2–6.6 km/s) is imaged in the lower crust, between 25 and 35 km depth, beneath Vancouver Island on the vertical cross sections AB, CD, EF, PQR, and MNO (Figures 6 and 8). The lateral extent of the low-velocity zone above the subducting plate is seen in the horizontal velocity slices at 27 and 33 km (Figure 9a and 9b). Cassidy and Ellis [1991] identified a crustal low-velocity zone between 20 and 26 km depth at station ALB (Figure 6a) from receiver function studies.

This broad low-velocity zone in the lower crust correlates with an approximately 5–8 km thick band of reflectivity in the lower crust beneath Vancouver Island identified from the 1984 Vancouver Island LITHOPROBE Vibroseis reflection lines [Clowes et al., 1987]. This band of reflective zone is generally referred to as the E reflectivity. Hyndman [1988] suggested that fluids released from the dehydration reactions occurring in the subducting slab can be trapped at this depth level and account for the seismic E reflectors. Calvert and Clowes [1990] suggested a shearing mechanism for these reflectors.

Earthquake hypocenters do not fall within this low-velocity zone (Figures 6 and 8), and it is inferred that any accumulated stress there is probably released by episodic, aseismic slow-slip events [Dragert et al., 2001]. Rogers and Dragert [2003] identified nonearthquake, tremor-like signals accompanying the episodic aseismic slip. The episodic tremors are distributed over a wide depth range of ~40 km with a peak at 25–35 km, and many of these events occur within or in the close vicinity to the E reflectors [Kao et al., 2005]. Approximately 50% of the episodic tremor events are located within or close to the E reflectors [Kao et al., 2005], whereas >90% of the local earthquakes tend to be located away from the reflectors [Calvert, 2004]. Kao et al. [2006] suggest that shear deformation and fluids may be closely related to the occurrence of episodic tremor. The low velocities imaged between 25 and 35 km depth beneath Vancouver Island coincides with the zone of postulated mechanisms of shearing, and presence of fluids inferred from EM measurements [Hyndman, 1988].

5.2. Olympic Subduction Complex

Eocene and younger clastic sedimentary rocks scraped off the subducting Juan de Fuca plate are accreted to North America to form a thick accretionary prism along the Cascadia subduction zone margin, mainly beneath the shelf and the continental slope [Tabor and Cady, 1978]. This accretionary prism has been uplifted, metamorphosed and exposed in the Olympic Peninsula [e.g., Brandon and Calderwood, 1990] (Figure 2). The low-velocity rocks of the Olympic Subduction Complex exhibit a sharp velocity contrast with the neighboring and overlying Crescent terrane rocks, as can be clearly observed on cross sections GH and IJ (Figure 7). A similar velocity contrast is also seen in laboratory studies; the Crescent terrane volcanic rocks exhibit a velocity contrast of up to 1.0 km/s with the low-velocity core rocks from the Olympic Peninsula [Brocher and Christensen, 2001]. From cross sections GH, IJ, and KL (Figure 7), we infer that the low-velocity rocks of the Olympic Subduction Complex underthrust the Crescent terrane and extend to at least 30 km depth. The Olympic Core rocks appear to extend downward to the top of the Juan de Fuca plate.

The eastward underthrusting of the low-velocity rocks of the Olympic Subduction Complex deforms the overlying Crescent Terrane and is probably the source of some seismicity in the depth range of 15 to 30 km (Figure 7b, model distance 50–100 km; Figure 8a, line QR model distance 20–130 km). However, earthquake mechanisms in the fore-arc [e.g., Wang et al., 1995] and GPS data [e.g., Hyndman et al., 2003] indicate northward compression of the fore-arc against the British Columbia buttress as the main source for the earthquakes in this region. The low-velocity rocks of the Olympic Subduction Complex themselves are inferred from our study to deform aseismically as evidenced by the lack of earthquakes within the low-velocity rocks.

5.3. Location of the Fore-arc Basins and Slab Dehydration

Locations of sedimentary basins outlined as velocity anomaly lows in the horizontal velocity slice at 3 km depth (Figure 10a) and as gravity anomaly lows (Figure 10b) show a strong correlation. The sedimentary basins are classified into two groups (Figure 10b); several larger basins in the Strait of Georgia and the Puget Sound, and several smaller basins in the Strait of Juan de Fuca. The sedimentary basins in the Strait of Georgia and the Strait of Juan de Fuca were previously detailed from 3-D tomographic P wave velocity models by Zelt et al. [2001] and Ramachandran et al. [2004]. Brocher et al. [2001], Tréhu et al. [2002], and Van Wagoner et al. [2002] mapped the locations of the sedimentary basins in Puget Sound from 3-D velocity models. The velocity model from the present study provides the first contiguous image of all the sedimentary basins in the Strait of Georgia, the Strait of Juan de Fuca, and the Puget Lowland (Figure 10a).

The thick red dashed line shown in Figure 10b marks the approximate position of the junction of the subducting Juan de Fuca crust and the fore-arc mantle wedge. East of this boundary, the basins in the Strait of Georgia and Puget Sound lie above the zone of inferred fore-arc mantle serpentinization. Dehydration, eclogitization, and densification of the slab crust [Peacock, 1993] and slab mantle dehydration and densification [e.g., Hacker et al., 2003] decrease the buoyancy of the oceanic plate with respect to the surrounding mantle [Rogers, 1983]. This phase change may be reflected as a small increase in the angle of subduction of young plates in the subduction zones of Cascadia, southern Chile, and the Nankai region of SW Japan [e.g., Rogers, 2002]. The change of subduction angle at such a shallow depth may foster fore-arc basin subsidence. Subcrustal erosion and cooling of the fore-arc continental crust may also contribute to fore-arc basin subsidence.

In contrast to the larger basins in the Strait of Georgia and Puget Sound, the small Clallam and Sequim basins in the Strait of Juan de Fuca lie in the synclinal depression formed in the Crescent Terrane. The synclinal depression was probably formed in the Crescent Terrane due to folding and faulting in the Strait of Juan de Fuca and Olympic Peninsula to accommodate the northward motion of the Coast Range block along the Coast Range Boundary fault [Snavely, 1987].

5.4. Juan de Fuca Plate Position

Our P wave velocity model is a significant improvement over previous velocity models because of the inclusion of a large set of earthquake picks from a broader distribution of permanent stations and a denser distribution of active source traveltimes. In the velocity model constructed by inverting first arrival traveltimes, the oceanic Moho is not expected to appear as a sharp boundary. Given the vertical cell size of 1.5 km, sharp velocity boundaries will be smeared over adjacent cells to either side in the vertical (depth) direction. However, the velocity model would be characterized by the gradient and the geometrical structure of the sharp boundary. This feature is utilized to infer the position of the oceanic Moho.

The velocity model images oceanic mantle rocks which are inferred to have velocities greater than ~7.6 km/s. This isocontour is assumed to define the top of the oceanic Moho. The top of the Juan de Fuca crust is drawn as a smooth surface approximately 7 km above the oceanic Moho, assuming an average oceanic crustal thickness. Using these constraints, the interpreted position of the top of the Juan de Fuca crust and mantle over a relatively large geographic region are shown on the cross sections in Figures 6, 7, and 8. The inferred position of the oceanic Moho beneath the Strait of Juan de Fuca and Olympic Peninsula is consistent with the depth estimates obtained by Tréhu et al. [2002] and Preston et al. [2003] through the analysis of first arrival traveltime data and wide-angle reflection data.

Estimates of crustal thickness and slab geometry from Figures 6, 7, and 8 match those of several previous workers, including Hyndman et al.'s [1990] interpretation that a short reflection segment at 10 s two-way time on LITHOPROBE line 84-01, referred to as the F reflection, represents the top of the Juan de Fuca crust. Our estimate for the depth to the oceanic crust agrees with that suggested by Drew and Clowes [1990] from refraction data modeling. These depth estimates were also substantiated by Calvert [2004] through analysis of reflection and refraction seismic data from Vancouver Island and adjoining waterways. Our estimates are also consistent with the depth estimates derived from receiver function analysis by Cassidy and Ellis [1993] and Cassidy [1995].

The agreement, however, is not universal. Nedimovic' et al. [2003] suggested that the position of the Juan de Fuca crust is geometrically below the E reflectors and is shallower than previous estimates by at least 6 km beneath Vancouver Island. If the depth to the top of the crust is brought shallower by 6 km, as suggested by Nedimovic' et al. [2003], the oceanic mantle would be represented by velocities in the range of ~7.2 km in our velocity model, which then would be indicative of extensive serpentinization of the oceanic mantle rocks. At this time, there is not enough evidence to suggest extensive oceanic mantle sepentinization in this region to account for a oceanic mantle velocity of ~7.2 km.

Through the analysis of P coda from teleseismic events, Nicholson et al. [2005] interpreted the position of the Juan de Fuca crust to coincide with the E reflection zone. Their estimates for the depth to the top of the Juan de Fuca crust beneath Vancouver Island are shallower than previous estimates by at least 10 km. It is increasingly difficult to correlate our velocity model with the position of the Juan de Fuca plate suggested by Nicholson et al. [2005] for the same reasons discussed for the model of Nedimovic' et al. [2003].

5.5. Juan de Fuca Plate Seismicity

It has been previously suggested that oceanic slab earthquakes at depths below ~35 km are induced by dehydration-embrittlement processes in the slab crust [e.g., Kirby et al., 1996] and in the slab mantle [e.g., Peacock, 2001; Hacker et al., 2003]. On the vertical cross sections (Figures 6 and 7), earthquakes occurring in the region with velocities in the range of 7.6–8.0 km/s are inferred to occur in the slab mantle. From relocated slab earthquakes beneath the Strait of Georgia, Cassidy and Waldhauser [2003] demonstrated that some of the slab earthquakes originated in the oceanic mantle. Beneath the Olympic Peninsula and Puget Sound, a significant number of in-slab earthquakes lie in regions with velocities between 7.6 and 8.0 km/s [Preston et al., 2003, supplemental online material, Table S1]. Almost fully eclogitized slab crust and weakly serpentinized upper mantle are expected to be in this velocity range (7.6 and 8.0 km/s).

Peacock [2001] proposed that seawater percolating along fault zones at the outer rise could hydrate and serpentinize the slab's upper mantle. This process may reduce the velocity of the oceanic mantle below the normal ~8.2 km/s on Figures 7, 8, and 9. As observed on the horizontal slice at 45 and 51 km depth (Figures 6, 7, 8, and 9d and 9e), the level of seismicity inferred in the slab mantle beneath the Washington margin is higher than the seismicity inferred along the British Columbia margin. This increased level of seismicity in the slab mantle beneath the Washington margin may reflect higher levels of oceanic mantle deserpentinization and fluid expulsion, facilitating seismic rupture.

The distribution of earthquakes used in this study (Figure 1, bottom left inset) is representative of the occurrence of earthquakes in the fore-arc upper crust and the Juan de Fuca slab [e.g., McCrory et al., 2004]. The in-slab earthquakes occurring below 35 km depth and inferred to be within the Juan de Fuca slab are summarized by depth range in Table 2. The earthquakes within a given depth range are further classified into three different groups according to the P wave velocity at the hypocenter: 6.8–7.2 km/s (slab upper crust with hydrous minerals), 7.2–7.8 km/s (mid to lower slab crust, partially eclogitized slab upper crust, and serpentinized upper mantle), and 7.8–8.2 km/s (slab mantle that may be partially serpentinized).

The crust and uppermost mantle of warm slabs dehydrate at shallow depths [e.g., Hacker et al., 2003; Peacock and Wang, 1999; Currie et al., 2002]. Thermal modeling for warm slabs like Juan de Fuca show that metamorphic reactions in the subducting oceanic crust can start as shallow as 40–50 km depth [Peacock et al., 2002]. Depending on the amount of mantle serpentinization, the hypocentral regions of in-slab earthquakes in the subducting mantle will have a velocity of 7.8–8.2 km/s. In contrast, the hypocentral regions of shallow earthquakes in the slab crust are expected to have a seismic velocity range of 6.8–7.2 km/s. Between 35–40 km depth in Table 3, the number of hypocentral regions having velocities between 6.8 and 7.2 km/s is lower than the number of hypocentral regions having velocities between 7.2 and 7.8 km/s. This finding suggests that more earthquakes occur in the slab's partially eclogitized crust than in unaltered crust. Hypocentral regions having velocities between 7.8 and 8.2 km/s velocity zone are inferred to occur within the slab's nearly fully eclogitized crust and/or partly serpentinized slab mantle. Between 40 and 50 km depths, significantly larger number of hypocenters occur in the region having a velocity range between 7.8 and 8.0 km/s (Table 3). This result suggests that these earthquakes occur close to the slab Moho. At 50–60 km depth interval, hypocenters in the 7.8–8.2 km/s velocity zone dominate the distribution, and below 60 km depth no hypocenters are observed in the velocity range of 6.8–7.8 km/s. Hence it is inferred that the in-slab earthquakes between 50 and 60 km depth occur in the nearly fully eclogitized oceanic crust and/or the oceanic upper mantle.



 

Figure 11. Location of the nine earthquakes (M > 4) inferred to originate close to the slab Moho (Table 4), shown by the stars. The depth to the top of the oceanic Moho is shown by the dashed line contours computed from the slab depth values (km) from McCrory et al. [2004] by adding an average oceanic crustal thickness of 7 km. Beneath southern Vancouver Island the slab depth and hence oceanic Moho depth contours are modified according to the results inferred from the present study.

Cassidy and Waldhauser [2003] showed that large earthquakes in the British Columbia margin do occur in the uppermost mantle, and there is little activity in the lower slab crust. Wada et al. [2004] showed that many events along the Nankai margin, SW Japan occur in the slab mantle. The recent, damaging earthquakes in the Cocos plate (Oaxaca, 1999, M 7.5), the Philippine Sea plate (Geiyo, 2001, M 6.7), and the Juan de Fuca plate (Nisqually, 2001, M 6.8) are inferred to have occurred close to the slab Moho [e.g., Wang et al., 2004]. Tréhu et al. [2002] inferred that the hypocenters of the damaging lower plate earthquakes in the Olympic Peninsula have hypocenters 0–2 km below the Moho of the subducting oceanic plate. Preston et al. [2003] inferred that earthquakes updip of the Moho's 45-km depth contour occur in the subducted oceanic mantle, and earthquakes located downdip of the Moho's 45 km depth contour occur primarily within the subducted crust. In our study, all nine earthquakes of M > 4 (Table 3) are located between 40 and 55 km depth (Table 4), and are characterized by a velocity of 7.8–8.1 km/s at the hypocenter. The seismic velocities at the hypocenters of these nine earthquakes indicate that these earthquakes most likely originated close to the subducting slab Moho. This observation is consistent with the slab depth contours (Figure 11) of McCrory et al. [2004].

Our finding that large in-slab earthquakes in the Juan de Fuca plate occur close to the slab Moho is also consistent with the hypothesis of Wang et al. [2004]. Wang et al. [2004] proposed that crustal densification in the upper few hundred meters of the oceanic slab crust will shatter the slab crust resulting in only smaller magnitude earthquakes, whereas large in-slab earthquakes have a tendency to occur deep inside the subducting slab. Dehydration along existing faults near the slab's lower crust and upper mantle can facilitate seismic rupture that can propagate for large distances, resulting in large in-slab earthquakes [Wang et al., 2004].

6. Summary

Our regional 3-D P wave velocity model presents the first regional image of the northern Cascadia subduction zone beneath SW British Columbia and NW Washington state. The megathrust zone above the Juan de Fuca plate is characterized by broad zone of low P wave velocities in the range of 6.4–6.6 km/s between 25 and 35 km depth and lies down dip of the megathrust locked zone. This low-velocity zone is virtually aseismic but lies within the region experiencing episodic tremor and slip. This region also coincides with the high electrical conductivity region mapped in previous magnetotelluric studies, a landward dipping band of seismic reflectors, and a low shear velocity layer modeled in receiver function studies. These low velocities are inferred to be due to either trapped fluids, highly sheared lower crustal rocks, and/or underthrust accretionary rocks.

A broad low-velocity zone in the fore-arc mantle having velocities between 7.2 and 7.6 km/s, imaged at depths of 37–45 km above the Juan de Fuca plate, is inferred to reflect a serpentinized upper mantle. This zone is similar in velocity to the low-velocity zones observed in the fore-arc mantle beneath central Japan, northern Costa Rica, eastern Aleutians, and the Andes. The younger sedimentary basins in the Strait of Georgia and Puget Sound are imaged above the zone of inferred fore-arc mantle serpentinization.

The large in-slab earthquakes (M > 4) occurring at depths between 40 and 55 km are inferred to originate close to the subducting slab Moho. Earthquakes occurring closer to the slab Moho are probably not limited in size by the thickness of the slab crust; the faulting could extend from the slab crust into the slab upper mantle with larger rupture surfaces than earthquakes that are restricted solely to the upper crust of the subducting slab.