Herb Dragert, Kelin Wang, Thomas S. James

PDF version [PDF, 1.4 Mb, viewer]

This article was originally published in Science magazine, Vol 292, May 2001.

Continuous Global Positioning System sites in southwestern British Columbia, Canada, and northwestern Washington state, USA, have been moving landward as a result of the locked state of the Cascadia subduction fault offshore. In the summer of 1999, a cluster of seven sites briefly reversed their direction of motion. No seismicity was associated with this event. The sudden displacements are best explained by ~2 centimeters of aseismic slip over a 50-kilometer- by-300-kilometer area on the subduction interface downdip from the seismogenic zone, a rupture equivalent to an earthquake of moment magnitude 6.7. This provides evidence that slip of the hotter, plastic part of the subduction interface, and hence stress loading of the megathrust earthquake zone, can occur in discrete pulses.

Great thrust earthquakes [moment magnitude (M_w > 8] repeatedly rupture the shallow (<25 km) portion of the Cascadia subduction interface (1, 2) where the oceanic Juan de Fuca plate descends beneath the North America plate (Fig. 1). Geodetic measurements over the past decade at sites on the Cascadia margin have confirmed that this seismogenic zone of the subduction fault is currently locked (3-9). Continuous motion of the converging plates produces tectonic loading of the locked segment, eventually leading to earthquake rupture. Downdip from the seismogenic zone, temperature-controlled rheology and friction allow smoother slip without producing earthquakes (10). If there is instability in the deep segment caused by time-varying rheology or friction, the aseismic slip may be episodic and could at times trigger an earthquake in the updip seismogenic zone (11-13). The mode of the aseismic slip below the seismogenic zone of subduction faults between great earthquakes has not been observationally constrained, although theory suggests that plastic instabilities in the pressure-temperature environment of the upper mantle and lower crust may give rise to transient enhancement of slip rates (14). Data from a contiguous set of seven continuous Global Positioning System (GPS) sites have now provided evidence for the occurrence of sudden aseismic slip over a large area of the deeper Cascadia subduction interface.

Analyses of GPS data are routinely carried out for the 14 continuous GPS sites identified in Fig. 1. These sites were established specifically for the study of crustal motions (15). Changes in the latitude, longitude, and height of sites relative to the GPS site at Penticton (DRAO) are estimated from daily data (16).

Fig. 1. Location of continuous GPS sites that are included in the routine analysis of GPS data carried out at the Geological Survey of Canada (GSC). Sites in Canada are operated and maintained by the GSC; U.S. sites, which form part of the PANGA (Pacific Northwest Geodetic Array) network, are operated by a consortium of university and government agencies. Bold (red) arrows show displacements (with respect to DRAO) due to the slip event. Error ellipses are double the 95% confidence limits derived from the formal regression errors of Table 1. Thin (black) arrows show 3- to 6-year average GPS motions with respect to DRAO (7). The two dashed lines show the nominal downdip limits of the locked and transition zones from the model of Flück et al. (20). Inset shows the approximate time interval of the transient signal at each site along a northwest-striking line. larger image [GIF, 65.8 kb, 1112 X 856, notice]

Total displacements are estimated by assuming a single step at the midpoint of the transient signal. Columns marked "Start" and "End" give day of year (DOY) 1999 for the transient signal's start and end dates estimated from graphs; dN, Nsig and dE, Esig are the estimated displacements and their errors for the North and East components, respectively.

To reduce further day-to-day scatter that is common to all sites of the network, daily residuals in individual components are computed for each site by removing linear trends, annual signals, and steps due to antenna setup changes and the slip event identified in this study, all estimated simultaneously by least-mean-squares regression. These daily residuals are averaged for all network sites, and these averages are subtracted from the raw time series (Fig. 2). The application of this regional filter reduces the means of the daily (rms) scatter in the north, east, and up components from 1.3, 1.4, and 4.7 mm to 0.8, 0.8, and 3.1 mm, respectively.

Fig. 2. Daily changes in longitude at ALBH with respect to DRAO. shows the raw time series from the Bernese 4.2 solutions. Linear trend (assumed constant) and the magnitude of the transient displacement were estimated from these (and the commensurate latitude) data. shows the filtered time series where network-averaged daily residuals have been removed. shows an expanded view of the transient displacement where five-point triangular smoothing has been carried out. Date of occurrence and the duration of the transient (shaded area) were estimated from these smoothed longitudinal data. Dotted red line shows linear trend. larger image [GIF, 18.6 kb, 672 X 806, notice]

Processing data with this precision has resulted in the identification of an unexpected episode of displacements at seven contiguous sites (Table 1 and Fig. 3). Total horizontal displacements, estimated from regression, ranged from 2 to 4 mm. Estimates of the time span for the displacements at individual sites ranged from 6 to 15 days (Table 1 and Fig. 2). The longer term northeastward motion of these sites, largest at outer coastal sites, is consistent with convergence between the Juan de Fuca and North America plates with no slip at the locked part of the subduction boundary that underlies the offshore continental slope. The newly detected transient motion is in the opposite direction of this northeastward motion (Fig. 1), and displacements are largest at ALBH and PGC5, sites that are located more than 100 km landward of the locked zone. The displacements attenuate rapidly to the east and south and less rapidly to the west and northwest. The time of their occurrence varies systematically, being earliest in the southeast and about 35 days later in the northwest region of detection, indicating a signal propagation parallel to the strike of the subducting slab at an equivalent speed of roughly 6 km per day.

Fig. 3. Filtered time series of daily variations in relative longitude with respect to DRAO Annual signals have also been removed. Red lines show the best-fitting linear trends, which are assumed constant before and after the transient. The vertical green bars indicate the earliest date of detection (day 230) of the transient. larger image [GIF, 42.3 kb, 1266 X 880, notice]

The fact that the transient motion is limited to a subset of contiguous sites, shows a systematic delay in arrival at these sites, and is not correlated with tropospheric parameters, rules out its generation as an artifact of data processing. Furthermore, the source of this transient deformation is unlikely to be within the continental crust. The background crustal stress field in the Cascadia forearc is characterized by compression parallel to the plate margin as a result of a secular northward motion of the forearc sliver (17) and/or oblique subduction (18). The motion of crustal faults in this stress environment should produce predominantly margin parallel displacements, but the observed displacements are nearly perpendicular to the plate margin. The direction of the displacements and the along-strike propagation lead us to conclude that slip occurred on the subduction interface. No displacement was detected at PABH, and the reduced, more northerly directed displacement at SEAT suggests that the slip zone pinches out just south of Seattle. No transients were observed at HOLB, located at the northern tip of Vancouver Island. The slip zone may end at the Nootka Fault, which underlies central Vancouver Island and marks the northern edge of the subducting Juan de Fuca plate (Fig. 1).

Fig. 4. Three-dimensional model of slip on the subduction interface. Dashed lines are depth contours of the interface. Slip direction is set constant at 235°, the direction of motion of the North America plate with respect to the Juan de Fuca plate. Dark shading indicates the plate interface area with full (2.1 cm) slip; lighter shading indicates area where slip tapers linearly from 2.1 cm to 0 cm updip. Panels, marked by the day of year 1999, show the total area of slip on the interface in three time slices and the commensurate evolution of the surface displacement vectors [broad ( yellow) 5 model; thin (red) with error ellipses 5 observed]. Day 240 is within the time interval of the GPS transient at PGC5 and ALBH, and their observed displacement vectors have been scaled, assuming a linear increase of the displacement with time. larger image [GIF, 69.1 kb, 660 X 1466, notice]

We use a three-dimensional model of faulting in an elastic half-space to constrain the geometry of the slip zone and to fit the observed displacements. The use of a purely elastic medium is justified because of the short duration of the transient event relative to the viscoelastic time constant constrained by postglacial rebound analysis (19). The geometry of the subducting slab compiled by Flu¨ck et al. (20) is used, and the rupture is allowed to propagate along the curved upper surface of the slab. The displacement at each GPS site is calculated by numerically integrating the contribution from every point on the entire fault surface using the Okada point-source solution (21). The small number of GPS sites and limited time resolution of the transient signal do not warrant a temporally continuous rupture model. Consequently, the evolution of the slip event is presented in three snapshots (Fig. 4). For each time step, the geometry of the slip area and the slip distribution, which are kept as simple as possible, are adjusted to fit the GPS observations.

The best estimate of the slip distribution along the plate interface (Fig. 4) consists of a full slip area with 2.1 cm of slip (dark shading) in the plate-convergence direction and an area where slip has been tapered linearly from the full 2.1 cm to 0 cm updip (light shading). This area of slip is downdip of the nominal locked and transition zones (Fig. 1). The 2.1-cm slip is equivalent to about half a year of plate convergence (~4 cm/year) and generates a maximum surface strain signal of the order of 20 X 10^-9 to 30 X 10^-9 (SWNE extension to the east and SW-NE compression to the west of ALBH). The downdip boundary of the slip zone is relatively sharp, as required by the rapid landward decrease in the observed GPS site displacements, whereas the updip boundary requires a more gradual transition from full slip to no slip. Ten days after its initial detection, the rupture is limited to an area south of Vancouver Island, with the full slip zone roughly bounded by the 30- and 40-km depth contours of the plate interface. The rupture then propagates (or steps) along strike to the northwest with a gradual decrease in width. Assuming a rigidity of 40 GPa, the total cumulative moment of the rupture is 1.35 X 10^-9 N.m, equivalent to an earthquake of M_w = 6.7. As the rupture propagates, the downdip limit of the slip zone is kept constant at a depth of about 42 km, where the temperature along the plate interface is about 500°C (22).

This slip model replicates the observations well. Refinements such as heterogeneous slip directions and magnitudes are possible, although they will not alter the main features of, nor the underlying physical processes reflected by, the model of Fig. 4. Observations of vertical displacements were not used in the derivation of the slip model because of their larger (3 to 4 mm) uncertainties, but they were generally consistent with the model predictions, which ranged from -3 to +3 mm.

No obvious seismic trigger could be identified for the silent rupture. In the 5-day period marking the beginning of the reversal of motion at SEAT, SEDR, and ALBH, the largest earthquake observed in the southern part of the model zone had a surface wave magnitude (Ms) of 2.8. During the slip, background seismicity was typical for both the overriding and downgoing plates, totaling about 80 regional events, none of which exceeded Ms 5 3.0. A moderate in-slab earthquake (M_s = 5.5, depth ~41 km) occurred 103 km WSW of Seattle 45 days before the slip event, and the long time separation makes a direct causal connection unlikely. The suggestion is that the convergent plate motion across the deeper plate interface varies with time. There can be longer periods ( possibly years) of slower sliding, punctuated by brief periods ( possibly days) of more rapid slip, brought on not by seismic triggers but by sudden changes in rheology or friction. Such episodic creep events may be common and characteristic of deep slab interface dynamics, but detection of the subtle signals with GPS demands densely-spaced continuous monitoring using stable monuments, combined with exacting data processing.

The time-varying coupling of plates across the deeper subduction interface can cause average annual motion of GPS sites to fluctuate with time. For example, the average velocity of ALBH relative to DRAO from 1993 to 1998 is 4.7 +/- 0.4 mm/year at N57.7°E +/- 3.0° (7), whereas the velocity estimate based on the linear trend observed before and after the slip displacement (Fig. 2) is in the same direction but almost twice as large (9.2 +/- 0.8 mm/year). The increase is consistent with stress accumulation on the deeper interface augmenting the deformation velocity due to the offshore locked thrust zone. This time-varying behavior must be considered when inferring crustal deformation from shorter spans of GPS data.

In other seismically active areas of the world, continuous strainmeter and GPS monitoring have detected a number of slow earthquakes occurring with slip rates over a range of time scales (23). Many observations suggest that the brittle, seismogenic part of a fault can slip aseismically (24-26), whereas evidence for slip events on the deeper part of faults in the plastic regime is more limited. A slip event in the western end of the Nankai subduction zone lasted for about a year in 1997 (23, 27), and the event occurred along the plate interface at about 40 km depth, similar to the Cascadia event. Extrapolating from a thermal model for central Nankai (28), the Nankai event likely occurred where the plate interface temperature is about 400° to 500°C and, therefore, is also downdip from the seismogenic zone. The 1-year duration of the Nankai event and the much shorter Cascadia event show that the deeper part of subduction faults can slip aseismically at a variety of time scales.

Silent deep-slip events could play a key role in the cumulative stress loading of the shallower seismogenic zone, each event bringing the locked zone closer to failure. It is conceivable that one of these slip events may keep propagating updip and evolve into a trigger mechanism for a great subduction thrust earthquake. The process of deep slip leading to a thrust earthquake is considered responsible for the 1960 (M_w = 9.5) Chilean earthquake (12) and the 1944 and 1946 (both M_w = 8.2) Nankai Trough earthquakes (13), and suggests that enhanced seismic hazard may accompany silent slip events. Denser continuous GPS networks spanning Puget Sound and Georgia Strait could characterize the timing and spatial distribution of episodic silent slip, thereby allowing realtime monitoring of the seismic potential of the subduction megathrust.

1. B. F. Atwater et al., Earthquake Spectra 11, 1 (1996).
2. K. Satake, K. Shimazaki, Y. Tsuji, K. Ueda, Nature 378, 246 (1996).
3. H. Dragert, R. D. Hyndman, Geophys. Res. Lett. 22, 755 (1995).
4. G. Khazaradze, A. Qamar, H. Dragert, Geophys. Res. Lett. 26, 3153 (1999).
5. M. M. Miller et al., Tectonics 20, 161 (2001).
6. H. Dragert, R. D. Hyndman, G. C. Rogers, K. Wang, J. Geophys. Res. 99, 653 (1994).
7. J. A. Henton, thesis, University of Victoria, British Columbia, Canada (2000).
8. R. M. McCaffrey, M. D. Long, C. Goldfinger, P. C. Zwick, Geophys. Res. Lett. 27, 3117 (2000).
9. M. H. Murray, M. Lisowski, Geophys. Res. Lett. 27, 3631 (2000).
10. R. D. Hyndman, K. Wang, J. Geophys. Res. 98, 2039 (1993).
11. W. Thatcher, Nature 299, 12 (1982).
12. A. T. Linde, P. G. Silver, Geophys. Res. Lett. 16, 1305 (1989).
13. A. T. Linde, I. S. Sacks, M. T. Gladwin, M. J. S. Johnston, P. G. Silver, Eos 79, F600 (1998).
14. G. Ranalli, H. H. Schloessin, Geophys. Monogr. Am. Geophys. Union 49, 55 (1989).

15. Twelve of the 14 sites analyzed use stable, geodeticquality monuments, either concrete piers anchored directly into bedrock, or deeply anchored drilled-and braced monuments (29). Antennas at SEAT and LIND are mounted on roofs of large buildings, but their data show noise characteristics similar to the other sites. All sites have well-documented histories of site activity, allowing the recognition of potential disruption of their time series due to instrumental changes.

16. The Bernese GPS Software Version 4.2 (30) is used for data analysis with the following strategy: DRAO is used as a fixed reference site; precise IGS (International GPS Service) satellite orbits are used and kept fixed; ionospheric-free phase solutions are used to determine relative positions of network sites; phase ambiguities are resolved and fixed to integer values; tropospheric zenith delay is estimated hourly with no a priori tropospheric model and using a dry Niell tropospheric mapping function; tropospheric gradients (tilt direction of the mapping function) are estimated every 6 hours; a nominal 10° cut-off elevation is used for satellites; and solid Earth tide, pole tide, and ocean loading corrections are applied.

17. R. E. Wells, C. S. Weaver, R. J. Blakely, Geology 26, 759 (1998).
18. K. Wang, Geophys. Res. Lett. 23, 2021 (1996).
19. T. S. James, J. J. Clague, K. Wang, I. Hutchinson, Quat. Sci. Rev. 19, 1527 (2000).
20. P. Flu¨ck, R. D. Hyndman, K. Wang, J. Geophys. Res. 102, 20539 (1997).
21. Y. Okada, Bull. Seismol. Soc. Am. 75, 1135 (1985).
22. K. Wang, T. Mulder, G. C. Rogers, R. D. Hyndman, J. Geophys. Res. 100, 12907 (1995).
23. S. Ozawa, M. Murakami, T. Tada, J. Geophys. Res. 106, 787 (2001).
24. K. Heki, S. Miyazaki, H. Tsuji, Nature 386, 595 (1997).
25. T. Sagiya, Eos 78, F165 (1997).
26. A. T. Linde, M. T. Gladwin, M. J. S. Johnston, R. L. Gwyther, R. G. Bilham, Nature 383, 65 (1996).
27. H. Hirose, K. Hirahara, F. Kimata, N. Fujii, S. Miyazaki, Geophys. Res. Lett. 26, 3237 (1999).
28. K. Wang, R. D. Hyndman, M. Yamano, Tectonophysics 248, 53 (1995).
29. J. O. Langbein, A F. Wyatt, A H. Johnson, A D. Hamann, P. Zimmer, Geophys. Res. Lett. 22, 3533 (1995).
30. G. Beutler et al., Bernese GPS Software Version 4.2 (Astronomical Institute, University of Berne, Berne, Switzerland, 2000).

31. We thank M. Schmidt and Y. Lu for their support in GPS network operations, and Central Washington University (M. Miller) and the University of Washington (A. Qamar) for providing GPS data from sites of the PANGA network, which was established with support from NSF. Supported by U.S. Geological Survey National Earthquake Hazards Reduction Program research grant 00HQGR0061; this paper is Geological Survey of Canada contribution no. 2001002. Copyright, Her Majesty the Queen in right of Canada (2001).