U.S. Department of the Interior
U.S. Geological Survey
SUBDUCTION ZONE AND CRUSTAL DYNAMICS OF WESTERN WASHINGTON: A TECTONIC MODEL FOR EARTHQUAKE HAZARDS EVALUATION
By Dal Stanley, Antonio Villaseñor, and Harley Benz
Denver Federal Center, Denver, CO 80225
Phone 303-273-8620 email dals@usgs.gov
Open-File Report 99-311 On-line Edition
This report is preliminary and has not been reviewed for conformity with U.S. Geological Survey editorial standards or with the North American Stratigraphic Code. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
1999
Web construction by Phil Powers, USGS, Central Region Hazards Team webmaster
Geological Setting
Gravity and Magnetic Maps of Washington
Magnetotelluric Profiles
Faults and Seismicity
Velocity Model Construction
Velocity Model 2D Cross-sections
Juan de Fuca Plate
Olympic Mts. Sedimentary Complex
Juan de Fuca Plate Arch
Seattle Basin and SWCC
Crustal Terrane Boundary
3D Visualization of Velocity Results
Sedimentary Basins and SWCC
General Results From 3D Visualization of Cross-sections
S-wave Model
Interpretive Geological Cross-sections
East-west Transect Interpretation
Subduction Thrust Frictional States
North-South Transect Interpretation
Comparison of Cascadia With Other Subduction Zones
Implications for Interplate and Crustal Seismic Hazards
Subduction Zone Source Ground Motions
Seismicity Patterns
Deformation Modeling
Model Construction
Model Traction Limits
Model Stress Results
Pn Anistropy
Plate Interactions
Surface Velocities and Fault Slip Rates
Stress Origins and Seismicity in the Puget Lowland
1100 yr B.P. Deformation
Appendix A-Resolution Tests
SUBDUCTION ZONE AND CRUSTAL DYNAMICS OF WESTERN WASHINGTON: A TECTONIC MODEL FOR EARTHQUAKE HAZARDS EVALUATION
By
Dal Stanley, Antonio Villaseñor, and Harley Benz
Abstract: The Cascadia subduction zone is extremely complex in the western Washington region, involving local deformation of the subducting Juan de Fuca plate and complicated block structures in the crust. It has been postulated that the Cascadia subduction zone could be the source for a large thrust earthquake, possibly as large as M9.0. Large intraplate earthquakes from within the subducting Juan de Fuca plate beneath the Puget Sound region have accounted for most of the energy release in this century and future such large earthquakes are expected. Added to these possible hazards is clear evidence for strong crustal deformation events in the Puget Sound region near faults such as the Seattle fault, which passes through the southern Seattle metropolitan area. In order to understand the nature of these individual earthquake sources and their possible interrelationship, we have conducted an extensive seismotectonic study of the region. We have employed P-wave velocity models developed using local earthquake tomography as a key tool in this research. Other information utilized includes geological, paleoseismic, gravity, magnetic, magnetotelluric, deformation, seismicity, focal mechanism and geodetic data. Neotectonic concepts were tested and augmented through use of anelastic (creep) deformation models based on thin-plate, finite-element techniques developed by Peter Bird, UCLA. These programs model anelastic strain rate, stress, and velocity fields for given rheological parameters, variable crust and lithosphere thicknesses, heat flow, and elevation. Known faults in western Washington and the main Cascadia subduction thrust were incorporated in the modeling process.
Significant results from the velocity models include delineation of a previously studied arch in the subducting Juan de Fuca plate. The axis of the arch is oriented in the direction of current subduction and asymmetrically deformed due to the effects of a northern buttress mapped in the velocity models. This buttress occurs under the North Cascades region of Washington and under southern Vancouver Island. We find that regional faults zones such as the Devils Mt. and Darrington zones follow the margin of this buttress and the Olympic-Wallowa lineament forms its southern boundary east of the Puget Lowland. Thick, high-velocity, lower-crustal rocks are interpreted to be a mafic/ultramafic wedge occuring just above the subduction thrust. This mafic wedge appears to be jointly deformed with the arch, suggesting strong coupling between the subducting plate and upper plate crust in the Puget Sound region at depths >30 km. Such tectonic coupling is possible if brittle-ductile transition temperatures for mafic/ultramafic rocks on both sides of the thrust are assumed.
The deformation models show that dominant north-south compression in the coast ranges of Washington and Oregon is controlled by a highly mafic crust and low heat flow, allowing efficient transmission of margin-parallel shear from Pacific plate interaction with North America. Complex stress patterns which curve around the Puget Sound region require a concentration of northwest-directed shear in the North Cascades of Washington. The preferred model shows that greatest horizontal shortening occurs across the Devils Mt. fault zone and the east end of the Seattle fault.
The northern Cascadia subduction zone along Oregon, Washington, and Vancouver Island (Figure 1) has been hypothesized as the potential location for a large megathrust earthquake (Heaton and Hartzell, 1986,1987; Weaver and Shedlock, 1996). Evidence for past great earthquakes on the subduction zone is extensive (Atwater et al., 1995) and includes signs of land-level changes, tsunamis (Satake et al., 1996), and liquefaction. In addition, geodetic (Mitchell et al., 1994), GPS, and strain (Savage et al., 1991) data provide evidence for recent stress buildup in the subduction zone. However, there is no historically recorded evidence for interplate earthquakes of significant magnitude (M>6), although native American legends tell of large pre-1800 earthquakes in Cascadia (Carver and Carver, 1996). Radiocarbon and tree-ring dating of drowned forests, assumed to be associated with major subduction zone earthquakes, show that the most recent large such event occurred about 1700 AD (Atwater et al., 1995). Satake et al. (1996) interpret that a M9.0 earthquake along the Cascadia subduction zone is a likely source of the large tsunami that struck in Japan during January, 1700 AD.
Although it is generally accepted that convergence between the Juan de Fuca plate (JDFP) and North American plate (NAP) occurs at a rate of about 30-46 mm/yr (Wilson, 1993), there is considerable uncertainty pertaining to the seismic or aseismic behavior of the megathrust and possible length and location of segments which might generate large earthquakes. Rogers et al. (1996) provide a review of the various models for the megathrust and possible earthquakes and resulting ground motions. Besides convergence rate, many other factors influence the location of source regions and moment release in great subduction zone earthquakes. Most analyses of probable source regions in the northern Cascadia subduction zone have been based upon heat flow models (Hyndman and Wang, 1993), deformation in the accretionary prism (Davis and Hyndman, 1989), and coastal uplift data (Dragert et al., 1994).
Frictional states along strike and downdip within the subduction thrust zone are key parameters required for quantification of seismic source regions. Frictional states are related to thrust zone temperatures that are unknown in the deeper parts of the subduction interface, but are also dependent upon deep fluid flow and lubricating minerals in the thrust. Subducted sediments may lubricate the thrust zone downdip from the accretionary prism, but there is little agreement on the actual amount of deeply subducted sediments and no accurate constraints on their pore-pressures or metamorphic grade, factors that determine seismic strength of the interface. Onshore uplift data has been interpreted to represent elastic deformation due to offshore subduction thrust coupling. However, without knowledge of offshore vertical motions, modeling of subduction zone dynamics from onshore uplift data alone is somewhat qualitative. Geodetic data indicate that active strain is occurring in the northeast subduction convergence direction (Savage et al., 1991), but focal mechanism and borehole breakout data (Zoback and Zoback, 1991) map a pattern of general north-south compression. Such complexities need to be better understood.
In this report we analyze the onshore part of the subduction zone and overlying crust in northern Cascadia using velocity models from seismic tomography, deformation models, and other data. We present arguments for crust-slab coupling at depths deeper than those previously assumed and discuss physical conditions which would allow this deep coupling. We show that joint deformation of the arch in the Juan de Fuca plate and mafic rocks just above the subduction thrust require strong coupling in the past. We discuss the impact of this interpreted deep coupling on earthquake hazards in the region. We also compare Cascadia with other subduction zones to stress the role of key subduction zone features, primarily utilizing a close analog from the western Mexico subduction zone.
The complexity of the velocity models and other data that we present requires a large number of color illustrations. Many of these illustrations are best viewed on high resolution color monitors, since color printers do not always adequately display shading parameters and other features. Adobe Acrobat (pdf) and postscript files for this report are contained in the files index.pdf (279Kb) and index.ps(936Kb), respectively. The Adobe Acrobat reader utility can be downloaded at http://www.adobe.com/products/acrobat/readstep.html. The full-resolution gif files for the illustrations are available as a UNIX tar file'figures.tar'. For hard copy versions of this report, only grayscale versions of the figures at medium resolution are available.
The northern Cascadia region is identified by the subducting Juan de Fuca plate offshore and the magmatic arc onshore (Figure 1). The Quaternary magmatic arc extends from northern California to British Columbia and includes large stratovolcanoes such as Mt. St. Helens, Mt. Hood, Mt. Shasta, Mt. Baker, and Mt. Rainier. An earlier Oligocene to Miocene volcanic arc is displaced somewhat westward, but overlaps the Quaternary arc. The Coast Ranges of Oregon and Washington are composed of voluminous basaltic flows of Eocene age that apparently formed in a nearshore pull-apart system, possibly in conjunction with a leaky transform (Wells et al., 1984; Babcock et al., 1992). These basalt flows (Figure 1) are known as the Siletz River volcanics in Oregon, Crescent Formation in Washington and Metchosin Formation on Vancouver Island. The extensive basalt terrane has been referred to as Siletzia (after Siletz River volcanics) by Irving (1979). The easternmost extent of Siletzia interpreted in Figure 1 was estimated using gravity and magnetic data (Finn, 1990). Pre-Tertiary rocks in Cascadia outcrop in the Klamath Mts. of northern California, in the Blue Mts. of Oregon, but most extensively in northcentral Washington, Vancouver Island, and in British Columbia. The trend of pre-Tertiary rocks in NE Oregon and SE Washington is obscured by thick (up to 10 km) Miocene volcanic rocks of the Columbia River basalt group. Several forearc basins filled with Tertiary marine sediments occur in western Washington and Oregon, formed on downwarps in the Crescent and Siletz River basement. In addition, thick nonmarine sediments of Tertiary age fill pull-apart and thrust-bounded basins in western Washington, most notably the Seattle and Tacoma basins.
There are numerous faults in western Oregon and Washington, but few appear to be seismically active (Stanley et al., 1996). The faults that appear to be most active (Figure 1) are the Portland Hills fault zone (PHFZ), St. Helens zone (SHZ) , western Rainier zone (WRZ), Seattle fault zone (SFZ), South Whidbey Island fault zone (SWF), Devils Mt. (DMF) and Darrington fault zones (DAR). In some figures of this report other than Figure 1 which show the DMF and DAR, we follow the locations as represented by Tabor (1994). These locations seem to follow deep velocity expression of these fault zones. However, the surface trace of the DMF in northern Puget Sound has recently been interpreted by Johnson et al. (1999b) from seismic reflection data to be as indicated in Figure 1. The Seattle and South Whidbey Island fault zones have been the focus of increased recent study because of the possible large impact on Seattle from a crustal earthquake on either of these zones (Johnson et al., 1994,1996). The Seattle fault zone passes through the southern part of the Seattle metropolitan area and may have been involved in one or more large earthquakes between 500 and 1700 years ago (Bucknam et al., 1992,1996). The evidence for these large earthquakes in the Seattle region includes uplifted shorelines, tsunamis, landslides, and turbidity currents in the southern Puget Lowland. The uplifted shoreline evidence for large earthquakes affecting the southern Puget Sound region between 500 and 1700 years ago can be narrowed to between 1000 and 1100 years age using radiocarbon dating and tree-ring correlation of tsunami and landslide deposits near Seattle, under the assumption that these deposits are coeval with shoreline uplifts (Bucknam et al., 1992). This range corresponds to that of one of three periods of large earthquakes determined from sand blows along the Washington coastline. Atwater (1992) found deeply sourced sand blows near the mouth of Copalis River in coastal Washington that occurred from 900-1300 years ago. Atwater (1992) and Bucknam et al., (1992) indicate that this latter evidence may be correlated with a crustal fault, such as the Seattle fault, because of the coincidence of timing with mapped uplift in southern Puget Sound.
Active tectonic forces are most evident in the Olympic Peninsula, where a rim of Crescent Formation has been deformed into a tight u-shape and thick Cenozoic marine sedimentary rocks (Figure 1) have been uplifted approximately 10 km (Brandon and Calderwood, 1990). Brandon and Calderwood (1990) relate the uplift to an arch in the underlying JDFP which was first mapped by Crosson and Owens (1987). They interpret that the arch was produced by oroclinal bending of the continental margin and underlying plate due to Basin and Range spreading.
Gravity and Magnetic Maps of Washington
Gravity and magnetic maps compiled for Cascadia (Finn, 1990;1998) have been very useful in outlining key geologic features in the region. Gravity data, in particular, have been useful in studying the distribution and thickness of the mafic crust of Siletzia. The mafic rocks appear as gravity highs (reds) in the map of Washington shown in Figure 2. The general picture obtained by viewing all of the gravity data for Cascadia is that the higher density rocks of Siletzia occupy a coherent block in the coast range of Oregon, but are separated into more discrete mafic centers in Washington. Vancouver Island appears as a large high because of the combination of Crescent Formation (or equivalent Metchosin Fm.) basalts, relatively high density Wrangellia terrane crust (Clowes et al., 1997), and an underlying mafic/ultramafic wedge. Sedimentary basins, such as the Seattle and Everett basins cause large gravity lows and the signal from these sedimentary rocks obscures the underlying Crescent Formation. A curving trend of high gravity values surrounding the Olympic Mts. outlines the rim of Crescent Formation which forms the backstop for imbrication of accretionary prism sedimentary rocks inside the rim. The rim is seen on the gravity data to turn inward on the core of the Olympic Mts. near the southwestern end (white line). This feature illustrates the strong encroachment of Siletzia blocks on the accretionary prism from the south, and is also observed in magnetotelluric (MT) models and seismic velocities discussed below.
Aeromagnetic maps for Cascadia also dramatically illustrate the distribution of highly magnetic Siletiza rocks, but are less useful in determining thicknesses and depth for the mafic rocks, because the magnetic data are very strongly affected by shallow rocks. Upward continuation and filtering of the aeromagnetic data can be used to grossly study deeper geologic features, but magnetic features deeper than about 1-2 km are obscured by the large signal caused by near-surface mafic rocks. One significant use of aeromagnetic data in Washington has been to investigate lineaments that can be interpreted in a seismotectonic manner. The aeromagnetic data used as a base in Figure 2 has been used by Finn and Stanley (1997) to point out a number of features that have seismotectonic significance. The most striking lineament in the aeromagnetic data is the long Olympic-Wallowa Lineament (OWL) defined by Raisz (1945), seen as a set of anomalies crossing the Cascades and Columbia Plateau region. Finn and Stanley (1997) have also pointed out a number of northeast-trending aeromagnetic and seismicity lineaments in western Washington that are aligned in the direction of JDFP subduction. The Seattle fault appears as an east-west set of banded anomalies along the south side of the corresponding gravity anomaly. These magnetic features extend slightly southeastward to connect with the OWL magnetic anomalies. Blakely et al. (1999) have studied a more detailed expression of the western part of the Seattle fault zone using high-resolution aeromagnetic surveys of the Puget Lowland, and the fault zone appears much as defined by marine seismic reflection profiles and previous aeromagnetic surveys.
Magnetotelluric (MT) surveys have also been very useful in sorting out the complex geology of western Washington. Stanley et al. (1984, 1987, 1992, 1994, 1996) describe MT results from the Cascades which outline tectonic features associated with hidden sedimentary rocks, intrusions, block boundaries, and seismicity zones. The locations of three regional MT profiles completed in our study of western Washington are shown in Figure 2. The MT data were modeled with the two-dimensional inversion program of Smith and Booker (1991). Because of the extreme complexity of the geology and somewhat unfavorable layout of the profiles with respect to the geologic strike, we used only transverse magnetic models (Figure 3) which minimize the three-dimensional effects in the data caused by the complexity of geology and profiles (Wannamaker et al., 1984). In addition to these three MT profiles, Aprea et al. (1998) completed an east-west MT profile through the center of the Olympic Mts. region and across the Puget Lowland into the Seattle Basin.
Our profile AA' (Figure 3a) crosses the southern and eastern rim of Crescent Formation and was designed to determine if Olympic Mts. sedimentary rocks occurred underneath the Crescent rim on the south. The western extent of this southern rim of Crescent Formation is evident just east of Quinalt Lake with resistivities of >800 ohm-m. These resistivities and similar values on profile BB' (Figure 3b) are high for Crescent Formation (Stanley et al., 1987) and probably reflect the absence of interbedded sedimentary rocks and alteration. Moderate resistivities (40-100 ohm-m) in the center part of profile probably represents the off-profile expression of core sedimentary rocks north of the profile. In the center of the Olympic Mts. region, Aprea et al. (1998) measured resistivities for these sedimentary rocks of 5-10 ohm-m so they represent a very strong resistivity contrast with the Crescent Formation. Eastward-dipping resistivities of >600 ohm-m are associated with Crescent Formation outcrops in the eastern part of the Olympic Mts. There is possible expression of a conductive region (~40 ohm-m) beneath the Hood Canal and Bremerton regions at depths of about 20 km, but this feature is very poorly constrained and probably due to 3D effects from the nearby Seattle Basin (just north and east of the east end of BB'). Aprea et al. (1998) used the same 2D modeling program that we used to model their data and also found a conductor at depths of about 20 km under Hood Canal (their conductor C4), but we attach little significance to this feature because of the very unfavorable orientation of their profile across the Puget Sound and Seattle Basin region. More robust results were derived by Aprea et al. (1998) for the Seattle basin sedimentary thickness, which is about 10 km in their model.
The model for profile BB' (Figure 3b) illustrates complexity in the sedimentary basins and underlying Crescent Formation basement. The only region of high resistivity Crescent Formation is on the north end, where the profile crosses the southern rim of Crescent Formation flanking the Olympic Mts. Sedimentary rocks in the Grays Harbor basin are up to 6 km thick and have resistivities of 1-10 ohm-m. There appears to be a thick section of underthrust sedimentary rocks underneath Crescent Formation in the Black Hills area,, possibly associated with present or older accretionary prisms.
On the profile CC' model (Figure 3c), the eastern margin of the Seattle basin is filled with 5-7 km of sedimentary rocks, in contrast to 10 km maximum thickness to the west of the profile (Aprea et al., 1998). According to Tabor (1994) the rocks in the central part profile CC', labeled 'North Cascades', are part of the western and eastern melange belt and those north of the Devils Mt. fault are in the Northwest Cascades system. Both of these rock systems consist of highly imbricated metasedimentary and metavolcanic rocks of varying grade, so it is not clear why the crust in the central part of the profile is much more resistive. However, it is possible that the high resistivities are related to a Cretaceous pluton, several of which occur in the North Cascades region (Tabor, 1994). There appears to be a strong resistivity contrast between pre-Tertiary rocks on the north and south sides of the Devils Mt. fault. North of the Devils Mt. fault, the Bellingham basin sediments of 5-7 km thickness are outlined with low resistivities. In addition, at depths of 15-20 km, there is a mid-lower crustal conductor which has also been observed in British Columbia and may be similar to a conductor in the crust of Oregon and Northern California at comparable depths. Stanley et al. (1990) interpreted that this crustal conductor was related to deep crustal fluids trapped below the brittle-ductile transition in a warming continental crust.
Crustal seismicity (<30 km depth) and regional faults are shown in Figure 4(a). Just south of the area covered in this figure, intense seismicity in the Mendocino triple junction region is related to deformation of the Gorda plate as it sheared under the North America plate by the colliding Pacific plate. The force of this collision is interpreted to be transmitted to the Juan de Fuca plate and the Coast Ranges of Oregon and Washington by Spence (1989) as a component of NS shear stress. Additional shear stresses in the region include a component of NW directed shear caused by extension in the Basin and Range province as interpreted by Zoback and Zoback (1991). Pezzopane and Weldon (1993) estimate that horizontal translations of 4-10 mm/yr are associated with this NW oriented shear.
As pointed out by Stanley et al. (1996), the tectonics of western Washington and Oregon fit the model of Beck et al. (1993) for transcurrent shear within an obliquely convergent margin where local stress directions are highly modified by block boundaries, buttresses and subduction zone bending. In Beck's model, a "sliver" of the continental plate margin becomes detached and moves parallel to the subduction zone. This margin "sliver" in Oregon and Washington is composed of the rigid, mafic Siletzia complex, which effectively transmits transverse (NS) compression. A buttress that resists northward translation of the margin sliver is formed by thick continental crust in the Vancouver Island and North Cascades region. Wang (1996) has analyzed the horizontal stresses of the Cascadia forearc sliver to illustrate the development of NS compression. In addition to this scenario of margin sliver dynamics, the Cascadia subduction zone acts as a window through which the active stresses generated by Pacific plate-North American plate interaction can be redistributed. The evidence for strong NS compression does not minimize the potential for subduction zone earthquakes, but simply indicates that this subduction zone contains a highly mobile margin.
No large earthquakes have been recorded on the JDFP thrust, but several M6-7 downdip, tensional earthquakes have occurred within the subducting plate. These intraslab earthquakes are indicated in Figure 4(b) along with other events at depths >30 km. With the exception of a M6.0 event off the SE end of Vancouver Island, all of these intraslab earthquakes have occurred in southern Puget Sound. They include the 1949 M7.1 earthquake which occurred near Olympia at a depth of 54 km, and a M6.5 earthquake which occurred near Tacoma in 1965 at a depth of 60 km. Focal mechanism solutions for the southern Puget Sound earthquakes and the one near Vancouver Island indicate normal faulting expected from slab pull forces occurring within the downdip parts of the subducting JDFP. There have been no earthquakes identified on the subduction thrust in the region offshore, and most of the small events that occur at oceanic plate depths beneath Washington have been related to trench-parallel bending stresses within the JDFP arch (Ma et al., 1996).
Benz et al. (1996) recently developed a local-earthquake, travel-time tomography method that we have applied in this study to map subsurface velocity structure in western Washington and southern Vancouver Island. Travel times between source and receivers are calculated using finite differences and simultaneous inversion for 3-D, P-velocity structure and earthquake locations is accomplished using a large matrix solver (LSQR) with smoothing constraints. Travel times are calculated using the technique of Podvin and Lecomte (1991), which estimates the wave-front (isochron) by solving the Eikonal equation across a finite-difference grid of rectangular, constant velocity cells. Standard formulations of the arrival time tomography problem require knowledge of the ray length within each cell sampled by a source-receiver pair, which is not explicitly known from finite-difference computations. Therefore, rays are found by backtracking along the perpendicular to the wave front (the steepest path) from source to receiver (Hole, 1992). Finite-difference travel-time calculations and ray-tracing are robust for the computation of direct, diffracted and refracted waves in complex velocity structures and extreme topography. S-wave phase data were used to improve earthquake relocation, assuming a constant Vp/Vs ratio of 1.77 obtained from a composite Wadati plot. The grid outlined in Figure 5 was used in the P-wave tomography. The grid consisted of 8x8x2 km cells (2 km vertical) with 35 columns (EW), 43 rows (NS), and 42 depths slices, for horizontal dimensions of 280x344 km. The top of the model grid was 4 km above sea level and the model extended to a depth of 80 km below sea level. In addition to this detailed model of western Washington, we developed a more regional model that covers much of the State of Washington. The grid for this model is shown in Figure 6. The cell size for this model was 10x10x2 km, but because of memory limitations, we had to cut the depth range off at 40 km. Results from this model proved to be useful for viewing the context of the more detailed model and examining large features.
Several selection criteria were used to provide the best phase data and coverage for the inversion. To obtain a reliable, well-resolved model it was necessary to select both well-recorded and well-distributed earthquakes with 10 or more P-wave readings and at least 2 S-wave readings (all recorded at less than 150 km separation). The addition of S-wave arrivals allowed more accurate relocation of earthquakes. All selected events were inside the station distribution (gap <180o). Additionally, to eliminate redundant data without affecting the coverage, earthquakes within a minimum distance of 1 km to other events have been eliminated, keeping the best recorded event. From 36,565 earthquakes in the PNSN catalog in the interval from 1970-1996, 3431 events were selected for the inversion process. This resulted in 61,906 arrival times recorded at the 87 seismic stations (Figure 5) in western Washington. These events also provided 16,972 S-wave arrival times that were later inverted to obtain a 3-D S-wave velocity model using the same grid as for the P-wave inversion. Resolution analyses for the cross-sections utilized in this paper are discussed in Appendix A.
Previous local earthquake tomographic studies have been done by Symons and Crosson (1998) of a smaller region than our study, with some similar results. Moran et al. (1996) have applied a similar method to a study of southwestern Washington and derive results compatible with those of our more regional study. Lees and Crosson (1990) performed a tomographic velocity study of much of western Washington with cubic cells of 5 km dimension. Their results were presented as slowness perturbations which were compared with local geology, but the subducting JDFP was not resolved in the inversion.
Velocity Model 2D Cross-sections
Selected vertical cross-sections through our three-dimensional P-wave velocity model are shown in Figures 7 and 8. The locations of the west-to-east cross-sections and south-to-north cross-sections are shown in Figure 5. The main features of the velocity cross-sections are listed below.
The subducting Juan de Fuca plate is clearly seen on all the west-east cross-sections except for row 36 (Figure 7a). We define the top of JDFP mantle by velocities of 7.5-7.75 km/s that dip down to the east. These high velocities and the gradients that encompass a narrow zone of earthquakes represent oceanic mantle rocks of ultramafic composition (peridotites). Above this well-defined, high-velocity expression of the Juan de Fuca mantle, 5-10 km of basalt flows and cummulates probably occur; thus, the actual top of the subducting plate may be at about 6.75 km/s velocities. The JDFP dips changes from about 10-12o to >35o at about -122.4o longitude, as indicated in Figure 7d and 7e, where deep earthquakes were available to map this change.
The use of S-wave phases combined with the P-wave velocity model greatly improved the relocation of earthquakes associated with the subducting plate. Most of these earthquakes are related to bending stresses in the plate, as determined by Ma et al. (1991) from focal mechanism solutions. There are no thrust events evident in the catalog for this part of the JDFP interface.
Very high velocities observed beneath western Vancouver Island (Figure 7a) have also been mapped with onshore-offshore seismic refraction profiles described by Spence et al. (1985). This high-velocity feature continues southward on the remaining cross-sections shown in Figure 7. Spence et al. (1985) interpreted that the high-velocity region has velocities of 7.7 km s-1 at depths of 16-30 km under southwestern Vancouver Island, and has an irregular bottom surface. Clowes et al. (1987) imaged the high-velocity feature with LITHOPROBE deep reflection data. The reflection image of the high-velocity zone is one of a flat-topped body beginning at depths of about 14-16 km and with a monotonically-dipping lower surface from 18 km on the southwest to about 33 km on the northeast and a total, across-strike width of 50 km. The geometry of the high-velocity zone interpreted from the refraction and reflection data matches our tomographic results (Figure 7a). Clowes et al (1987) discussed the nature of the high-velocity feature under Vancouver Island using laboratory observations of rock velocities and tectonic information. They rejected the possibility that the high-velocities were related to Crescent Formation because they could clearly see an intervening reflective region between known Crescent Formation and the deeper high-velocity feature. They also rejected the notion that the high-velocity section is composed of metamorphosed Olympic core sedimentary rocks because lab measurements on similar rocks produced velocities no higher than 6.5 km s-1, much lower than the 7.7 km s-1 required. Their two remaining alternatives for the origin of the deep, high-velocity feature both assumed mafic, and perhaps ultramafic rocks added to the base of the accretionary wedge. One scenario involved accreted mafic and ultramafic rocks added through a specific episode of oceanic slab accretion (rather than continued underthrusting), while the alternative relied on a more continuous, stacked addition of oceanic slab components. Inspection of our velocity models indicates that the high velocities occur through a wide region in the lower crust, but becoming shallower the southern Vancouver Island region and with apparent tectonic thickening in the area of the JDFP arch. Symons and Crosson (1998) referred to these high velocity rocks as a 'mantle wedge' in their tomographic velocity studies of the Puget Sound region.
The seismic refraction models of Miller et al. (1997) include velocities of 7-7.4 km/s in the lower crust. Their north-south refraction profile was located east of Puget Sound and did not cross the thickest part of the high-velocity, lower-crustal rocks shown in Figures 7 and 8 where our velocities are in the range from 7.0 to 7.75 km/s. However, Miller et al. (1997) recognized the anomalous nature of the lower crust in western Washington and also found that it was highly reflective in wide-angle recordings. Upper mantle velocities cited by Miller et. al. (1997) for western Washington, including their data, range from 7.6 to 7.85 km/s. Our tomographic results suggest a value of about 7.75-7.8 km/s east of Puget Sound. The fact that lower crustal velocities are very close to upper mantle velocities makes it difficult to estimate the most probable lithology and origin for the anomalous lower crust, including the wedge that shallows to the north under Vancouver Island. We initially believed from looking at east-west velocity profiles like that of Figure 7d and others which suggest a wedge tip, that the high-velocity units consisted of upturned continental mantle wedge, or a slice of older oceanic lithosphere as envisioned by Clowes et al (1987). However, gaps and the general complex geometry of the high-velocity section tends to disfavor the rocks as part of a mantle wedge of sheared oceanic lithosphere. We lean to the interpretation that the anomalous lower crust consists of voluminous mafic and ultramafic cumulates which formed in the large rift system that generated the Crescent Formation. Velocities in mafic rocks such as garnet granulites or garnet gabbros are 7-7.3 km/s and ultramafic rocks such as mafic eclogite, dunite, and pyroxenite are 7.4-7.9 km/s (Christensen, 1979). Babcock et al. (1992) find no surface evidence for ultramafic rocks in western Washington or Oregon to match with the voluminous basaltic rocks that they interpret were associated with a rifted margin, with the exception of some fairly rare cumulate layers in the basalts. However, the rift phase may have included a large ponding of ultramafic melt from the mantle, which then partially solidified and differentiated to generate the basalts seen in such volume at the surface. This solidified mafic/ultramafic body at the base of the crust may explain the high lower-crustal velocities in western Washington. We prefer this interpretation to that of an upturned section of mantle wedge.
Olympic Mts. Sedimentary Complex
Thick sedimentary rocks in the Olympic Mts. are imaged in detail in the west-east cross-sections (Figure 7cde). In Figure 7c, the Olympic sedimentary complex appears to be up to 30 km thick and thrust beneath higher-velocity Crescent Formation basalts for a distance of about 30-40 km. These sedimentary units probably extend all the way to the top of the JDFP. The north-south cross-sections (Figure 8) show a view of the Olympic sedimentary complex from the orthogonal direction. The tight compression of this sedimentary complex is exhibited in Figure 8a to Figure 8c as one progresses west to east. The sedimentary sections on Figure 8b and 8c near the eastern margin of the Olympic Mts. are overlain by Crescent Formation. On Figure 8c, a section of Olympic sedimentary rocks appear to be shoved down to lower crustal depths (30 km). The velocity sections show that the thickest, lowest velocity part of the Olympic sedimentary complex trends northeast along the axis of the JDFP arch. This may illustrate a critical taper effect, where shallower dip along the JDFP arch axis requires increased imbrication of the accretionary sediments.
The arch in the subducting plate is evident on the north-south cross-sections (Figure 8a-d). It has an amplitude of about 10 km, spans the region from 47°N to 48.5°N and is well imaged to at least 122.4°W (Figure 8d). We derived depth contours (Figure 4b) for the top of the subducting Juan de Fuca plate using our P-wave velocity model and relocated earthquakes. We picked the top of the oceanic mantle by well-constrained velocity horizons at 7.75 km/s. Where this velocity horizon was not so clear, but the seismicity pattern still showed the arch, we used a level near the bottom of the earthquake band which corresponds to the 7.75 km/s velocity level elsewhere.
Our contours for 30, 40, 50, and 60 km (Figure 4b) outline the arch identified by Crosson and Owens (1987) and Weaver and Baker (1988). The JDFP plate is found at depths of 50-60 km in the Puget Sound region and 15-20 km near the coastline. It has been assumed that the depth to the plate is 80-100 km beneath the volcanic arc, but seismicity is too limited to accurately determine this depth. However, teleseismic tomography by Bostock and Vandecar (1994) and Humphreys and Dueker (1994) both show the JDFP at these depths below the volcanic arc in Washington State. The most significant part of our plate contour map is the northeast trend of the arch axis, which can also be seen in velocity isosurfaces presented in a later section. This detail is a new finding from the shape of the plate contours shown by Crosson and Owens (1987) and Weaver and Baker (1988), but the teleseismic tomography of Bostock and Vandecar (1994) also show a NE trend to the deeper part of the JDFP arch.
Low velocities associated with the Seattle Basin extend to ~10 km depth near 47.6oN (Figure 8d). The Seattle fault occurs on the south margin of the basin. Low-velocities occur in the southern Washington Cascades conductor (SWCC) found with magnetotelluric (MT) surveys (Stanley et al., 1992) on the north-south cross-section of Figure 8e at about 46.9°N. This section of low velocity rocks correspond to conductive regions delineated in the MT models of Stanley et al. (1996). Mt. Rainier is located just east of the boundary of the SWCC (Stanley et al., 1996).
The terrane boundary between Eocene to Quaternary oceanic components (mainly Crescent Formation) on the west and pre-Tertiary components of the continental plate on the east is imaged on several of the west-east cross-sections (Figure 7c-e). In the mid-crust at depths of 15-25 km, the terrane boundary appears to occur at about 122.25°W to 122.35°W, based on the velocity structure and tight, vertical seismicity concentrations (circled areas near 122.3o W). These concentrations may represent dextral slip along the deep contact of Siletzia with pre-Tertiary basement on the east. Such a terrane boundary location fits a previous interpretation of the contact from magnetic and gravity data (Finn, 1990). A sharp increase in plate dip occurs at about the same longitude as the terrane boundary.
3D Visualization of Velocity Results
The cross-sections of Figure 7 and 8 are two-dimensional (2D) displays of the velocity and seismicity information from our models. However, many of the details of the velocity model can only be appreciated with a three-dimensional visualization method. The 2D displays were made with commercial display software for X-windows (output to Postscript). Public domain Generic Map Tools (GMT, Wessel and Smith, 1991) software was also used to make similar 2D cross-sections and many of the map illustrations in this report. For 3D visualization we used a different commercial software package for X-windows display. The 3D software used has only minimal transparency capability for overlaying different parameter fields, but because of its ease of use, it proved to be a valuable tool. In this section we discuss selected results from 3D visualization of the velocity model. Unfortunately, the need to keep the 3D display grid of the same dimensions as the model prevented us from plotting all earthquakes on the 3D cross-sections, but details of seismicity can best be viewed with successive 2D sections, as presented in the 'Seismicity Patterns' portion of this report.
One of the more important results from the 3D visualization is a clear image of the JDFP structure under western Washington shown in Figure 9. The plot represents a P-wave velocity isosurface taken at 7.5 km/s and viewed from the southeast. The coastline is only crudely plotted at the surface level because the sample points were limited to the grid dimensions (35x43). The perspective view causes apparent shifting of the coastline to the northwest. Some indication of the surface projection of geographic features on the isosurface can be obtained by dropping a vertical from the geographic point to the isosurface. This is true of the other plots that we present in this section, since they are all viewed with the same perspective parameters and from the same direction.
The JDFP arch is the smooth part of the isosurface which gently dips under the Puget Sound region. The black line outlines the simplest part of the JDFP arch, but the expression of the plate surface continues to the south (left), where blocks of Crescent Formation have been pushed against it. The topography of the mafic wedge outlined in the 2D sections of Figures 7 and 8 is defined in this isosurface view, and there is an obvious 'gap' between high-velocity wedge sections on the north and the south.
An overhead view (looking to the northeast) of the 7.5 km/s isosurface is shown in Figure 10. The red lines are 'streamlines' calculated by the display program and represent the direction of local maximum dip. The dip lines show that the axis of the JDFP arch dips to the northeast, almost exactly in the direction of JDFP convergence with North America (white arrow). The JDFP arch dips rather gently to the southeast, but on the north it is truncated again the Vancouver Island region near the Devils Mt. fault (DMF). The JDFP is being forced against the subduction zone buttress formed by thick, pre-Tertiary crust under Vancouver Island and the North Cascades of Washington. This buttress forms a backstop for the north side the arch and for north-south compression of mafic blocks in western Washington.
A good view of the northern buttress can be obtained from an isosurface generated with the more regional model (Figure 11). The isosurface from the regional model was taken at 7.0 km/s, so very little of the JDFP arch is delineated. Additional roughness on the mafic wedge topography in this isosurface is caused by viewing these mafic rocks at a higher level. Some Crescent Formation basement rocks are lumped with the deeper wedge units in this shallower isosurface. The dashed white line denotes the extent of the mafic wedge units. The mafic wedge units are significantly shallower under southern Vancouver Island. The smooth, raised region of the isosurface in the area just north and east of the DMF is the more uniform-velocity, pre-Tertiary crust of the Vancouver Island-North Cascades buttress region. Note that the DMF and DAR closely follow this boundary, suggesting that these thrust faults are controlled by the southern boundary of the buttressing crust (at least in the Cascades). In the maps of this report which show the DMF and DAR, we follow the positions used by Tabor (1994) because these seem to match the deep expression of these faults in the velocity models. However, Johnson et al. (1999b) has found evidence from marine seismic data that the location of the DMF in the near-surface is 10 km or more south of the deep location at depths of 20-30 km. This is to be expected since the DMF is a north-dipping thrust. The boundary of the buttress region extends in a linear fashion southeast along the OWL for some distance. This implies that the shear zone associated with the OWL is located along the southern boundary of the pre-Tertiary buttress. Other pre-Tertiary sections occur outside the buttress region, but the crust apparently is not the coherent rigid anchor block formed by pre-Tertiary units north of the OWL and the DMF, and east of the DAR. The velocity expression of this boundary in the Cascades region also coincides with a NW-trending heat flow low to be discussed in a later section.
The tilted, perspective view of the 7.5 km/s isosurface (Figure 9) forms a convenient base for display of perspective views of velocity cross-sections that illustrate key tectonic features of western Washington. The first 3D cross-sectional view illustrates the key sedimentary basins in western Washington (Figure 12). This display consists of the set of sections for velocity model columns 16 through 30 and outlines low-velocity sedimentary rocks. As indicated in the 2D view (Figure 8d), low-velocity rocks in the Seattle basin extend to about 10 km depth. The use of multiple cross-sections in the 3D view of Figure 12 shows an offset in basin sediments that occurs along a dextral slip zone in Puget Sound, which has moved the western margin of the Seattle basin north of the eastern part (Johnson et al. 1999a). The offset in depth to the base of the Quaternary section interpreted by Johnson from industry seismic reflection data is shown in Figure 13. The offset is quite similar to that imaged in the tomographic velocities, even though the latter sample a much deeper section (4-10 km depths).
The Tacoma and Everett basins are also evident in the 3D display of Figure 12, but the largest low-velocity feature is associated with the SWCC. Low-velocities occupy the region of the conductive crustal rocks mapped with MT surveys (white lined area), except at the location of two large intrusions. The larger of the intrusions occurs between Mt. St. Helens and Mt. Adams and the other is just south of Mt. Rainier. The velocity model shows a velocity contrast of >1 km /sec in the rocks within the SWCC conductive anomaly (apart from the intrusive bodies), primarily in the upper 5-7 km. In our models, velocities range from 4.75 to 6.0 km/s for the SWCC. These are only slightly higher than the velocities of the Eocene to Quaternary sedimentary rocks of the Seattle basin (4.5-5.75 km/s). This velocity range is compatible with Eocene and possible Cretaceous marine sedimentary rocks as favored by Stanley et al. (1992, 1994), especially in the southern Cascades region where known Eocene sedimentary rocks have been moderately metamorphosed by high heat flows in the Cascades and contain extensive coal beds and interlayered volcanic flows (Stanley, et al., 1994). The range of velocities from our tomographic model is also for compatible with sedimentary velocities documented by Stanley et al. (1994) in the SWCC based on well-logs, seismic-reflection stacking velocities, and compaction curves. We expect that deeper in the conductivity anomaly (which Stanley et al., 1992 interpret extends to >10 km depth), velocity contrasts with surrounding rocks are small.
Previous seismic refraction surveys revealed no details of velocity structure for the SWCC. Miller et al. (1997) modeled a long north-south refraction profile completed just west of the Cascade volcanic belt in western Washington which crossed through the long axis of the SWCC. Their ray-trace model included a upper crustal layer with velocity gradients from 4.0 to 5.5 km/s and a mid-crustal layer with velocity gradients from 6.1 to 7.1 km/s. No details of the SWCC were evident in the model. One factor that may have prevented recognition of velocity variations with the SWCC may be the high-velocity volcanic flows (>6.0 km/s) in the SWCC region documented by Stanley et al. (1994). These flows cause a velocity reversal to the lower-velocity sedimentary rocks beneath them, creating a blind zone which negates the ability of refraction methods to determine velocities in the sedimentary section. Our tomographic models were not affected by the high-velocities in the same way, although a slight bias may be present.
The most extensive active seismic experiment through the region has been described by Parsons et al. (1998). This experiment extended from the Washington coastline into the Columbia Plateau and involved 1500 seismometers and 17 large explosive sources. The data were interpreted through both standard refraction and wide-angle reflection techniques. Although Parsons et al. (1998) found evidence for a velocity contrast at the SHZ, no details were forthcoming on the velocity structure of the SWCC; this failure stems from the same reasons noted for the north-south survey of Miller et al. (1997).
In a local earthquake velocity model derived for the greater Mt. Rainier region, Moran (1997) found both high and low velocities associated with the SWCC. These variations were due to the effects of the surface volcanic flows and intrusive bodies shown in Figure 12, which Moran did not recognize. However, part of the problem in recognizing the velocity structure of the SWCC is visualization of its 3D signature. The nature of the SWCC and features such as the offset in the Seattle basin are relatively clear when visualized in the manner of Figure 12, but when only horizontal depth slices are used (as by Moran, 1997), the detail sometimes overwhelms these features. It is obvious to us that the southern Washington Cascades has anomalously conductive and low-velocity upper crust, especially in the upper 5-7 km. There is a almost a direct correlation between the region of high electrical conductivities in the SWCC and seismic P-wave velocity at these depths, but the correlation is more difficult to find at greater depths. This difficulty may be due to the effects of interlayered volcanic flows and coal beds in the sedimentary rocks and compaction/metamorphism which increase the velocities significantly. Neither of these effects decrease the electrical conductivity, and the latter can increase conductivities (Raab et al., 1998) by producing connected zones of graphite and sulfides. Also, electrical conductivities in sedimentary rocks vary from 10 to 0.001 mho/m (four orders of magnitude), whereas P-wave velocities vary only from 3.0 to 6.5 km/s (much less than one order of magnitude). Thus, it is not surprising that the electrical signature of a thick section of Tertiary marine and nonmarine sedimentary rocks with interbedded volcanic flows and coal beds appears much more intense than the acoustic signature. Deep Vibroseis reflection data shown by Stanley et al. (1992, 1994) show the character of the deep sedimentary section (down to 8-10 km) better than either the refraction models or our tomography.
General Results From 3D Visualization of Cross-sections
A suite of east-west 3D cross-sections is presented in Figures 13-25. The north-south cross-section at column 8 is retained on all the images. The column 8 section extends across the southern end of Vancouver Island. In the web pages for these image, a clicker is provided on the left side to navigate through the images by going forward or backwards; alternately a map selection utility is also provided. The high velocity feature beneath Vancouver Island (Figure 14) is the shallow, northern terminus of the mafic wedge feature that is concentrated beneath the Puget Sound region. On the row 14 east-west section (Figure 14), SWCC low-velocities dip to the east beneath Mt. Rainier, as shown also in MT profile interpretations by Stanley et al. (1996). The low velocity rocks are intruded by higher velocity Tertiary and Quaternary intrusive bodies. A better view of the intrusion near Mt. Rainier was developed from the more detailed tomographic model of the Mt. Rainier region (Figure 15). The intrusion in Figure 15 seems to originate in a northwest trending dike system at depth. The row 16 section (Figure 16) outlines remnants of the SWCC north of Mt. Rainier, but the most significant feature is the complicated deep structure on the JDFP near the southern margin of the most massive part of the mafic wedge. The mechanical resistance provided by the mafic wedge may cause the JDFP to take a locally steep dip.
Figure 17 (row 18) outlines part of the Tacoma basin as low velocities near Olympia. Accretionary prism sedimentary rocks in the southern part of the Olympic Peninsula where the profile transect crosses have thicknesses of about 10-15 km and are mildly thrust downward to the east. These sedimentary rocks have maximum thicknesses inside the rim of Crescent Formation and over the axis of the JDFP arch. The row 20 section (Figure 18) outlines the southwestern extent of this rim (labeled Crescent Fm.) noted in gravity data (Figure 2) and MT models (Figure 3a). Row 22 (Figure 19) shows continued thickening and eastward underthrusting of the prism sedimentary rocks. Details of the backstop formed by the Crescent Formation rim show how the complexities of this backstop probably affect the imbrication of prism sedimentary rocks.. Note that the Crescent Formation is buttressed by the mafic wedge, otherwise there would be no resistance to eastward spreading of the prism sedimentary rocks and the thin Crescent Formation rim. Row 23 (Figure 20) extends through the southern part of the Seattle basin and row 24 (Figure 21) crosses the maximum thickness of the basin sedimentary rocks (about 10 km). Row 25 (Figure 22) continues to image the Seattle basin. This latter cross-section also illustrates how the accretionary prism sedimentary rocks thicken over the axial part of the JDFP arch. Row 26 (Figure 23) indicates an important tectonic behavior in which increased eastward thrusting of prism rocks is accompanied by eastward translation of Crescent Formation backstop through the 'gap' in the mafic wedge. Since buttressing of the Crescent Formation does not occur in this gap, there is increased underthrusting and eastward translation. This is even more evident on row 27 (Figure 24) where lower-velocity, younger prism sedimentary units near the coast fill in behind this increased eastward transport. Row 28 (Figure 25) shows continuation of this movement of accretionary prism and Crescent Fm. units through the mafic wedge gap. Rows 29 (Figure 26) and 30 (Figure 27) illustrate the effects of return of mafic wedge buttressing on the Crescent Formation and accretionary prism complex. Rows 32 (Figure 28), 33 (Figure 29), and 34 (Figure 30) show the emergence of the high-velocity rocks of the mafic wedge and Crescent Formation to their shallowest depths beneath southern Vancouver Island.
For viewing velocity information in north-south cross-sections, we step from west to east starting with column 10 (Figure 31). Figure 31 shows the JDFP arch, as well as the start of compression and eastward thickening of accretionary prism rocks. The bumpiness on the southwestern flank of the JDFP arch is related to Crescent Formation buildup or discrete blocks. This Crescent buildup is better illustrated on column 12 (Figure 32), where it connects to shallow Crescent Formation mapped in MT profiles (Figure 3a) and gravity data (Figure 2) well west of outcrops of Crescent. This more westerly extent of Crescent Formation is significant because it turns inward on the core of the Olympic Mts. (Figure 2), illustrating the dramatic nature of north-south compression. This geometry of the Crescent rim partially explains tectonic features mapped by McCrory (1997) north of Grays Harbor and just south of the hidden Crescent buildup. These tectonic features include nine ridges which show a progressive clockwise shift in orientation from south to north. According to McCrory, 'the southernmost ridges trend east-west and the northern ones trend NW-SE. Farther north, Quaternary thrust fault and fold trends shift to an even more northerly orientation (N10oW) (McCrory, 1996). The progressive shift in structural orientation suggests that the study area serves as a transitional zone between north-south-directed compression to the south and east-west-directed compression to the north. This shift occurs adjacent to a major boundary between crustal blocks of widely differing rheologies.' We agree and suggest that the pattern occurs precisely around the southern margin of the hidden Crescent Formation buildup.
Continuing with column 14 (Figure 33), the reverse nature of the Crescent-prism contact on the north side of the Olympic Mts. is illustrated. This section crosses east of the shallowest part of the hidden Crescent on column 12, so the contact is not so dramatic on the south side. The column 16 (Figure 34) section crosses the Peninsula just west of Hood Canal and the thick prism sedimentary complex is now significantly narrowed and located just over the axis of the JDFP arch. On column 18 (Figure 35) extreme compression of the prism is occurring, with Crescent Formation thrust over the top of the sedimentary rocks. The tightly compressed prism units are being forced into the mafic wedge gap on column 20 (Figure 36) and 21 (Figure 37) and almost closed off by surrounding Crescent Formation. The Seattle basin has begun to appear over the gap in the mafic wedge. The downwarping which formed the Seattle basin probably was assisted by a better facility for compression of the deep crustal units filling the wedge gap. The thickest part of the Seattle basin sediments occur over the downdip, eastern part of the wedge gap, as shown in the sections for columns 22-25 (Figures 38-41). Details of the Everett and Bellingham basins are shown on column 24. The SWCC and interior intrusions are imaged on columns 23-26 (Figs 39-42).
We have also inverted for S-wave velocities to study their utility in seismotectonic analysis of western Washington. The starting S-wave model was derived from the starting P-wave model using a Vp /Vs of 1.77. Some raypath coverage was lost for the S-wave model, especially in the southern and eastern parts of the model due to a lack of high quality S-wave arrival times. A cross-section of S-wave model velocities from row 23 (Figure 5) is shown in Figure 43 compared to a P-wave velocity section. Differences in the S-wave and P-wave models are most conveniently viewed as variations in Vp /Vs ratios (Figure 43cd).
As summarized by Christensen (1996) for deep crustal rocks, Vp /Vs and Poisson's ratio variations are primarily related to SiO2 content and alteration state. In general, felsic rocks in the deep crust have Vp /Vs ratios closer to the crustal average of 1.77 which we assumed in our study. However, mafic rocks generally have Vp /Vs ratios from about 1.81 to 1.85. One of the more variable Vp /Vs ratios comes from metamorphosed basalts (Christensen, 1996). Zeolite facies basalts have very high Vp /Vs ratios of 1.861 at large pressures (1 GPa or 30 km) and greenschist facies basalt have much lower Vp /Vs ratios of 1.766 at the same pressures. Such a large range may be important in Vp /Vs data analysis for western Washington, where basalts dominate the crustal composition.
Laboratory measurements by Kern (1982) at 6 kbar (20 km) confining pressure and 20o C indicate Vp /Vs ratios of about 1.8 for basalt, 1.69-1.73 for ultramafic rocks, 1.67-1.7 for granite, and 1.48-1.52 for quartzite samples. Kern (1982) found a nearly linear temperature coefficient in Vp and Vs for the plutonic and metamorphic rock samples he measured up to a temperature of about 500o C and 6 kbar confining pressure. The results for basalts were more complex, exhibiting very small changes in velocity to around 350o C where an abrupt, complicated decrease in both Vp and Vs occurred. Kern (1982) interpreted that this discontinuous velocity decrease is caused by dehydration reactions of zeolites and chlorites in the basalt samples. The dehydration reactions produce solid-fluid systems that gave rise to formation of microcracks and widening of old cracks in the basalt samples. Kern's basalt samples showed that the discontinuous velocity decrease also represented a decrease in Vp /Vs of 3-5%. The crack systems appear to affect Vp more than Vs in the measured samples.
Eberhart-Phillips et al. (1989) found from sandstone samples that clay content had a large effect on Vp /Vs at room temperatures and pressures of 0.4 kbar. Sandstones with no clay content and porosities of 6-8 percent had Vp /Vs ratios of about 1.56, whereas other samples with high clay content had Vp /Vs values of up to 1.8. This result can be explained by dramatic decreases of Vs in clay-rich rocks. Thus, expected Vp /Vs ratios in western Washington could vary over a relatively large range with the occurrence of large percentages of basaltic crust and thick sedimentary rocks with unknown percentages of clay.
The Vp /Vs cross-sections (Figure 43cd) show that Olympic core sedimentary rocks and Seattle basin sediments have low Vp /Vs values of 1.5-1.7. In the upper part of the Seattle basin and in other Quaternary sedimentary strata in the depth range 0-1 km, we expect that Vp /Vs ratios will be much higher due to low S-wave velocities in water-rich, porous sediments. However, these depths are above the coverage of our velocity models. The subducting Juan de Fuca plate crust and mantle are both indicated as having Vp /Vs values of 1.77-1.82, acceptable values for basaltic and gabbroic rocks. However, we caution that the deeper parts of the S-wave model has much lower resolution than the P-wave model. We have ruled out much significance to the Vp /Vs results below depths of about 20-25 km. The shallow part of the Crescent Formation along the east side of the Olympic Peninsula (Figure 43b) and in the Black Hills area (Figure 43cd) have the highest Vp /Vs ratios. These sections of near-surface basalt show up more clearly on the Vp /Vs data than in either the P-wave or S-wave velocity sections. This suggests that the basalt is relatively unfractured in these areas and possibly of low-metamorphic grade as required by Christensen's (1996) and Kern's (1982) lab measurements. Glassley (1974) has identified two separate volcanic members in the Crescent Formation found in the eastern Olympic Mts. The two volcanic members were interpreted by Glassley to have unrelated magma suites. The eastern, structurally lower member has compositions like those in inter-arc-basin tholeites. The upper, more effusive, member has compositions very similar to Columbia River Group basalts and Babcock (1992) places both members in a tectonic model which involves forearc rifting. However, it is clear that the section of Crescent Formation (most likely the lower member) interpreted from the Vp model under the Puget Sound region has highly variable Vp /Vs values in the depths near its base where most earthquakes occur. This high variability may be an indication of a the complex fracturing and possibly fluids which generates the denser seismicity shown on the cross-section; however, more detailed models of this region and a thorough evaluation of P and S-wave model uncertainties in this area will be required to more fully investigate this possibility. At this point we present the S-wave results primarily for completeness and to illustrate the type of S-wave information that can be obtained for western Washington.
Interpretive Geological Cross-sections
East-west Transect Interpretation
A geologic interpretation of velocity structures through the center of the Olympic Mountains-Puget Sound region is shown in the east-west section of Figure 44(a). This geologic interpretation was primarily developed using the east-west velocity cross-section of Figure 7c (row 23), although guided by other nearby sections. The velocity models provide no information about the offshore region, so our geological sketch has been augmented using marine deep reflection and refraction data from the Flueh et al. (1998) profile, SO12, acquired offshore from the Olympic Mountains. In particular, the section of Crescent Formation offshore (Figure 44a) was interpreted by Flueh et al. (1998) and likely connects to the southwestern end of Crescent rim which we have identified from our velocity models, gravity data, and MT models. Subduction zone dip offshore is assumed to be about 5-7o and about 10-12o between the coastline and longitude -122.4o, as indicated in the velocity model. East of -122.4o longitude, the plate bends to a dip of about 30-40o. This change in dip is probably caused by the effects of phase changes in the deep part of downgoing slab as previously inferred by Stanley (1984). Slab pull forces from the higher density, lower part of the downgoing slab are exhibited in the large downdip tensional earthquakes of 1949 and 1965 which occurred within slab, just below the change in dip (Figure 44a). Spence (1987) has pointed out that in many of the world's subduction zones, slab pull forces cause downdip tensional earthquakes just downdip from such fundamental changes in plate dip. Typically the region just updip from the change in dip is locked during initiation of these tensional events. If the subduction thrust were not locked, then slab pull forces would propagate further updip and tensional earthquakes would occur nearer the trench. We interpret that the location of the relatively large intraslab, tensional earthquakes is strong circumstantial evidence that the subduction thrust is locked in the region just west of -122.4o longitude, where mafic rocks are in contact along the thrust.
The T-axes for both the 1949 and 1965 earthquakes are downdip, as is the other intraslab event east of the southeast end of Vancouver Island (Figure 45) near Pender Island (Rogers, 1983). However, as derived by Baker and Langston (1987) the focal mechanism for the 1949 earthquake is complicated. Baker and Langston found distinct pulses in the waveforms which they assumed to be successive source events along a single fault. Under these assumptions, source directivity models indicated that a east-west (±15o) oriented fault dipping at 45o (±15o) with left-lateral slip was the source of the earthquake. Since the JDFP dips at 35-40o, the fault plane within the plate would have moved as a left-lateral, reverse fault. However, they noted that this determination is dependent on one reference station (WES) and that a combination of radiation effects from rupture complexities could mimic the source directivity assumed. We prefer the alternate choice of fault planes, or the plane directed N14oE with 75o dip and -140o rake. The T-axis is still downdip and this geometry for the fault plane is closer to the trend of constant depth contours (Figure 4b) of the JDFP. The dip and rake of this fault plane is similar to that expected for a slab-pull generated, tensional earthquake.
Kirby (et al., 1996, 1998) has hypothesized that the gabbro-eclogite phase transition plays the main role in producing large tensional intraslab earthquakes such as those in 1949 and 1965 under Puget Sound. The phase transition reduces volume by 13% and can produce shrinkage stresses (Kirby et al., 1996) that produce downdip tension. He argues (oral comm., 1999) that plate driving forces alone do not produce enough stress to cause intraslab earthquakes. Wang et al. (1997) has analyzed the force balances of ridge push, slab pull, and the forces that resist these plate driving forces. Ridge push is resisted by basal drag in the asthenosphere and slab pull is resisted by similar viscous drag, but at steeper angles of attack. These forces are delicately balanced near the bend in the subducting slab and can be highly concentrated in bending moments. The forces can produce tensional earthquakes only in zones of weakness, since the strength of the oceanic slab is probably on the order of 200-600 MPa (Denlinger, 1992) under moderate strain rates. There is considerable evidence that normal-faulting, intermediate-depth, intraslab earthquakes are routinely produced above the depths where the phase transition occurs (Isacks and Molner, 1971). For instance, the Makran M6.5 tensional earthquake of April 18, 1983 may be a close analogue to the Puget Sound events in 1949 and 1965. The Makran quake occurred at a depth of about 65 km just downdip from the shallow-to-steep bend in the subducting plate and had downdip T-axes. The subducting plate age is 65-95 Ma in the Makran region; thus, the gabbro-eclogite phase transition occurs at depths far below the 1983 tensional earthquake (Hacker, 1996). Bending stresses near the shallow-to-steep bend were probably involved in faulting of the oceanic slab, but the driving force for the bend and for the stress drops were related to slab pull forces, which represent the downdip integration of all densification effects in the slab. Another classical example of a large, shallow normal-faulting earthquake in a subducting slab is the 1977 M7.9 Sumba event, which Spence (1986) has related to slab pull forces. This earthquake ruptured the upper 28 km of the slab lithosphere along a strike-length of 200 km, as defined by numerous aftershocks. Aftershocks were not recorded after the 1949 and 1965 Puget Sound intraslab earthquakes, but a larger event spanning a longer strike length might produce a number of aftershocks.
Subduction Thrust Frictional States
It is difficult to quantitatively determine the frictional state of the subduction thrust. Most recent analyses of the Cascadia subduction zone have been presented in papers such as those of Davis et al. (1990) and Hyndman and Wang (1993, 1995). These authors assumed that any seismic locked zone must be located offshore where the subduction thrust is less than 20 km deep. Their analysis was based upon transition to ductile behavior of metamorphosed sedimentary rocks at temperatures between 350o and 450o C.
Tichelaar and Ruff (1993) surveyed world subduction zones for temperature characteristics of the maximum depth of interplate seismic coupling. Using heat flow measurements as constraints, they modeled the frictional shear stress and, in turn, calculated the temperature at the down-dip edge of the seismogenic zone. In their analysis of world subduction zones, they found the global average frictional shear stress to be:
(a) 14±5 MPa, assuming stress is constant downdip along the seismogenic zone;
(b) 30±5 MPa, assuming a constant friction coefficient, with stress increasing linearly through the downdip seismogenic zone.
Under assumption (a) of constant stress, they found that a single critical temperature of about 250o C fit the subduction zone data. However, for assumption (b) of constant friction coefficient, they found a bimodal distribution of seismogenic edge temperatures with peaks at 400oC and 550o C. They inferred that a temperature of about 400oC might be characteristic of the maximum depth of subduction thrust earthquakes when the upper plate was composed of silicic rocks and about 550o C if the upper plate contained mafic rocks.
In the Ryukyu-Kyushu arc and Japan trench, Kao and Chen (1991) found that large to moderate-sized earthquakes initiated in the depth range from 20-30 km occurred at interplate thrust temperatures calculated to be 300-600o C, higher than those assumed to limit thrust earthquakes by Hyndman and Wang (1993,1995). Kao and Chen (1991) interpret that the accretionary wedge is too weak to accumulate seismogenic strain and generally view the seismogenic region as occurring relatively deeper than do Hyndman and Wang (1993,1995).
Regarding the assumed thermal regime along the subduction thrust, in Figure 44a we show assumptions from Hyndman and Wang (1993) concerning the location of the 350 and 450o C isotherms as they intersect the subduction thrust, but add error bars to indicate minimum uncertainties in the actual location of these isotherms at the interface. We agree with Tichelaar and Ruff (1993) who cautioned that the thermal regime of the subduction interface cannot be determined with the same confidence as in other tectonic settings. Major sources of uncertainty where the plate interface is shallow offshore include the thermal effects of channelized upflow of fluids (depicted in arrows on Figure 44a) along the decollement and upward along faults (Davis et al., 1990; Shi and Wang, 1994). Under the onshore region, where the interface is deeper, thermal advection from updip flows of fluids due to metamorphism in the oceanic plate and thrust zone can shift the isotherm depths at a rate of 0.5-1 km/Ma (Peacock, 1987). Advection from metamorphism and sediment dewatering causes isotherms updip from this process to be shifted to shallower levels, but due to cooling at larger depths, lower isotherms will be shifted downdip. Frictional heating and endothermic/exothermic mineral reactions are other factors which are poorly known in the Cascadia subduction zone. In addition, thermal conductivities and fluid processes within a complex continental plate must be considered as sources of uncertainty in the vertical positioning of the critical isotherms. In a review of all these factors, we infer that ±2.5 km errors in the position of the isotherms portrayed by Hyndman and Wang (1993) are minimum uncertainties. For a subduction dip of 10o, a 5 km uncertainty translates to a difference of 29 km downdip in where the critical isotherms intersect the interplate thrust.
To support their view of the offshore locked part of the subduction thrust Hyndman and Wang (1993, 1995) modeled vertical uplift (from leveling and geodetic measurements) and horizontal velocity data (from GPS). They fit these two data sets using dislocation models for offshore Vancouver Island with a locked zone of 60 km downdip width and a region of transition from stick-slip to free-slip also of 60 km downdip width. However, as pointed out by McCaffrey (1997), the vertical uplift data can just as easily be modeled with a model having underplating near the coastline that represents only 10-20% of the downdip width assumed by Hyndman and Wang (1995) to be locked. Nevertheless, horizontal velocities of about 7-10 mm/yr observed on southern Vancouver Island and near Neah Bay on the Olympic Peninsula (Dragert and Hyndman, 1995) require models with a significant region of offshore thrust zone coupling, and possibly extending onshore (Flück et al., 1997). Recent GPS results from Goldfinger et al. (1999) have been modeled with a wider locked zone in central Oregon than assumed by Hyndman and Wang (1995). Goldfinger et al. (1999) find that their GPS results require a locked region extending beneath the Oregon Coast Range down to depths of 40 km or more. This is in contrast to earlier assumptions that there is poor coupling in the Oregon part of the megathrust. Their data also show margin-parallel motion diagnostic of northward translation of a margin sliver. A significant difference between the Goldfinger et al. (1999) and Hyndman and Wang (1995) dislocation models was the incorporation of a freely-slipping zone from the trench to 50 km landward in the Golderfinger et al. models, as we also interpret in Figure 44a. Qamar and Khazaradze (1999) find from their dislocation modeling of Washington GPS data that the combined locked and transition zone must be 200 km wide (downdip), extending well past the coastline. This result is identical to that found by Hyndman and Wang (1995) for a Washington transect across the subduction zone. Our anelastic strain models discussed in a later section are compatible with having the subduction thrust locked in the whole region where the subduction thrust dips at 10-12o (from the coast to Puget Sound). Thus, current models for subduction zone locking vary significantly between investigators.
Keeping in mind uncertainties in the position of the key isotherms, we show our interpretation of frictional states for the upper part of the subduction thrust. The thrust is assumed to be freely slipping from the trench to approximately 50 km landward, due to the effects of high-pore pressure fluids in the toe of the accretionary prism and to the presence of layered clays such as illites and smectites, which lubricate the thrust zone (Vrolijk, 1990). Layered clays will have been metamorphosed into a more shear-resistant species at temperatures of >200o C. The region between ~200 and 350o C can provide significant frictional, or stick-slip forces along the thrust zone, unless extensive high-pore pressure fluids exist at the decollement. This latter question has not been resolved for Cascadia, although evidence from many world subduction zones with thick, young accretionary prisms suggest that such high pressure fluid regions are common. If this possibility is correct, then the downdip extent of the locked region would be considerably reduced and pushed downdip due to depression of critical isotherms from advection along the decollement. The relationship of deep decollement fluid flow to frictional states along the subduction thrust is poorly understood, but may fit the model of seismic pumping developed by Sibson et al. (1988) for high-angle thrust faults. Without adequate constraints, however, our model assumes that there is stick-slip behavior along a significant part of the thrust down to about 20 km depth.
Downdip from the offshore stick-slip region, at depths >20 km, temperatures for silicic minerals in the thrust are above the brittle-ductile transition. In this region there may be a transition from stick-slip to free-slip, similar to that portrayed by Hyndman and Wang (1995). However, in the velocity results from our modeling and from the refraction models of Flueh et al. (1998), there is a significant wedge of Crescent Formation on top of the JDFP that extends from offshore to onshore as depicted in Figure 44a. This mafic section in contact with the JDFP would mean that the thrust zone was not above the ductile transition temperature, unless there were intervening subducted sediments. Flueh et al. (1998) show that all sedimentary rocks in the accretionary prism are scraped off at the outboard contact with the Crescent Formation units near the coastline. In addition, Flueh et al. (1998) find that the dip in the JDFP offshore is significantly shallower than previous assumed by Hyndman and Wang (1995); thus, a larger part of plate boundary lies within the seismogenic zone as outlined by Hyndman and Wang (1993, 1995). These factors place large uncertainties in the actual downdip extent of the seismogenic zone postulated by Hyndman and Wang (1993, 1995). The region of the thrust where Olympic core rocks rest on the downgoing JDFP may be freely slipping or transitional, since the silicic core rocks would be above the brittle-ductile transition, but where the mafic wedge under Puget Sound occurs in the upper plate of the thrust, temperatures are not high enough to make this contact ductile. Thus, the deepest part of the thrust just updip from the slab bend may be locked due to the properties of this contact. Below the deep locked zone, a transition to free-slip must occur rapidly near the change in slab dip, as densification of the subducting slab causes it to be negatively buoyant and pull away from the upper plate.
Temperatures in the region of the thrust overlain by the mafic wedge are below the ductility transition temperature for mafic rocks, so stick-slip behavior would be expected unless significant fluids or serpentinites lubricate the thrust zone. Serpentinite minerals are stable to at least 550o C and may exist in any ultramafic portions of the subducting oceanic crust. In addition, dehydration of the oceanic crust may produce 1-2% of water that is free to move across the subduction thrust and cause retrograde metamorphism of peridotites in the mafic/ultramafic wedge. Both the subducting oceanic basalts and overlying mafic wedge rocks would be in the greenschist to upper amphibolite facies, with the main hydrous minerals being amphibole, epidote, and chlorite (Tatsumi, 1989). Water released from the oceanic crust probably is rapidly taken up by the upper plate mafic/ultramafic rocks and a buffering, or negative feedback, process involving the local P-T path may take place to maintain the greatest stability. When fluid is released due to prograde metamorphism of the subducting crust, temperature near the thrust is reduced. Conversely, when this fluid is absorbed in retrograde reactions in the upper plate, exothermic heat production occurs to raise the temperature. Too much water release will cool the upper part of the oceanic plate and retard prograde metamorphism. The subduction velocity also enters into the buffering process, coupled with reaction rates. The dewatering of the lower plate and formation of serpentines and other retrograde minerals in the upper plate is probably a very inhomogeneous process along the thrust region. Although the role of dewatering of the slab is often discussed, geodetic and thermal constraints could be used to infer that the total effect on the upper plate is small. Water entering the mafic/ultramafic upper plate in the Puget Sound region would cause retrograde metamorphic reactions to produce serpentine at temperatures below 600o C. This reaction produces heat and a volume increase of 16% in the ultramafic rock of the upper plate. Extensive fluid flow into the upper plate would raise heat flow and cause positive elevation changes, neither of which are observed in the Puget Sound region (Figure 44a). In addition, if the mafic wedge were completely hydrated, seismic velocites would be much lower than we measure, as indicated in recent petrophysical/metamorphic models developed for Cascadia by Bradley Hacker, U.C. Santa Barbara (writ. comm., 1999).
It might be assumed that water-bearing greenschist/upper amphibolite minerals, primarily serpentinites, could act to lubricate the thrust and prevent large shear stresses. Moore et al. (1986) have done extensive laboratory studies on the behavoir of serpentinite gouges under varying sliding velocites, temperatures, and pressures. They found, that above 400o C and pressures to 2.5 kbar, the serpentinite behaved in stick-slip fashion and supported higher stresses at 600o C. Sliding velocities of 4.8 um/s and 4.8 x 10-2 um/s were used and the lower velocity supported higher stress drops. Both of these velocites are much higher than Cascadia subduction rates, therefore we expect that serpentinites in the hypothesized deep locked zone (depths 30-40 km, temperatures 450-600o C) will behave in stick-slip fashion. In addition, Moore et al. (1997) measured coefficient of friction (COF) in water-rich gouges made from the serpentinite minerals chrysotile, lizardite, and antigorite at temperatures up to 200o C. Chrysotile formed the weakest gouge with a COF of about 0.2. Lizardite and antigorite gouges had COFs ranging from 0.4 to 0.52. All three minerals were stronger at higher temperatures because free water between mineral layers had been driven off (Diane Moore, oral comm., 1999). Talc has also been cited as a mineral which reduces friction along the subduction thrust (Hyndman et al., 1997) and could occur within the stability field of the JDFP at depths of 30-40 km (Bose and Navrotsky, 1998). However, laboratory data from Summers and Byerlee (1977) show that talc under high confining pressures (to 6.98 Kbar) has a COF nearly as high as antigorite and lizardite (about 0.4). Talc will not generally be present unless there is substantial free silica present, as in subducted sedimentary fluids (Bose and Navrotsky, 1998). Talc reacts with forsterite to form antigorite and enstatite (Bose and Navrotsky, 1998) at depths of the mafic/ultramafic wedge; this reaction may reduce the amount of talc above the thrust if there is abundant forsterite in the ultramafic rocks of the mafic wedge and in the oceanic lithosphere. In addition, at temperatures of 400o and 600o C, Moore et al. (1986) found that serpentine minerals under pressure and variable sliding rates converted to forsteritic olivine, with minor evidence for talc.
If high-pore pressure fluids do not exist and serpentinite minerals have relatively high COF's at large depths and temperatures, then another mechanism might make the deep locked zone implied in Figure 44a a freely slipping zone. This mechanism involves ductility induced during slow straining of ultramafic rocks undergoing fluid-driven metamorphism. Although prograde reactions strengthen the rocks in the downgoing slab, retrograde reactions soften the upper plate ultramafic rocks. As suggested by David Lockner, USGS (oral comm., 1999) laboratory data show that such rocks may deform somewhat ductility under small strain rates during retrograde reactions. There is some evidence that the front of the mantle wedge undergoes such flow deformation, or dragging, downward along the top of the subducting slab (Tatsumi, 1989). The toe of the mantle wedge must be slowly circulating to maintain high temperatures in the arc region (Peacock, 1993). Such slow creep may be a key anelastic process that is involved in north-south compression in the coast ranges. For the deep part of the subduction thrust, it may mean that the implied deep locked zone is creeping due to ductility that is related not to brittle transition temperatures, but to retrograde softening under small strain rates. Thus, depending upon the water supply from the downgoing slab and strain rates, the implied deep locked zone could act somewhat elastically and be seismogenic, or be slowly creeping, analogous to portions of the San Andreas crustal fault zone. The obvious paucity of small earthquakes within the mafic/ultramafic wedge (Figures 7, 44a) could be interpreted to mean that this fluid-driven ductility extends somewhat uniformly thoughout the wedge. However, the lack of small earthquakes in the wedge may also mean that it is a very strong body, and not subject to the same bending moment strains as the underlying JDFP and the overlying Crescent Formation blocks. Without new large earthquakes, the constraints from geodetic data provide the best hope that we have to decipher the elastic/anelastic behavoir of the deep part of the subduction thrust under western Washington (and Oregon).
Dynamics of formation and maintenance of the JDFP arch may support the view that the deep part of the thrust is coupled. The formation of the plate arch has been interpreted by Brandon and Calderwood (1990) to have been caused by rotation of continental lithosphere about a pivot point in the Puget Sound area due to Basin and Range extension. If they are correct and the plate arch was caused by forcing the subducting JDFP into a curving continental margin of the North American plate, then there must have been strong coupling forces between the two plates during the time of development of the arch. Cascadia convergence velocities are almost as large now as during the time of formation of the arch (middle Miocene) and there is ample evidence for continued Basin and Range extension and north-south compression in the Coast Ranges. If Brandon and Calderwood (1990) are correct about the mechanics of formation of the arch, then we interpret from these boundary conditions that relatively strong coupling exists around the area of the plate arch. Otherwise, crustal motions would not have influenced subducting plate geometry.
The arch is also likely maintained by membrane stresses within the subducting JDFP as it is forced or pulled into the concave continental margin. Chiao (1991) has modeled the geometry of the arch that would be formed during forcing or pulling a viscous sheet into a margin like that of northern Cascadia. The modeled arch looks quite similar to that under the Olympic Mts. region as viewed in the contours of Figure 4b. Chiao (1991) points out that his modeled arch has a shallower subduction dip along its axis, as observed in our tomographic results. The shallower dip along the arch axis provides greater rigidity to the subducting JDFP. During periods when the plate is moving most rapidly under the concave margin, the arch probably becomes more prominent as the need to provide greater rigidity is met. During this period, north-south compression forces crustal blocks again the steeping sides of the arch. When JDFP-NA convergence velocities slow (due to either increased subduction thrust locking or decreased spreading rates at the Juan de Fuca ridge) contact with margin crustal blocks under north-south compression prevents collapse of the arch and this becomes a period of maximum traction forces along the deep thrust. North-south translation of mafic crustal blocks into the arch area can be viewed as causing asperities in the Puget Sound region, similar to the asperities caused by partial subduction of seamounts in several world subduction zones (Scholtz and Small, 1997).
Another factor which may maintain the arch is the draping effect caused by mineral phase changes in the JDFP below depths of about 50 km. The phase change from gabbro to eclogite causes a decrease in volume for dry minerals of about 17% (Ahrens and Schubert, 1975). Stanley (1984) inferred that this phase transition probably caused the dramatic change in dip in the JDFP under western Washington, as well as a large residual gravity high located over the dip change. In the curving part of the downgoing JDFP arch, the volume reduction acts like a tightening string in a draping skirt. Peacock and Wang (1998) state that the phase change may be beginning at ~45 km depth and be completed by 80-100 km depth under Cascadia. This interpretation is also supported by new phase diagrams for Cascadia developed by Bradley Hacker (writ. comm., 1999).
Other areas of northern Cascadia with mafic Siletzia crust outside the study area, such as in western Oregon, might be expected to have strong deep coupling between the two major plates. The primary difference between western Washington and Oregon is that the JDFP appears to dip at a steeper angle beneath western Oregon than below western Washington (Michaelson and Weaver, 1986). In a review of Cascadia seismic profile results, Trehu et al. (1994) infer that thick Siletzia crust in the Oregon coast range causes the steeper dip in the subducting plate. Trehu et al. (1999) now interpret high-reflectivity zones in seismic data from the Oregon coast range as showing direct contact between the JDFP and overlying Siletzia mafic crust at depths of 35-40 km. The recent GPS results of Goldfinger et al. (1999) have been modeled to indicate that subduction thrust coupling extends to depths of about 40 km under the coast range of Oregon. These two results seem compatible with our interpretation of a mafic-mafic locked zone under Puget Sound.
Magnetotelluric models have been advanced that seem to contradict the results of Trehu et al. (1999) and Goldfinger et al. (1999). EMSLAB magnetotelluric models Wannamaker et al., (1989) were interpreted to show an east-dipping conductive zone under western Oregon that they inferred to be associated with subducted sediments between the mafic lower crust of the continental plate and the subducting JDFP. MT models of data from across Vancouver Island (Kurtz et al., 1990) were interpreted to indicate a similar conductive region to that postulated by Wannamaker et al. (1989). However, as noted by Jiracek et al. (1989), both of these models were derived using forward modeling with two-dimensional (2D) programs and were not generalized inversions of the data. Jiracek et al. (1989) modeled the same data set as Wannamaker et al. (1989) using a 2D Backus-Gilbert inversion and found little evidence for an east-dipping conductive zone above the JDFP. In addition, Jiracek et al. (1989) noted the evidence for three-dimensional effects in the data which require caution in using 2D modeling methods for the EMSLAB data. Aprea et al. (1998) have used a 2D inversion method to model MT data across the Olympic Mts.-Seattle area of Washington and although they found a restricted mid-crustal conductor (10-20 km depth) under the Seattle Basin region and area just to the west, there was no inferred conductor just above the JDFP. Our analysis of the Aprea et al. (1998) data suggests to us this conductive zone (their conductor C4) is likely an artifact of applying 2D methods to the highly 3D geometry of the Seattle Basin and Olympic Mt. region.
North-South Transect Interpretation
A north-south geological interpretation is shown in Figure 44b. This interpretation was compiled by combining features on three closely spaced velocity cross-sections at columns 18,19, and 20. The structure changes rapidly in an east-west direction between these slices, but they represent the most critical data required to understand the JDFP-continental crust dynamics in the region. In addition to features noted in velocity model columns 18-20, the location of the Everett and Bellingham basins from columns 23-25 are shown on the interpretation. Also, characterization of the mixed Crescent Formation basement and interleaved thick Eocene-Miocene sedimentary rocks in southwestern Washington as interpreted from MT models is indicated, since the velocity model did not sort out this complexity. The Black Hills block of Crescent Formation is shown as indicated by the Vp/Vs results and from gravity data (Figure 2). This block is bounded on the north by the NW-SE oriented 'Legislature fault', a term informally used for the probable reverse fault required by very steep gravity gradients there (Gower, et al., 1985).
Mid-crustal velocities along the transect and in most of western Washington are in the range of 6.5-7 km/s, which we identify as Crescent Formation without interleaved sedimentary rocks. Lower-crustal velocities are 7-7.75 km/s, a high value that we assign to the mafic/ultramafic cummulates inferred in the bottom of the rift system associated with formation of the Crescent Formation. We interpret that this high-velocity section has been tectonically thickened over the axial region of the JDFP arch, as suggested by the velocity model. The gap between the northern and southern portions of the mafic wedge may have developed from separate magma bodies within the initial rift. North-south translation in the margin sliver has continued to close this gap, leading to tectonic thickening as the northern and southern parts of the mafic/ultramafic complex come into contact. The northeastern nose of the thick accretionary prism sediments in the Olympic Peninsula are being compressed in the gap and surrounding Crescent Formation also must conform to the geometry of the gap. We infer that the formation of the Seattle basin is partially related to the greater compliance for downwarping over the gap as the thick block of Crescent Formation south of the Seattle fault is driven northward. A light blue pattern is assigned to the upper part of this thick block of Crescent Formation in Figure 44b, because it appears to have significantly lower velocities than the main Crescent Formation below. Similarly, we used a light blue pattern for the Black Hills on the interpretative section, because it was well resolved only in the Vp/Vs results.
We infer that separate blocks of Crescent Formation may occur just over the gap, under the Seattle basin, as suggested by looking at the rapidly changing geometry of Crescent type velocities in this area. A contact between these blocks or a fault in one of them may explain a narrow zone of earthquakes (Figure 47 and 48) in the basin basement that extends northeastward from the Seattle fault offset (Johnson et al., 1999a) across the southern part of Seattle. These earthquakes are part of two swarms studied by Yelin (1982). He located one swarm 10 km east of Seattle, in Lake Washington, trending NE, at depths of about 20 km. A second swarm was located under Seattle, also trending NE, and located at depths of 10-12 km. Relocation of selected earthquakes in these swarms with our 3D velocity migrated them to a tight northeast zone through Seattle. These earthquakes are also shown in seismicity cross-sections of figures 51 (row 22 and 24) and 52 (column 24).Yelin (1982) identified focal mechanisms for several of these earthquakes with ENE and NE P-axes, which he interpreted was direct evidence for coupling of the upper crust to the subducting JDFP. We agree with this interpretation and see other evidence for seismicity lineaments that reflect JDFP northeast-directed shear stresses. The most noticeable of these is indicated in Figure 48 where a northeast trend of earthquakes at depths of 10-20 km cuts across southern Hood Canal toward the west half of the Seattle fault zone.
Earthquakes from column 19 are shown plotted on the interpretive section. Focal mechanism for the two large intraslab earthquakes (nos. 5 and 6, Figure 44b) and three important crustal earthquakes near the transect (nos. 1,2, and 3) are plotted on the cross-section. T-axes for the intraslab earthquakes (5 and 6) are oriented downdip in the steeply dipping part of the slab (Figure 44a). This T-axis orientation likely reflects slab pull and slab bending moment stresses. Alternately, Kirby (1998) has suggested that phase changes may be primarily responsible for the downdip T-axes. Kirby et al. (1996) calculated that densification from mafic oceanic crust being transformed to eclogite places the subducting oceanic crust in deviatoric tension and the oceanic mantle in compression. Hacker (1996) points out that differences in reaction rates between fine-grained basalt and coarse-grained gabbro means that the Kirby et al. (1996) calculations of deviatoric stress in the oceanic crust are oversimplified. In addition, since the large intraslab events appear to have occurred within the JDFP mantle (Figure 7e) at velocities of over 8.0 km/s and not in the oceanic crust, the Kirby et al. (1996) explanation may not be valid for these earthquakes.
Earthquakes in the JDFP crust are occurring, but appear to be small (white dots, Figure 44b). Ma et al (1996) have done an extensive study of focal mechanisms for 42 intraslab earthquakes, including deriving principal stress directions using the method of Gephart and Forsyth (1984). Ma et al (1996) found that the T-axes for intraslab events (including a high-quality subset of eight events with magnitudes >M3.0) were in agreement with that expected from bending stresses due to known arching of the JDFP. If these events had been connected primarily with densification phase changes the T-axes would have been oriented downdip, but in reality they were primarily oriented along contours of the arch. We contend that the downdip, gradational metamorphic front associated with this phase change is only one of the factors to be considered in understanding the localization of large intraslab earthquakes in Cascadia. There were three and possibly four such events located in southern Puget Sound (Figure 4b) and most of the intraslab energy release has occurred in a very small part of the JDFP. Southern Puget Sound occurs near a double flexure point for the JDFP, with dip changing from about 10o to 35o and just off the south part of the nose of the arch. In addition, there appear to be even smaller-scale wrinkles in the slab, as indicated in Figures 9 and 10, which could greatly concentrate bending stresses (Chapple and Forsyth, 1979). The southern Puget Sound intraslab earthquakes also occurred just downdip from the thickest part of the mafic wedge, which we infer may be partially locked to the JDFP. We interpret that a combination of bending moment stresses directions, a locally strong locked thrust updip, integrated slab pull forces from the deep slab, and possibly localized phase change effects all contributed to the concentration of large intraslab earthquake in southern Puget Sound.
Related to the issue of gabbro-ecologite transition as a mechanism for producing intraslab earthquakes, we find it difficult to determine the downdip extent of the oceanic crustal layer in our velocity models. Down to depths of 40 km and to longitudes of about -123o the velocity models (Figure 7) show a structure compatible with having about 7 km of oceanic crust with velocities of about 7 km/s. At larger depths downdip, we do not detect a consistent low-velocity layer which could be construed as oceanic crust. Velocity reversals under small portions of the mafic wedge (Figure 7d) are assumed to be due to very high velocities in the middle of the wedge (>7.75 km/s). The gabbro-eclogite transition should raise oceanic crustal velocities from about 7.0 to more than 8.0 km/s (Gubbins et al., 1994), and where this transition has occurred to completion the oceanic crust would not be detectable. Thus, the two intraslab events in 1949 and 1965 could have occurred in eclogitic oceanic crust, but this would place the phase transition completion zone very high in the JDFP. The eclogite transition would be complete at depths of 40-50 km, which seems very unlikely (Peacock and Wang, 1998). Difficulties in mapping a deep oceanic crustal layer before it has undergone the eclogite phase transition have been reviewed by Gubbins et al. (1994). He points out that most slabs show up in tomographic studies as high-velocity features, with a notable exception found from SW Japan, where a waveguide 10 km thick and velocity of 7.0 km/s was interpreted to persist down to only 50-60 km (Oda et al., 1990) and then vanish. The evidence for the waveguide came from inferred trapped mode energy. Oda et al. (1990) note that the waveguide could have continued down deeper than 50-60 km, but that earthquakes at these depths were located in the underlying slab mantle and did not produce waveguide energy. The alternative, preferred by Kirby (1996), is that the waveguide was terminated by the gabbro-eclogite transition.
The three crustal focal mechanisms shown are important for tectonic analysis of the region
near Seattle. The Robinson Point earthquake (no. 3) has been studied by Dewberry and Crosson (1996), who found
that it was a high-stress-drop event with a reverse faulting mainshock. They mentioned possible correlation with
the Seattle Fault and also that the earthquake might have initiated near the brittle-ductile transition. We infer
in Figure 44b that the Robinson Point earthquake initiated near the base of the Crescent
Fm., possibly at a decollement as envisioned by Pratt et al. (1997). However, it is unlikely that it was related
to the brittle-ductile transition, since the rocks are all mafic and temperatures probably only about 400o
C. The Port Orchard earthquake (no. 2) is next most significant earthquake near the transect. Yelin and Crosson
(1982) studied this earthquake and we used the fault plane from their focal mechanism determination that trends
N78oW, and dips 60o to the south, which seems to fit the pattern of other crustal earthquakes
in Puget Sound better than the alternative (N18oW strike, 80oE dip). This earthquake may
have been associated with a decollement at the base of the Crescent Formation or on a steeply dipping boundary
fault between blocks. The final crustal earthquake included on the transect interpretation is the June, 1997 Bremerton
M4.9 (no. 1) earthquake. Oregon State catalog moment tensors for this event had one nodal plane similar to that
of the Port Orchard event, striking N85oW and dipping 65o S. Art Frankel (writ. comm.,1999),
USGS, has also computed a focal mechanism for the Bremerton earthquake using local array data and obtains a result
similar to the Oregon State solution, but with the focal depth at 6 km. This earthquake is commonly assumed to
have occurred on the Seattle fault zone, an assumption that appears to be a realistic interpretation.
Comparison of Cascadia With Other Subduction Zones
As noted by Rogers et al. (1996), one of the closest analogs to the Cascadia subduction zone in regard to plate age, subduction rate, and size may be the Rivera-Cocos plate system off western Mexico. The ages of the Rivera and Cocos plates are less than 10 Ma (Klitgord and Mammerickx, 1982). Subduction rates vary from about 2.3 to 6.8 cm/y for the Mexican portions of Rivera and Cocos plates (Eissler and McNally, 1984; Pardo and Suarez, 1995). These age and subduction rate parameters are similar to those of the JDFP.
The Cocos plate is much larger than the JDFP, extending under the region from Manzanillo, Mexico to Colombia, and thus may not be a particularly good detailed analogue for the JDFP. However, the Rivera plate is very similar in geometry to the JDFP. Plate velocities are similar and there is a transition from oblique convergence to mostly strike-slip involving small microplates and stepping transforms. The main difference in the two analogous plates is the lack of a thick accretionary prism and plate arch for the Rivera plate.
The Rivera plate, as discussed by Eissler and McNally (1984), has low levels of oceanic plate seismicity, similar to Cascadia. Several large (M7-8) interplate earthquakes associated with the larger Cocos plate to the south have occurred. Aftershocks of several of these large earthquakes indicate that stress drops are occurring in a region of the subduction interface from 10 to 25 km in depth (Singh, 1988). The mainshocks generally originated in the deeper parts of this region. Because of the similar ages and subduction rates, it might be expected that this seismogenic depth zonation would be similar to that of Cascadia. However, the Rivera-Cocos system has a very narrow accretionary prism with only 0.6 km of sediments. The apparent shallow seismogenic region is likely due to the absence of a thick, wet, accretionary prism required to lubricate the subduction zone and force the updip seismogenic edge to greater depths (Byrne et al., 1988).
Tichelaar and Ruff (1993) and Byrne et al. (1988) also note an apparent shallow limit to aftershocks of several large subduction earthquakes in western Mexico. However, this pattern may be an oversimplification, based on recent seismicity analysis and focal mechanism solutions by Pardo and Suarez (1995). They interpret that the maximum depth of seismogenic coupling for the Rivera plate is about 40 km, based upon focal mechanisms of earthquakes with hypocenters accurately located using both local and teleseismic recordings. Pardo and Suarez (1995) infer that the Rivera plate has a deeper coupling extent and could generate larger earthquakes than the Cocos plate. In the Guerrero seismic gap within the Cocos plate system (roughly between Zihuatanejo and Acapulco), Kostoglodov et al. (1996) found concentrated zones of low-level, thrust seismicity at depths of 10-22 km on the subducting plate interface, merging into a deeper zone of seismicity that extends to 40 km depth. The events were relocated using a tomographic velocity model and appear to us to be the highest quality seismicity sections available for the Mexican subduction zone. The seismicity patterns outline a section of the subducting plate dipping at about 5-100 from the trench to the coast, and a more steeply dipping interface (27-320 ) down to 40 km depth, where another inflection in plate dip occurs. Focal mechanism solutions indicate thrust motion along the plate interface in both of the seismogenic regions above 40 km depth. The exceptions were two normal faulting, intraslab mechanisms near the two abrupt changes in plate dip, which fit the concepts of Spence (1987) for the relationship between locked regions, plate dip changes, and tensional intraplate earthquakes. It appears that most of the large western Mexican earthquakes (Singh, 1988) initiated in the upper, shallower dipping part of the subduction interface, but the 1979 Petatlin M7.6 earthquake (Valdes et al., 1982) ruptured the deeper, more steeply dipping part of the interface, fitting the pattern of recent seismicity noted by Kostoglodov et al. (1996). The aftershock pattern of the 1979 Petatlin earthquake overlaps that of the lesser (M7.5) of the two great Michoacan 1985 earthquakes (UNAM, 1986). The fact that the Mexican subduction interface is seismogenic to depths of over 40 km has important implications for northern Cascadia. We infer that the deeper Mexican interplate region of seismicity is possible because the subducting oceanic lithosphere is in contact with mafic lower crust or mantle wedge in the upper plate at these depths, as we believe may be the case in the Puget Sound and Vancouver Island region.
Another applicable study from the Mexican subduction zone comes from McNally et al. (1986), who found that aftershocks consisting of downdip-tension, intraslab earthquakes preceded the 1973 Colima (M7.6), 1979 Petatlan (M7.6), and 1985 Michoacan (M8.1) interplate earthquakes during the 2.5-4 years before the main interplate events. Intraslab downdip, tensional earthquakes were not observed in the few years following the mainshock, interplate events. This suggests that the thrust was locked until intraslab tension build up to the point where the traction forces were overcome and both interplate and intraslab forces were readjusted. This analysis follows that of Spence (1987) for the dynamics of slab pull and locked interplate thrust regions. The 1949 and 1965 intraslab events in southern Puget Sound indicate that the deep part of the thrust is still locked and other large intraslab and subduction thrust earthquakes are expected.
Hyndman et al. (1997) noted that some island arc subduction zones apparently have two seismogenic zones downdip along the interplate thrust. The deeper of the two zones occurs at depths of 40-50 km (Pacheco et al., 1993). Hyndman et al. (1997) suggest that the upper seismic zone occurs above depths where hydrated serpentinite causes stable sliding, rather than stick-slip behavoir. Similarly, they infer that the deeper seismic zone may occur where serpentinite begin to be dehydrated and supports stick-slip behavoir, in the range 400-600o C, as documented in lab studies by Moore et al., 1986. This is very close to our interpretation for western Washington, except that the deeper zone occurs at depths of 30-40 km because of the higher temperatures along the thrust zone in Cascadia.
Implications for Interplate and Crustal Seismic Hazards
Subduction Zone Source Ground Motions
Our hypothesis of a deep source region for interplate earthquakes has important implications for ground motion hazards in the highly populated Puget Sound region. The model presented by Hyndman and Wang (1993,1995) has the nearest part of the locked region approximately 200 km from Seattle. The analysis of Spence (1987), Tichelaar and Ruff (1993) and Byrne et al. (1988) indicates that nucleation of large interplate thrust earthquakes is most typically in the deeper portions of the locked zone, with slip and aftershocks propagating updip from the source region. This means that a large interplate earthquake in the model of Hyndman and Wang (1993,1995) would largely involve the part of the megathrust 200-300 km from Seattle. In the model that we propose, the source region also may include the interplate thrust at depths from 30 to 40 km (Figure 44a) where mafic /ultramafic rocks in the upper plate are interpreted to be in contact with oceanic crust. This contact occurs at distances of 50-80 km from Seattle. This postulated deeply coupled p
