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

 


 

Frame Version of Contents

CONTENTS

Abstract

Introduction

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

Mafic Wedge

Olympic Mts. Sedimentary Complex

Juan de Fuca Plate Arch

Seattle Basin and SWCC

Crustal Terrane Boundary

3D Visualization of Velocity Results

Isosurfaces

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

Summary

Acknowledgments

References

Appendix A-Resolution Tests

Figure Captions

Illustrations


SUBDUCTION ZONE AND CRUSTAL DYNAMICS OF WESTERN WASHINGTON: A TECTONIC MODEL FOR EARTHQUAKE HAZARDS EVALUATION

By

Dal Stanley, Antonio Villaseñor, and Harley Benz

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.

Introduction

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.


Geological Setting

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 Profiles

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.


Faults and Seismicity

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).


Velocity Model Construction

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.

Juan de Fuca Plate

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.

Mafic Wedge

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.

Juan de Fuca Plate Arch

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.

Seattle Basin and SWCC

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).

Crustal Terrane Boundary

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.

Isosurfaces

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.

Sedimentary Basins and SWCC

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).

S-wave Model

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 part of the megathrust wraps around the nose of the plate arch, as outlined by the view of the mafic wedge shown in the velocity isosurfaces of Figures 9 and 11. The coupled region can be extended further northwestward along Vancouver Island based upon seismic and gravity evidence for mafic rocks in the upper plate (Clowes et al., 1997). The northern boundary of the onshore, locked part of the subduction thrust is interpreted to extend to the onshore projection of the Nootka fault zone, where two large crustal earthquakes in 1946 and 1957 (Figure 4a) occurred. Hyndman et al. (1979) interpreted that these large events were due to stress coupling of the subducting plate to the Vancouver Island crust, therefore we extend our onshore locked region as far as the projection of the Nootka fault (Figure 4a). This postulated deep locked zone of the megathrust extends for over 250 km around the curving subduction zone. If the mafic crust in the Oregon coast range is also locked to the JDFP as suggested by Goldfinger et al. (1999) and Trehu et al. (1999), then the area of concern is much larger.

A megathrust source model that includes the deep locked zone has not been modeled directly, but some indications of the relative ground motions that could be expected from the deeper source can be derived from an analysis of simulated ground motion by Cohee et al. (1991). They estimated the amplitude and duration of strong motions from hypothetical M8 subduction zone earthquakes in northern Cascadia by assuming that the source zone could be characterized by individual asperities. Primary S-wave and Moho postcritical reflections were calculated at horizontal distances from 0 to 150 km using a 2D velocity structure derived from refraction models (Taber and Lewis, 1986) for western Washington. Scattering and attenuation were empirically modeled using corrected accelerograms from M7 earthquakes at Michoacan, Mexico and Valparaiso, Chile. The technique was validated by modeling M8 earthquakes from Michoacan and Valparaiso. Site conditions were also included in the simulations. Peak horizontal ground accelerations (PGA) were calculated assuming that 60% of the slip occurred on either shallow (18.9 km), middle (26.5 km), or deep (34.2 km) asperities with a strike length of 150 km in the Seattle region. For Puget Sound, the relationships between PGA and distance to the fault were highly variable for the shallow, middle, and deep asperities, but the picture was somewhat simplified by plotting PGA with respect to distance to the asperity. In this case, a simple relationship between PGA and distance to the asperity was evident. For the Hyndman and Wang model (1993,1995), the distance from Seattle to the asperity is slightly beyond the plot limits, but at the maximum distance of 150 km used by Cohee et al. (1991), PGA was about 0.04g for a bedrock site. If the asperity caused by the deeply locked thrust in Figure 44a corresponds to a depth of 50 km and horizontal distance of 50 km from Seattle (actual distance of 70.7 km), then the PGA from Cohee et al. (1991) would be about 0.16g for a bedrock site, a factor of four larger than the Hyndman and Wang model. Cohee et al. (1991) also computed simulated pseudo-velocity (PSV) spectra for their asperity models and found that the horizontal response spectra for modeled interplate events was similar in shape to that of the 1949 M7.1 and 1965 M6.5 intraslab earthquakes in southern Puget Sound. For a M8 subduction zone earthquake, the estimated velocities between 1 and 10 Hz were 2 times larger than the M7.1 intraslab earthquake and the estimated duration of strong motion was 45-60 sec versus 10-15 for the M7.1 event. The simulations by Cohee et al. (1991) are useful to point out the dangers in assumptions about the interplate source zone being well offshore.

Seismicity Patterns

With a background provided by the preceding discussions of velocity structure and tectonics, we provide more details of seismicity patterns for completeness and additional insights. We start with a set of maps of seismicity for four depth ranges. The map of Figure 46 shows all of the earthquakes from the Pacific Northwest Seismic Network, including those relocated using the tomographic velocity model. The general pattern of seismicity and tectonic features is outlined, but depth slicing is needed to isolate some of the patterns. In Figure 47 earthquakes at depth from 0-10 km are plotted. For this depth range, dense shallow seismicity occurs along the SHZ and WRZ, near Seattle, and along the DAR. There is a distinct zone of earthquakes trending NE from the southern third of the WRZ. Numerous patches of shallow seismicity also occur in the central Columbia Basin region, especially along east-west faults within the folded basaltic rocks that occupy this region. The dense patch of seismicity near the city of Wenatchee occurs along the intersection of NW trending faults in the North Cascades with the margin of the Columbia Basin.

Earthquakes in the depth range from 10-20 km are shown in Figure 48. The trend of the SHZ and WRZ is still outlined by earthquakes in this depth range, and there are interesting stepover features from the SHZ to the WRZ. The most evident stepover extends NE across Mt. St. Helens and then across the southern end of the WRZ where the shallower (0-10 km) earthquakes delineate a NE trend. This zone actually cuts across the SHZ near Mt. St. Helens from the southwest, suggesting a fault or contact whose intersection with the SHZ may be important in controlling the location and magmatic access of Mt. St. Helens. The other, more diffuse stepover pattern extends from the northern part of the SHZ to toward the WRZ at the latitude of Mt. Rainier. There is a dramatic increase in the number of earthquakes in this depth range in the Puget Sound region surrounding Seattle, extending eastward to the margin of the North Cascades. In addition, earthquakes are concentrated along the DMF-DAR.

For depths of 20-30 km (Figure 49), seismicity is largely concentrated in the Puget Lowland region, centered on Seattle. Except for a band of earthquakes near the eastern margin of the Puget Lowland, the event concentration appears to be largely controlled by stress near the base of the key block of Crescent Formation portrayed in the interpretation of Figure 44b. The stress focus is likely due to the large block of Crescent Formation being driven along the contours of the plate arch by north-south compression and northeast convergence. This block may be tightly coupled to the underlying mafic wedge complex and JDFP below that. These earthquakes involve some subtle NE banding south of the Seattle fault and NW banding north of the Seattle fault. This subtle banding may be evidence for shearing in the directions indicated by focal mechanisms, with P-axes trending NE in southern Puget Sound and NW in northern Puget Sound.

For depths of 30-60 km (Figure 50), earthquakes are primarily within the upper part of the subducting Juan de Fuca plate and reflect bending stresses within the arching part of the plate under the Olympic Peninsula and Puget Sound region. The events outline the nose of the arch, just updip from where the dip in the subducting plate changes from about 12 degrees to more than 35 degrees (Figure 44a). Earthquakes from Seattle to Olympia are co-located with overlying crustal earthquakes and underlying intraslab events in southern Puget Sound. We interpret this co-location as evidence that the upper crust, lower crust, and JDFP are tightly coupled in southern Puget Sound and possibly elsewhere around the arch.

To continue investigation of seismicity patterns, we used a set of east-west (in row direction from velocity model) cross-sections of seismicity shown in Figure 51. Earthquakes are plotted on top of P-wave velocity contours. A few details of these sections are worthy of note:

The western Rainier zone (WRZ) is evident on row 10, possibly as an east-dipping reverse fault. The WRZ seismicity is at a maximum on row 12, which crosses the Mt. Rainier area. Seismicity under Mt. Rainier is largely due to gravitational effects from the massive volcanic edifice and rock and ice falls. The WRZ linear seismicity and earthquakes under Mt. Rainier continue on row 14. The final expression of the WRZ may be on row 16, where earthquakes form a diffuse NS zone northwest of Mt. Rainier continuing into the Puget Lowland (Figure 48, 10-20 km depth slice). On row 18, no seismicity is associated with the Tacoma basin, but events near the top of the Juan de Fuca plate mantle (~7.75 km/s) are prominent. In addition, earthquakes concentrated at depths of 15-25 km are common west of -122o longitude. As suggested in the previous sections, these events may be associated with the base of the Crescent Fm, but are probably not related to the brittle-ductile transition because the rocks are mafic. However, brittle-ductile control may be exhibited in shallowing of events to the east of the terrane boundary, as indicated in rows 18 and 20, where the crust is silicic (pre-Tertiary metamorphic rocks) and heat flow increases near the Cascades volcanoes. Seismicity in the Crescent Formation basement is more intense at row 20. The terrane boundary between Crescent Formation on the west, juxtaposed with Cretaceous and older metasedimentary basement on the east, is somewhat easier to define from velocity and seismicity characteristics. The velocity expression is subtle, consisting mainly of a depression of 5 km in the 7 km/s contour. The seismicity expression is somewhat clearer, with an eastern end to the depth- banded seismicity in the Crescent Fm., followed by a slight gap in seismicity east of that point, and then a transition to seismicity which is distributed both vertically and horizontally in the pre-Tertiary crust. Continuity of the seismicity near the base of Crescent Formation is most evident on rows 22, 24, and 26. The terrane boundary seems well defined on these sections, exhibiting the characteristics just noted for row 20. There is a distinct northeast trend of earthquakes recorded in the basement of the Seattle basin, which have been studied by Yelin (1982) as noted earlier. These events may be located on a northeast trending fault within a Crescent Formation block or between blocks as portrayed in the interpretation of Figure 44b. The terrane boundary is expressed as a narrow, vertical zone of seismicity at depths from 22 to 30 km on rows 22 and 24.

One significant feature of rows 26, 28, and 30 is the north-trending, east-dipping fault that is evident on the velocity model at -121.8oW. The May 3, 1996 M5.3 Duvall earthquake is interpreted to have occurred on this reverse fault, as indicated by the moment tensor solution for nodal plane 1 (strike 350o, dip 52oE, rake 78o). Historic earthquakes apart from the Duvall event and its aftershocks continue on this fault zone to row 28 and possibly to row 30. Except for the events associated with this reverse fault, seismicity patterns on rows 28, 30, and 32 become increasingly more diffuse as the more intense seismicity in southern and central Puget Sound diminishes.

North-south, or column direction cross-sections of earthquakes are shown in Figure 52. On column 8 only minor seismicity has been recorded near the top of the Juan de Fuca plate. A few more events occur on column 10, mainly near the top of the Juan de Fuca plate. Continued increases in seismicity near the top of subducting plate are present on the column 12 section, with very little difference on column 14. A well-defined outline of the top of the Juan de Fuca plate from the seismicity is evident on column 16, and there is an increase in upper crustal seismicity. Dramatic increases in crustal seismicity in the same region as basin formation occur on column 18 and most of the seismicity is inferred to occur near the base of the Crescent Formation (at a decollement ??) where it is in contact with the underlying mafic wedge. Maximum crustal seismicity occurs where Crescent Formation fills the gap in the mafic wedge (from column 18 to 20). Continuation of extensive crustal seismicity is seen on column 22, but somewhat diminished near base of Crescent Formation under the Seattle basin region. The column 24 section crosses the SHZ (tight vertical pattern at about 46.2 degrees latitude) and a northeast trending stepover towards the WRZ. A tight pattern of events in the basement of the Seattle basin is the northeast trending alignment of earthquakes that passes through southern Seattle and studied by Yelin (1982). Some deep crustal seismicity occurs near the Devils Mt. fault zone (DMF) on columns 22 and 24 (and possibly 26), which separates slivers of pre-Tertiary crust on the south from the more massive pre-Tertiary crust of a northern buttress region under southern Vancouver Island and the North Cascades. For column 26, east of the terrane boundary in the Puget Lowland, crustal seismicity in pre-Tertiary rocks spreads vertically throughout the crust between 47 and 48 degrees latitude. On column 26 the sliver of pre-Tertiary crust between the SWIF and the DMF-DAR is seismically quiet, but some deep seismicity occurs near the DMF-DAR. Diffuse expressions of the SHZ and WRZ can be noted on the column 26 plot. By column 30, coherent patterns in the seismicity become more difficult to detect.

 Deformation Modeling

Model Construction

A key part of understanding the tectonics and earthquake hazards of western Washington involves modeling of deformation and stress data for the area. The western Washington part of the Cascadia subduction zone is complex, with the JDFP arch plunging in the direction of subduction dip and crustal blocks showing strong evidence for north-south compression. To further complicate matters, P-axes (Mulder, 1995) from earthquake focal mechanisms in southwestern British Columbia, northern Puget Sound, and the North Cascades region trend northwest (Figure 53). The pattern of stress apparently swings from NE, through NS, and then to NW across the Puget Sound region. This rapid change in stress direction reflects a large amount of work performed during anelastic deformation of the lithosphere. If the lithosphere behaved strictly elastically, then the work would become potential energy. Work done for a unit volume in the principal shear stress direction is equal to the product of the stress and strain rate in that direction. This anelastic work over a long time period is converted to either heat or dislocation energy release (earthquakes). When combined with modeled and geological estimates of crustal shortening, an analysis of work functions could be conceivably be applied to model total crustal energy release from earthquakes in the Puget Sound region.

To study the stress complexities and other important features, we employed thin-plate, finite-element models of the Cascadia subduction zone using programs developed by Peter Bird, UCLA (Bird et al., 1996). These finite-element programs predict anelastic strain rate, stress, and velocity under the assumptions of a given crustal and mantle rheology. Other inputs include crust and mantle thicknesses, heat flow, and elevation for each node point. Deformation is by frictional sliding or nonlinear dislocation creep (whichever requires less shear stress). The program neglects all elastic strain accumulation and release, but solves for time-averaged velocities, fault slip rates, anelastic strain rates, and stresses. Rheological assumptions can be varied, as can the frictional coefficients for crustal faults and internal frictional limits for the continuum elements. Low-angle subduction thrusts are assigned traction limits which are independent of other faults. Velocity boundary conditions are specified around the edges of the finite element grid. We used Bird's flat-earth program, 'PLATES', for our study, under the belief that earth curvature does not affect the pattern of stress in the region. For larger regions, Bird has developed a spherical earth approximation, called 'SHELLS'.

The finite-element grid consists of triangular, six-node isoparametric elements representing the continuum for both crust and mantle layers (Figure 54). Faults are specified as separate isoparametric elements of infinitesimal width with three nodes in common to adjacent triangular continuum elements (Bird and Baumgardner, 1984). The faults are assigned a lower coefficient of friction than the continuum, which can be studied parametrically. Bird (1996) has supplied programs for constructing the finite element meshs and inclusion of faults, as well as for display of the computed results in Postscript form (Bird, 1996). However, we found it more convenient to develop our own interactive mesh generation and display software. For mesh construction we utilized the public domain, 2D-grid-generation program 'TRIANGLE' written by Jonathon Skewchuk (1996) of Carnegie Mellon University (http://www.cs.cmu.edu/~quake/triangle.html). This program is very flexible and has its own X-windows viewer for the constructed grid. It allows specification of linear subdomains, which can denoted as fault traces. The program will automatically take care of finer gridding required along the faults (subdomains). In addition, regions of the grid can be specified to have different triangle sizes. In order to employ this useful program we developed an interface to the input and node numbering specifications for 'PLATES'. The interface program was written using commercial display software and allows interactive design of the grid region and fault traces in X-windows. The fault traces are selected from previously digitized and plotted fault locations; a few fault points are selected with the mouse to determine grid density along the fault. The meshing is then handled through 'TRIANGLE'. Another interface program that we developed takes the outputs from Bird's 'PLATES' program and plot the results, along with the grid, and other parameters as needed. This program also can be used to interactively input and modify velocity boundary conditions, as well as sample elevation and heat flow data files for new grids.

Velocity boundary conditions were based on Juan de Fuca velocities relative to North America from Wilson (1993); Pacific plate motions relative to North America were taken from DeMets et al. (1990). On the continental interior part of the model, initial velocity conditions were taken from Basin and Range shear values estimated from very-long-baseline-interferometry (VLBI) measurements (Angus and Gordon, 1991), GPS data (Prescott et al., 1998), and fault shortening calculations (Pezzopane and Weldon, 1993). In the final or best models, we adjusted the interior values as necessary to model stress parameters. Topography was taken from the ETOPO5 digital elevation data set available from NOAA (http://www.ngdc.noaa.gov/mgg/global/global.html). Heat flow values (Figure 54) were obtained from Dave Blackwell's data used for the heat flow map of North America (Blackwell and Steele, 1992). Crustal thickness values were selected from various seismic refraction profiles (Schultz and Crosson, 1996; Miller et al., 1997; Parsons et al., 1998) and from our own velocity models. These ranged from 12 km for oceanic crust to over 40 km in the Cascades part of the model where the crust is thickest (Schultz and Crosson, 1996). Models were tested both with constant mantle thickness and with isostatically adjusted thicknesses. In the latter procedure, as suggested by Bird (1996), mantle thickness were adjusted to obtain an isothermal surface at the base of the lithosphere. This involved an iterative calculation which was not always satisfying for our data, because of extremely high heat flow contrasts across the model. Exaggerated strain anomalies were found in important areas such as near Mt. Rainier and in the Oregon Cascades where there are strong heat flow gradients. We believe these are not realistic and reverted to constant mantle thicknesses of 40 km under the continental part of the model and 22 km under the oceanic part. Combined with the rheological parameters used, the model reflects mainly variation in crustal parameters which are better known than mantle compositions and thicknesses.

Parameters kept fixed in all of the models were crust(mantle) densities of 2800(3330) kg/m-3 (at 0oK); volumetric thermal expansion coefficients of 2.4x10-5(3.1x10-5 )oK-1; thermal conductivities of 3.0(4.06) W/m-1 oK-1; radioactive heat productions of 4.55x10-7(3.2x10-8) W/m-3; asthenosphere temperature of 13000 C and density of 3168 kg/m-3. Power law (n=3), thermally controlled, dislocation creep was assumed to dictate maximum possible shear stress, , according to the given rheology as defined by Kirby (1983):

 


Where es is shear strain rate, z is depth, and T is absolute temperature. The initial value of the constants A, B, and C for crust(mantle) were A=2.3x109(5.4x104) Pa/s1/3; B=4000(18,314)oK; C=0(0.17)oK/m-1. The maximum shear stress at any point (plasticity limit) was 500 MPa. It was quickly discovered that we had to assume a more mafic composition to the crust in order to simulate the known stress and velocity parameters. This is required because of the dominant control of north-south compression by thick mafic crust in the Oregon and Washington coast ranges. The best value for rheology involved setting B=9000oK, a value between that of diabase and quartz diorite (Kirby, 1983), which may be similar to that of the basalts found in the coast ranges. This effectively deepened the brittle-ductile transition over the region. Although not realistic for eastern Oregon and other areas, this was the only solution possible apart from modifying the program to allow specific rheology at each node point. This modification is a worthy goal, but not yet implemented.

Crustal faults (Figure 55) were assumed to be either vertical strike-slip faults or reverse faults dipping at 45o . The JDF subduction thrust was assumed to dip at 12o, although tests were run at different dips, since dip under western Washington changes as one moves inland and from north to south. A selected set of borehole breakout data (Zoback and Zoback, 1991), focal mechanisms from earthquakes, fold axis directions, and dike systems as indicator of regional stress were simulated in the models (Figure 53). In addition, we attempted to fit recently obtained GPS horizontal motion data. Very little geological shortening data is available for crustal faults in Washington, although we will subsequently discuss some data from Johnson et al. (1999a) for the Seattle fault and recent estimates for the Devils Mt. fault by Johnson (oral comm., 1999).

More than 500 separate runs were made with different grids, boundary conditions, and physical parameters. Many runs were initially made with only five iterations to see if the changes were reasonable, and if so, then the computation was continued to 50 iterations. An acceptable velocity error for convergence was specified at 1E-10. For simple models convergence was reached before 50 iterations, but in the more complex final grids, convergence to the specified error was not reached. However, for interpretive purposes, the models converged to undiscernible stress, strain, and velocity differences by 10-20 iterations. The models were calculated on a Sun UltraSparc 10 (300 MHz) and required about 40 minutes per run for 50 iterations. Model stresses represent integrated vertical stress referenced to the standard density model of a mid-ocean rise (Bird, 1996).

Model Traction Limits

We found that the best fit to stress patterns and surface velocities was obtained by using a subduction thrust traction limit of 10-30 MPa, with anything in this range being acceptable for simulating stress patterns. However, to best fit the GPS horizontal velocities, we needed to increase thrust traction limits to 30 MPa in order to have high enough velocities in the coast areas of Washington. The range of 10-30 MPa for thrust tractions compares well with that found by Tichelaar and Ruff (1993) for world subduction zones. For constant stress assumptions, they found tractions of 14±5 MPa and for assumption of constant friction coefficient they obtained tractions of 30±5 MPa.

North-south shear from Pacific plate coupling to North America at the Mendocino and Queen Charlotte transforms and the NE shear from JDFP subduction are delicately balanced with NW shear from extension in the Basin and Range province. The model-determined range of fault friction for the San Andreas fault, Mendocino fracture, and Queen Charlotte fault in the model was found to be .11 to .2, with the preferred value being about .17. A similar value was used by Bird in studies in Alaska and southern California (Bird, 1996; Bird and Kong, 1994).

A regional deformation model encompassing most of the western U. S. (including Cascadia) And British Columbia has been completed by Geist (1998). Geist used Bird's 'LARAMY' program (Bird, 1989) which differs from 'PLATES' in that time-stepping is used and temperature fields are computed at each time step to take into account thermal diffusion and shear strain. This program is used for study of continental scale deformation through long periods of geologic time, whereas 'PLATES' is more useful for studying neotectonic stress and fault-slip patterns. Geist (1998) found a preferred value for subduction zone traction limit of 15 MPa, within our preferred range, suggesting long term stability in this order of traction limit. He also inferred that, if values of less than 5 MPa were used, then Basin and Range extension would occur throughout the Pacific Northwest due to interior shear from Pacific-North America interaction. We observe the same effect, especially when the crust is made more silicic in rheology than our assumption for the coast ranges. Bird (1978;1996) similarly found subduction shear limits of 10-20 MPa in the western Pacific and Alaskan subduction zones. Our limits and those of Geist (1998) for strike-slip fault friction in Cascadia compare well with low values found for southern California and Alaskan strike-slip faults by Bird and Kong (1994) and Bird (1996) respectively. Low fault shear stresses are likely related to the effects of fluids and lubricating minerals in the upper, brittle parts of the fault zones (Byrne et al., 1988). Below the brittle-ductile transition, mineral plasticity can explain low shear stresses.

Model Stress Results

The best fit to stress patterns is shown for model r1_6_99n in Figure 56. The velocity boundary conditions used are indicated along the margins of the grid. Model stress directions are oriented north-south in the coast ranges of Oregon and Washington as suggested by borehole breakouts and earthquake focal mechanisms. The cold, rigid, mafic crust of the coast ranges acts as a stress guide to transmit margin-parallel shear from the Mendocino triple junction. Borehole breakouts in the northern Great Valley and near the Mendocino triple junction are oriented northeast, but our stress model does not indicate this orientation except in the region where Oregon, California, and Nevada meet. The northeast stress direction is dominant in eastern Oregon where ductile flow occurs across the region of high heat flow (Figure 54). This is partially related to the stress free boundary created by this high-heat-flow region, but primarily from loss of the coast range stress guide that channels NS compression.

In southwestern Washington (Figure 56), model stress directions turn northeastward as Juan de Fuca subduction stresses overcome the margin parallel shear. This change in direction is partially assisted by high heat flow in the southern Washington Cascades and low heat flow just north of Mt. Rainier (Figure 56). North-south compressive principal stress continues across the Olympic Peninsula in westernmost Washington. The modeled northeast trend of stress across the Cascades ends abruptly near the OWL, where the heat flow low just north of Mt. Rainier provides stronger crust and principal stress directions trend northwest. This heat flow boundary near the OWL also coincides with the P-wave velocity expression of the OWL on the 7.0 km/s isosurface of Figure 11. P-axes from focal mechanism determinations in north-central Washington and southwestern British Columbia are oriented dominantly northwest, as demonstrated by Mulder (1995). To produce this pattern of stress directions, we were forced to use boundary conditions of 3-10 mm/yr north of the OWL directed parallel to Pacific plate motion (N50oW), similar to those of northern California and northwestern Nevada. Velocities of near zero were required across the mesh boundary in eastern Oregon. We interpret that interior shear in the Basin and Range produced by Pacific plate-North America interaction is most effectively concentrated near the Mendocino triple junction and the OWL because of low heat flows in these regions. The more ductile lithosphere of eastern Oregon does not respond rapidly to the NW shear motion, but over longer geologic periods, should express the same effect, as demonstrated by Geist (1996).

The pattern of stress vectors shown in Figure 56 would change radically under alteration of two main assumptions. By reducing JDF thrust tractions to the <5 MPa value, maximum principal stress in western Oregon trends almost true north (Figure 57), and in the Cascades of Oregon and eastward, trends more northerly than for higher traction limits. This effect cause a more classical margin-sliver stress pattern (Wang, 1996) than for higher traction limits. The balance point for this transition is relatively sharp, and the pattern of Figure 57 begins suddenly at tractions of less than 5 MPa. The other assumption in our model, of relatively mafic crust, heavily influences the stress pattern in eastern Oregon and elsewhere. For instance, if the rheology parameter B=9000oK assumed for the preferred model is reduced to B=6000oK, then extension is the principal stress in most of eastern Oregon and Washington, and compressive principal stresses in southwestern Washington do not show clear a NE orientation as observed from focal mechanisms.

Pn Anistropy

An interesting comparison to the stress pattern in our preferred model comes from Pn anisotropy calculations by Hearn (1996). Pn rays travel along the crust mantle-boundary and sample the upper mantle and crust. Pn travel times are affected not only by lateral variations in crust and mantle velocity, but also by significant amounts of laterally varying seismic anisotropy. When the earth has seismic velocity anisotropy, energy traveling in the principal anisotropy direction arrives earlier than that traveling in the minor direction. The anisotropy can be caused by alignment of certain minerals, such as olivine, along the maximum stress. Sets of faults and geologic boundaries can also cause anisotropy Hearn (1996) reformulated a tomography algorithm to include lateral variations in both velocity and horizontal anisotropy. Using this tomographic method, Hearn found Pn anisotropy in the western U.S. ranged up to 7.5%. We show a subset of this data for Cascadia in Figure 58. His results compared favorably with SKS measurements (Silver and Chan, 1991; Bostock and Cassidy, 1995) of shear-wave splitting in the western U. S. The SKS method uses polarized shear-wave, core phases to measure anisotropy beneath a seismic station. Hearn interprets from his results and that of SKS splitting measurements that the fast direction of Pn anisotropy occurs in the maximum compressive stress direction and is most likely related to present-day mantle deformation. This interpretation seems to be compatible with our stress model for Cascadia, with good correlation in many of the details of our stress directions and Hearn's Pn anisotropies.

Northeast-directed compressive stress is indicated in most of Oregon from Pn anisotropies, except in the coast range, the anisotropy azimuths are north-south oriented. In southwestern Washington, anisotropy azimuths turn NE across the Cascades and southern Puget Sound, but remain north-south oriented in the coast range, just as in our deformation model. Very large values of NW oriented anisotropies start near the OWL and continue northward across the Columbia Plateau, as suggested by our stress model and focal mechanism determinations in the region north of the OWL. Very low values of anisotropy occur in high heat flow regions of east-central Oregon, suggesting a very uniform, ductile crust and upper mantle. A change from NW oriented to NE oriented anisotropy occurs in northeastern California, similar to the change in stress directions from our model. Two main areas of disagreement can be noted. The first is the large anisotropy in a NE direction just north and east of the Mendocino triple junction. This direction of maximum compression also seems to be required by borehole breakout and focal mechanism determinations (figure 53). Our stress model does not match these two data sets well, but could be corrected if we locked the southern part of the JDFP with very high traction forces, or made the Pacific plate thrust under the Mendocino triple junction in a NE direction. Both of these options may be valid. As partially tested in our modeling exercises, neither of these options would change the stress pattern significantly in the western Washington region, where we focused our attention in this study. Another area of disagreement between Hearn's anisotropies and our stress model is NW oriented anisotropies near southern Vancouver Island, and NE oriented values east of southern Vancouver Island. This results does not agree with our stress model, with SKS splitting results from Bostock and Cassidy (1995) or with NW trending focal mechanism determination of P-axes (Mulder, 1995) east of Vancouver Island. We suspect that these anisotropies may be poorly determined because of the location of this region near the corner of the tomographic grid. However, the general agreement between Hearn's (1995) Pn anisotropies and our stress model, compared to other observations, suggest that mantle stress directions are broadly similar to those reflected in the crust and are controlled by the interactions of the JDF-North America-Pacific plate system as described earlier. In addition, crustal rheology and heat flow are the main factors that modify the stress patterns caused by this interaction. Especially significant in the Pn results is the strong anisotropy observed near the OWL.

Plate Interactions

The evidence from our stress modeling and Pn anisotropies point to a major region of plate interaction shear near the OWL. The most important elements of Pacific-North America-Juan de Fuca plate interactions are summarized in Figure 59. Primary extension in the Basin and Range province is caused by the shearing of the Pacific plate against North America. The OWL is co-linear with the Queen Charlotte fault and forms part of the part of the northern exit zone for Pacific plate induced shear to the transform system. The more plastic region of eastern Oregon deforms uniformly in the central area between the projections of the two transform boundaries defined by the Blanco fracture zone and the Queen Charlotte fault. Shear exits interior North America in a complex manner toward the Queen Charlotte transform by coupling with the thick continental crust of the North Cascades. Other parts of this northern exit region encompassing the OWL may be defined by Eastern Snake River plain boundary faults to the southeast (Figure 59), and Devils Mt.-Darrington-San Juan faults to the northwest. This pattern of NW shear explains the NW trending focal mechanisms north of the OWL and could explain the large (~M7.2) crustal earthquake of 1872 in the North Cascades (Malone and Bor, 1975). Significantly, the Milton-Freewater M5.0 earthquake of July, 1936 and a M5.4 event in April, 1945 (Figure 59) were located near the OWL and probably reflect the shear concentration discussed above.

There might appear to be a problem with our velocity boundary conditions required near the OWL to reproduce NW trending stress directions in the North Cascades. As noted in a tectonic analysis by Hooper and Conrey (1989), evidence for crustal thinning on the south side of the OWL relative to the north suggests right-lateral motion. The OWL represents a fundamental crust boundary, probably a suture zone, with the north side anchored in thick continental of the North Cascades on the NW and the Idaho batholith on the SE (Hooper and Conrey, 1989). Thick Miocene basalt flows form a rigid plate, with the maximum thickness centered in a ponded region (Pasco Basin), just south of the OWL. We interpret that this shallow crust rigid plate is being pushed along like a hockey puck over the relatively thick layer of sedimentary rocks (2-5 km thick) that occur beneath the basalts (Hooper and Conrey, 1989), accentuating the NW shear caused by Basin and Range spreading. On the north side of the OWL, the thick crustal anchor of the North Cascades means that near-surface translation probably does not occur in very noticeable fashion. Instead, the NW directed shear causes slow creep to be distributed over the whole of the North Cascades. This type of strain partitioning due to JDF-North America-Pacific plate interaction means that large crustal earthquakes such as the 1872 North Cascades earthquake are more likely to occur in the North Cascades than south of the OWL. In summary, small-plate behavior in the Columbia Plateau region suggests right-lateral slip, but regional stresses north of the OWL require a region of left-lateral motion with respect to the Puget Sound area along the northeast part of the OWL, the DAR, and the DMF. This left-lateral shear can also be better visualized if the region between the OWL and the Coast Range plutonic complex of British Columbia (CRPC, Figure 59) is considered to be highly transpressive, such that the broad right-lateral shear zone induced by Pacific-NA interaction develops interior block-specific, left-lateral shear.

A southern exit zone for the JDF-North America-Pacific plate interaction probably occurs in northern California where shear along the eastern Sierra Nevada block connects to the Mendocino triple junction. This shear zone may coincide with a Lake Almanor-Tahoe seismicity zone and fault belt (Unruh et al., 1996) which is co-linear with the Blanco transform. Incipient fractures and seismicity trends suggest to Hudnut (1998) that the Blanco transform will eventually form the southern boundary of the Juan de Fuca plate. This shear province is probably a continuation of the Fish Lake-Furnace Creek fault zones and related fractures along the east side of the Sierra Nevada (Reheis and Dixon, 1996). The Brothers fault zone in central Oregon may be a shear locus midway between the northern and southern boundaries, bounding the Mesozoic-cored Blue Mountains block of northeastern Oregon (Pezzopane and Weldon, 1993). Northwest stress is concentrated near the southern and northern margins of the shear region and tends to zero in the higher temperature center, as suggested by our models and Pn anisotropy data (Figure 58).

Principal stress in the center of the Basin and Range province is northwest-oriented extension. This is clearly demonstrated by the simple 'PLATES' model of Figure 60, where input heat flow, surface topography, and crustal thicknesses were uniform. The only controlling parameters are crustal rheology and transform boundaries and their motions. The model illustrates that northwest oriented extension should be occurring within the Basin and Range irrespective of heat flows and topography. This implies that the high heat flows are related to lithospheric response to the shear effects directly, modulated by other mechanisms. The local variations in stress directions in our detailed model are caused by heat flow and interior velocity variations not accounted for in the simple, regional model. The two main deviations in the detailed model (Figure 56) from the simple stress picture of Figure 60 are the northeast trending compressive stress directions in northern California and eastern Oregon, and northwest directed principal stress in north-central Washington. Apart from these subprovince variations the regional picture of stress should look similar to that of Figure 60.

Surface Velocities and Fault Slip Rates

Horizontal nodal velocities for the mesh in the best fitting anelastic model are shown in Figure 61. These velocities represent the effects of creep over several seismic cycles, and not the interseismic elastic effects that typically are measured with GPS instrumentation. However, if there is viscous creep along the subduction thrust, or if anelastic creep in the upper plate is large, then the GPS results may reflect such effects. Significantly, the fit to GPS velocities measured at station FTS1 (Qamar et al., 1998) is excellent, both in direction and magnitude. The model results illustrate that the anelastic signal from NS compression is essential to modeling GPS azimuths in the coast ranges. Previous dislocation models (Flück et al., 1997; Qamar et al., 1998) did not incorporate this interaction and do not closely model the GPS azimuths at stations SEAT, RPT1, and FTS1. Model velocities at the two Puget Sound GPS stations (SEAT and RPT1) are relatively similar in direction to the GPS results, but too large in magnitude for the short term GPS measurements. This correspondence between the surface velocities from our anelastic model and the GPS results suggests that there may be a large anelastic component to the GPS velocities. For the northern Washington stations, NEAH and WHD1, both directions and magnitude are off much more than the other locations because the anelastic signal is much smaller.

Fault slip components from the best model are shown in Figure 62. Red arrows represent the calculated horizontal shortening across the faults and the black arrows are strike-slip components. For the modeled reverse faults, the uplift component is the same as the horizontal shortening, since we assumed all reverse faults dip at 45o. Relatively small shortening (<0.5mm/yr) occurs on the bounding thrust fault in the accretionary prism in the Olympic Mts. (OLMT), with the largest values being on the north side of the curving thrust. The west end of the Seattle fault zone has low shortening rates of <0.3 mm/yr, but the east end of the fault zone has somewhat higher rates of ~0.6 mm/yr. These rates are comparable to Quaternary shortening rates across the Seattle fault estimated by Johnson et al. (1999a) using seismic reflection images of folded Quaternary strata

The largest shortening rates in western Washington occur across the Devils Mt. fault, which is advantageously oriented to accommodate the maximum amount of combined north-south compression and northeast JDFP shortening. Shortening rates range from 0.7 to 2.3 mm/yr from west to east along the Devils Mt. fault elements and sinistral strike-slip occurs at rates of 0.8 to 1.5 mm/yr. S. Y. Johnson (oral comm., 1999) has made initial estimates of horizontal shortening across the Devils Mt. fault on north Whidbey Island and on one well profile and obtains values of ~0.20 mm/yr. Johnson has also mapped a large number of faults somewhat concordant with the Devils Mt. fault in marine seismic data. These faults suggest left-lateral transpression and shortening distributed over tens of kilometers. Geologic shortening may be distributed across the relative wide fault system, whereas our model assumes all of the shortening occurs on the single model fault. The Darrington fault has shortening rates of 1.9 mm/yr and sinistral strike-slip of 1.7 mm/yr.

The South Whidbey Island fault was modeled with very small slip rates, reflecting the fact that it occurs in the zone between the OLMT and DMF and is shorter than these two regional faults. Sinistral strike-slip along the OWL ranges from 0.6 to 2.0 mm/yr. An interesting result from the model is the relative large northeast shortening across the WRZ just south of Mt. Rainier, at about 1.1 mm/yr. The localization of this slip feature is probably related to the higher plasticity of the crust near the heat flow anomaly at Mt. Rainier (Figure 55). This slip feature correlates with a NE trending zone of earthquakes that steps over from the SHZ and across the WRZ (Figure 48). The SHZ has a small, purely dextral slip rate of 0.4 mm/yr in the center of the fault near Mt. St. Helens. The Portland Hills fault zone (PHFZ) has a dextral slip rate of 0.25 mm/yr and negligible shortening.

Stress Origins and Seismicity in the Puget Lowland

A final view of stress patterns and seismicity in the Puget Sound region is provided by the compilation in Figure 63. Stress vectors from our model are shown in this illustration, along with relocated earthquakes (red dots) at depths from 15 to 30 km in the crust. The outline of a key block of Crescent Formation (dashed black line) south of the Seattle fault was determined using several lines of evidence. The northern edge is evident in the seismic velocity, gravity, and magnetic data. The western margin is defined in Vp/Vs results and in MT models. The southern margin is constrained by gravity and magnetic data and by Vp/Vs results. The eastern margin is defined by the distinctive seismicity patterns and by the slight shift in velocity contours in vertical cross-sections (Figure 7). The smaller block of Crescent Formation north of the Seattle fault is outlined by inferences from the seismicity patterns, from magnetic and gravity data, and from the velocity model. We also show a small, separate block at the northeast corner of the larger block (near Seattle), based upon the tight NE trend of earthquakes in the basement near Seattle. Tom Brocher (USGS, pers. comm.) has noted a basement structure under the Seattle basin in deep reflection data, which may be related to this small corner block.

Seismicity (red dots) in the large block of Crescent Formation is extensive, but the region outside the block along the western, southern, and eastern margins is seismically very quiet. Slightly east of the key Crescent block, near the eastern margin of the Puget Lowland, and north of the NW end of the OWL, small earthquakes are again numerous. We have tentatively connected the structure represented by the OWL to the east end of the Seattle fault zone (queried short segment) because of the favorable alignment of topographic features and magnetic anomalies, along with the change in seismicity north of this alignment. The Seattle fault zone and the OWL were connected in our deformation models, but the resulting northwest shear was focused on the much longer DMF-DAR system which was also connected to the OWL. In reality, northwest-directed shear is probably distributed along a set of faults north of the OWL, including the Seattle, Whidbey Island, and Devils Mt. fault zones. Seismicity in pre-Tertiary crust north of the OWL is diffuse and largely concentrated where stress vectors turn from NE to NW, reflecting a large amount of tectonic work. The seismicity pattern in the key block of Crescent south of the Seattle fault is located in the lower part of the block, as indicated in the cross-sections of Figure 51 and 52. The earthquakes are more clustered than in the pre-Tertiary crust to the east and may be indicative of shear fracturing of the brittle mafic block along ENE and NNW directions (dashed lines), as would be expected from the maximum compressive stress directions.

This key block of Crescent Formation may be tightly coupled to the mafic wedge beneath and the subducting JDFP below that, therefore it picks up the NE subduction stress very effectively. It also must turn the corner around the nose of the JDFP arch (outlined in green contours, Figure 63). This geometry and coupling may help explain anomalous GPS velocities in the Puget Lowland observed by Qamar and Khazaradze (1999) that are difficult to model either in our anelastic models or their dislocation models. We have previously used a 'train-wreck' analogy to portray the complexities of the mafic block dynamics in western Washington. In this analogy, mafic blocks of Siletzia in the western Oregon and Washington Coast Range behave somewhat like a train pushing a set of boxcars. The rear, pusher engine is in the area of the Mendocino triple junction and the thick, continuous mafic rocks of the Oregon coast range couples the force to less continuous mafic crust in western Washington. Sedimentary basins between blocks of Crescent in southwestern Washington act as viscous shock absorbers.. The front end of the train is just north of Seattle and the cars have been jammed into the bend in the subduction zone and JDFP plate arch. Stress is being progressively applied to these front cars in the Seattle region. Although this analogy is useful in portraying the complexities, our current analysis focuses attention on the 'boxcar' or key block of Crescent Fm. just south of the Seattle fault. The stress dynamics seem to be mainly related to local, vertically distributed tractions between the key Crescent block, underlying mafic wedge, and the JDFP which concentrates NS compressive stress.

Larger earthquakes in the Puget Lowland region are noted by the solid white dots in Figure 63. The 1949 and 1965 earthquakes (red magnitude labels) are the two large intraslab events shown on the cross-sections of Figures 44a and 44b. The 1939 earthquake may either be an intraslab event or deep crustal earthquake; no detailed analysis has been done of this event or the large event in 1946. Malone and Bor (1979) report depths of 40 km for the 1939 and larger 1946 events. The M6.1 earthquake in 1946 is of considerable interest because of two M5 events during the same week (2/15/46). The two M5 events in February, 1946 occurred 8 hours and 6 days after the M6.1 1946 event, but are clearly not aftershocks. They may be related to structures along the southern margin of the key Crescent Formation block, since gravity data (Figure 2) indicate that the linear boundary denoted as the 'Legislature fault' is not too far from the line formed by the three earthquakes in February, 1946. The pattern formed by the large and medium earthquakes in southern Puget Sound, including events such as the 1997 Bremerton and 1995 Robinson Point earthquakes illustrates that this part of western Washington has the highest seismic energy release. The only other medium-size earthquakes are the 1945 M5.7 event near the NW end of the OWL and the 1996 Duvall M5.3 earthquake.

A M5.4 earthquake in 1954 was noted by Malone and Bor (1979) with a focal depth of 50 km, so it may be an intraslab event similar to those nearby in 1949 and 1965. Kirby (1998) has inferred that these intraslab events are related to the gabbro-eclogite phase transition. However, as discussed in a previous section, we find difficulty in this interpretation because of the evidence from the velocity structure that the 1949 and 1965 earthquakes occurred at depths where velocities are > 8.0 km/s, indicative of the events having occurred within the slab mantle where the phase transition is not applicable. We interpret that the intraslab earthquakes in southern Puget Sound are localized by complex bending strains in the subducting slab. Slab contours (green) are shown in Figure 63. The contour intervals from 70 km and deeper are dashed because they are inferred from the steep dip part of the plate, arch structure, and deep plate contours of Bostock and Vandecar (1994). However, inflections in the 60 and 70 km contours are documented in velocity cross-sections such as those in Figures 16-18. There are two primary components of curvature on the slab: (1) horizontal radiality from the arch structure and (2) shallow-to-steep curvature in a vertical plane as indicated by the velocity cross-sections. The tighter-radius inflection south (downdip) of the 1949 and 1965 intraslab events would focus bending moments in the area of these earthquakes (Chapple and Forsyth, 1979), possibly leading to normal faulting in the slab. The contact of the JDFP and southern part of the mafic wedge, as viewed in Figure 9, may cause additional bending moments in this locality and be the cause of the tighter-radius inflection. In summary, the southern Puget Sound region seems to be an area of complex and tightly-focused stresses.

1100 yr B.P. Deformation

We suggest that there may be a relationship between uplift observed in southern Puget Sound (Bucknam et al., 1992) and the hypothesized deep interplate source region. Bucknam et al. (1996) has found two broad areas of sudden uplift in southern Puget Sound which occurred around 1000-1100 years ago. The first area is along the southwestern side of the Seattle fault and a second area of about 800 km2 occurs 30-70 km south of the Seattle fault. They point out that the size and location of the second uplift area indicate that it was not caused by slip on the Seattle fault. Sherrod (1998, 1999) has found extensive evidence for land submergence in southern Puget Sound which occurred between 1140±80 and 1010±50 yr B.P. radiocarbon dates. Sherrod (1998) speculates that two or more faults were involved in producing large earthquakes near this time interval, with the 'Legislature fault' discussed earlier (Figure 44b) being a likely candidate for one of them.

Logan et al. (1998) summarized evidence for large earthquakes in the Olympic Mts.-Puget Sound region consisting of landslide-dammed lakes and submarine landslides. Their radiocarbon dates are shown in Figure 64, along with the outline of a key block of Crescent Fm. between the 'Legislature' and Seattle fault zones. The landslide-dammed lakes were caused by large rock slides, indicative of significant earthquakes (R. L. Schuster, oral. comm., 1999) and the submarine slides, which were mainly in Lake Washington, were assumed to be also related to large earthquakes. The tight time clustering of radiocarbon dates for the ~1100 yr B.P. earthquake and broad geographic distribution of these sites in southern Puget Sound and eastern Olympic Mts. is striking. We suggest that the grouping of earthquake indicators, including evidence for crustal deformation, for the 1100 yr B.P. earthquakes should be interpreted as a major single crustal deformation event, and not as a sequence of unrelated fault ruptures. Specifically, we infer that there was nearly simultaneous deformation and faulting related to the key block of Crescent Formation outlined in Figures 50 and 54 that involved several boundary faults, including the 'Legislature' fault, Seattle fault, and small faults in the eastern Olympic Mts. that may be connected in an unknown fashion.

Evidence cited earlier from Atwater (1992) for deep sand blows near the Washington coastline at the mouth of the Copalis River in the interval 900-1300 years ago correlates with the timing of earthquake deformation in Puget sound. Meyers et al. (1996) have also found earthquake-subsidence events near the mouth of the Columbia River, one of which is radiocarbon dated at 1110±60 yr B.P. These earthquake indicators near the Washington coastline correlate with the Puget Sound indicators dated at ~1100 yr B.P., suggesting a common cause. Atwater (1992) documents other sites of sand volcanoes that may be related to large earthquakes during the 900-1300 year interval along the northwest coastline of the Olympic Peninsula, again suggesting a possible large subduction zone earthquake about this time. Supporting evidence for a large Cascadia subduction earthquake during this interval also is cited by Carver et al. (1996) who interpret that a tsunami deposit near Crescent City, California dated at 670-1173 years ago was due to a single Cascadia event. We interpret that the best explanation for the widespread evidence for large earthquakes at ~1100 yr B.P. involves a subduction thrust earthquake initiated in the deep locked region under Puget Sound. Such an event would have created large ground motions in the Puget Sound region. If the earthquake involved only 100-150 km strike length of the thrust with a M~8.0 in the Puget Sound region, and possibly did not propagate to the shallowest part of the thrust offshore, it would not have generated a great tsunami similar to the one documented for the 1700 A. D. subduction zone earthquake (Satake et al., 1996). The proposed earthquake may have destabilized the frictional points in the overlying the key block of Crescent Fm., causing release of both JDFP generated stress and that arising from margin parallel shear. Deformation in the key block of Crescent Formation and surrounding blocks would have been complex, as observed in paleoseismicity results, dependent upon the block thicknesses, strength, and local state of stress. Large crustal earthquakes might have occurred over a period of several hundred years, as the stresses reoriented to the loss of the deep traction force. Remote triggering of earthquakes within the Puget Sound region may have also played a factor in the paleoseismic observations. An example of large crustal fault motion connected with a subduction thrust earthquake comes from the 1964 Alaska great earthquake (M8.4). The largest uplifts were documented along the Patton Bay fault that runs through Montague Island, with offsets of over 10 meters measured (Plafker, 1969). Although the Seattle fault is not oriented parallel to the subduction zone, release of underlying traction forces at the plate boundary might be the trigger for the key Crescent block on the south side of the Seattle fault to jump to the north (because of stored north-south compressive stress). The cause of triggering of the deep locked zone could be seismic or aseismic slip on the offshore locked zone, a buildup of high-pore pressures in the thrust, or the occurrence of a large intraslab earthquake.

Summary

We have conducted a comprehensive review of geological and geophysical models and related information in order to provide useful seismotectonic details of the western Washington region. Gravity and magnetic maps, combined with magnetotelluric models, provide significant information on the nature of crustal lithology and structure, but the use of velocity models derived from local earthquake tomography has provided the most encompassing set of data for this analysis. The P-wave velocity models and a relocated earthquake data set provide a range of tectonic information for the northern Cascadia subduction zone in western Washington. Accurate images of the subducting Juan de Fuca plate have been obtained with these results, and structures in the subduction zone, such as the arch in the oceanic plate are better understood as a result. Details of a northern buttress which forms a backstop for deformation of the plate arch and for crustal blocks moving northward have also been provided by the velocity models. The velocity models outline key blocks of Crescent Formation, dynamics of north-south compression, unmapped faults, and sedimentary basin structures. In addition, the velocity models were used to map a thick, high-velocity region of the lower crust that we relate to a mafic/ultramafic wedge. This mafic/ultramafic wedge is inferred to have developed in the rift system which also generated the Crescent Formation and other coast range basalts.

In order to better understand stress patterns in Cascadia, we have employed finite-element, thin-plate modeling methods to compute anelastic responses of the lithosphere to various boundary conditions. These boundary conditions involved assumptions about rheology, horizontal velocities, heat flow, and crust/mantle thicknesses. Our results show that the combination of a highly mafic crust and low heat flow in western Oregon and Washington cause this region to act as a stress guide for margin-parallel shear. The models easily replicate transition from NS compression to NE compression in the southern Washington Cascades. However, to model indicated NW-oriented compression in the North Cascades and northern Columbia Plateau region, we were required to use boundary motions of 3-10 mm/yr north of the Olympic-Wallowa lineament. This finding may provide a better understanding for large crustal earthquakes in the region, such as the 1872 M7.2 quake in the North Cascades.

We have reviewed evidence for the downdip extent of the seismogenic part of the subduction thrust in northern Cascadia and propose a revised model that takes advantage of new information on upper-plate lithology from our velocity models and those of Flueh et al. (1998). Seismic reflection data from Flueh et al. (1998) suggest that all accretionary prism sediments are offscraped west of the shoreline. Using our velocity models, we interpret that the deep part of the subduction thrust under the Puget Sound region has the oceanic plate in contact with mafic/ultramafic rocks in the upper plate. Since this contact is at temperatures below the brittle-ductile transition for mafic rocks, we infer that the thrust may be locked down to 40 km depth, with a possible transitional region updip. If this assumption is correct, then ground motions from a subduction thrust earthquake sourced at these depths would be much greater in the Puget Sound region that one sourced from an offshore part of the thrust, as preferred by Hyndman and Wang (1993, 1995). Our stress models and dislocation models from others, based on new GPS data, may require a deep locked region of the subduction zone.

Crustal seismicity in the southern Puget Sound region appears to be controlled by a key block of Crescent Formation occurring just south of the Seattle fault. Earthquakes in this block are concentrated near its base and appear to reflect shear fracturing compatible with the stress directions derived from our finite-element models. We suggest that the complex pattern of deformation in southern Puget Sound which has been dated at about 1100 yr B.P. is related to destabilization of stress in this block and nearby blocks, possibly triggered by a thrust earthquake in the deep source zone. The concept of a subduction earthquake during this interval is supported by paleoseismic evidence from the coastal area of Washington. Such an event probably would have ruptured only the 150-200 km strike length of the deep thrust in the Puget Lowland region and would not have been as large (~M8) as the event of 1700 AD (~M9) which apparently ruptured the entire offshore region of the subduction thrust. The other major hazard in the Puget Sound region is from large intraslab earthquakes like the M7.1 April 1949 event. This event and other intraslab earthquakes in southern Puget Sound are interpreted to be related primarly to complex curvatures which focus bending moment stresses within the subducting slab. If this inference is correct, then future large intraslab earthquakes are to be expected in southern Puget Sound.

Acknowledgments

The authors with to thank the staff of the Pacific Northwest Seismograph Network for providing the arrival time data used in the tomographic inversion. The authors also wish to thank Tom Hearn for supplying results from his Pn anisotropy studies and Dave Blackwell for providing heat flow data from his files. Bradley Hacker provided early results from his metamorphic facies and velocity model for Cascadia, which proved to be very useful. Many of the figures in this report were made with Generic Mapping Tools (GMT) from Wessel and Smith (1991).


References

Bird, Peter, 1996, Computer simulations of Alaskan neotectonics: Tectonics, 15, 225-236.

Kirby, S. H., 1983, Rheology of the lithosphere: Reviews of Geophysics, 21, 1458-1487.

Peacock, S. M., 1987, Thermal effects of metamorphic fluids in subduction zones: Geology, 15, 1057-1060.

Peacock, S. M., 1996, Thermal and petrologic structure of subduction zones: in Subduction Top to Bottom, Geophysical Monograph 96, AGU, 119-133.

Peacock, S. M., 1996, Large-scale hydration of the lithosphere above subducting slabs: Chemical Geology, 108, 49-59.

Raisz, E., 1945, The Olympic-Wallowa lineament: American Journal of Science, 243-A, 479-485.

Schultz, Andreas P., and Crosson, Robert S., 1996, Seismic velocity structure across the central Washington Cascade Range from refraction interpretation with earthquake sources: Journal of Geophysical Research, B, Solid Earth and Planets, 101 (12), p. 27,899-27,915.

Vrolijk, P., 1990, On the mechanical role of smectite in subduction zones: Geology, 18, 703-707.

Wannamaker, P. E., Hohman, G. W., and Ward, S., 1984, Magnetotelluric responses of three-dimensional bodies in layered earths: Geophysics, 49, 1517-1533.


Appendix A


Resolution Studies For Western Washington Tomography

Details of the inversion method used in this study are contained in Benz et al. (1996), but in this appendix we discuss specific resolution tests for our model. The resolving power of the method used to obtain the velocity structure from the arrival time information was tested empirically by mapping the ray coverage (defined here as the cumulative ray-length inside each cell, expressed in km), and by the reconstruction of synthetic checkerboard models, a method commonly used in earthquake tomography. The checkerboard test involved constructing a synthetic model, consisting of checkers of 40x40x10 km with alternating velocity perturbations ±10% over a homogeneous half space of 6.5 km/s. Synthetic travel times were calculated for the checkerboard model with the same distribution of sources and receivers used in the inversion. This synthetic data set was inverted using the same 8x8x2 km grid and identical smoothing weights to those employed in the actual inversion. The smoothing constraints are described in (Benz et al., 1996). The resolution tests provide
confidence in our geologic interpretation of this paper regarding important structural features.
Typical results for the tests are indicated in Figure A-1 for an east-west vertical slice for row 23 (see Figure 7d as well). Below the P-wave velocity model (Figure A-1a) used as part of the geologic interpretation, we show the ray-path coverage (Figure A-1b). The resolution should be best in the parts of the model where ray coverage is greatest. However, more important than total ray length through a cell is the azimuthal sampling of the cell. For some cells, the total number of ray paths may be large, but they may all be traveling in a similar direction (e.g. toward one receiver). Better velocity resolution is obtained when the rays have a wide distribution of azimuths (to various receivers). This is the case for upper crust below the Puget lowlands down to 30 km, due to the large number of earthquakes and their distribution. Another well-sampled region is the top of the subducting Juan de Fuca plate, due both to the presence of intraplate seismicity and the large number of Pn rays travelling through the refracting surface with a wide variation in azimuth. This effect can be observed in the coverage plot (Figure A-1b).
The synthetic checkerboard model is shown in Figure A-1c, and the reconstruction in Figure A-1d. An excellent recovery of the shape and amplitude of the anomalies can be observed below the Puget lowlands down to 30 km, coincident with the region of high coverage shown in Fig A-1b. In the rest of the model, the shape of the anomalies is in general well recovered, except for some regions in the lower crust, where smearing due to rays traveling in similar directions is observed. The synthetic checkerboard model does not contain a refractor or a surface with a large velocity gradient (comparable to the top of the Juan de Fuca plate). Therefore there are no head waves travelling through that interface, and the resolution derived from the checkerboard tests for the Juan de Fuca plate underestimates the actual model resolution.
Checkerboard tests show that velocity anomalies with wavelengths on the order of 40 km in horizontal and 10 km in vertical are well recovered for most of the model volume. In our discussion we have only considered large velocity anomalies, located in well sampled regions of the model, providing confidence on the interpretations derived from our results.
We have also included resolution results for four horizontal slices with test model depths starting at 0, 10, 20, and 30 km respectively. These results show that horizontal resolution of the deep crustal features that we use in our geologic interpretation is very good. In addition to these simple checkerboard tests, we are currently developing a more meaningful resolution study that utilized a subducting plate and upper plate velocity structure simplified from the inferred geology.


Figure Captions

Figure 1--Geologic and tectonic index map of northern Cascadia. Geology adapted from Snavely and Wells (1996). Secondary thrust faults (dotted) in the Olympic Mts. adapted from Brandon et al. (1998). Cross pattern is mafic Coast Range rocks of Siletz River, Crescent, and Tillamook Volcanics (informally known as the Siletzia terrane after Irving, 1979). Closed area of Siletzia in southwestern Olympic Mts. (just below M in Mountains) inferred from gravity data, magnetotelluric models, and tomographic velocity results. Dotted line in western Oregon and Washington represents probable extent of Siletzia in the subsurface as determined from gravity and magnetic data. Asterisks are Cascade volcanoes: B=Mt. Baker; G=Glacier Peak; R=Mt. Rainier; SH=Mt. St. Helens; A=Mt. Adams; H=Mt. Hood; J=Mt. Jefferson. Key faults and slip zones in the region are: OWL=Olympic-Wallowa lineament; PHFZ=Portland Hills Fault zone; SHZ=St. Helens zone; WRZ=Western Rainier zone; SFZ=Seattle fault zone; SWF=South Whidbey Island fault; DMF=Devils Mt. fault zone; DAR=Darrington fault zone. The region covered by our tomographic velocity model is indicated by the rectangle.



Figure 2--This image consists of colored, isostatic gravity values superimposed on grayscale, aeromagnetic total field values (from Finn et al., 1998). Scale at right is in milligals. High values (reds) are most commonly related to thick mafic rocks. White dashed lines are locations of magnetotelluric profiles. Solid white lines are regional faults. SFZ=Seattle fault zone;WIF=South Whidbey Island fault; LFZ='Legislature' fault zone; SHZ=St. Helens fault zone; WRZ=western Rainier zone; PHFZ=Portland Hills fault zone; OWL=Olympic-Wallowa lineament; DMF=Devils Mt. fault zone; DAR=Darrington fault zone; OLMT=Olympic Mt. thrust (informal to this report); LRT=Leech River thrust; SB=Seattle basin; EB=Everett basin; BB=Bellingham basin: GB=Grays Harbor basin; TB=Tacoma basin.

Figure 3-Two-dimensional models of magnetotelluric profiles. Locations of the profiles are shown in Figure 2. Transverse magnetic mode data were used for the models.

Figure 4(a)--Crustal earthquakes with magnitudes >M4.0 and depths <30 km in northern Cascadia from National Earthquake Information Center (NEIC) database. Selected fault zones indicated as in Figure 1. OWL=Olympic-Wallowa lineament of Raisz (1945). Published focal mechanisms for earthquakes larger than M6.0 in the region shown are from Hyndman et al. (1979), Baker and Langston (1987), Langston and Blum (1977), and Ludwin et al., (1992). There is no focal mechanism determination for the large 1874 earthquake in northern Washington. Contours of depth (in km) to subducting Juan de Fuca plate mantle determined from velocity modeling in this report are shown as dotted lines.

Figure 4(b)--Earthquakes with magnitude >M4.0 and at depths from >30 km, from NEIC data base. Focal mechanisms and fault abbreviations as in Figure 2(a). Two deep intraslab earthquake in Puget Sound are also shown as solid dots, with no focal mechanisms available.

Figure 5--Model grid coverage for tomographic inversion. Triangles are seismic stations of the Pacific Northwest Seismic Network used in our tomographic study. Circles are earthquakes used in the inversion. Open triangles are volcanoes. The model grid used for P-velocity inversions is shown by the light grid, with rows and column intervals labelled. Rows and columns shown as cross-sections in Figures 6 and 7 are indicated by the bolder lines. These lines are indicated in the center of the cells, used for calculation of latitude and longitude for the cross-sections. The location of the four intraslab earthquakes in southern Puget Sound are shown to allow close registration with the velocity model.

Figure 6--Grid outlines for regional tomographic model and for detailed model of Mt. Rainier. Solid black rectangle is the grid for western Washington model (shown in Figure 5). Large dotted rectangle is the grid outline for a more regional model and red rectangle is the grid outline for a detailed model of The Mt. Rainier region.

Figure 7--Vertical cross-sections of P-wave velocity model along west-east direction. Row numbers (from Figure 5) are at the lower left of each panel. Velocity scale indicated on the right side of figure. Vertical scales in depth (kilometers) below sea level and horizontal scale in longitude. Earthquakes (dots) plotted are those relocated with the velocity model. Black dots are earthquake above 30 km depth and white dots are those at or below 30 km in depth. Only sampled cells (those traversed by rays) in the model grid are assigned colors from the velocity scale, but contours extend to edge of model, and the starting model can be deduced from the values at the extreme left or right edges of the plot. Areas enclosed by white, dashed circles (c, d, and e) are earthquake patterns inferred to be associated with the terrane boundary discussed in the text. Thin white line is inferred position of JDFP mantle from the velocity model and seismicity.

Figure 8--Vertical cross-sections of P-wave velocity model along column (north-south) direction. Features and construction of sections similar to that of Figure 7. Column numbers from the model grid of Figure 5 noted as ‘ncol=' below the sections. The arch in the subducting JDFP as defined by the top of JDFP mantle is denoted by the curving white lines.

Figure 9--Isosurface of P-wave velocity model for value of 7.5 km/s. View is to the northwest. In all of the 3D plots of the following figures, the X direction is north, Y is east, and Z is down.

Figure 10--Overhead view of 7.5 km/s isosurface. Red lines are streamlines, or directions of maximum dip, calculated by the visualization program.

Figure 11--Isosurface of P-wave velocities from the large, regional model (Figure 6) for value of 7.0 km/s. View is to the northwest. Blue line is outline of Washington State. Dashed line bounds the mafic wedge.

Figure 12--Multiple vertical-slice display of cross-sections for model rows 14-26. Sedimentary basins (SB=Seattle basin; EB=Everett basin; BB=Bellingham basin) in western Washington and SWCC (southern Washington Cascades conductor) are noted in the plot.

Figure 13--Map of elevation of base of Quaternary from industry seismic reflection and USGS marine reflection data. Unpublished results from Sam Johnson, USGS, 1999. X and Y axes are in meters.

Figure 14--Row 14 velocity section. See map of Figure 5 for location.
Figure 15-View of pluton near Mt. Rainier from the detailed model of Mt. Rainier area (grid outline in Figure 6).
Figure 16--Row 16 velocity section. See map of Figure 5 for location.
Figure 17--Row 18 velocity section. See map of Figure 5 for location.
Figure 18--Row 20 velocity section. See map of Figure 5 for location.
Figure 19--Row 22 velocity section. See map of Figure 5 for location.
Figure 20--Row 23 velocity section. See map of Figure 5 for location.
Figure 21--Row 24 velocity section. See map of Figure 5 for location.
Figure 22--Row 25 velocity section. See map of Figure 5 for location.
Figure 23--Row 26 velocity section. See map of Figure 5 for location.
Figure 24--Row 27 velocity section. See map of Figure 5 for location.
Figure 25--Row 28 velocity section. See map of Figure 5 for location.
Figure 26--Row 29 velocity section. See map of Figure 5 for location.
Figure 27--Row 30 velocity section. See map of Figure 5 for location.
Figure 28--Row 32 velocity section. See map of Figure 5 for location.
Figure 29--Row 33 velocity section. See map of Figure 5 for location.
Figure 30--Row 34 velocity section. See map of Figure 5 for location.
Figure 31--Column 10 velocity section. See map of Figure 5 for location.
Figure 32--Column 12 velocity section. See map of Figure 5 for location.
Figure 33--Column 14 velocity section. See map of Figure 5 for location.
Figure 34--Column 16 velocity section. See map of Figure 5 for location.
Figure 35--Column 18 velocity section. See map of Figure 5 for location.
Figure 36--Column 20 velocity section. See map of Figure 5 for location.
Figure 37--Column 21 velocity section. See map of Figure 5 for location.
Figure 38--Column 22 velocity section. See map of Figure 5 for location.
Figure 39--Column 23 velocity section. See map of Figure 5 for location.
Figure 40--Column 24 velocity section. See map of Figure 5 for location.
Figure 41--Column 25 velocity section. See map of Figure 5 for location.
Figure 42--Column 26 velocity section. See map of Figure 5 for location.

Figure 43--Combined results from P-wave and S-wave tomographic models of western Washington viewed as Vp/Vs values. (a) P-wave velocity section for model row 23; (b) S-wave velocity section for model row 23; (c)Vp/Vs ratios along model row 23; ( d) north-south Vp/Vs cross-section at model column 14.

Figure 44(a)-Geologic interpretation of east-west velocity model through the Seattle region, coresponding to profile AA' of Figure 4a and based on row 23 and nearby velocity sections. JdFc=Juan de Fuca crust; JdFm=Juan de Fuca mantle. Estimated error bars for isotherms are indicated near intersection of isotherms and JDFP. Red dots are earthquakes. Frictional states along the subduction thrust are indicated by line patterns in the legend (lower left in figure). Bold arrows denote upflow of fluids along the subduction thrust. Offshore geology constructed from seismic reflection and refraction results of Flueh et al. (1998). (b) Geologic interpretation of north-south velocity model through the Puget Sound region, east of Seattle, using velocity cross-sections at columns 16 and 18 as a primary guide. Black dots and white dots are earthquakes. Focal mechanisms for six important earthquakes are shown (table below figure). Vertical velocities are from Holdahl et al. (1989).

Figure 45--Selected focal mechanisms for some of the largest earthquakes in Cascadia. Also shown are solid white circles indicating locations of 1936 Milton-Freewater M5.4 and April 29, 1945 M5.0 earthquake in the Cascades which may be associated with the OWL. Small numbers by JDF plate earthquakes are hypocentral depths. Large numbers are occurrence years of events when noted.

Figure 46-This map shows all of the earthquakes (from the Pacific Northwest Seismic Network) for 1970-1995, including those relocated using our tomographic velocity model. The colors represent depth ranges for the earthquakes with red : 0-10 km, green:10-20 km, blue:20-30 km, black: 30-60 km. PHFZ=Portland Hills fault zone, SHZ=St. Helens zone, WRZ=western Rainier zone, SFZ=Seattle fault zone, SWF=South Whidbey Island fault, LRT=Leech River thrust, DMF=Devils Mt. fault zone, DAR=Darrington fault zone, OWL=Olympic-Wallowa lineament.

Figure 47--This image shows earthquakes occuring at depth from 0 to 10 km, with the blue dots representing events relocated with the tomographic velocity model and red dots are events as located in the PNSN (Pacific Northwest Seismic Network) catalog.

Figure 48--Earthquakes from 10-20 km depth are shown in this plot, with blue dots representing relocated earthquakes and red dots being the relocated events.

Figure 49--Earthquakes from depths of 20-30 km shown on this plot. Blue dots are relocated earthquakes and red dots are unrelocated events.

Figure 50--Earthquakes from depth of 30-60 km. Blue dots are relocated earthquakes and red dots are unrelocated events.

Figure 51--East-west sections of seismicity at rows 10-32, indexed by two. Black circles are earthquakes above 30 km depth and white circles are earthquakes below =>30 km depth. Contours from the velocity model are plotted, and selected velocity cell values are shaded. The latitude of each column and the row number are shown below the plots.

Figure 52--South-north sections of seismicity at columns 8-30, indexed by two. Black circles are earthquakes above 30 km depth and white circles are earthquakes below =>30 km depth. Contours from the velocity model are plotted, and selected velocity cell values are shaded. The column number and longitude for each section is shown below the plot.

Figure 53--Stress and deformation indicators. Borehole breakout direction and some focal mechanisms were taken from Zoback and Zoback (1991). Remaining focal mechanisms are from Mulder (1995) and Pacific Northwest Seismic Network (Qamar and Ludwin, 1992). Dikes in the Mt. Rainier are adapted from T. Sisson (1998, writ. comm.); Fold axis directions inferred from Goldfinger et al. (1991) and McCrory (1996,1997). GPS velocities from Qamar et al. (1998).
Figure 54--Heat flow values used in models. Data provided by Dave Blackwell, SMU, from data base used to generate heat flow map of North America (Blackwell and Steele, 1992)

Figure 55--Finite-element grid used for 'plates' program to model stress, strain, and velocity parameters. Dip of faults are coded as in the legend below the grid. Faults modeled are OLMT (Olympic Mt. thrust), DMF-DAR (Devils Mt.-Darrington fault zone), WIF (Whidbey Island fault), SFZ (Seattle fault zone), OWL (Olympic-Wallowa lineament), WRZ (western Rainier zone), SHZ (St. Helens zone), GRZ (Goat Rocks zone from Stanley et al., 1996), PHFZ (Portland Hills fault zone).

Figure 56--Best fitting stress vectors. Red arrows represent compressive stress and blue arrows (where shown) are extensional stress. Bold black arrows are velocity boundary conditions on edges of model grid. Yellow dotted line is edge of model grid.

Figure 57--Stress vectors for model with low traction forces on subduction zone. Red arrows represent compressive stress and blue arrows (where shown) are extensional stress. Bold black arrows are velocity boundary conditions on edges of model grid. Yellow dotted line is edge of model grid.
Figure 58--Pn anisotropy directions from Hearn, 1996. Anisotropy amplitudes range from 0.0 to 0.34 km/s.

Figure 59--Plate interactions in western North America. OWL=Olympic-Wallowa lineament, SJF=San Juan fault, DMF=Devils Mt. fault, BFZ=Brothers fault zone, THAG=Tahoe-Almanor graben, MTJ=Mendocino triple junction, FLFV=Fish Lake-Furnace Valley fault zone, BLF=Blanco fracture zone, QCF=Queen Charlotte fault. Three large crustal earthquakes discussed in the text are noted by the white circles: 1872 M7.2, July 1936 M6.1, and April, 1945 M5.7 (Malone and Bor, 1979).

Figure 60--Finite-element stress model of western U. S. using the 'PLATES' program. Blue lines are indicators of maximum horizontal stress as compressive stress and red lines are where maximum horizontal stress is extension.

Figure 61--Best fitting surface velocities from 'plates' models. Velocites in mm/yr (see scale). GPS stations (triangles) and velocities from Qamar et al. (1998). Blue lines are faults noted in Figure 55.

Figure 62--Fault slip rates from 'plates' model. Slip rates are in mm/yr. Horizontal shortening noted by red arrows and strike-slip motions by black arrows. Large red arrow is motion scale (10 mm/yr).

Figure 63--Earthquakes at depths from 15-30 km and outline of key block of Crescent Formation. (dashed line). Red arrows are maximum horizontal stress (compression) directions from our finite-element models.

Figure 64--Evidence for ~1100 yr B.P. earthquake from radiocarbon data. From Logan et al. (1998), Bucknam et al. (1992), and Karlin and Abella (1992). Dashed line is outline of key block of Crescent Formation. Red arrows are maximum horizontal stress (compression) directions from our anelastic models.


Figure A1-(a) velocity slice for row 23 (as in Figure 5c); (b) cummulative ray path length in each model cell for actual inversion;(c) checkerboard model for resolution tests, velocity contrasts are ±10% and 50x50x10 km in dimension;(d) resolution results from checkerboard test, in percent P-wave perturbation.
USA.gov logo