Viscoelastic-cycle model of interseismic deformation in the northwestern United States

convergence of the block and of the de Fuca mechanical that


I N T RO D U C T I O N
In this study, we construct a kinematic model of crustal deformation in the Pacific Northwest. The modelling methodology incorporates deformation from different sources: time-dependent post-seismic viscoelastic relaxation, time-independent (cycle-averaged) relaxation, and fault creep, and it is constrained by GPS and fault data. Our results allow us to estimate the partitioning of slip at the Mendocino Triple Junction and western Washington, evaluate the coseismic slip distribution of the 1700 Cascadia-megathrust earthquake, and evaluate the effects of lateral variations in crustal rigidity and mantle viscosity. The results illuminate the physical mechanisms shaping the interseismic crustal velocity field and highlight the complexity of rationalizing the velocity field with a combination of well constrained and poorly constrained dislocation sources, as seen through the filter of the crust and mantle rheology.
The western North America Plate boundary zone ( Fig. 1) accommodates right-lateral shear along the ∼1300-km-long San Andreas fault (SAF) system along the Pacific-North America Plate boundary and east-west convergence along the ∼1200-km-long Cascadia megathrust (CSF) along the Juan de Fuca-North America Plate boundary. The plate boundary zone generally extends several hundred kilometres into North America through the Eastern California Shear Zone (ECSZ), the extensional system of the Basin and Range, and the backarc fold and thrust belt of eastern Washington. Rationalizing the deformation in the rapidly deforming SAF and CSF, the slowly deforming plate interior, and the transitions among them present special challenges. It is possible to address many of them using the detailed images of the crustal velocity field now available.
The crustal velocity field shown in Figs 2 and 3 consists of 3650 horizontal velocity vectors derived over two decades of observation using the Global Positioning System (GPS). Although there is considerable redundancy, particularly in southern California, the interseismic crustal velocity field is imaged almost everywhere to a resolution of 50 km or less, with densest coverage near major fault zones at a resolution of a few km. This velocity field and the associated strain field (Fig. 4) span spatial scales of 10 1 to 10 3 km. At the shortest scale of ∼10 km, it informs us on the degree of strain accumulation around fault zones, which are related to longterm slip rates and elastic plate thickness (Savage & Burford 1973;Savage & Prescott 1978). At the largest scale of ∼10 3 km, it informs us on the long-term slip rates of the major fault zones in the plate boundary zone (Fig. 1) and the major faults that bound the western North American Plate (e.g. the Queen Charlotte transform and the Gulf of California transform and spreading system) and the forces driving slow deformation of the plate interior (Elsasser 1969). At intermediate scales of 10 1 -10 3 km, it informs us on the flow pattern produced by repeated slip events on faults, which are related to long-term slip rates, effective elastic plate thickness (and lateral variations thereof), and the viscosity structure of the ductile lower crust and upper mantle (Savage & Prescott 1978;Thatcher & Rundle 1979;Cohen 1982;Thatcher 1983;Dixon et al. 2002Dixon et al. , 2003Meade & Hager 2004).
As GPS velocity fields within continents have become more comprehensive, models of crustal deformation have been correspondingly refined. Continuum models were initially advanced to explain the long-wavelength component of the crustal velocity field in terms of a viscous lithosphere with a (typically large) viscosity designed to mimic distributed deformation in continents (England & McKenzie 1982;England & Molnar 1997;Flesch et al. 2000). The recognition that major faults tend to localize deformation and regions between the faults tend to behave rigidly at nearly all timescales led to the development of block models which explicitly account for fault locking effects (Hashimoto & Jackson 1993;Thatcher 2003). At the same time, the observation of strong time dependence in the crustal deformation following large (M 7) earthquakes compelled the development of viscoelastic models to simulate stress re-adjustment of the lower crust and mantle following these events (Nur & Mavko 1974;Savage & Prescott 1978;Thatcher & Rundle 1979Cohen 1982). More detailed sampling of the postearthquake deformation field, particularly after the 1992 Landers and 1999 Hector Mine, California, earthquakes, allowed the refinement of viscoelastic models in terms of viscoelastic structure as well as rheology (Maxwellian versus transient; linear versus nonlinear) (Pollitz et al. 2000(Pollitz et al. , 2001Pollitz 2003a;Freed & Bürgmann 2004;Freed et al. 2007;Johnson et al. 2006). The viscoelastic structures derived from post-earthquake geodetic studies in the western US agree well with those derived from palaeo-shoreline analyses over longer time periods (e.g. Bürgmann & Dresen 2008;Thatcher & Pollitz 2008;Hammond et al. 2009). These models are consistent with the expected reduction in rock strength with increasing temperature. They also place block models in a broader context because block models are theoretically an end-member case of viscoelastic models in the limit of high sublithosphere (asthenosphere) viscosity (Savage 1983;. The utility of the viscoelastic-cycle model extends beyond explaining timedependent crustal deformation following large earthquakes, and it is generally applicable to the interseismic velocity field (e.g. Thatcher & Rundle 1979;Thatcher 1983;Savage & Lisowski 1998;Savage 2000;Dixon et al. 2003;Hetland & Hager 2003;Johnson et al. 2006;Smith & Sandwell 2006;Pollitz et al. 2008;Hearn et al. 2009).

671
The purpose of this study is threefold.
(1) To update and expand the GPS velocity field of Pollitz et al. (2008) with additional observations, drawing on newly available sources. This includes crustal velocities at sites of the Plate Boundary Observatory (PBO) which benefit from longer observation times than were available in the earlier study. In doing so we more than triple the amount of data considered by Pollitz et al. (2008).
(2) To provide a better treatment of the deformation sources in the Pacific Northwest, which were explored incompletely by Pollitz et al. (2008). In this study, we revise the kinematic parameters of sources along the Cascadia megathrust, benefitting from the constraints on net relative motions deduced by McCaffrey et al. (2007), add sources associated with the Explorer plate boundary, and incorporate numerous onland faults which contribute minor, but substantial, signals to the observed crustal velocity field.
(3) To test whether the viscoelastic-cycle approach of Pollitz et al. (2008) is capable of rationalizing the expanded data set, which captures more detail over a greater breadth of tectonic domains than the earlier data set.
The actively deforming northwestern US includes the Mendocino triple junction, where the style of faulting changes from northwestdirected shear to east-west shortening along the Cascadia megathrust; northwestern California and western Oregon and Washington, which record this east-west shortening; western Washington, which is a diverse deformation zone accommodating the northward motion of the Oregon Coast block and northwestward convergence of the Juan de Fuca Plate. The regional deformation has been rationalized with a block model , which accounts well for the budget of slip along the major tectonic boundaries and numerous minor fault zones. Here we use the model of Pollitz et al. (2008) in which interseismic crustal deformation is rationalized as the product of viscoelastic relaxation from past slip events on identified faults and steady distributed deformation in the regions surrounding the faults, with further allowance for the effects of lateral variations in depth-averaged rigidity and creeping faults. The use of a viscoelastic model necessitates consideration of GPS data over a very large domain because of the long-range effects of viscoelastic deformation cycles, hence the use of a velocity field spanning the western US (Figs 2 and 3). This velocity field affords the opportunity to also examine the nature of strain accumulation in the continental interior and first-order characteristics of the viscoelastic rheology.
The USGS campaign measurements are described in numerous prior publications (Savage et al. , 1999a(Savage et al. ,b, 2001aThatcher et al. 1999;Prescott et al. 2001;Svarc et al. 2002a,b;Savage et al. 2004;Hammond & Thatcher 2004). They are generally conducted at intervals of 3-4 yr, and the associated velocity field is a composite of such measurements conducted between 1991 and 2007. The Payne et al. (2008) Payne et al. (2008), for a total of 3650.
In general, all contributing data sets have been processed in slightly different realizations of fixed North America. Any possible disparities in reference frames are corrected to first order by referring each individual data set to the PBO data set at common sites using a Helmert transformation. Each Helmert transformation is parametrized by the three Cartesian components of an Euler vector, and the derived transformation is applied to all velocity vectors of the data set. The non-PBO velocity vectors shown in Figs 2 and 3 have already been corrected in this manner. Fig. 5 shows the mean misfit between the PBO velocity field and three selected velocity fields prior to and after alignment. The rotation prescribed by each respective Helmert transformation reduces the initial misfit substantially.

T H E O R E T I C A L F R A M E W O R K
The driving forces of earthquakes originate in the convective system of Earth's interior, which generate toroidal flow fields and basal tractions on the lithosphere (e.g. Humphreys & Coblentz 2007). From a kinematic viewpoint, however, the effects of long-lived driving forces are embraced by the history of earthquake faulting extending indefinitely into the past (Savage & Prescott 1978). In this context, crustal deformation is envisaged to be driven by the post-earthquake relaxation of a ductile lower crust and mantle underlying an elastic upper crust. The treatment allows for both cyclic and steady-state dislocation sources in the upper crust. We employ eq. (1) of Pollitz et al. (2008) v inst (r) = n n d 3 r m(r ) : (1) Here V refers to the volume of the lithosphere, which is assumed to be populated with discrete fault surfaces, and G (d) (r, r , t) is the response of the viscoelastic system at point r and time t to a unit dislocation source applied at point r and time 0. Quoting Pollitz et al. (2008), the terms of eq. (1) represent: (1) Viscoelastic relaxation from all known/estimated past major regional earthquakes. Letting n define the nth (discrete) fault surface, fault geometry and slip of these events are represented through the moment-release rate densityṁ(r ) at points r ∈ n . Time of last event and recurrence interval on nth fault are t n and T n , respectively.
(2) Interseismic-cycle averaged velocity produced by viscoelastic relaxation from moment-release rate density on faults m .
(3) Interseismic-cycle averaged velocity produced from moment release on dislocations distributed throughout the remaining volume.
(4) Secular deformation arising from steady creep at points r ∈ cr corresponding to creeping fault surfaces.
(5) The effects of lateral heterogeneity in shear modulus δμ(r ) and bulk modulus δκ(r ) at points r ∈ V .
Note that deformation associated with the second and third categories is time-independent, and that the moment rate density tensoṙ m (fault) (r ) is proportional to the slip rateṡ m on fault surface m . Under the assumptions of time-independence and negligible depth variation in crustal strain rate, the moment release rateṁ (δμ) (r ) associated with deformation of the fifth category depends on δμ(r ), δκ(r ), and the observed strain rate field using eq. (8) of Pollitz et al. (2008).
We employ a time-dependent model for larger-magnitude sources for which sufficient information (i.e. slip estimates and date of occurrence) is available, and a time-independent model for remaining sources. To implement the time-dependent model, we use the viscoelastic structure shown in fig. 7 of Pollitz et al. (2008). It consists of a 20-km-thick elastic lithosphere underlain by viscoelastic lower crust and mantle. The lower crust is assumed linear Maxwellian with a viscosity of 2 × 10 19 Pa s, and the upper mantle is a Burgers body (e.g. Pollitz 2003a) with transient and steady state viscosities of 5 × 10 17 and 10 19 Pa s, respectively.
The time-independent model (category 2 source in eq. 1) is a valid approximation to a time-dependent model provided that the material relaxation time of the asthenosphere is larger than the mean recurrence interval of the fault (e.g. Savage & Prescott 1978). With a relaxation time of 5 yr for the steady-state component of the mantle rheology and typical recurrence intervals greater than 100 yr, the interseismic velocity would vary considerably with time and thus the approximation is not strictly valid. On the other hand, if such faults are loaded not by relaxation bur rather by steady slip at depth at the long-term slip rate, then the employed term in eq. (1) is applicable. Moreover, if a deforming zone is occupied by several faults, each at a different stage into its respective seismic cycle, then the resultant interseismic velocity field will statistically tend to that predicted by the time-independent model (Pollitz 2003b).
The time-dependent sources are summarized in table 3 of Pollitz et al. (2008) and shown in Fig. 6. These include the Cascadia megathrust which ruptured in 1700 and parts of the San Andreas fault system which ruptured in 1906 (northern California) and 1857 (southern California), and they are modelled with deformation of category 1. We considered in addition several large earthquakes which occurred around the MTJ, including  Fig. 3, resolved in the direction parallel to the profile, that is, in the N67 • E direction on profile CD and in the N30 • E direction on profile AB. Velocity is averaged over those GPS sites within 20 and 40 km of the profile for profiles AB and CD, respectively, and within a ±15 km distance from the plotted point in the direction parallel to the profile. Solid curves are observed velocity with interpolation (i.e. Fig. 2); long dashed curves are observed velocity after correction for the effects of the viscoelastic cycle on the megathrust using the Cascadia slip rates presented in Table 2; short dashed curves are the residual velocity field (i.e. Fig. 10).
No claim to original US Government works, GJI, 181, 665-696 Journal compilation C 2010 RAS the 1992 April 25 M w = 7.1 and 1994 September 1 M w = 7.0 Cape Mendocino earthquakes. The post-earthquake relaxation signals from any of the largest events (based on the assumed viscoelastic model) do not exceed 0.5 mm yr −1 in the year 2000 and diminish with time, and they are therefore considered negligible compared with other deformation sources.
All other discrete faulting sources are modelled as timeindependent deformation of categories 2 and 4. In order to account for distributed faulting and effective lateral variations in elastic parameters, we append deformation of categories 3 and 5.
Fault #86 is meant to accommodate NE-SW horizontal shortening within the Olympic Mountains. This region has rapid short-term uplift (Savage et al. 1991), and the crustal velocity field obtained after correcting for elastic locking effects on the megathrust has several mm yr −1 residual horizontal contraction parallel to the plate convergence direction (Mazzotti et al. 2002). This is true for the present model as well, which exhibits about 14 mm yr −1 net shortening across the Olympic Peninsula, about 8 mm yr −1 of which remains after correction for the effects of the viscoelastic cycle on the megathrust (profile CD of Fig. 7). Fault #86 is a crude approximation to a likely broad zone of shortening and uplift of the region, which may be frontally accreting and thus concentrated in the western Olympic Peninsula (Pazzaglia & Brandon 2001). Fault #93 represents a collection of thrust faults mapped on the Washington coast between Grays Harbor and the Quillayute River (near 47 • N) accommodating northward shortening (McCrory 1996;McCrory et al. 2002;McCaffrey et al. 2007). The onland structures extend considerably offshore to the shelf part of the accretionary margin. The chosen dislocation plane represents the western portion of a diffuse boundary between the northward-translating OC block and the North American Plate. Faults #87 and 88 follow well-defined lineaments of seismicity. Fault #87 is the Hood Canal-Discovery Bay fault, which is expressed chiefly through sharp geophysical discontinuities and appears to deform Quaternary deposits (Haug 1998;Johnson et al. 1999). Fault #88 may accommodate rightlateral strike slip based on focal mechanisms .
The Tacoma and Seattle fault zones (faults #89-90) are part of a broad zone of N-S contraction between the OC block and Vancouver Island and the Coast Range of British Columbia, which is considered to behave as a rigid backstop. This zone accommodates about 4 mm yr −1 N30 • E contraction, of which about 3 mm yr −1 remains after removal of the viscoelastic cycles on the megathrust (profile AB of Fig. 7). We model the zone using the south-dipping floor thrust and north-dipping back thrust of Brocher et al. (2004) which bound the Seattle Uplift between the Tacoma and Seattle fault zones. This overall zone of convergence may include other active faults such as the Olympia fault . Slip inferred on the idealized Tacoma and Seattle faults may represent, in part, slip on other components of this plate boundary zone. The amount of modelled slip on the summed faults #89 and 90 (which are dipping 45 • ) is constrained to be 5 mm yr −1 , which corresponds to summed horizontal slip rate of 3.5 mm yr −1 .
With ω given by the appropriate angular velocity vector from eqs (2) to (4), the long-term relative motion rate and azimuth at a point r representative of the footwall of a fault segment on the megathrust is theṅ The slip ratesṡ m on the megathrust provided by eq. (5) Since, for a given megathrust fault segment, the two parameters under the constraint in eq. (6) have only one degree of freedom, we may reduceṡ seismic m andṡ aseismic m to one free parameter by supposing that (1) the fault 'creeps' at the rateṡ m and (2) the fault has an additional slip rateṡ seismic m that is released in periodic seismic events and is compounded by steady creep at the rate −ṡ seismic m .
With this simplification, model parameters includeṡ seismic m on the megathrust andṡ m on other faults. We assume an average recurrence interval of Cascadia megathrust earthquakes of T = 500 yr (Adams 1990;Atwater & Hemphill-Haley 1997). This treatment is an improvement with respect to the slip model of Pollitz et al. (2008) (their fig. 16) because that study did not account for the aseismic component of the total slip budget on the megathrustṡ aseismic m . It has been noted by Leonard et al. (2004) that slip in the 1700 Cascadia earthquake may have exceeded the slip accumulated in an average ∼500-800-yr average earthquake cycle. In trial models we have allowed for an 'excess' slip above that which would be accumulated in a 500 yr earthquake cycle and found that such slip No claim to original US Government works, GJI, 181, 665-696 Journal compilation C 2010 RAS is not needed to fit the interseismic velocity field, nor does the introduction of such slip improve the fit.

I N V E R S I O N F O R D E F O R M AT I O N PA R A M E T E R S
For a given viscoelastic structure the interseismic velocity field depends upon a combination of slip ratesṡ m (for faults implemented with the time-independent relationship, i.e. category 2 and 3 sources) and characteristic slip values s m for faults implemented with the time-dependent relationship (i.e. category 1 sources). In all models, table 3 source faults-excluding the Cascadia megathrust faults (#30-45) and 'VE' faults-are implemented with a timeindependent relationship between slip rate and crustal velocity. Among these, source faults #11, 14-15 and 86 are implemented as category 4 sources and the remainder as category 2 sources. The Cascadia megathrust faults (#30-45) and 'VE' faults listed in table 3 of Pollitz et al. (2008) are implemented with a time-dependent relationship between recurring earthquake slip and subsequent crustal velocity, that is, category 1 sources. The time-dependent crustal deformation field depends on the time elapsed since the last major event and the recurrence interval, which are given in that table.
The distributed moment releaseṁ (V) (r ) (category 3 source) and lateral variations in rigidity δμ(r ) (category 5 source) are expanded in terms of a set of smooth basis functions covering portions of the western US . The expansion coefficients,ṡ seismic m , and s m are estimated by linear inversion of the GPS velocity field using the methodology given in section 4 of Pollitz et al. (2008). Several of the slip rates are fixed at independently determined values in areas where the GPS data is deemed insufficient to estimate a slip rate .

R E S U LT S
In all models, the slip on the Cascadia megathrust is partitioned into seismic and aseismic slip subject to the relationship and inequalities in eq. (6). On a given segment, in instances whereṡ seismic m consistently exceeds these bounds in trial inversions (i.e. slip rate less than zero or greater than the long-term slip rateṡ m ),ṡ seismic m is held fixed at the endpoint value (zero orṡ m ).
In the MTJ area (Fig. 8), both right-lateral slip and reverse slip are allowed on the Little Salmon, Russ, and Mad River faults (faults #2, 2a and 3). In trial inversions, we find that non-negative slip rates are realized for the dip-slip component of the Mad River fault zone and strike-slip components of the Little Salmon and Russ faults. Consequently, the dip-slip component of the Little Salmon and Russ faults and the strike-slip component of the Mad River fault are fixed at zero. Inversions without restriction on the strike slip of the Little Salmon fault and Russ fault or reverse slip of the Mad River fault yield ≥10 mm yr −1 corresponding rates. Therefore, we assign 10 mm yr −1 of strike slip on the Little Salmon fault and Russ fault and 10 mm yr −1 reverse slip on the Mad River fault. All assigned rates are larger than the 4-5 mm yr −1 geological slip rate estimated for the Mad River fault (McCrory 2000), but it reflects the tendency of inversions to distribute >30 mm yr −1 strike slip on the entire system of faults. It is possible to vary the distribution of right-lateral slip on the various faults, and trade-offs with the right-lateral slip rate of the Russ fault are explored in Section 7.2.
We refer to the resulting model, consisting of slip-rate values and distributions of δμ and vertically integratedṁ (V ) as the Western US (WUS) deformation model. Inverted slip-rate parameters are listed in Table 2, and Fig. 9 shows the resulting distribution of δμ. The distribution of vertically integratedṁ (V ) is similar to that presented and discussed in section 8.3 of Pollitz et al. (2008). Fig. 10 shows the residual fit of the model to the GPS data set. Root-mean-square residuals are 2.1 mm yr −1 .
The pattern of observed and modelled strain rate fields are shown in Fig. 11, where the strain field is depicted as the second invariant of the strain rate field I 2 (Jaeger & Cook 1984). The WUS deformation model replicates the observed strain field in most areas, including the SAF system, ECSZ, and Walker Lane, with a small residual strain rate field (Fig. 11c).

Lateral Variations in Rigidity
Lateral variations in rigidity (Fig. 9) are largest along parts of the SAF, ECSZ and ISB. Low rigidity values could be interpreted equally well as steady slip along the fault zone (averaged over the depth extent of the elastic layer) or reduced effective elastic plate thickness. The former interpretation is supported by the correlation of these zones with seismicity ( fig. 3 of Pollitz et al. 2008), while the latter interpretation is supported by the correlation of these zones with elevated surface heat flow (e.g. Humphreys et al. 2003) and low seismic velocities in the uppermost mantle (Pollitz 2008). Pollitz et al. (2008) notes that the pattern may be indicative of Rigidity perturbations are also associated with elevated strain rates ( fig. 4 and table 6 of Puskas & Smith 2009), negative δμ tending to amplify the predicted strain rate regardless of tectonic regime. That is, relatively low rigidity along the ECSZ and Walker Lane tends Figure 12. Model GPS velocity field shown for four separate components: (a) that associated with major-fault earthquake cycles (i.e. those of Fig. 6), including steady slip on parts of the SAF and Cascadia megathrust, (b) lateral variations in rigidity and distributed deformation, (c) ridges and transform faults associated with the Gorda, Juan de Fuca, and Explorer plates, (d) all other faulting sources, (e) an additional rotation to a fixed North America reference frame (eq. 8 of Pollitz et al. 2008) and (f) the total model velocity field, which is the sum of the components in (a)-(e) and is directly comparable with the observed velocity field in (g, which is subsampled from the velocity field in Fig. 2). Note that for visual clarity, the plotted velocity fields greatly subsample the actual velocity fields.
to amplify right-lateral horizontal strain along NW-trending axes, whereas along the ISB it tends to amplify E-W normal faulting with a small amount of counter-clockwise rotation (Fig. 4). This association suggests that the mechanism of laterally varying rigidity is generally needed to produce elevated geodetic strain rates that may not be produced by other means (e.g. post-earthquake relaxation), and it highlights the utility of interseismic strain rates for illuminating this mechanism (Chéry 2008). If steady slip is the causative mechanism in the rigidity pattern (Fig. 9), then such slip is inferred along essentially two continuous zones beginning in the southern ECSZ and continuing along both the ECSZ (to southern Oregon) and the ISB; a gap in δμ along the ISB in northern Arizona is due to the paucity of data in that region.
To illustrate the composition of the net strain rate field in the western US, the model velocity field field is divided into four  separate components [Figs 12a-d-those associated with major earthquake cycles, lateral variations in rigidity plus distributed deformation, i.e. that arising from theṁ (δμ) andṁ (V ) terms of eq. 1], oceanic ridges and transform faults, and all other sources. The total model velocity field in Fig. 12(f) is composed of these components plus an additional rotation (Fig. 12e), and it replicates the observed velocity field (Fig. 12g). The corresponding model strain rate field is shown in Fig. 13. It is noteworthy that combined lateral No claim to original US Government works, GJI, 181, 665-696 Journal compilation C 2010 RAS Figure 14. (a) Observed strain rate field derived from the GPS velocity field (Fig. 3), represented by the amplitudes and directions of the principal strain rate axes (thick and thin line segments denoting a principal contractile or tensile strain rate axis, respectively) and rotation rate (indicated by color shading). (b) Residual strain rate field derived from the residual velocity field (Fig. 10).
variations in rigidity and distributed deformation (Fig. 13b) account for ∼10 mm yr −1 net motion of the west coast with respect to the plate interior (Fig. 12). In terms of strain, they contribute a substantial signal to I 2 in the ECSZ, Walker Lane, and part of the ISB. When these components are omitted, rms residuals increase from 2.1 to 3.3 mm yr −1 (i.e. the residual variance of the velocity field increases by a factor of 2.5) with correspondingly greater I 2 (Fig. 11d). An F-test indicates that the inclusion of these components is significant at nearly the 100 per cent confidence level, similar to the finding of Pollitz et al. (2008).

Mendocino triple junction
The observed strain rate field (Fig. 14a) illustrates the transition from predominantly right-lateral shear on the NNW-trending SAF system to east-west compression along the megathrust. The transition occurs near the northern terminations of the Maacama and Bartlett Springs faults and northward through the Little Salmon and Mad River faults, where contractile principal strain axes encourage dip-slip. In the plate interior where effects of the megathrust are diminished, the deformation consists largely of positive (clockwise) rotation, consistent with the rotation of the Oregon Coast block. Fig. 14(b) shows the residual strain field. The WUS deformation model explains the primary patterns in the observed strain field, including most of the contractile strain and rotation around the MTJ and Humboldt Bay area faults. The slip rates on the northern extensions of the SAF, Maacama and Bartlett Springs faults (e.g. the Russ, Little Salmon and Mad River faults) carry implications for the partitioning of slip at the complex Mendocino triple junction. The WUS deformation model (Table 2) is the case where the Russ fault accommodates 10 mm yr −1 strike-slip motion. Inversions that include its slip rate as an unknown generally yield up to 20 mm yr −1 strike slip motion. We test the trade-off of estimated slip rates with respect to the Russ-fault slip rate by inverting for slip rates on the Little Salmon and Mad River faults for trial values of the Russ-fault slip rate between 0 and 20 mm yr −1 (Fig. 15). The only non-negative slip rates obtained are for strike slip on the Little Salmon fault and reverse slip on the Mad River fault. These and similar inversions without constraints on the slip values consistently yield ∼25 mm yr −1 reverse slip on the Mad River fault, which we deem too large because such a large value would allow little remaining shortening to occur in the megathrust; we therefore fix its reverse slip at 10 mm yr −1 in these inversions. When strike slip on the Russ fault is factored in, the results in Fig. 15   this magnitude are applicable, then the slip may be distributed on numerous NW-trending structures (Fig. 8).
In the WUS deformation model (Russ-fault slip rate fixed at 10 mm yr −1 ), the Mad River fault zone accommodates most of the dip-slip on crustal faults of the southern Gorda convergent margin, and this implies that a significant portion of the net Gorda-CC convergence may be accommodated by crustal shortening in the Humboldt Bay region. This pattern is consistent with the observed velocity and strain fields after removing the contributions of the GDZ and Juan de Fuca and Explorer plates on the respective models (Fig. 16). The plotted velocity field removes the contributions of all dislocations associated with boundaries lying on one or more of these plates, that is, faults #30-46 and 81-84, plus a rigid rotation (determined in the inversion) that adjusts the net model velocity field to a fixed North America reference frame (section 5.3 of Pollitz et al. 2008). Since almost all of the remaining dislocations after this correction are associated with the Pacific-North America Plate boundary, the velocity field shown in Fig. 16 is effectively with respect to a SAF-neutral reference frame, that is, an average of velocity fields that would be obtained in fixed Pacific and fixed North America reference frames, respectively. Contractile NE-SW strain rates of ∼3 × 10 −7 yr −1 are seen around the Mad River fault are explained by dip-slip on the that fault. Right-lateral shear strain rates evaluated on the Little Salmon fault are 1.5 × 10 −7 yr −1 , leading to rates of ∼15 mm yr −1 strike slip on this fault when inverted without an upper limit. There are substantial shear strain rates ∼2 × 10 −7 yr −1 evaluated on the Russ fault, which leads to 10 mm yr −1 rates of estimated right-lateral strike slip. Combined strike slip and dip-slip on these models is sufficient to explain most of the observed crustal deformation around Humboldt Bay, as evidenced by the small residual tensor strain and rotation (Fig. 14).
The preferred slip rate of ∼10 mm yr −1 on the Mad River fault zone is robust with respect to alternative seismic and aseismic slip rates on the megathrust, for example, in tests with variableṡ m anḋ s seismic m on megathrust segment #30b. This GPS-derived rate is larger than geological estimates indicating ∼ 4 − 5 mm yr −1 across the fault zone (McCrory 2000). Much of the additional reverse slip may be accommodated in the offshore region where numerous thrust structures are mapped (Carver 1992;Burger et al. 2002;Gulick et al. 2002).

Cascadia megathrust
The coseisemic slip distribution of the 1700 Cascadia earthquake may be estimated for megathrust fault #m as whereṡ m is the seismic component of slip on the Cascadia megathrust for segments #30-45 (Table 1), the recurrence interval T is 500 yr. The resulting slip distribution is shown in Fig. 17. The seismic component of slip on the Cascadia megathrust is given byṡ seismic m for faults #30-45 and 81. The relatively low values at latitude ∼43.7-46 • N indicate relatively little accumulation of seismic strain off Oregon, in good agreement with independent inferences (Mitchell et al. 1994;Wang et al. 2003;Pollitz et al. 2008;Burgette et al. 2009). The formal uncertainties inṡ seismic m all lie between 0.5 and 1 mm yr −1 , translating into corresponding uncertainties in 1700 slip of 0.25-0.5 m. Due to trade-offs with other model parameters, particularly viscoelastic structure, this slip distribution is highly uncertain. A robust feature of the model, however, is the large coseismic slip demanded for faults #42 and 43, without which the large interseismic contraction observed in most of northwestern Washington (e.g. Fig. 3) cannot be matched. Tests also indicate that end-member slip distributions yield relatively poor fits to the GPS data set. For example, ifṡ seismic m are constrained to equal the relative plate motion rateṡ m (i.e. the megathrust is fully locked), then average root-mean-square misfit increases from 2.1 to 2.5 mm yr −1 . If this fully locked model is regarded as a reference model, then an F-test indicates that the preferred slip distribution (i.e. Table 2, with variable ratioṡ seismic m /ṡ m ) is significant at nearly 100 per cent confidence.
It is noteworthy that coseismic slip is concentrated in the upper half of the slab interface south of 44 • N and in the lower half of the slab interface between 44 • N and 47 • N. This could reflect first-order differences in the age of the subducted slab (e.g. Wilson 2003), the volume of sediment being carried with the descending slab, and the properties of the overriding Silezia Plate between 44 • N and 47 • N (Wells et al. 1998). Based on the position of the coast relative to the upper and lower halves of the slab interface, the coseismic slip model would predict substantial coseismic coastal subsidence from southern Oregon (Cape Blanco) to northern Washington. Coastal subsidence predicted by the 1700 coseismic slip model is shown in Fig. 18, where it is compared with estimates of coseismic subsidence off the coasts of Washington, Oregon and Vancouver Island (Leonard et al. 2004). However, predicted subsidence is about twice as great as observed subsidence. This reflects the fact that most of the coast is located directly above the downdip edge of modelled slip. Predicted subsidence is very sensitive to the location of the lower edge of a number of segments. Between 44 • and 46 • N, actual slip (on faults #36, 38, 40) need not extend as deep as 20 km. Similarly, slip off Cape Blanco (fault #35) need not be restricted to the upper 10 km. We construct a revised coseismic slip model in which the slip on faults #36, 38 and 40 is moved updip by 20, 15 and 12 km, respectively, and slip on fault #35 is moved downdip by 5 km; the width of the slipping regions and slip values are unchanged. In this revised model (dashed curve in Fig. 18, the predicted subsidence agrees roughly with the observation. Moreover, these revisions have negligible impact on predicted interseismic motions. In detail, the interpretation of coastal subsidence measurements is hampered by uncertainties in the corrections made of raw measurements (global Figure 16. Left-hand panel: velocity field obtained after correcting the observed GPS velocity field (Fig. 3) for the effect of deformation associated with all GDZ, Juan de Fuca, and Explorer plate boundaries. The sources that contribute to the correction are faults #30-46 and 81 of Table 1. Right-hand panel: strain rate fields corresponding to the plotted velocity fields, represented by the amplitudes and directions of the principal strain rate axes (thick and thin line segments denoting a principal contractile or tensile strain rate axis, respectively) and rotation rate (indicated by color shading). It is derived from the velocity field using the velocity-gradient determination method described in appendix A of Pollitz & Vergnolle (2006). sea level rise; tectonically induced interseismic uplift; postglacial rebound- Satake et al. 2003;Leonard et al. 2004) and in the relationship between 1700 coseismic slip and intersesimic velocity using an idealized rheology with poorly constrained viscosity values.
The WUS deformation model does not incorporate interseismic vertical strain rates derived from combined levelling profiles and tidal records (e.g. Mitchell et al. 1994;Burgette et al. 2009). These data may be able to refine the location of the locked-transition zone boundary on the Cascadia subduction fault in future models as vertical land-level changes occur above the locked portion of the fault.

Mantle viscosity around the Cascadia megatrust
The contribution of post-1700 relaxation to the crustal velocity field is such as to produce strong westward velocity at the coast early in the cycle (with corresponding E-W extension) and moderate eastward velocity in the middle to late part of the cycle (with corresponding E-W contraction) ( fig. 19 of Pollitz et al. 2008). The present contribution is such as to produce E-W contraction south of 47 • N and SW-NE contraction north of 47 • N (Fig. 4). The WUS deformation model, which uses the laterally homogeneous structure ( fig. 7 of Pollitz et al. 2008, mantle viscosity of 10 19 Pa s) and source parameters of Table 2, however, underpredicts the amount of eastward motion at the coast between 44 • and 47 • N by ∼3-4 mm yr −1 (Fig. 19) and consequently the amount of E-W contraction at the coast (Fig. 20). The post-1700 relaxation is sensitive to the mantle viscosity, and an increase in viscosity produces larger coastal contraction by virtue of using an infinite series of relaxation from past (assumed periodic) 1700-type events, as prescribed by the first term of eq. (1). We tested the effect of increasing the mantle viscosity for the time-dependent response functions of faults #36-41 off Oregon and southern Washington, using the original viscoelastic structure for all other time-dependent sources, that is, those of table 3 of Pollitz et al. (2008). Modest changes in mantle viscosity yield almost identical slip rates on Cascadia megathrust segments but change the amount of predicted coastal convergence. A factor of three increase in mantle viscosity suffices to replicate observed E-W motions (Fig. 19) and consequently the observed coastal contraction (Fig. 20). This indicates that the effective mantle viscosity, at least near the Oregon and southern Washington sections of the Cascadia megathrust, is larger than elsewhere in the western US, likely because of the presence of the downgoing slab.

Western Washington
Reverse faulting within the Olympic Peninsula and eastern Puget Lowland zones explains observed contraction within these regions. The observed contraction, after correction for the effects of the viscoelastic cycle on the megathrust, is about 8 mm yr −1 N67 • E within the Olympic Peninsula and 3 mm yr −1 N30 • E around the eastern Puget Lowland (Fig. 7). The residual contraction is negligible (Fig. 7). Within the Olympic Mountains this is accomplished with reverse faulting on a fictitious structure localized in the western Olympic Peninsula. It is a proxy for a likely continuous crustal flow field associated with the accretionary wedge (Pazzaglia & Brandon 2001). The horizontal shortening of 10 mm yr −1 across fault #86 (representing the Olympic peninsula) needed to match the observed shortening considerably exceeds the ∼4 mm yr −1 shortening inferred by Pazzaglia & Brandon (2001) based on stream terrace analysis; it correspondingly implies that about 50 per cent of the observed instantaneous shortening across the Olympic peninsula is due to long-term shortening of the Olympic peninsula, and the remaining 50 per cent due to the effect of the slab. Our result is larger than the ∼20 per cent attributable to long-term shortening stated by  fig. 10a). Thus the amount of long-term shortening inferred by Pazzaglia & Brandon (2001) is about mid-way between the shortening used in the model of McCaffrey et al. (2007) and the present study. This highlights a clear difference between the block-model and viscoelastic-cycle approaches, which we attribute to the time dependence of interseismic velocity resulting from viscoelastic cycles on major faults, including the Cascadia megathrust. The disparity between the two approaches could be reduced if the velocity field were corrected for time-dependent effects, which is feasible in the block modelling approach (Puskas & Smith 2009), or if the western Olympic Mountains structure were simply included in the catalogue of faults used in the block model. In the realm of the present approach, the long-term shortening rate in the Olympic peninsula trades off with local viscosity structure, and it is conceivable that a lower long-term shortening rate is admissible provided that local mantle viscosity is higher, as was noted for southern Washington and Oregon in Section 7.4.
In the eastern Puget Lowland the observed contraction is accommodated with a net 5 mm yr −1 slip on the 45 • -dipping Tacoma/Seattle faults, amounting to 3 mm yr −1 horizontal slip rate. The inferred 4 mm yr −1 on the Seattle fault is considerably larger than the ∼1 mm yr −1 Holocene slip rate estimated on the Seattle fault zone (Johnson et al. 1999). This may reflect reverse slip on additional faults over a broader region that is being projected onto the Seattle fault in our model, for example, additional faults around Puget Sound (Johnson et al. 1996 or the Olympia fault . The ∼3 mm yr −1 geodetic shortening rate is consistent with the inference of ∼36 m horizontal shortening of the region during the Holocene, based on palaeoseismic investigations of 11 faults in the Eastern Puget Sound (Sherrod et al. 2008b).
About 3 mm yr −1 left-lateral slip is inferred on the Hood Canal fault (Table 2). This is consistent with left-lateral slip on the nearby Saddle Mountain East fault (Wilson et al. 1979) and Canyon River fault (Walsh & Logan 2007). This zone is interpreted to be a shear zone caught between the Olympic Mountains and the Figure 19. Left-hand panel: observed velocity field and (middle and right) residual velocity fields using a mantle viscosity of 10 19 Pa s and 3 × 10 19 Pa s for the post-1700 response of megathrust segments #36-41. The residual velocity field with mantle viscosity 10 19 Pa s exhibits ∼3-4 mm yr −1 unpredicted eastward motion at coastal sites between 44 • and 47 • N, leading to an underprediction of coastal contractile strain rate (Fig. 20). This is remedied by increasing the mantle viscosity to 3 × 10 19 Pa s. northward-translating East Puget Lowland Witter et al. 2008). The 3 mm yr −1 inferred for the Hood Canal fault may be distributed on several left-lateral faults on the eastern border of the Olympic Mountains.
Slip rates of 2 mm yr −1 given in Table 2 for the both Devils Mountain fault and southern Whidbey fault are much smaller than the values that the inversion would formally yield for these faults (12.2 ± 1 and 6.9 ± 1 mm yr −1 , respectively. The formally inverted values may result from neglect of systematic signals in this area, for example, several slowly slipping other minor faults  or an inaccurate viscoelastic model. We deemed these values unrealistically large given the palaeoseismic history of these Figure 20. Magnitude of east-west contractile strain using the observed velocity field (Fig. 19) and as predicted by models derived using the original viscoelastic structure ( fig. 7 of Pollitz et al. 2008, mantle viscosity of 10 19 Pa s) and a revised viscoelastic structure (mantle viscosity of 3 × 10 19 Pa s for the post-1700 response of megathrust segments #36-41; mantle viscosity of 10 19 Pa s for all other time-dependent sources listed in table 3 of Pollitz et al. 2008). The contractile strain rate is evaluated at longitude 123.6 • W using the method in section 5.3 of Pollitz et al. (2008). faults, which suggest on the order of 1 mm yr −1 slip rates (Johnson et al. 1996Sherrod et al. 2008a) and repeated the inversion assigning the given values. The left-lateral transpression on the Devils Mountain fault (Table 2) may be distributed on several faults within the eastern Strait of Juan de Fuca (Johnson et al. , 2004. Similarly, reverse slip on the southern Whidbey fault may be distributed on several subparallel strands (Johnson et al. 1996;Sherrod et al. 2008a).

C O N C L U S I O N S
Using a new compilation of regional GPS data spanning the western US, we construct a viscoelastic-cycle model of interseismic crustal deformation. This model (named the WUS deformation model) depends on the viscoelastic stratification of the lithosphere and upper mantle, fault slip rates, and slip of selected past large earthquakes. It is driven by time-independent relaxation (i.e. deformation rates averaged over all seismic cycles) for numerous minor faults and time-dependent (i.e. viscoelastic) relaxation for major faults, and it uses numerous onland faults and includes the offshore GDZ and the Juan de Fuca and Explorer plates.
The detailed data set and interpretation method reveal new details of selected regions. Models require significant amounts of strike slip continuing from the San Andreas fault system northward of the MTJ to major Humboldt Bay faults (Russ and Little Salmon faults), and part of the GDZ convergence is accommodated by dip-slip on the Mad River fault and likely other nearby thrust structures. Mantle relaxation from the 1700 Cascadia earthquake and similar-sized prior events shapes the velocity fields near the coasts of northern California, Oregon, Washington, and Vancouver Island. For an assumed viscoelastic structure, the interseismic velocity field yields a rough estimate of the amount of coseismic slip in the 1700 earthquake, increasing to 20 m off northern Washington. SW-NE contraction in western Washington demands additional sources of faulting in the upper plate. Based on the amount of contraction remaining after correcting for the viscoelastic cycle on the megathrust, two main sources of such faulting are identified as the Tacoma/Seattle fault zones and an inferred zone of steady slip within the Olympic peninsula. These and other minor sources of faulting in western Washington accommodate the large SW-NE stresses imparted by the descending slab off the Olympic peninsula and by the northwardtranslating Oregon Coast block.

A C K N O W L E D G M E N T S
This study benefited from the high-quality geodetic data from the PBO provided by UNAVCO. The paper was improved by the critical comments of two anonymous reviewers and the Associate Editor. We thank Jim Savage and Wayne Thatcher for their comments on a preliminary draft.   Fig. B1 illustrates the relative velocities among several plates at a latitude of about 41 • N. It highlights the fact that the eastern GDZ, particularly the southeast GDZ, moves several mm yr −1 southeastward relative to the southwestern GDZ because of combined north-south contraction and east-west extension within the plate as deduced by mapping of magnetic lineations (Wilson 1989;Gulick et al. 2001;Chayton et al. 2004) and seismic focal mechanisms (Fig. 8). We represent the internal GDZ kinematics with a distribution of left-lateral strike-slip faulting in the internal deformation zone, amounting to 10 mm yr −1 north-south shortening and 10 mm yr −1 of east-west extension across the zone. Hence the relative motion between the southeast GDZ and northern California is up to 10 mm yr −1 larger than would be expected for the relative motion between the southwestern GDZ and northern California. The azimuth of relative motion on the megathrust, however, varies little with latitude in this kinematic model. This is a consequence of variable motion of the northern California coast (Wilson 1989). The northward component of this motion diminishes with increasing latitude because of rightlateral slip along faults north of the MTJ. The northward component of motion along the northern California coast and the adjacent GDZ covary such as to produce an azimuth of relative motion ∼68 • for slippage along megathrust segments #30, 30a, 31, 31a, 32 and 33. Inland northern California in Fig. B1 refers to the area east of the SAF-Maacama-Bartlett Springs faults and has similar motion to the OC block, though the latter may differ slightly because strike slip from the Eastern California Shear Zone leads into southern Oregon and may terminate at the megathrust .

A P P E N D I X B : V E L O C I T Y R E L AT I O N S H I P S A RO U N D M E N D O C I N O T R I P L E J U N C T I O N
Half-spreading rates on the Gorda Ridge are 27.5 mm yr −1 in the north and 14 mm yr −1 in the south (Wilson 1989). Based on the spreading rate of the southern Gorda Ridge and the 17 mm yr −1 east-west extension across the GDZ, the slip rate on the Mendocino transform is prescribed as 28 mm yr −1 west of 126 • W (where the deformation zone intersects the transform fault) and 45 mm yr −1 east of 126 • W. Figure B1. Velocity relationships at latitude ∼41 • N among several plates bordering the Mendocino triple junction-North America and Pacific plates, Gorda deformation zone, and northern California points (represented with variable northward velocity to account for relative convergence between the Sierra Nevada block and the Oregon Coast block to its north. Numerals refer to fault numbers indexed in Table 2 of main text. Southeast GDZ velocities are uncertain but have greater southward and eastward velocity components relative to the southwest GDZ because of combined north-south shortening and east-west extension of the GDZ (Wilson 1989).