- Email: [email protected]
Interseismic locking of the plate interface in the northern Cascadia subduction zone, inferred from inversion of GPS data
Earth and Planetary Science Letters 231 (2005) 239 – 247 www.elsevier.com/locate/epsl
Interseismic locking of the plate interface in the northern Cascadia subduction zone, inferred from inversion of GPS data Shoichi Yoshiokaa,*, Kelin Wangb,c, Stephane Mazzottib a
Department of Earth and Planetary Sciences, Kyushu University, Hakozaki 6-10-1, Higashi ward, Fukuoka 812-8581, Japan b Geological Survey of Canada, Pacific Geoscience Centre, 9860 West Saanich Rd, Sidney, B.C., Canada V8L 4B2 c School of Earth and Ocean Sciences, University of Victoria, Victoria, B.C., Canada Received 29 April 2004; received in revised form 6 November 2004; accepted 10 December 2004 Editor: R.D. van der Hilst
Abstract We inverted GPS velocities from 20 continuous and 53 campaign sites in the northern Cascadia subduction zone using a Bayesian inverse method to estimate the locking state of the plate interface. The results are consistent with previous estimates based on thermal arguments and forward modeling. They suggest that the completely locked segment of the plate interface is offshore and that the degree of locking gradually decreases landward. The very gradual transition from full locking to full slip approximates the effect of stress relaxation that is not included by the elastic model assumed for the inversion. D 2005 Elsevier B.V. All rights reserved. Keywords: Cascadia subduction zone; ABIC; GPS data inversion; interseismic deformation
1. Introduction At the Cascadia subduction zone, the Juan de Fuca (JDF) plate subducts beneath the North America (NA) plate. The last great earthquake at this margin occurred in 1700 with an estimated moment magnitude of 9 [1]. Most elastic dislocation models constrained by GPS (Global Positioning System) measurements of crustal deformation show that the * Corresponding author. Tel.: +81 92 642 2646; fax: +81 92 642 2684. E-mail address: [email protected] (S. Yoshioka). 0012-821X/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.epsl.2004.12.018
megathrust plate boundary is fully locked, and that the locked zone is mostly offshore where the subducting plate is in contact with the large accretionary prism [2–4]. In addition to interseismic deformation, the southern Cascadia forearc is subject to a long-term motion in the north–northwest direction [5]. Uncertainties in the kinematics of this forearc motion permit models that allow the locked zone to slip at a fraction of the plate convergence rate (incomplete locking) in southern Cascadia [6]. The northernmost Cascadia forearc is little affected by this forearc motion, and the interpretation of GPS data is simpler. However, all published models for this
240
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
region are of the trial-and-error type. Because of the non-uniqueness of the inverse problem of inferring slip distribution on the plate interface from limited surface deformation measurements, different slip models may provide similarly good fit to GPS observations. One way to address this problem is to carry out a formal inversion [7,8]. With clearly stated model constraints and assumptions, an inversion will yield a best estimate of the slip pattern and provide an error estimate that quantifies the probability of other, less preferred, models. In this paper, we report the results of an inversion study in northern Cascadia.
2. Plate convergence and GPS data The four open arrows in Fig. 1 along the accretionary-prism deformation front represent the convergence velocities of the Juan de Fuca (JDF) plate
relative to the North America (NA) plate. Following Mazzotti et al. [9], we have used an Euler pole at 35.288N/116.378W with a rotation rate of 1.5048/ m.y., based on the Pacific-NA Euler vector of DeMets and Dixon [10] and the most recent (0.4 Ma) JDFPacific Euler vector of Wilson [11]. Within the map area, the average convergence direction is N55.38E and the convergence rate increases northward from 38 to 47 mm/yr. Shown in the same figure are GPS velocities with respect to stable North America as defined in the ITRF2000 reference frame [12]. Velocities at the northernmost 20 continuously monitoring sites of the Western Canada Deformation Array (WCDA) and the Pacific Northwest Geodetic Array (PANGA) were included in our analysis (solid arrows). The continuous GPS data were processed for the period of January 1996–August 2002. The most recently installed stations had been recording for about 3 yr.
100 km
Vancouver Island 50˚N
No
ot
ka
F.
Explorer plate
nF
1σ=2 mm
Washington
atio
10 mm/yr
orm
Def
48˚N
t
ron
46˚N
Juan de Fuca plate 128˚W
North America plate 126˚W
124˚W
122˚W
Fig. 1. Tectonic setting and GPS observations. Triangles denote locations of volcanoes. The barbed line represents the location of the deformation front. The four open arrows are plate motion vectors of the Juan de Fuca plate relative to the North America plate. Velocities at continuous GPS sites are shown as black arrows and campaign sites as gray arrows. Nootka F.—Nootka fault zone.
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
For campaign-style GPS data (gray arrows), we only used measurements made on Vancouver Island. Three geodetic networks each having 15–20 sites were surveyed twice during the 1991–1994 and 1996– 1999 periods. Typical surveys consisted of 2–6 occupations of 5–24 h for each site, repeated after 5–6 yr [9,13]. GPS data from various sources were processed using a consistent double-difference technique with common reference station DRAO (49.32258N, 119.62508W) held fixed. Uncertainties in continuous GPS velocities are based on a frequency-dependent noise model [14]. Uncertainties in the campaign velocities are based on error budget considerations [9]. Uncertainties in the final GPS velocities, shown as error ellipses in Fig. 1, also include those incurred during the transformation of the velocity field into the stable North America reference. Further details of the GPS data have been discussed in Mazzotti et al. [9]. The most preeminent feature of the GPS velocities is a northeastward motion, decreasing from ~15–20 mm/yr along the west coast to 2–3 mm/yr and less on the mainland (Fig. 1). The landward decrease reflects strain accumulation as a result of the locking of the plate interface of the subduction zone. Full locking (maximum slip deficit rate) of the shallow seismogenic part of the plate interface causes the largest strain accumulation, and full slipping (zero slip deficit) causes no strain. The velocity pattern is episodically and briefly interrupted by slow slip events that occur downdip from the seismogenic zone [15]. The asso-
(a)
(b)
241
ciated signal in the GPS time series was removed by assuming that these transient events do not contribute to the longer-term deformation of the upper plate [9].
3. Inversion method In an elastic dislocation model for interseismic deformation, slip deficit on the plate interface is represented by fault slip in a reverse direction [16], usually called the back-slip rate. The slip deficit is parameterized as follows in our inversion. A model region that is roughly rectangular in plan view is used to define the segment on which slip deficit is to be determined. The 3-D curved plate interface in the region is divided into 75 smaller segments (Fig. 2). A uniform back-slip rate is assumed for each curved segment. Surface velocities due to the uniform slip on each segment are determined by numerical integration using a point source dislocation solution [17]. The slip deficit vectors for the 35 segments, expressed in parameter vector a, is related to GPS velocities, expressed in data vector d, by d=Ha, where matrix H depends on the dislocation theory and how the slip is parameterized. Parameters a are determined from data d using an inverse method based on a Bayesian Information Criterion (ABIC) [18,19]. Uncertainties in the GPS velocities are assumed to have a Gaussian probability density function (pdf) with a covariance matrix r 2E, where the scaling factor r 2 is to be determined. For this analysis, we assume
(c)
(d)
50˚N
48˚N
46˚N
128˚W
124˚W
128˚W
124˚W
128˚W
124˚W
128˚W
124˚W
Fig. 2. (a–c) Three model regions with slightly different orientations used for the inversion. For each model, all GPS data were inverted to infer slip deficit distribution, assumed to be uniform on each of the 75 segments. Dashed lines are contours of the plate interface according to McCrory et al. [20] at 20 km intervals. Open triangles denote locations of GPS stations. (d) The final model constructed by piecing together segments from the northernmost three rows of (a), central row of (b), and southernmost three rows of (c).
242
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
that data uncertainties are uncorrelated, and hence E is diagonal. A similar a priori pdf is assumed for the 35 slip deficit vectors, but in place of the inverse of a covariance matrix is a matrix G/q 2 that controls the smoothness of the slip deficit, where q 2 is an unknown scaling factor. The general form of the G matrix is given in Yabuki and Matsu’ura [19]. It can be shown that minimizing the following objective function yields the best estimate for a, S ðaÞ ¼ ðd HaÞT E1 ðd HaÞ þ a2 aT Ga
constraints, and therefore the trade off between data fit and the smoothness of slip deficit distribution. A large a 2 gives a smoother slip distribution at the expense of a poorer data fit. A very small a 2 implies a conventional least squares method and may give an excellent fit but a rugged slip distribution. The optimal value of a 2 is simultaneously determined with a by the inversion algorithm using the ABIC. The determination of r 2 and the a posteriori covariance of a and other details of the method are explained by Yabuki and Matsu’ura [19]. There is no parameter that needs to be adjusted in an ad hoc fashion. Because the present method requires a straight southwestern boundary (in plan-view), the curvature
ð1Þ
The hyperparameter a 2=r 2/q 2 quantifies the relative weights of the a priori data variances vs. smoothness
0
-10
0
-8
-60
-40
-20
100 km
50˚N
80
0
0
0
-60
-4
-2
48˚N
46˚N
40 mm/yr
128˚W
126˚W
124˚W
122˚W
Fig. 3. Slip deficit vectors determined by the inversion. The black, gray, and open arrows represent results on the 75 fault segments for the models shown in Fig. 2a, b, and c, respectively. Dashed lines are contours of the plate interface according to McCrory et al. [20] at 20 km intervals.
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
of the deformation front presents a problem. Each of the three options shown in Fig. 2a, b, and c will cut off some shallow area of the plate interface. To overcome this problem, we decided to perform inversion on the three individual model regions on the 3-D curved plate interface recently proposed by McCrory et al. [20] in Fig. 2 a, b, and c but use results from relevant segments of each model to piece together a final model (Fig. 2d). In Fig. 3, we show the slip deficits determined using all three different model regions. The slip vectors from these models are mutually consistent for overlapping areas. All inversions fit the GPS data equally well (not shown for individual models). Therefore, the abovementioned procedure introduces little error. The northernmost row of the fault segments in Fig. 2a is located in a different tectonic regime where the Explorer plate subducts beneath NA at a slower rate than JDF. No GPS data from that tectonic regime are used in our analysis, but these segments are included such that the model
243
boundary is set farther away from the main region of interest to minimize potential artifacts.
4. Results The slip deficits determined by the inversion are shown in Fig. 4a, with the error ellipses indicating the 95% confidence interval (2r). The error ellipses, based on the a posteriori covariance matrix of the slip vectors [21], quantify the errors in the slip vectors and are a formal description of the model sensitivity and resolution. A larger error ellipse results from poorer data constraints. As explained above, the northernmost five vectors should be interpreted with caution. The southernmost row of vectors cannot be very well constrained by the sparse GPS observations used for the southern region (Fig. 2), as reflected in their larger error ellipses.
(a)
(b)
100 km
100 km
0
00
-1
, A
-8
0
0
-6
-4
0
0
00
-2
-1
-8
0
0
-6
0
-4
-2
50˚N
A 48˚N mm/yr
46˚N
40
40
20
20
0
0
40 mm/yr 128˚W
, B
mm/yr
126˚W
B
40 mm/yr 124˚W
122˚W
128˚W
126˚W
124˚W
122˚W
Fig. 4. (a) Slip deficit vector distribution (open arrows), which is uniform within each fault segment, pieced together from three individual inversions as illustrated in Fig. 2d. The solid arrows are plate motion vectors. (b) Real slip rate vectors (arrows) calculated from the slip deficits of (a), shown as the motion of the subducting plate relative to the upper plate. Zero slip means full locking. Lines A–AV and B–BV indicate locations of the profiles shown in Fig. 6. In both (a) and (b), the magnitude of the slip rates is contoured with discrete grayscale at 10 mm/yr intervals, and the depth of the plate interface with dashed lines at 20 km intervals.
244
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
Surface velocities predicted by the model fit the GPS data very well, and almost all of them are within the 1-r data error ellipses of Fig. 1. Results for the model of Fig. 2a are shown in Fig. 5 as an example. The root root-mean-square difference between predicted and observed velocities is 0.8 mm/yr, and the reduced v 2 fit is 1.0. The calculated surface velocities predicted by the models in Fig. 2b and c are almost the same as those shown in Fig. 5. In Fig. 4b, the slip deficit has been translated into actual slip rate. The vectors represent the motion of the subducting plate relative to the upper plate and have been obtained by subtracting the back-slip rates of Fig. 4a from (steady) plate convergence rates defined by the Euler pole described in Section 2. Important features of the inversion results are summarized and discussed as follows: (1)
The slip rate vectors of the most seaward segments off Washington are all within error
(2)
ellipses (Fig. 4b), indicating that within model errors this offshore segment is fully locked. This is consistent with the thermal arguments by Hyndman and Wang [22] and lends support to all previous trial-and-error based dislocation models for northern Cascadia. The degree of locking very gradually decreases landward, as shown by the increasing size of the slip rate vectors (Fig. 4b). From the fully locked region, the slip rates decrease northeastward. The pattern is in qualitative agreement with the recent 3-D dislocation model CAS3D2 [4] that used a more gradual decrease than earlier models constrained by mostly pre-GPS data [23]. Comparison with CAS3D-2 along profiles A–AV and B–BV (Fig. 4b) is made in Fig. 6. The unphysical negative, i.e., normalfaulting type, slip seen in the most seaward segment of profile B–BV is an artifact of the inversion, reflecting the diminished resolving power of the GPS data for the far offshore part
50˚N
48˚N
obs. 10 mm/yr 1σ =2 mm cal. 10 mm/yr
100 km 128˚W
126˚W
124˚W
122˚W
Fig. 5. Observed (solid arrows) and calculated (open arrows) horizontal velocities predicted by the model shown in Fig. 2a. The thick black and thin gray solid arrows represent permanent and campaign data, respectively.
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
40
70
30
60
20
50
10
40
0
30
-10
20
-20
10
-30
0
Slip Rate (mm/yr)
B
Depth (km)
A'
B'
40
70
30
60
20
50
10
40
0
30
-10
20
-20
10
-30
Depth (km)
Slip Rate (mm/yr)
A
245
0 0
50
100
150
200
250
Distance From Deformation Front (km) Fig. 6. Comparison of slip rates along two profiles shown in Fig. 4b with the CAS3D-2 model (gray line). Dashed lines represent 2r errors for the inversion results. The dotted lines represent depths to the plate interface along the cross sections.
(3)
(4)
of the plate interface. The lack of resolution is reflected in the large error bars (and error ellipses in Fig. 4) of the slip values in this area. This is a common problem in defining offshore fault motion/locking from land-based geodetic observations. The slip rate vectors in the southernmost region of the model are less northerly than plate convergence (Fig. 4b), or equivalently, the slip deficit vectors are more northerly than plate convergence (Fig. 4a). We think this is due to the effect of the forearc motion mostly south of our study area (see Section 1). GPS stations in northern Washington are slightly affected by this motion and move in a more northerly direction than plate convergence [24]. This effect is translated into a more northerly slip deficit. Along the updip end of the megathrust, the degree of locking is fairly uniform south of 488N but becomes smaller further north (Fig.
4b). There does not seem to be a lack of data constraint in that area. The age of the subducting plate decreases northward, and the higher temperature of the younger incoming plate may result in aseismic slip of the plate interface all the way to the seafloor deformation front. Alternatively, the locked zone may be extremely narrow and far offshore because of the higher temperature but unresolved by land GPS data and our model parameterization. In addition, the subduction of the Explorer plate just north of the Nootka fault (Fig. 1) is much slower [25], and may affect deformation on land just south of the fault. In cross-section view along profile A–AV (Fig. 6), the landward increase in slip rate towards full plate convergence rate is more gradual than assumed in the CAS3D-2 model. This result is consistent with findings by Mazzotti et al. [9], who related this behavior to possible long-term upper plate deformation.
246
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
5. Discussions One technical point that deserves consideration is the smoothness constraint for the inversion. The results will be smoother if correlation between parameters is introduced, such that neighboring segments have a tendency to have similar slip rates. The G matrix, which is composed of the second derivative of slip distribution with respect to two coordinates taken on the curved plate interface, introduces cross correlation only for neighboring segments, and the effect on smoothness is very small. The inversion results will also be smoother if individual parameters are more tightly constrained relative to data variance. This is determined by the hyperparameter a 2. Although the optimal values (0.00798, 0.00714, and 0.00409 for the models shown in Fig. 2a, b, and c, respectively) of a 2 are automatically determined using ABIC, it is useful to investigate the sensitivity of our results for changes in a 2. For this purpose, we conducted inversion tests by fixing a 2 values without using the ABIC. Using an a 2 value smaller than the optimal value by a factor of five resulted in better data fit but more rapid spatial changes in both slip directions and rates than in Fig. 4a. Since the data include observation errors, smaller misfit does not necessarily indicate better solution. The effects on the misfit and smoothness were the opposite if we increased the a 2 value by a factor of five. However, the difference between the results of both these tests and the results of Fig. 4a was less than 20% in either rates or directions. The relatively low sensitivity of our inversion results to the smoothing constraint can be attributed to the high density and internal consistency of the GPS data. When interpreting the results, it is important to remember the limitations of the elastic half-space model assumed for the inversion. The upper mantle of active margins is viscoelastic, with a viscosity commonly believed to be around 1019 Pa s (see review by Wang [26]). Coseismic deformation and some short-term transients such as the episodic silent slips [15] can be described using elastic models. For interseismic deformation, especially at a time as long as 300 yr after a great earthquake, the effect of viscoelastic stress relaxation may become important. In elastic models, effects of stress relaxation on surface deformation have to be attributed to equiv-
alent slip [27]. For this reason, the zone of downdip transition from full locking to full slip is called the beffective transition zoneQ [4]. It may represent some actual creep along the plate interface, but it also accounts for the effect of stress relaxation in the upper mantle. An equivalent fault slip distribution implicitly represents a rheological structure and the integrated effect of stress relaxation since the last earthquake, but we are generally not in the position to determine uniquely the rheology and deformation history from the limited geodetic data at Cascadia [26]. The modern geodetic observations and in a sense the elastic model provide a bsnapshotQ of an evolving deformation field.
6. Conclusions We have inverted GPS data in the northern Cascadia using an ABIC-based inversion algorithm to infer the locking state of the plate interface in the subduction zone. The GPS data include velocities derived at 20 continuously monitoring stations in British Columbia and northern Washington and 53 campaign sites on Vancouver Island. The inversion results support the previous conclusion based on thermal arguments and trial-anderror dislocation models that the shallower plate interface is fully locked and that the locked zone is mostly offshore. The results are also consistent with recent models in that the downdip transition from fully locked to full slip takes place over a very wide area, but this transition should be understood to be mostly an equivalent description of the effect of stress relaxation that is not directly considered by the elastic model. The inversion also reveals some other complications, such as a tapering of the degree of locking near the northern end of the subduction zone. Northward motion of the Cascadia forearc takes place mostly south of our study area, but it may still slightly affect the southernmost segments of our model.
Acknowledgments Continuous GPS data used in this study are from the Western Canada Deformation Array (WCDA) and northern stations of the Pacific Northwest Geodetic
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
Array (PANGA). Campaign GPS surveys were carried out by the Geodetic Survey Division of Geomatics Canada. We thank T. Yabuki for sharing his original inversion source code, P. Wessel and W. H. F. Smith for making available the Generic Mapping Tool (GMT), with which the figures for this paper were created, and E. Hearn, E. Hetland, H. Dragert, and the editor, R. van der Hilst for their valuable comments. The work was initiated when SY was a visiting scientist at the Pacific Geoscience Centre in 2001 supported by Japan Society for the Promotion of Science. Geological Survey of Canada contribution 2004228.
References [1] K. Satake, K. Wang, B.F. Atwater, Fault slip and seismic moment of the 1700 Cascadia earthquake inferred from Japanese tsunami descriptions, J. Geophys. Res. 108 (B11) (2003) 2535, doi:10.1029/2003JB002521. [2] G. Khazaradze, A. Qamar, H. Dragert, Tectonic deformation in western Washington from continuous GPS measurements, Geophys. Res. Lett. 26 (1999) 3153 – 3156. [3] M. Miller, D.J. Johnson, C.M. Rubin, H. Dragert, K. Wang, A. Qamar, C. Goldfinger, GPS-determination of along-strike variation in Cascadia margin kinematics: implications for relative plate motion, subduction zone coupling, and permanent deformation, Tectonics 20 (2001) 161 – 176. [4] K. Wang, R.E. Wells, S. Mazzotti, R.D. Hyndman, T. Sagiya, A revised dislocation model of interseismic deformation of the Cascadia subduction zone, J. Geophys. Res. 108 (2003), doi:10.1029/2001JB001227. [5] R.E. Wells, R.W. Simpson, Northward migration of the Cascadia forearc in the northwestern U.S. and implications for subduction deformation, Earth Planets Space 53 (2001) 275 – 283. [6] R. McCaffrey, Crustal block rotations and plate coupling, in: S. Stein, J. Freymueller (Eds.), Plate Boundary Zones, AGU Geodynamics Series, vol. 30, 2002, pp. 101 – 122. [7] S. Yoshioka, T. Yabuki, T. Sagiya, T. Tada, M. Matsu’ura, Interplate coupling and relative plate motion in the Tokai district, central Japan, deduced from geodetic data inversion using ABIC, Geophys. J. Int. 113 (1993) 607 – 621. [8] S. Yoshioka, T. Yabuki, T. Sagiya, T. Tada, M. Matsu’ura, Interplate coupling in the Kanto district, central Japan, deduced from geodetic data inversion and its tectonic implications, Tectonophysics 229 (1994) 181 – 200. [9] S. Mazzotti, H. Dragert, J. Henton, M. Schmidt, R. Hyndman, T. James, Y. Lu, M. Craymer, Current tectonics of northern Cascadia from a decade of GPS measurements, J. Geophys. Res. 108 (B12) (2003) 2554, doi:10.1029/2003JB002653.
247
[10] C. DeMets, T.H. Dixon, New kinematic models for the Pacific-North America motion from 3 Ma to present: I. Evidence for steady motion and biases in the NUVEL-1A model, Geophys. Res. Lett. 26 (1999) 1921 – 1924. [11] D.S. Wilson, Confidence intervals for motion and deformation of the Juan de Fuca plate, J. Geophys. Res. 98 (1993) 16053 – 16071. [12] Z. Altamimi, P. Sillard, C. Boucher, ITRF2000: a new release of the international terrestrial reference frame for earth science applications, J. Geophys. Res. 107 (B10) (2002) 2214, doi:10.1029/2001JB000561. [13] J.A. Henton, GPS Studies of Crustal Deformation in the Northern Cascadia Subduction Zone, PhD thesis, University of Victoria, British Columbia (2000). [14] A. Mao, C.G.A. Harrison, T. Dixon, Noise in GPS coordinate time series, J. Geophys. Res. 104 (1999) 2797 – 2816. [15] H. Dragert, K. Wang, T.S. James, A silent slip event on the deeper Cascadia subduction interface, Science 292 (2001) 1525 – 1528. [16] J.C. Savage, A dislocation model of strain accumulation and release at a subduction zone, J. Geophys. Res. 88 (1983) 4984 – 4996. [17] T. Maruyama, Statical elastic dislocations in an infinite and semi-infinite medium, Bull. Earthq. Res. Inst. Univ. Tokyo 42 (1964) 289 – 368. [18] H. Akaike, Likelihood and the Bayes procedure, in: J.M. Bernardo, M.H. DeGroot, D.V. Lindley, A.F.M. Smith (Eds.), Bayesian Statistics, University Press, Valencia, Spain, 1980, pp. 143 – 166. [19] T. Yabuki, M. Matsu’ura, Geodetic data inversion using a Bayesian information criterion for spatial distribution of fault slip, Geophys. J. Int. 109 (1992) 363 – 375. [20] P.A. McCrory, J.L. Blair, D.H. Oppenheimer, S.R. Walter, Depth to the Juan de Fuca slab beneath the Cascadia subduction margin: a 3-D model for sorting earthquakes, U.S. Geological Survey Digital Data Series, 1 CD-ROM (2004). [21] D.D. Jackson, M. Matsu’ura, A Bayesian approach to nonlinear inversion, J. Geophys. Res. 90 (1985) 581 – 591. [22] R.D. Hyndman, K. Wang, Thermal constraints on the zone of major thrust earthquake failure: the Cascadia subduction zone, J. Geophys. Res. 98 (1993) 2039 – 2060. [23] P. Flqck, R.D. Hyndman, K. Wang, Three-dimensional dislocation model for great earthquakes of the Cascadia subduction zone, J. Geophys. Res. 102 (1997) 20539 – 20550. [24] S. Mazzotti, H. Dragert, R.D. Hyndman, M.M. Miller, J.A. Henton, GPS deformation in a region of high crustal seismicity: N. Cascadia forearc, Earth Planet. Sci. Lett. 198 (2002) 41 – 48. [25] R. Riddihough, Recent movements of the Juan de Fuca plate system, J. Geophys. Res. 89 (1984) 6980 – 6994. [26] K. Wang, Elastic and viscoelastic models for subduction earthquake cycles, in: T. Dixon (Ed.), The Seismogenic Zone Experiment, in press. [27] J.C. Savage, Viscoelastic-coupling model for the earthquake cycle driven from below, J. Geophys. Res. 105 (B11) (2000) 25525 – 25532.
Interseismic locking of the plate interface in the northern Cascadia subduction zone, inferred from inversion of GPS data Shoichi Yoshiokaa,*, Kelin Wangb,c, Stephane Mazzottib a
Department of Earth and Planetary Sciences, Kyushu University, Hakozaki 6-10-1, Higashi ward, Fukuoka 812-8581, Japan b Geological Survey of Canada, Pacific Geoscience Centre, 9860 West Saanich Rd, Sidney, B.C., Canada V8L 4B2 c School of Earth and Ocean Sciences, University of Victoria, Victoria, B.C., Canada Received 29 April 2004; received in revised form 6 November 2004; accepted 10 December 2004 Editor: R.D. van der Hilst
Abstract We inverted GPS velocities from 20 continuous and 53 campaign sites in the northern Cascadia subduction zone using a Bayesian inverse method to estimate the locking state of the plate interface. The results are consistent with previous estimates based on thermal arguments and forward modeling. They suggest that the completely locked segment of the plate interface is offshore and that the degree of locking gradually decreases landward. The very gradual transition from full locking to full slip approximates the effect of stress relaxation that is not included by the elastic model assumed for the inversion. D 2005 Elsevier B.V. All rights reserved. Keywords: Cascadia subduction zone; ABIC; GPS data inversion; interseismic deformation
1. Introduction At the Cascadia subduction zone, the Juan de Fuca (JDF) plate subducts beneath the North America (NA) plate. The last great earthquake at this margin occurred in 1700 with an estimated moment magnitude of 9 [1]. Most elastic dislocation models constrained by GPS (Global Positioning System) measurements of crustal deformation show that the * Corresponding author. Tel.: +81 92 642 2646; fax: +81 92 642 2684. E-mail address: [email protected] (S. Yoshioka). 0012-821X/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.epsl.2004.12.018
megathrust plate boundary is fully locked, and that the locked zone is mostly offshore where the subducting plate is in contact with the large accretionary prism [2–4]. In addition to interseismic deformation, the southern Cascadia forearc is subject to a long-term motion in the north–northwest direction [5]. Uncertainties in the kinematics of this forearc motion permit models that allow the locked zone to slip at a fraction of the plate convergence rate (incomplete locking) in southern Cascadia [6]. The northernmost Cascadia forearc is little affected by this forearc motion, and the interpretation of GPS data is simpler. However, all published models for this
240
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
region are of the trial-and-error type. Because of the non-uniqueness of the inverse problem of inferring slip distribution on the plate interface from limited surface deformation measurements, different slip models may provide similarly good fit to GPS observations. One way to address this problem is to carry out a formal inversion [7,8]. With clearly stated model constraints and assumptions, an inversion will yield a best estimate of the slip pattern and provide an error estimate that quantifies the probability of other, less preferred, models. In this paper, we report the results of an inversion study in northern Cascadia.
2. Plate convergence and GPS data The four open arrows in Fig. 1 along the accretionary-prism deformation front represent the convergence velocities of the Juan de Fuca (JDF) plate
relative to the North America (NA) plate. Following Mazzotti et al. [9], we have used an Euler pole at 35.288N/116.378W with a rotation rate of 1.5048/ m.y., based on the Pacific-NA Euler vector of DeMets and Dixon [10] and the most recent (0.4 Ma) JDFPacific Euler vector of Wilson [11]. Within the map area, the average convergence direction is N55.38E and the convergence rate increases northward from 38 to 47 mm/yr. Shown in the same figure are GPS velocities with respect to stable North America as defined in the ITRF2000 reference frame [12]. Velocities at the northernmost 20 continuously monitoring sites of the Western Canada Deformation Array (WCDA) and the Pacific Northwest Geodetic Array (PANGA) were included in our analysis (solid arrows). The continuous GPS data were processed for the period of January 1996–August 2002. The most recently installed stations had been recording for about 3 yr.
100 km
Vancouver Island 50˚N
No
ot
ka
F.
Explorer plate
nF
1σ=2 mm
Washington
atio
10 mm/yr
orm
Def
48˚N
t
ron
46˚N
Juan de Fuca plate 128˚W
North America plate 126˚W
124˚W
122˚W
Fig. 1. Tectonic setting and GPS observations. Triangles denote locations of volcanoes. The barbed line represents the location of the deformation front. The four open arrows are plate motion vectors of the Juan de Fuca plate relative to the North America plate. Velocities at continuous GPS sites are shown as black arrows and campaign sites as gray arrows. Nootka F.—Nootka fault zone.
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
For campaign-style GPS data (gray arrows), we only used measurements made on Vancouver Island. Three geodetic networks each having 15–20 sites were surveyed twice during the 1991–1994 and 1996– 1999 periods. Typical surveys consisted of 2–6 occupations of 5–24 h for each site, repeated after 5–6 yr [9,13]. GPS data from various sources were processed using a consistent double-difference technique with common reference station DRAO (49.32258N, 119.62508W) held fixed. Uncertainties in continuous GPS velocities are based on a frequency-dependent noise model [14]. Uncertainties in the campaign velocities are based on error budget considerations [9]. Uncertainties in the final GPS velocities, shown as error ellipses in Fig. 1, also include those incurred during the transformation of the velocity field into the stable North America reference. Further details of the GPS data have been discussed in Mazzotti et al. [9]. The most preeminent feature of the GPS velocities is a northeastward motion, decreasing from ~15–20 mm/yr along the west coast to 2–3 mm/yr and less on the mainland (Fig. 1). The landward decrease reflects strain accumulation as a result of the locking of the plate interface of the subduction zone. Full locking (maximum slip deficit rate) of the shallow seismogenic part of the plate interface causes the largest strain accumulation, and full slipping (zero slip deficit) causes no strain. The velocity pattern is episodically and briefly interrupted by slow slip events that occur downdip from the seismogenic zone [15]. The asso-
(a)
(b)
241
ciated signal in the GPS time series was removed by assuming that these transient events do not contribute to the longer-term deformation of the upper plate [9].
3. Inversion method In an elastic dislocation model for interseismic deformation, slip deficit on the plate interface is represented by fault slip in a reverse direction [16], usually called the back-slip rate. The slip deficit is parameterized as follows in our inversion. A model region that is roughly rectangular in plan view is used to define the segment on which slip deficit is to be determined. The 3-D curved plate interface in the region is divided into 75 smaller segments (Fig. 2). A uniform back-slip rate is assumed for each curved segment. Surface velocities due to the uniform slip on each segment are determined by numerical integration using a point source dislocation solution [17]. The slip deficit vectors for the 35 segments, expressed in parameter vector a, is related to GPS velocities, expressed in data vector d, by d=Ha, where matrix H depends on the dislocation theory and how the slip is parameterized. Parameters a are determined from data d using an inverse method based on a Bayesian Information Criterion (ABIC) [18,19]. Uncertainties in the GPS velocities are assumed to have a Gaussian probability density function (pdf) with a covariance matrix r 2E, where the scaling factor r 2 is to be determined. For this analysis, we assume
(c)
(d)
50˚N
48˚N
46˚N
128˚W
124˚W
128˚W
124˚W
128˚W
124˚W
128˚W
124˚W
Fig. 2. (a–c) Three model regions with slightly different orientations used for the inversion. For each model, all GPS data were inverted to infer slip deficit distribution, assumed to be uniform on each of the 75 segments. Dashed lines are contours of the plate interface according to McCrory et al. [20] at 20 km intervals. Open triangles denote locations of GPS stations. (d) The final model constructed by piecing together segments from the northernmost three rows of (a), central row of (b), and southernmost three rows of (c).
242
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
that data uncertainties are uncorrelated, and hence E is diagonal. A similar a priori pdf is assumed for the 35 slip deficit vectors, but in place of the inverse of a covariance matrix is a matrix G/q 2 that controls the smoothness of the slip deficit, where q 2 is an unknown scaling factor. The general form of the G matrix is given in Yabuki and Matsu’ura [19]. It can be shown that minimizing the following objective function yields the best estimate for a, S ðaÞ ¼ ðd HaÞT E1 ðd HaÞ þ a2 aT Ga
constraints, and therefore the trade off between data fit and the smoothness of slip deficit distribution. A large a 2 gives a smoother slip distribution at the expense of a poorer data fit. A very small a 2 implies a conventional least squares method and may give an excellent fit but a rugged slip distribution. The optimal value of a 2 is simultaneously determined with a by the inversion algorithm using the ABIC. The determination of r 2 and the a posteriori covariance of a and other details of the method are explained by Yabuki and Matsu’ura [19]. There is no parameter that needs to be adjusted in an ad hoc fashion. Because the present method requires a straight southwestern boundary (in plan-view), the curvature
ð1Þ
The hyperparameter a 2=r 2/q 2 quantifies the relative weights of the a priori data variances vs. smoothness
0
-10
0
-8
-60
-40
-20
100 km
50˚N
80
0
0
0
-60
-4
-2
48˚N
46˚N
40 mm/yr
128˚W
126˚W
124˚W
122˚W
Fig. 3. Slip deficit vectors determined by the inversion. The black, gray, and open arrows represent results on the 75 fault segments for the models shown in Fig. 2a, b, and c, respectively. Dashed lines are contours of the plate interface according to McCrory et al. [20] at 20 km intervals.
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
of the deformation front presents a problem. Each of the three options shown in Fig. 2a, b, and c will cut off some shallow area of the plate interface. To overcome this problem, we decided to perform inversion on the three individual model regions on the 3-D curved plate interface recently proposed by McCrory et al. [20] in Fig. 2 a, b, and c but use results from relevant segments of each model to piece together a final model (Fig. 2d). In Fig. 3, we show the slip deficits determined using all three different model regions. The slip vectors from these models are mutually consistent for overlapping areas. All inversions fit the GPS data equally well (not shown for individual models). Therefore, the abovementioned procedure introduces little error. The northernmost row of the fault segments in Fig. 2a is located in a different tectonic regime where the Explorer plate subducts beneath NA at a slower rate than JDF. No GPS data from that tectonic regime are used in our analysis, but these segments are included such that the model
243
boundary is set farther away from the main region of interest to minimize potential artifacts.
4. Results The slip deficits determined by the inversion are shown in Fig. 4a, with the error ellipses indicating the 95% confidence interval (2r). The error ellipses, based on the a posteriori covariance matrix of the slip vectors [21], quantify the errors in the slip vectors and are a formal description of the model sensitivity and resolution. A larger error ellipse results from poorer data constraints. As explained above, the northernmost five vectors should be interpreted with caution. The southernmost row of vectors cannot be very well constrained by the sparse GPS observations used for the southern region (Fig. 2), as reflected in their larger error ellipses.
(a)
(b)
100 km
100 km
0
00
-1
, A
-8
0
0
-6
-4
0
0
00
-2
-1
-8
0
0
-6
0
-4
-2
50˚N
A 48˚N mm/yr
46˚N
40
40
20
20
0
0
40 mm/yr 128˚W
, B
mm/yr
126˚W
B
40 mm/yr 124˚W
122˚W
128˚W
126˚W
124˚W
122˚W
Fig. 4. (a) Slip deficit vector distribution (open arrows), which is uniform within each fault segment, pieced together from three individual inversions as illustrated in Fig. 2d. The solid arrows are plate motion vectors. (b) Real slip rate vectors (arrows) calculated from the slip deficits of (a), shown as the motion of the subducting plate relative to the upper plate. Zero slip means full locking. Lines A–AV and B–BV indicate locations of the profiles shown in Fig. 6. In both (a) and (b), the magnitude of the slip rates is contoured with discrete grayscale at 10 mm/yr intervals, and the depth of the plate interface with dashed lines at 20 km intervals.
244
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
Surface velocities predicted by the model fit the GPS data very well, and almost all of them are within the 1-r data error ellipses of Fig. 1. Results for the model of Fig. 2a are shown in Fig. 5 as an example. The root root-mean-square difference between predicted and observed velocities is 0.8 mm/yr, and the reduced v 2 fit is 1.0. The calculated surface velocities predicted by the models in Fig. 2b and c are almost the same as those shown in Fig. 5. In Fig. 4b, the slip deficit has been translated into actual slip rate. The vectors represent the motion of the subducting plate relative to the upper plate and have been obtained by subtracting the back-slip rates of Fig. 4a from (steady) plate convergence rates defined by the Euler pole described in Section 2. Important features of the inversion results are summarized and discussed as follows: (1)
The slip rate vectors of the most seaward segments off Washington are all within error
(2)
ellipses (Fig. 4b), indicating that within model errors this offshore segment is fully locked. This is consistent with the thermal arguments by Hyndman and Wang [22] and lends support to all previous trial-and-error based dislocation models for northern Cascadia. The degree of locking very gradually decreases landward, as shown by the increasing size of the slip rate vectors (Fig. 4b). From the fully locked region, the slip rates decrease northeastward. The pattern is in qualitative agreement with the recent 3-D dislocation model CAS3D2 [4] that used a more gradual decrease than earlier models constrained by mostly pre-GPS data [23]. Comparison with CAS3D-2 along profiles A–AV and B–BV (Fig. 4b) is made in Fig. 6. The unphysical negative, i.e., normalfaulting type, slip seen in the most seaward segment of profile B–BV is an artifact of the inversion, reflecting the diminished resolving power of the GPS data for the far offshore part
50˚N
48˚N
obs. 10 mm/yr 1σ =2 mm cal. 10 mm/yr
100 km 128˚W
126˚W
124˚W
122˚W
Fig. 5. Observed (solid arrows) and calculated (open arrows) horizontal velocities predicted by the model shown in Fig. 2a. The thick black and thin gray solid arrows represent permanent and campaign data, respectively.
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
40
70
30
60
20
50
10
40
0
30
-10
20
-20
10
-30
0
Slip Rate (mm/yr)
B
Depth (km)
A'
B'
40
70
30
60
20
50
10
40
0
30
-10
20
-20
10
-30
Depth (km)
Slip Rate (mm/yr)
A
245
0 0
50
100
150
200
250
Distance From Deformation Front (km) Fig. 6. Comparison of slip rates along two profiles shown in Fig. 4b with the CAS3D-2 model (gray line). Dashed lines represent 2r errors for the inversion results. The dotted lines represent depths to the plate interface along the cross sections.
(3)
(4)
of the plate interface. The lack of resolution is reflected in the large error bars (and error ellipses in Fig. 4) of the slip values in this area. This is a common problem in defining offshore fault motion/locking from land-based geodetic observations. The slip rate vectors in the southernmost region of the model are less northerly than plate convergence (Fig. 4b), or equivalently, the slip deficit vectors are more northerly than plate convergence (Fig. 4a). We think this is due to the effect of the forearc motion mostly south of our study area (see Section 1). GPS stations in northern Washington are slightly affected by this motion and move in a more northerly direction than plate convergence [24]. This effect is translated into a more northerly slip deficit. Along the updip end of the megathrust, the degree of locking is fairly uniform south of 488N but becomes smaller further north (Fig.
4b). There does not seem to be a lack of data constraint in that area. The age of the subducting plate decreases northward, and the higher temperature of the younger incoming plate may result in aseismic slip of the plate interface all the way to the seafloor deformation front. Alternatively, the locked zone may be extremely narrow and far offshore because of the higher temperature but unresolved by land GPS data and our model parameterization. In addition, the subduction of the Explorer plate just north of the Nootka fault (Fig. 1) is much slower [25], and may affect deformation on land just south of the fault. In cross-section view along profile A–AV (Fig. 6), the landward increase in slip rate towards full plate convergence rate is more gradual than assumed in the CAS3D-2 model. This result is consistent with findings by Mazzotti et al. [9], who related this behavior to possible long-term upper plate deformation.
246
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
5. Discussions One technical point that deserves consideration is the smoothness constraint for the inversion. The results will be smoother if correlation between parameters is introduced, such that neighboring segments have a tendency to have similar slip rates. The G matrix, which is composed of the second derivative of slip distribution with respect to two coordinates taken on the curved plate interface, introduces cross correlation only for neighboring segments, and the effect on smoothness is very small. The inversion results will also be smoother if individual parameters are more tightly constrained relative to data variance. This is determined by the hyperparameter a 2. Although the optimal values (0.00798, 0.00714, and 0.00409 for the models shown in Fig. 2a, b, and c, respectively) of a 2 are automatically determined using ABIC, it is useful to investigate the sensitivity of our results for changes in a 2. For this purpose, we conducted inversion tests by fixing a 2 values without using the ABIC. Using an a 2 value smaller than the optimal value by a factor of five resulted in better data fit but more rapid spatial changes in both slip directions and rates than in Fig. 4a. Since the data include observation errors, smaller misfit does not necessarily indicate better solution. The effects on the misfit and smoothness were the opposite if we increased the a 2 value by a factor of five. However, the difference between the results of both these tests and the results of Fig. 4a was less than 20% in either rates or directions. The relatively low sensitivity of our inversion results to the smoothing constraint can be attributed to the high density and internal consistency of the GPS data. When interpreting the results, it is important to remember the limitations of the elastic half-space model assumed for the inversion. The upper mantle of active margins is viscoelastic, with a viscosity commonly believed to be around 1019 Pa s (see review by Wang [26]). Coseismic deformation and some short-term transients such as the episodic silent slips [15] can be described using elastic models. For interseismic deformation, especially at a time as long as 300 yr after a great earthquake, the effect of viscoelastic stress relaxation may become important. In elastic models, effects of stress relaxation on surface deformation have to be attributed to equiv-
alent slip [27]. For this reason, the zone of downdip transition from full locking to full slip is called the beffective transition zoneQ [4]. It may represent some actual creep along the plate interface, but it also accounts for the effect of stress relaxation in the upper mantle. An equivalent fault slip distribution implicitly represents a rheological structure and the integrated effect of stress relaxation since the last earthquake, but we are generally not in the position to determine uniquely the rheology and deformation history from the limited geodetic data at Cascadia [26]. The modern geodetic observations and in a sense the elastic model provide a bsnapshotQ of an evolving deformation field.
6. Conclusions We have inverted GPS data in the northern Cascadia using an ABIC-based inversion algorithm to infer the locking state of the plate interface in the subduction zone. The GPS data include velocities derived at 20 continuously monitoring stations in British Columbia and northern Washington and 53 campaign sites on Vancouver Island. The inversion results support the previous conclusion based on thermal arguments and trial-anderror dislocation models that the shallower plate interface is fully locked and that the locked zone is mostly offshore. The results are also consistent with recent models in that the downdip transition from fully locked to full slip takes place over a very wide area, but this transition should be understood to be mostly an equivalent description of the effect of stress relaxation that is not directly considered by the elastic model. The inversion also reveals some other complications, such as a tapering of the degree of locking near the northern end of the subduction zone. Northward motion of the Cascadia forearc takes place mostly south of our study area, but it may still slightly affect the southernmost segments of our model.
Acknowledgments Continuous GPS data used in this study are from the Western Canada Deformation Array (WCDA) and northern stations of the Pacific Northwest Geodetic
S. Yoshioka et al. / Earth and Planetary Science Letters 231 (2005) 239–247
Array (PANGA). Campaign GPS surveys were carried out by the Geodetic Survey Division of Geomatics Canada. We thank T. Yabuki for sharing his original inversion source code, P. Wessel and W. H. F. Smith for making available the Generic Mapping Tool (GMT), with which the figures for this paper were created, and E. Hearn, E. Hetland, H. Dragert, and the editor, R. van der Hilst for their valuable comments. The work was initiated when SY was a visiting scientist at the Pacific Geoscience Centre in 2001 supported by Japan Society for the Promotion of Science. Geological Survey of Canada contribution 2004228.
References [1] K. Satake, K. Wang, B.F. Atwater, Fault slip and seismic moment of the 1700 Cascadia earthquake inferred from Japanese tsunami descriptions, J. Geophys. Res. 108 (B11) (2003) 2535, doi:10.1029/2003JB002521. [2] G. Khazaradze, A. Qamar, H. Dragert, Tectonic deformation in western Washington from continuous GPS measurements, Geophys. Res. Lett. 26 (1999) 3153 – 3156. [3] M. Miller, D.J. Johnson, C.M. Rubin, H. Dragert, K. Wang, A. Qamar, C. Goldfinger, GPS-determination of along-strike variation in Cascadia margin kinematics: implications for relative plate motion, subduction zone coupling, and permanent deformation, Tectonics 20 (2001) 161 – 176. [4] K. Wang, R.E. Wells, S. Mazzotti, R.D. Hyndman, T. Sagiya, A revised dislocation model of interseismic deformation of the Cascadia subduction zone, J. Geophys. Res. 108 (2003), doi:10.1029/2001JB001227. [5] R.E. Wells, R.W. Simpson, Northward migration of the Cascadia forearc in the northwestern U.S. and implications for subduction deformation, Earth Planets Space 53 (2001) 275 – 283. [6] R. McCaffrey, Crustal block rotations and plate coupling, in: S. Stein, J. Freymueller (Eds.), Plate Boundary Zones, AGU Geodynamics Series, vol. 30, 2002, pp. 101 – 122. [7] S. Yoshioka, T. Yabuki, T. Sagiya, T. Tada, M. Matsu’ura, Interplate coupling and relative plate motion in the Tokai district, central Japan, deduced from geodetic data inversion using ABIC, Geophys. J. Int. 113 (1993) 607 – 621. [8] S. Yoshioka, T. Yabuki, T. Sagiya, T. Tada, M. Matsu’ura, Interplate coupling in the Kanto district, central Japan, deduced from geodetic data inversion and its tectonic implications, Tectonophysics 229 (1994) 181 – 200. [9] S. Mazzotti, H. Dragert, J. Henton, M. Schmidt, R. Hyndman, T. James, Y. Lu, M. Craymer, Current tectonics of northern Cascadia from a decade of GPS measurements, J. Geophys. Res. 108 (B12) (2003) 2554, doi:10.1029/2003JB002653.
247
[10] C. DeMets, T.H. Dixon, New kinematic models for the Pacific-North America motion from 3 Ma to present: I. Evidence for steady motion and biases in the NUVEL-1A model, Geophys. Res. Lett. 26 (1999) 1921 – 1924. [11] D.S. Wilson, Confidence intervals for motion and deformation of the Juan de Fuca plate, J. Geophys. Res. 98 (1993) 16053 – 16071. [12] Z. Altamimi, P. Sillard, C. Boucher, ITRF2000: a new release of the international terrestrial reference frame for earth science applications, J. Geophys. Res. 107 (B10) (2002) 2214, doi:10.1029/2001JB000561. [13] J.A. Henton, GPS Studies of Crustal Deformation in the Northern Cascadia Subduction Zone, PhD thesis, University of Victoria, British Columbia (2000). [14] A. Mao, C.G.A. Harrison, T. Dixon, Noise in GPS coordinate time series, J. Geophys. Res. 104 (1999) 2797 – 2816. [15] H. Dragert, K. Wang, T.S. James, A silent slip event on the deeper Cascadia subduction interface, Science 292 (2001) 1525 – 1528. [16] J.C. Savage, A dislocation model of strain accumulation and release at a subduction zone, J. Geophys. Res. 88 (1983) 4984 – 4996. [17] T. Maruyama, Statical elastic dislocations in an infinite and semi-infinite medium, Bull. Earthq. Res. Inst. Univ. Tokyo 42 (1964) 289 – 368. [18] H. Akaike, Likelihood and the Bayes procedure, in: J.M. Bernardo, M.H. DeGroot, D.V. Lindley, A.F.M. Smith (Eds.), Bayesian Statistics, University Press, Valencia, Spain, 1980, pp. 143 – 166. [19] T. Yabuki, M. Matsu’ura, Geodetic data inversion using a Bayesian information criterion for spatial distribution of fault slip, Geophys. J. Int. 109 (1992) 363 – 375. [20] P.A. McCrory, J.L. Blair, D.H. Oppenheimer, S.R. Walter, Depth to the Juan de Fuca slab beneath the Cascadia subduction margin: a 3-D model for sorting earthquakes, U.S. Geological Survey Digital Data Series, 1 CD-ROM (2004). [21] D.D. Jackson, M. Matsu’ura, A Bayesian approach to nonlinear inversion, J. Geophys. Res. 90 (1985) 581 – 591. [22] R.D. Hyndman, K. Wang, Thermal constraints on the zone of major thrust earthquake failure: the Cascadia subduction zone, J. Geophys. Res. 98 (1993) 2039 – 2060. [23] P. Flqck, R.D. Hyndman, K. Wang, Three-dimensional dislocation model for great earthquakes of the Cascadia subduction zone, J. Geophys. Res. 102 (1997) 20539 – 20550. [24] S. Mazzotti, H. Dragert, R.D. Hyndman, M.M. Miller, J.A. Henton, GPS deformation in a region of high crustal seismicity: N. Cascadia forearc, Earth Planet. Sci. Lett. 198 (2002) 41 – 48. [25] R. Riddihough, Recent movements of the Juan de Fuca plate system, J. Geophys. Res. 89 (1984) 6980 – 6994. [26] K. Wang, Elastic and viscoelastic models for subduction earthquake cycles, in: T. Dixon (Ed.), The Seismogenic Zone Experiment, in press. [27] J.C. Savage, Viscoelastic-coupling model for the earthquake cycle driven from below, J. Geophys. Res. 105 (B11) (2000) 25525 – 25532.