TY - JOUR
AU - Evans,, E.L.
AB - Summary GPS velocity fields in the Western US have been interpreted with various physical models of the lithosphere-asthenosphere system: (1) time-independent block models; (2) time-dependent viscoelastic-cycle models, where deformation is driven by viscoelastic relaxation of the lower crust and upper mantle from past faulting events; (3) viscoelastic block models, a time-dependent variation of the block model. All three models are generally driven by a combination of loading on locked faults and (aseismic) fault creep. Here we construct viscoelastic block models and viscoelastic-cycle models for the Western US, focusing on the Pacific Northwest and the earthquake cycle on the Cascadia megathrust. In the viscoelastic block model, the western US is divided into blocks selected from an initial set of 137 microplates using the method of Total Variation Regularization, allowing potential trade-offs between faulting and megathrust coupling to be determined algorithmically from GPS observations. Fault geometry, slip rate, and locking rates (i.e. the locking fraction times the long term slip rate) are estimated simultaneously within the TVR block model. For a range of mantle asthenosphere viscosity (4.4 × 1018 to 3.6 × 1020 Pa s) we find that fault locking on the megathrust is concentrated in the uppermost 20 km in depth, and a locking rate contour line of 30 mm yr−1 extends deepest beneath the Olympic Peninsula, characteristics similar to previous time-independent block model results. These results are corroborated by viscoelastic-cycle modelling. The average locking rate required to fit the GPS velocity field depends on mantle viscosity, being higher the lower the viscosity. Moreover, for viscosity ≲ 1020 Pa s, the amount of inferred locking is higher than that obtained using a time-independent block model. This suggests that time-dependent models for a range of admissible viscosity structures could refine our knowledge of the locking distribution and its epistemic uncertainty. Space geodetic surveys, Transient deformation, Dynamics of lithosphere and mantle 1 INTRODUCTION The Cascadia megathrust of the Pacific Northwest (Fig. 1) is recognized to be seismogenic (e.g. Adams 1990; Atwater & Hemphill-Haley 1997) and accumulating strain (Savage et al. 1981, 1991). The amount and distribution of fault locking on the megathrust is critical for understanding earthquake hazards (Petersen et al. 2014), the potential for damaging tsunamis, and the nature of the transition to aseismic slip or episodic tremor and slip at greater depth (Wech et al. 2009). The amount of Cascadia megathrust fault locking has been quantified using block models (e.g. McCaffrey et al. 2007, 2013; Schmalzle et al. 2014; Evans et al. 2015) and viscoelastic-cycle models which account for earthquake-cycle effects (e.g. Pollitz et al. 2010). The former models interseismic deformation as time-independent, while the latter models interseismic deformation as time-dependent because it is driven by viscoelastic relaxation of the ductile lower crust and upper mantle. Although the forward geophysical models used in the two approaches are very different, the two approaches are equivalent for self-consistent fault geometries, the block model being an endmember case of the viscoelastic-cycle model in the limit of small ‘Savage parameter’ (Hetland & Hager 2005) (ratio of earthquake recurrence interval to asthenosphere material relaxation time). This has been demonstrated by Savage & Prescott (1978) for the case of a very long strike-slip fault in an elastic layer over a viscoelastic substrate and Johnson & Fukada (2010) for more general fault networks that comprise a block model; Johnson & Fukada (2010) refer to their resulting time-dependent formulation as a viscoelastic block model. Viscoelastic block models are a natural extension of the Savage & Prescott (1978) viscoelastic-cycle model for time-dependent deformation where crustal blocks are enforced to behave on average rigidly over long time periods. They have been implemented in northern California (Johnson & Fukada 2010) and southern California (Chuang & Johnson 2011; Hearn et al. 2013) to rationalize GPS crustal velocity fields in the presence of likely time-dependent deformation arising primarily from viscoelastic mantle asthenosphere relaxation. They have thus far not been applied to crustal deformation fields observed in the Pacific Northwest despite ample evidence that mantle relaxation is an important part of the earthquake cycle in megathrust settings (Wang et al. 2012). In order to put the Pacific Northwest deformation sources into a wider context, we use the Western US GPS velocity field to infer the plate-scale sources of deformation. We employ methods that are closely related to both the viscoelastic cycle models of Pollitz et al. (2008, 2010) and the block models of Evans et al. (2015), but differ in several respects. First, we show how the (time-dependent) viscoelastic block model emerges from the Pollitz et al. (2010) viscoelastic-cycle model, from which follows the equivalence between the time-independent block model and viscoelastic-cycle model in the limit of small Savage parameter. We then use the observed GPS velocity field to explore two classes of time-dependent crustal deformation models in the Pacific Northwest. The first is the viscoelastic block model, which requires the decomposition of the crust into distinct blocks. The second is the more general viscoelastic-cycle model, which does not enforce long-term rigidity of any specific crustal domains and hence does not require the definition of blocks. In order to realize earthquake-cycle effects in the viscoelastic block model, time-dependent earthquake-cycle deformation above the average interseismic rates are implemented as a correction to the observed velocity field; the corrected velocity field is then interpreted with a conventional block model following the methods of Evans et al. (2015) which are recapitulated below. In the viscoelastic cycle model, earthquake-cycle effects are realized directly through evaluation of viscoelastic relaxation from an idealized earthquake history of the model faults. Using the Western US models, we shall focus on resulting estimates of fault locking on the Cascadia megathrust. It is important to note that the interseismic velocity field governed by the viscoelastic models is time-dependent, but with a snapshot of this velocity field we may infer properties of the megathrust that are assumed time-independent, that is, the degree of coupling. For simplicity we shall limit consideration of time-dependent earthquake-cycle effects to those contributed by the multi-century period Cascadia earthquake cycle. A set of models with varying mantle asthenosphere viscosity shall highlight the difference between time-dependent and time-independent models, as well as illuminate the epistemic uncertainty in megathrust locking. 2 VISCOELASTIC BLOCK MODEL 2.1 Model prescription In this section, we build upon basic theory of viscoelastic-cycle deformation presented in earlier studies to derive a formalism for the viscoelastic block model. Following Pollitz et al. (2008), crustal deformation is envisaged to be the result of post-earthquake relaxation of a ductile lower crust and mantle underlying an elastic upper crust, with further allowance for fault creep. This general prescription for interseismic velocity is here shown to reduce to a sum of block motions, ‘backslip’-driven motions, and a correction term that accounts for variations in interseismic velocity over a seismic cycle. Let Γn define the nth (discrete) fault surface, assumed to be within an elastic volume (i.e. the upper crust). Slip events on Γn are assumed to be identical and periodic with recurrence interval Tn. Denoting the time of the last event with tn, this implies repeating identical slip events at an infinite succession of past times tn − jTn, j ≥ 0. The fault geometry and slip of these events are represented through the moment-release density m(r΄) at points r΄ ∈ Γn. This is complemented by steady creep on fault surfaces Γcr represented through the moment-release rate density |$\mathbf { \dot{m}}^{\rm (cr)} (\mathbf { r^{\prime }})$| at points |$\mathbf { r^{\prime }}\in \Gamma _{\text{cr}}$|⁠, also assumed to be within the elastic upper crust. Then interseismic velocity at a point r and time t (0 < t − tn < Tn) is \begin{eqnarray} \mathbf {\rm v}^{\rm inst}(\mathbf {r}, t ) &=& \mathbf {\rm v}^{\rm cycle}(\mathbf {r}, t ) \, + \, \mathbf {\rm v}^{\rm creep}(\mathbf {r}) \nonumber\\ &=& \sum _n \int _{\Gamma _n} \text{d}^2 \mathbf { r^{\prime }}\, \mathbf {m}(\mathbf { r^{\prime }}) : \sum _{j\ge 0} \mathbf {\dot{G}}^{(d)} (\mathbf {r},\mathbf { r^{\prime }},t - t_n + j T_n) \nonumber \\ &&+\, \int _{\Gamma _{\text{cr}}} d^2 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}^{\rm (cr)}(\mathbf { r^{\prime }}) : \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},\infty ). \end{eqnarray} (1) Here the colon denotes tensor double dot product, and G(d)(r, r΄, t) is the displacement response of the viscoelastic system at point r and time t to a unit dislocation source applied at point r΄ and time 0. In terms of tensor components, with xp denoting the pth receiver Cartesian coordinate and |$x^{\prime }_q$| the qth source Cartesian coordinate, |$G_{kpq}^{(d)}$| is meant to represent the xk −component displacement resulting from a unit moment tensor component mpq, and it is related to the Green's function G commonly defined in seismic source theory via |$G_{kpq}^{(d)} = \partial _{x^{\prime }_q} G_{kp} \, \hat{x}_k$|⁠, where Gkp is the xk −component displacement response to a unit point force in the xp direction (e.g. eq. 3.20 of Aki & Richards 1980). Then |$\mathbf {m} : \mathbf {\dot{G}}^{(d)}$| denotes contraction over p and q to yield the xk −component of velocity. The first term of eq. (1) represents viscoelastic relaxation from infinite sequences of past earthquakes in idealized seismic cycles, as described above. The second term represents steady creep. The Green's function in the limit of complete relaxation is appropriate for this term as it represents a limit of combined static displacement and relaxation of a moment release m(cr) per time interval T as T → 0 and |$\mathbf {m}^{\rm (cr)} / T \rightarrow \mathbf { \dot{m}}^{\rm (cr)}(\mathbf { r^{\prime }})$|⁠, that is, \begin{eqnarray} \mathbf {\rm v}^{\rm creep}(\mathbf {r}, t ) &=& \lim _{T\rightarrow 0} \int _{\Gamma _{\text{cr}}} \text{d}^2 \mathbf { r^{\prime }}\, \frac{\mathbf {m}^{\rm (cr)}(\mathbf { r^{\prime }})}{T} :\nonumber\\ &&\quad \left[ \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},0^+) \, + \, \sum _{j\ge 0} \mathbf {\dot{G}}^{(d)} (\mathbf {r},\mathbf { r^{\prime }}, j T) \times T \right] \end{eqnarray} (2) which reduces to the second term of eq. (1). The intersesimic time average of the first term in eq. (1) is \begin{eqnarray} \mathbf {\bar{\rm v}}^{\rm cycle}( \mathbf {r}) &=& \sum _n \int _{\Gamma _n} \text{d}^2 \mathbf { r^{\prime }}\, \mathbf {m}(\mathbf { r^{\prime }}) :\! \left[ \sum _{j\ge 0} \frac{1}{T_n} \int _{0^{+}}^{T_n} \, \mathbf {\dot{G}}^{(d)} (\mathbf {r},\mathbf { r^{\prime }},t {+} j T_n) \text{d}t \!\right] \, \nonumber \\ &=& \sum _n \int _{\Gamma _n} \text{d}^2 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}(\mathbf { r^{\prime }}) : \left[ \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},\infty ) - \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},0^+) \right] \nonumber\\ \end{eqnarray} (3) where the average moment release rate is \begin{equation} \mathbf { \dot{m}}(\mathbf { r^{\prime }}) = \frac{ \mathbf {m}(\mathbf { r^{\prime }})}{T_n} \, \, \, ( \mathbf { r^{\prime }}\in \Gamma _n) \end{equation} (4)|$\mathbf {\bar{\rm v}}^{\rm cycle}( \mathbf {r})$| in eq. (3) is the ‘expected interseismic velocity’ described by Pollitz (2003). Each term represents the contribution of the nth fault to the interseismic deformation field, which may be interpreted as an average over the viscoelastic cycle of that fault or the time-dependent contribution to vcycle in the limit of infinite viscosity of the relaxing components of the earth model. In the latter case, let η represent the viscosities of these relaxing components and add this to the argument of the Green's function where necessary. Then |$\mathbf {\bar{\rm v}}^{\rm cycle}( \mathbf {r})$| can be written in the alternative form \begin{equation} \mathbf {\bar{\rm v}}^{\rm cycle}( \mathbf {r}) = \sum _n \int _{\Gamma _n} \text{d}^2 \mathbf { r^{\prime }}\, \mathbf {m}(\mathbf { r^{\prime }}) : \left[ \lim _{\mathbf {\eta } \rightarrow \infty } \sum _{j\ge 0} \mathbf {\dot{G}}^{(d)} (\mathbf {r},\mathbf { r^{\prime }}, j T_n ; {\mathbf {\eta }} ) \right] \end{equation} (5) where limη → ∞ denotes the limit as the viscosities tend to arbitrarily large values. Substituting eqs (3) and (5) into eq. (1) yields \begin{eqnarray} \mathbf {\rm v}^{\rm inst}(\mathbf {r}, t) &=& \delta \mathbf {\rm v}^{\rm cycle}(\mathbf {r}, t) + \sum _n \int _{\Gamma _n} \text{d}^2 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}(\mathbf { r^{\prime }}) :\nonumber\\ && \left[ \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},\infty ) - \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},0^+) \right] \nonumber \\ &&\quad+\, \int _{\Gamma _{\text{cr}}} \text{d}^2 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}^{\rm (cr)}(\mathbf { r^{\prime }}) : \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},\infty ) \end{eqnarray} (6) where \begin{eqnarray} \delta \mathbf {\rm v}^{\rm cycle}(\mathbf {r}, t ) &=& \sum _n \int _{\Gamma _n} \text{d}^2 \mathbf { r^{\prime }}\, \mathbf {m}(\mathbf { r^{\prime }}) : \left[ \sum _{j\ge 0} \mathbf {\dot{G}}^{(d)} (\mathbf {r},\mathbf { r^{\prime }},t - t_n + j T_n ; {\mathbf {\eta }} )\right. \nonumber\\ &&\left.- \, \lim _{\mathbf {\eta } \rightarrow \infty } \sum _{j\ge 0} \mathbf {\dot{G}}^{(d)} (\mathbf {r},\mathbf { r^{\prime }}, j T_n ; {\mathbf {\eta }} ) \right] \end{eqnarray} (7) If the fault network consisting of the combined sets {Γn} and {Γcr} bound crustal blocks that behave rigidly when averaged over long time intervals, then the second and third terms of eq. (6) are readily interpreted in terms of a combination of rigid block motions and backslip along the bounding faults. Let the volume occupied by block l be denoted by Vl. Then for |$\hat{\mathbf {r}}\in V_l$| and a ‘self-consistent’ prescription of long-term slip rates on the faults, \begin{eqnarray} \mathbf {\rm v}^{\rm block}(\mathbf {r}) &=& \sum _n \int _{\Gamma _n} \text{d}^2 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}(\mathbf { r^{\prime }}) : \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},\infty ) \nonumber\\ &&\quad+ \, \int _{\Gamma _{\text{cr}}} \text{d}^2 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}^{\rm (cr)}(\mathbf { r^{\prime }}) : \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},\infty ) = {\mathbf {\omega }}_l\mathbf {\times } \hat{\mathbf {r}} \end{eqnarray} (8) where ωl is the Euler vector of the lth block, and \begin{equation} \mathbf {\rm v}^{\rm backslip}(\mathbf {r}) = \sum _n \int _{\Gamma _n} d^2 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}(\mathbf { r^{\prime }}) : \left[ - \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},0^+) \right] \end{equation} (9) represents the deformation from backslip on the fault system. The simplest interpretation of eq. (8) is that it represents the motions of a set of thin elastic plates overlying a fluid, behaving like blocks of wood floating on water, provided that the prescribed relative motions between each plate are consistent with their long-term rigid rotation. Note that the collection of creeping-fault surfaces {Γcr} are envisaged to overlap partially or completely with the periodically-rupturing fault surfaces {Γn} ; where they overlap, the moment rates |$\mathbf { \dot{m}}(\mathbf { r^{\prime }})$| and |$\mathbf { \dot{m}}^{\rm (cr)}(\mathbf { r^{\prime }})$| must sum to a net moment rate consistent with the long-term slip rate on the fault. The backslip given by eq. (9) would then include only those fault surfaces that release strain episodically, that is, |$\mathbf { \dot{m}}(\mathbf { r^{\prime }})$| in that equation excludes any slip rate due to steady creep. From eqs (6) to (9), the instantaneous interseismic velocity is \begin{equation} \mathbf {\rm v}^{\rm inst}(\mathbf {r}, t ) = \delta \mathbf {\rm v}^{\rm cycle}(\mathbf {r}, t) \, + \, \mathbf {\rm v}^{\rm block}(\mathbf {r}) \, + \, \mathbf {\rm v}^{\rm backslip}(\mathbf {r}). \end{equation} (10) This is equivalent to the viscoelastic block model of Johnson & Fukada (2010), and the implied stressing rate is consistent with the related formulation of Devries & Meade (2016). Eqs (7) and (10) confirm that block model formulations that involve only the vblock(r) and vbackslip(r) terms are equivalent to a physically-based model driven by viscoelastic earthquake cycles in the limit of small Savage parameter, that is, very short earthquake recurrence time(s) and/or large viscosity of the relaxing components, as noted in previous studies (McClusky et al. 2001; Pollitz 2003; Meade & Hager 2005). The δvcycle term represents time-dependent viscoelastic-cycle effects. Such effects are already part of the interseismic motions, and this is embodied in the prescription given by eq. (7) in which δvcycle is that part of the viscoelastic relaxation above that given by the high-viscosity limit. Hearn et al. (2013) refer to δvcycle as the ‘ghost transient’, that part of the relaxation over and above that which is implicitly part of the block-motion solution. 2.2 Implementation The calculation of δvcycle in the form of eq. (7) is done using the viscoelastic-gravitational normal mode approach (Pollitz 1997). Note that the second term of eq. (7) involves the evaluation of an infinite sum of velocity viscoelastic-response functions |$\mathbf {\dot{G}}^{(d)}$| evaluated in the high-viscosity limit. This is readily accomplished using the normal mode approach, as the infinite summation and high-viscosity limit may be evaluated analytically for each contributing mode. In practice, δvcycle can be implemented on a fault-by-fault basis when knowledge of the earthquake history is sufficient for this purpose. Following Hearn et al. (2013), we use the δvcycle term—evaluated for selected faults —as a correction to the observed interseismic velocity field vobs(r) and then interpret the remainder in terms of a block model. Faults may include an unknown amount of fractional locking f which varies from 0 for creeping to 1 for completely locked. If |$\mathbf { \dot{m}}^{\rm total}(\mathbf { r^{\prime }})$| is the total moment density rate, then \begin{eqnarray} \mathbf { \dot{m}}(\mathbf { r^{\prime }}) &=& f (\mathbf { r^{\prime }}) \times \mathbf { \dot{m}}^{\rm total}(\mathbf { r^{\prime }}) \nonumber \\ \mathbf { \dot{m}}^{\rm (cr)}(\mathbf { r^{\prime }}) &=& ( 1 - f (\mathbf { r^{\prime }}) ) \times \mathbf { \dot{m}}^{\rm total}(\mathbf { r^{\prime }}) . \end{eqnarray} (11) Starting with trial values of the collection of Euler vectors {ωl}, we first prescribe the long-term slip on a fault bounding blocks l and k. The long-term vector velocity of the block l side relative to the block k side at point r on the fault surface Γlk between these two blocks is \begin{equation} \mathbf {s}_{lk} = ( {\mathbf {\omega }}_l- {\mathbf {\omega }}_k) \times \mathbf {r}. \end{equation} (12) Assume that the hanging wall of the fault is on block l, and let |$\hat{n}$|⁠, |$\hat{d}$|⁠, and |$\hat{t}$| denote the normal, down-dip, and along-strike vectors of the footwall block. Then the moment tensor density associated with seismic slip for r΄ ∈ Γlk is (eq. 3.21 of Aki & Richards 1980) \begin{eqnarray} \mathbf { \dot{m}}(\mathbf { r^{\prime }}) &=& \lbrace ( \mathbf {s}_{lk} \cdot \hat{n}) [(\lambda + 2 \mu) \hat{n} \hat{n} + \lambda \hat{d} \hat{d} + \lambda \hat{t} \hat{t}] \nonumber\\ &&\quad+\, ( \mathbf {s}_{lk} \cdot \hat{d}) \mu [\hat{n} \hat{d} + \hat{d} \hat{n}] + ( \mathbf {s}_{lk} \cdot \hat{t}) \mu [\hat{n} \hat{t} + \hat{t} \hat{n}] \rbrace \nonumber \\ && \times ( 1 - f_{lk} ) \end{eqnarray} (13) where λ and μ are Láme parameters and flk is the fraction of slip being accommodated seismically. The |$( \mathbf {s}_{lk} \cdot \hat{d})$| and |$( \mathbf {s}_{lk} \cdot \hat{t})$| terms of eq. (13) represent slip rate in the down-dip and along-strike directions, respectively, while the |$( \mathbf {s}_{lk} \cdot \hat{n} )$| term represents tensile opening (or contraction) across the fault. Tensile opening is problematic to deal with in block models because long-term motion of this type is unsustainable for a non-vertical fault. Since tensile opening does not typically occur in tectonic deformation, we choose to set the |$( \mathbf {s}_{lk} \cdot \hat{n} )$| term to zero for the Cascadia megathrust. This practice is identical to that of other authors (e.g. Savage 1983; McCaffrey et al. 2007). 3 VISCOELASTIC CYCLE MODEL The foregoing section shows how different forms of the block model (time-dependent with δvcycle(r, t) ≠ 0 or time-independent with δvcycle(r, t) = 0) result from the viscoelastic cycle model provided that the distribution of long-term slip rates permits eq. (8) to represent long-term motions as rigid block rotations. In the absence of such a self-consistent fault network, we may implement the viscoelastic cycle model without a strict requirement of long-term rigidity within any defined crustal volume. We retain here four of the five contributions originally presented by Pollitz et al. (2008): \begin{eqnarray} \mathbf {\rm v}^{\rm inst}(\mathbf {r}, t ) &=& \sum _n \int _{\Gamma _n} \text{d}^2 \mathbf { r^{\prime }}\mathbf {m}(\mathbf { r^{\prime }}) : \left[ \sum _{j\ge 0} \mathbf {\dot{G}}^{(d)} (\mathbf {r},\mathbf { r^{\prime }},t - t_n + j T_n) \right] \nonumber \\ &&+\, \sum _m \int _{\Gamma _m} \!\text{d}^2 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}^{\rm }(\mathbf { r^{\prime }})\! :\! \left[ \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},\infty ) - \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},0^+) \right] \nonumber \\ &&+\, \int _{\Gamma _{\text{cr}}} \text{d}^2 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}^{\rm (cr)}(\mathbf { r^{\prime }}) : \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},\infty ) \nonumber \\ &&+\, \int _{V} \text{d}^3 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}^{( \delta \mu )}(\mathbf { r^{\prime }}) : \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},\infty ) \end{eqnarray} (14) Here V refers to the volume of the lithosphere, which is assumed to be populated with discrete fault surfaces. The first and third terms are those of eq. (1) and represent viscoelastic relaxation from all known/estimated past major regional earthquakes on a set of faults Γn and secular deformation arising from steady creep on a set of faults Γcr, respectively. The second term represents the interseismic average of viscoelastic relaxation from a set of faults Γm resulting from an average moment release rate |$\mathbf { \dot{m}}$|⁠; this prescription is suitable for those faults where the earthquake history is unknown. The fourth term represents the effects of lateral heterogeneity in shear modulus δμ(r΄) and (implicitly) bulk modulus δκ(r΄) at points r΄ ∈ V; |$\mathbf { \dot{m}}^{( \delta \mu )}(\mathbf { r^{\prime }})$| depends on δμ(r΄) and δκ(r΄) and the observed strain rate through eq. (8) of Pollitz et al. (2008). That is, \begin{equation} \mathbf { \dot{m}}^{( \delta \mu )}(\mathbf { r^{\prime }}) = - 2 \, \delta \mu (\mathbf { r^{\prime }}) \, \mathbf {\dot{\sf {D}}} (\mathbf { r^{\prime }}) - \delta \kappa ( \mathbf { r^{\prime }}) \, \nabla \cdot \mathbf {\rm v}^{\rm inst}(\mathbf { r^{\prime }}) \, \mathbf {I} \end{equation} (15) where |$\mathbf {\dot{\sf {D}}} (\mathbf { r^{\prime }})$| is the deviatoric strain rate tensor, ∇ · vinst(r΄) is the dilatational strain rate, and I is the identity tensor. If the strain rate field is assumed constant as a function of depth, then |$\mathbf {\dot{\sf {D}}} (\mathbf { r^{\prime }})$| and ∇ · vinst(r΄) are completely determined by the observed strain rate field at earth's surface. Additional assumptions of depth-independent δμ(r΄) and a scaling relation (⁠|$\delta \kappa = \frac{5}{3} \delta \mu$|⁠) permits this term to be represented in terms of laterally variable δμ averaged over the elastic thickness (Pollitz et al. 2008). 4 DATA SET We use observations of the crustal velocity field compiled from numerous GPS data sets. These data span collectively the entire ∼2500 km long and ∼1500 km wide plate boundary zone of the western US (Fig. 1). As summarized by Evans et al. (2015), GPS velocity fields are contributed by McClusky et al. (2001), Shen et al. (2003), Hammond & Thatcher (2005), Williams et al. (2006), McCaffrey et al. (2007) and Loveless & Meade (2011). These are appended with GPS velocities from the Plate Boundary Observatory (PBO) (http://pboweb.unavco.org) and USGS campaign measurements. The USGS campaign measurements are described in numerous prior publications (Thatcher et al. 1999; Savage et al. 1998, 1999a,b, 2001a,b; 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 2011. The McCaffrey et al. (2007) data set is a compilation of various continuous and campaign GPS measurements, including USGS campaign measurements conducted up to 2006; the measurements in this compilation span the time interval 1991–2006. The PBO measurements were initiated in 2004 and represent up to 7 yr of continuous measurements. After removal of reference stations and outliers, we employ altogether 1896 velocity vectors. In general, all contributing data sets have been processed in slightly different realizations of fixed North America. Velocity fields are referred to a common reference frame by aligning the velocity fields at collocated sites using a six-parameter transformation (Loveless & Meade 2011). The resulting velocity field relative to a fixed North America is shown in Fig. 2, which is subsampled from the full 1896 velocity vectors shown in Supporting Information Fig. S1. Figure 1. Open in new tabDownload slide Tectonic provinces, M > 3 earthquakes (black circles), major plate boundary faults and deformation zones (highlighted in white) in the Pacific Northwest and wider western United States. Deformation zones include CNSB (Central Nevada seismic belt), ECSZ (Eastern California Shear Zone), WLSB (Walker Lane seismic belt) and ISB (Intermountain seismic belt). Other abbreviations are: GH (Grays Harbor), MTJ (Mendocino triple junction), PL (Puget Lowland), YSP (Yellowstone Plateau) and ESRP (eastern Snake River Plain). Modified from Puskas & Smith (2009). Figure 1. Open in new tabDownload slide Tectonic provinces, M > 3 earthquakes (black circles), major plate boundary faults and deformation zones (highlighted in white) in the Pacific Northwest and wider western United States. Deformation zones include CNSB (Central Nevada seismic belt), ECSZ (Eastern California Shear Zone), WLSB (Walker Lane seismic belt) and ISB (Intermountain seismic belt). Other abbreviations are: GH (Grays Harbor), MTJ (Mendocino triple junction), PL (Puget Lowland), YSP (Yellowstone Plateau) and ESRP (eastern Snake River Plain). Modified from Puskas & Smith (2009). Figure 2. Open in new tabDownload slide GPS data used in this study subsampled from the original 1896 velocity vectors (Supporting Information Fig. S1) for visual clarity. Superimposed are the Juan de Fuca plate boundaries and slab depth contours of McCrory et al. (2006) in intervals of 10 km. Figure 2. Open in new tabDownload slide GPS data used in this study subsampled from the original 1896 velocity vectors (Supporting Information Fig. S1) for visual clarity. Superimposed are the Juan de Fuca plate boundaries and slab depth contours of McCrory et al. (2006) in intervals of 10 km. 5 APPLICATION TO THE NORTHWESTERN US 5.1 Viscoelastic structure In order to evaluate the viscoelastic-cycle correction δvcycle(r, t) in the viscoelastic block model, as well as the various components of deformation in the viscoelastic-cycle model, we require viscoelastic Green's functions G(d)(r, r΄, 0+), G(d)(r, r΄, ∞), and |$\mathbf {\dot{G}}^{(d)} (\mathbf {r},\mathbf { r^{\prime }}, t)$|⁠. These are calculated using a combination of static deformation response functions (Pollitz 1996) and post-earthquake response functions (Pollitz 1997) using a suitable viscoelastic model. As a reference structure we choose the ‘northeast domain’ structure determined by Pollitz (2015) for the Mojave Desert region (Fig. 3). Derived from post-earthquake observations following the M7.1 1999 Hector Mine, California, earthquake, it represents the Basin and Range side of the region (Fig. 1). Although the heat flow northeast of the Hector Mine rupture is higher than that southwest of the rupture, the effective mantle asthenosphere viscosity is higher, possibly reflecting the joint effects of strain rate and mantle water content (Pollitz 2015). We shall explore several viscoelastic structures that have either higher or lower viscosities than this reference model. Figure 3. Open in new tabDownload slide Panels (a) and (b) show depth-dependent transient and steady-state viscosity, respectively, for the reference viscoelastic structure. The lower crust is assumed to be a Maxwell solid. Panels (c ) and (d) show the depth-dependent elastic bulk and shear modulus, respectively. Figure 3. Open in new tabDownload slide Panels (a) and (b) show depth-dependent transient and steady-state viscosity, respectively, for the reference viscoelastic structure. The lower crust is assumed to be a Maxwell solid. Panels (c ) and (d) show the depth-dependent elastic bulk and shear modulus, respectively. Along the Cascadia megathrust, dislocation sources may extend deeper than the base of the upper elastic layer, and these sources require a special treatment. Let |$r_{\rm elast}^{\prime }$| be the radius at the base of the elastic layer (in this case, at 20 km depth —Fig. 3). Then |$\mathbf {\dot{G}}^{(d)} (\mathbf {r},\mathbf { r^{\prime }}, t)$| (the integrand of the first contributing term to vinst(r)) for sources r΄ deeper than |$r_{\rm elast}^{\prime }$| would generally require the evaluation of post-seismic relaxation for sources that may be within a relaxing medium. In order to avoid the complexities of time-dependent deformation from such sources, we implement a time-independent formulation for the contribution of sources with |$r^{\prime } <r_{\rm elast}^{\prime }$| . In the viscoelastic-cycle formulation the contribution of these sources to vinst(r) takes the form \begin{equation} \mathbf {\bar{\rm v}}^{\rm cycle}( \mathbf {r}) = - \sum _n \int _{\Gamma _n| r^{\prime } <r_{\rm elast}^{\prime }} d^2 \mathbf { r^{\prime }}\, \mathbf { \dot{m}}(\mathbf { r^{\prime }}) : \left[ \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},0^+) \right] \end{equation} (16) which follows from eq. (3) making use of the property \begin{equation} \mathbf {G}^{(d)}(\mathbf {r},\mathbf { r^{\prime }},\infty ) = 0 \, \, \, {\rm for} \, \, \, r^{\prime } <r_{\rm elast}^{\prime } \end{equation} (17) The associated viscoelastic-cycle correction in the block model is simply \begin{equation} \delta \mathbf {\rm v}^{\rm cycle}(\mathbf {r}, t) = 0 . \end{equation} (18) 5.2 Choice of time-dependent versus time-independent earthquake cycle deformation In order to estimate the earthquake cycle correction δvcycle(r) needed in the viscoelastic block model, as well as specify explicitly ‘time-dependent’ contributions to the viscoelastic cycle model (i.e. those governed by the first term of eq. 14), it is necessary to provide estimates of slip rate, recurrence interval, and time of last rupture on selected faults. Pollitz et al. (2008) and Pollitz et al. (2010) implemented earthquake cycles with time dependence for numerous faults in the Western US, for example, different sections of the San Andreas fault, faults in the Eastern California Shear Zone, and faults in the Basin and Range (Fig. 1; fig. 6 of Pollitz et al.2010; table 3 of Pollitz et al.2008). For simplicity, we choose to include time-dependent cycle effects for only the Cascadia megathrust. In practice, this means that for the viscoelastic block model, only the Cascadia megathrust is included in the fault set used to calculate δvcycle(r) via eq. (7). We assume an average recurrence interval of Cascadia megathrust earthquakes of T = 500 yr (Adams 1990; Atwater & Hemphill-Haley 1997). We further assume that slip events on the Cascadia megathrust have taken place with a moment tensor density distribution |$\mathbf { \dot{m}}(\mathbf {r}) \times T$| at successive times T, 2T, 3T, …, back in time from the last event in the year 1700. Here |$\mathbf { \dot{m}}(\mathbf {r})$| is the moment tensor density rate associated with fault locking that is related to the slip rate components and fractional locking factor via eq. (13). 5.3 Viscoelastic block model inversion The Western US is initially divided into 137 blocks following Evans et al. (2015). These blocks span the entire plate boundary zone from the Salton Trough to the northern Cascadia subduction zone. The Cascadia magathrust is represented by the McCrory et al. (2006) slab-depth contours. The slab-depth contours define a set of 1916 triangular elements that accommodate shear dislocations to represent the descent of the Juan de Fuca slab beneath the overriding continental plate. Total Variation Regularization (TVR) is used to minimize the model misfit to a viscoelastic-cycle corrected GPS velocity field vobs(r) − δvcycle(r) simultaneously with grouping together blocks that have similar Euler vectors (Evans et al. 2015). This grouping determines the fault geometry, that is, the set of boundaries between blocks that are retained after similar blocks are grouped. This is performed by optimizing an objective function that includes a weighting parameter λ for dissimilarity of a set of Euler vectors {ωl} , that is, eq. (4) of Evans et al. (2015). The number of unique ωl, that is, the number of kinematically distinct blocks retained in a TVR optimization, usually decreases with increasing λ (Evans et al. 2015). We define the locking rate on a fault patch as the effective slip rate which contributes to elastic strain accumulation. The locking rate on the megathrust consists of a strike-slip component |$f( \mathbf {s} \cdot \hat{t} )$| and dip-slip component |$f ( \mathbf {s} \cdot \hat{d} )$|⁠, where f is the locking fraction and s is the relative plate velocity vector. The objective function includes a smoothness constraint on f (i.e. eq. 3 of Evans et al. 2015) but no bounds, that is, the constraint 0 ≤ f ≤ 1 is not strictly enforced (Evans et al. 2015). For a given viscoelastic structure, the corrected GPS velocity field is realized by employing successive estimates of δvcycle(r), which requires an iterative procedure. With |$\mathbf { \dot{m}}$| prescribed by eq. (13) for the current collection of Euler vectors {ωl} , eqs 7 and 9 show that the δvcycle(r) and vbackslip(r) components of model velocity vinst(r) depend on the set of unknown parameters {flk} , while from eq. (8) the vblock(r) component of vinst(r) depends upon the set of unknown Euler vectors {ωl} . Starting with initial values, a corrected GPS velocity field vobs(r) − δvcycle(r) can be constructed and fit to the model velocity field with TVR by solving for an updated set of parameters {flk} and {ωl} . This forms the basis for another iteration to solve for these parameters. 5.4 Viscoelastic-cycle model inversion For given viscoelastic structure and recurrence intervals of faults in Γn, the unknowns in the viscoelastic cycle model are generally the distribution of long-term locking (or seismic slip) rates embodied in |$\mathbf { \dot{m}}(\mathbf { r^{\prime }})$| (or m(r΄) through eq. (4)), creep (or aseismic slip) rates embodied in |$\mathbf { \dot{m}}^{\rm (cr)}(\mathbf { r^{\prime }})$|⁠, and the lateral distribution of δμ. The Cascadia megathrust is implemented via the first term of eq. (14), and for all other faults we implement time-independent (i.e. cycle-averaged) motions prescribed by the second or third term of eq. (14). For the latter set of faults we fix slip rates at the preferred slip rate values in table 2 of Pollitz et al. (2010). This reduces the set of unknowns to the distribution of |$\mathbf { \dot{m}}(\mathbf { r^{\prime }})$| on the Cascadia megathrust and the lateral distribution of δμ. On the Cascadia megathrust, we employ previously estimated angular velocity vectors to constrain the total slip rate. This stabilizes our procedure for determining the cycle correction δvcycle(r, t) and permits the effect of viscoelasticity to be assessed more directly. The total slip rate south of the mid-Olympic Peninsula (Fig. 1) is prescribed by the Oregon Coast (OC) to Juan de Fuca (JdF) plate angular velocity vector \begin{equation} \mathbf {\omega }_{\rm OC-JdF} = (12.70^{\circ }{\rm N}, -110.46^{\circ }{\rm E}, 0.554^{\circ }\,{\rm Myr^{-1}}). \end{equation} (19) The total slip rate north of the mid-Olympic Peninsula is prescribed by the Olympics (Olym) to Juan de Fuca (JdF) plate angular velocity vector \begin{equation} \mathbf {\omega }_{\rm Olym-JdF} = (33.43^{\circ }{\rm N}, -116.17^{\circ }{\rm E}, 1.389^{\circ }\,{\rm Myr^{-1}}). \end{equation} (20) The OC-JdF motion above is determined as the difference between the OC-NA angular velocity vector of (45.16°N, −119.01°E, −1.019° Myr−1) and the Juan de Fuca to North America (NA) angular velocity vector of (33.92°N, −115.34°E, −1.513° Myr−1) given in table 2 of McCaffrey et al. (2007). Similarly the Olym-Jdf motion is determined as the difference between the Olym-NA and JdF-NA angular velocity vectors given in table 2 of McCaffrey et al. (2007). These motions determine the vector velocities slk to be implemented in eq. (13). Together with eqs (14) and (4), this prescribes the forward model of instantaneous crustal velocity in terms of the distribution of locking fraction f on the slab interface, in addition to the other unknowns. The strain rate field is needed in order to relate lateral variations in shear modulus δμ or bulk modulus δκ to equivalent sources of deformation via eq. (15). It is obtained from the GPS velocity field using the method of appendix A of Pollitz & Vergnolle (2006), which is based on the method of Shen et al. (1996). The obtained strain rate pattern is very similar to that presented in fig. 4 of Pollitz et al. (2010). Figure 4. Open in new tabDownload slide Distribution of locking rate on the Cascadia megathrust, obtained from viscoelastic block model inversion with λ = 1600 using three values of ηasth: (a) 4.4 × 1018 Pa s, (b) 4.0 × 1019 Pa s and (c) 3.6 × 1020 Pa s. Superimposed are the Juan de Fuca plate boundaries and slab depth contours of McCrory et al. (2006) in intervals of 10 km. Figure 4. Open in new tabDownload slide Distribution of locking rate on the Cascadia megathrust, obtained from viscoelastic block model inversion with λ = 1600 using three values of ηasth: (a) 4.4 × 1018 Pa s, (b) 4.0 × 1019 Pa s and (c) 3.6 × 1020 Pa s. Superimposed are the Juan de Fuca plate boundaries and slab depth contours of McCrory et al. (2006) in intervals of 10 km. We use simulated annealing to infer the distributions of f and δμ subject to the constraint 0 < f < 1, with no constraint imposed on δμ as it may be either positive or negative. We impose no smoothing on these distributions but impose a small amount of squared amplitude damping on δμ. For a set of trial values of f and δμ, the only other free parameters—three independent rotation parameters that align vinst(r, t) optimally with a fixed North America reference frame (section 5.1 of Pollitz et al.2008)—are determined through least squares inversion. 6 RESULTS 6.1 Viscoelastic block model We derive a range of viscoelastic models from the reference model by scaling all viscosities by a constant factor; the resulting viscoelastic models can then be parameterized with just the viscosity in one depth interval, for example, the mantle asthenosphere viscosity ηasth, equal to steady-state viscosity η1 in the 43–220 km depth interval (Fig. 3b). For a given value of ηasth, inversions are performed for trial values of λ until a solution is attained with 30 blocks chosen by TVR. Two iterations of the procedure described in Section 5.3 are generally sufficient to converge on a TVR solution. For each of three values of ηasth, the resulting distributions of locking rate (i.e. the locking fraction times the total slip rate f|s|) obtained from the TVR inversion are shown in Fig. 4. For each case, the associated root-mean-square (RMS) misfit is 1.9 mm yr−1. Locking is generally maximal in the upper 20 km depth along the megathrust, and most sections deeper than 25 km have very small f and hence are continually slipping at close to the relative plate motion rate. The viscoelastic block models fit the data well, as illustrated for the case ηash = 4.0 × 1019 Pa s in the Pacific Northwest (Fig. 5) and other regions (Supporting Information Figs S2–S4). Figure 5. Open in new tabDownload slide Fit of viscoelastic block model to observed horizontal GPS velocity field in the Pacific Northwest using the reference viscoelastic model (ηasth = 4.0 × 1019 Pa s, Fig. 2). Figure 5. Open in new tabDownload slide Fit of viscoelastic block model to observed horizontal GPS velocity field in the Pacific Northwest using the reference viscoelastic model (ηasth = 4.0 × 1019 Pa s, Fig. 2). RMS residuals of the viscoelastic block models generally decreases with increasing number of blocks retained in the inversion (Fig. 6) regardless of viscosity. The residuals also depend upon viscosity, and the RMS pattern could form the basis for evaluating the trial viscoelastic structures. However, the pattern in Fig. 6 shows negligible differences in RMS among the different ηasth values for a fixed number of blocks. Although the viscoelastic block model for the lowest-viscosity case (ηasth = 4.4 × 1018 Pa s) performs worst for a number of block models, the model performance is essentially identical for all considered viscosities for the case of 30 retained blocks. Although 30 blocks are retained when inferring each of the locking distributions shown in Fig. 4, the block geometries differ slightly for different ηasth, for example, in Nevada (Fig. 4), depending on viscosity. We explored the alternative of fixing the model geometry at the 30 reference-model blocks defined in Evans et al. (2015). This results in negligible differences in inferred megathrust locking pattern. Figure 6. Open in new tabDownload slide RMS residual velocities of the time-dependent and uncorrected (time-independent) block models as a function of the number of blocks retained in the inversion. Figure 6. Open in new tabDownload slide RMS residual velocities of the time-dependent and uncorrected (time-independent) block models as a function of the number of blocks retained in the inversion. Using the distributions of locking factor obtained from a block model inversion for a given mantle asthenosphere viscosity ηasth, we may extract the cycle-velocity correction δvcycle(r) as well as the final corrected velocity field vobs(r) − δvcycle(r). These velocity fields are shown in Fig. 7 and Supporting Information Fig. S5, respectively, for the same trial values of mantle viscosity. Figure 7. Open in new tabDownload slide Cycle-velocity correction δvcycle(r) at GPS sites obtained using final distribution of Cascadia locking rate from block model inversions, three values of ηasth: (a) 4.4 × 1018 Pa s, (b) 4.0 × 1019 Pa s and (c) 3.6 × 1020 Pa s. Figure 7. Open in new tabDownload slide Cycle-velocity correction δvcycle(r) at GPS sites obtained using final distribution of Cascadia locking rate from block model inversions, three values of ηasth: (a) 4.4 × 1018 Pa s, (b) 4.0 × 1019 Pa s and (c) 3.6 × 1020 Pa s. 6.2 Viscoelastic cycle model For the same trial viscoelastic structures, we infer the distributions of Cascadia megathrust locking fraction f and perturbation in shear modulus δμ averaged vertically over the elastic plate thickness. The resulting locking rate and δμ distributions are shown in Figs 8 and 9 and Supporting Information Fig. S6, respectively. The RMS misfit of all models is 2.2 mm yr−1. The Cascadia megathrust is locked primarily at depth less than 25 km along nearly the entire ∼1300 km along-arc distance. This is similar to the inferred backslip distribution using the viscoelastic block model in section 6.1, as well as the backslip distributions obtained by McCaffrey et al. (2013) (their fig. 6) and Schmalzle et al. (2014) (their fig. 6). Figure 8. Open in new tabDownload slide Distribution of locking rate on the Cascadia megathrust (f × the relative plate motion rate |s|) obtained from the viscoelastic-cycle model optimization for three values of ηasth: (a) 4.4 × 1018 Pa s, (b) 4.0 × 1019 Pa s and (c) 3.6 × 1020 Pa s. Superimposed are major plate boundaries and the slab contours of McCrory et al. (2006) in intervals of 10 km. Figure 8. Open in new tabDownload slide Distribution of locking rate on the Cascadia megathrust (f × the relative plate motion rate |s|) obtained from the viscoelastic-cycle model optimization for three values of ηasth: (a) 4.4 × 1018 Pa s, (b) 4.0 × 1019 Pa s and (c) 3.6 × 1020 Pa s. Superimposed are major plate boundaries and the slab contours of McCrory et al. (2006) in intervals of 10 km. Figure 9. Open in new tabDownload slide Distribution of shear modulus perturbation averaged vertically over the elastic plate thickness obtained from the viscoelastic-cycle model optimization for ηasth = 4 × 1019Pa s. Faults from the USGS Quaternary faults database are superimposed. Figure 9. Open in new tabDownload slide Distribution of shear modulus perturbation averaged vertically over the elastic plate thickness obtained from the viscoelastic-cycle model optimization for ηasth = 4 × 1019Pa s. Faults from the USGS Quaternary faults database are superimposed. To illustrate the composition of the net velocity field in the western US, the model velocity field obtained on the reference model is divided into four separate components (Fig. 10)—those associated with the Cascadia earthquake cycles, an additional rotation to align the model velocity field with fixed North America, oceanic ridges and transform faults, all other faulting sources, and lateral variations in rigidity (i.e. that arising from the |$\mathbf { \dot{m}}^{( \delta \mu )}$| term of eq. 14). The net prediction of the viscoelastic model velocity on the reference model is the sum of those in Fig. 10. The fit of this model velocity field to the observed GPS velocity field is shown in Fig. 11 for the Pacific Northwest and Supporting Information Figs S7–S9 for other regions in the Western US. Fig. S10 summarizes this fit for the entire western US and includes the distribution of locking rate on model faults. It illustrates the dominance of the San Andreas fault in southern and central California and the patterns of distributed faulting in northern California, the ECSZ, and Basin and Range. Figure 10. Open in new tabDownload slide Viscoelastic-cycle model velocity field on the reference viscoelastic model shown for four separate components: (a) that associated with the Cascadia megathrust earthquake cycle, (b) an additional rotation to a fixed North America reference frame (eq. 8 of Pollitz et al. 2008), (c) ridges and transform faults associated with the Gorda, Juan de Fuca and Explorer plates, (d) all other faulting sources (table 2 of Pollitz et al. 2010) and (e) lateral variations in rigidity. The plotted velocity fields are subsampled for visual clarity. Superimposed are the model faults in green. Figure 10. Open in new tabDownload slide Viscoelastic-cycle model velocity field on the reference viscoelastic model shown for four separate components: (a) that associated with the Cascadia megathrust earthquake cycle, (b) an additional rotation to a fixed North America reference frame (eq. 8 of Pollitz et al. 2008), (c) ridges and transform faults associated with the Gorda, Juan de Fuca and Explorer plates, (d) all other faulting sources (table 2 of Pollitz et al. 2010) and (e) lateral variations in rigidity. The plotted velocity fields are subsampled for visual clarity. Superimposed are the model faults in green. Figure 10. Open in new tabDownload slide (Continued.) Figure 10. Open in new tabDownload slide (Continued.) Figure 11. Open in new tabDownload slide Fit of viscoelastic cycle model to observed horizontal GPS velocity field in the Pacific Northwest using the reference viscoelastic model (Fig. 2). Figure 11. Open in new tabDownload slide Fit of viscoelastic cycle model to observed horizontal GPS velocity field in the Pacific Northwest using the reference viscoelastic model (Fig. 2). There is practically no variation in inferred δμ pattern with ηasth (Supporting Information Fig. S6). For the case ηash = 4 × 1019 Pa s, as noted by Pollitz et al. (2010), the pattern of δμ in Fig. 9 may be equally well interpreted as either real variations in the seismic (i.e. high frequency) elastic rigidity, variations in elastic plate thickness in the absence of variations in elastic moduli within the elastic plate, active distributed deformation, or a combination of the three. The equivalence with active distributed deformation follows from the fact that steady deformation on faults prescribed by the |$\mathbf { \dot{m}}^{\rm (cr)}(\mathbf { r^{\prime }})$| term of eq. (14) can be equated with the volumetric sources |$\mathbf { \dot{m}}^{( \delta \mu )}(\mathbf { r^{\prime }})$| for a suitably dense distribution of faults. The hypothesis that the δμ variations represent primarily variations in effective elastic plate thickness is supported by the correlation of the inferred δμ pattern with mantle seismic velocity (e.g. Lin et al. 2009; Pollitz & Snoke 2010), negative δμ being correlated with relatively low seismic velocity that is presumably of thermal origin. The alternative interpretation of distributed deformation may be applicable to areas such as the ECSZ where long-term deformation is accommodated on a great number of fault strands which are not represented in the employed fault model, that is, the faults in Fig. 10(d). Identification of active distributed deformation is important in the context of seismic hazard analysis, as ‘off-fault’ deformation is thought to be an important component of the budget of strain accumulation and release in California (Field et al. 2013). It is noteworthy that lateral variations in rigidity account for ∼10 mm yr−1 net motion of the west coast with respect to the plate interior (Fig. 10e). 7 DISCUSSION 7.1 Distribution of locking on Cascadia megathrust Referring to Figs 4 and 8, the viscoelastic block model results and viscoelastic-cycle model results share consistent features. First, all estimates of locking rate are maximal in the upper 20 km depthwise regardless of ηasth. Next, the 30 mm yr−1 contour line of backslip rate tends to be deepest under the Olympic Peninsula. The distribution of locking rate on the Cascadia megathrust depends systematically upon the viscoelastic model. The average amount of locking in the upper 20 km tends to be greater the lower the viscosity. From eq. (7), the time-independent model corresponds to cycle velocity correction δvcycle(r) = 0 and is the high-viscosity limit of the time-dependent model. The time-dependent models have non-zero δvcycle(r). This correction, shown in Fig. 7 is generally directed ENE, that is, landward, so that the corrected velocity fields tend to have larger oceanward velocity than the uncorrected velocity field (Fig. 2). At first glance this appears inconsistent with the inferred greater amount of locking for ηasth = 4.4 × 1018 Pa s and 4.0 × 1019 Pa s relative to the time-independent block model. However, the strain field associated with δvcycle(r) is strongly WSW-ENE tension along the Washington and Oregon coast, so that the corrected velocity fields have more WSW-ENE contraction than the uncorrected velocity fields. Relative to the time-independent block model, this time-dependent velocity and strain pattern must be accommodated with a greater amount of shallow locking. This is verified by consideration of along-strike average locking rate in both shallow (22.5–0 km) and deeper (45–22.5 km) depth ranges (Figs 12a and c). The amount of shallow locking is greater at lower viscosity in order to produce greater interseismic contractile strain along the coast, but the amount of deeper locking is lower at lower viscosity in order to produce overall greater oceanward interseismic motions (Figs 12b and d). Figure 12. Open in new tabDownload slide Locking rate on the Cascadia megathrust averaged along strike for two depth ranges (a,c) or as a continuous function of depth (b,d) and for three values of ηasth. (a) and (b) are results from the viscoelastic block model, (c) and (d) results from the viscoelastic cycle model. In panels (a) and (c) the high-viscosity limit is plotted at viscosity 2 × 1021 Pa s for visual purposes. Figure 12. Open in new tabDownload slide Locking rate on the Cascadia megathrust averaged along strike for two depth ranges (a,c) or as a continuous function of depth (b,d) and for three values of ηasth. (a) and (b) are results from the viscoelastic block model, (c) and (d) results from the viscoelastic cycle model. In panels (a) and (c) the high-viscosity limit is plotted at viscosity 2 × 1021 Pa s for visual purposes. In the viscoelastic block model, a slight reversal in the trend of shallow (22.5–0 km) locking rate versus viscosity occurs at ηasth ≳ 3 × 1020 Pa s (Fig. 12a). This transition corresponds to a Maxwell relaxation time of 270 yr, which is of the same order as the recurrence interval (500 yr) and the elapsed time since the last event (300 yr). For ηasth ≲ 1020 Pa s, the viscoelastic block model-derived average locking rate is greater than that obtained from the time-independent (i.e. uncorrected) block model (Fig. 12a). However, the shallow average locking rate at any finite ηasth viscoelastic cycle model exceeds the average in the high-viscosity limit (Fig. 12b). This difference likely reflects the different methodologies used to derive the locking distributions. At shallower (22.5–0 km) depths, average locking rate along the Cascadia megathrust obtained in the viscoelastic block model (Fig. 12a) is about twice that obtained in the viscoelastic cycle model (Fig. 12c). This partially reflects differences in methodology, but we also find variability in viscoelastic block model locking rates depending on the constraints imposed on the backslip optimization. The results presented here use Laplacian smoothing but otherwise no constraint on the shape of the locking profile. An alternative viscoelastic block model further imposes slip at the full relative Juan de Fuca to Oregon Coast Block plate rate at the trench axis, and at the same level of misfit it results in similarly large shallow locking rates but smaller deeper (45–22.5 km) locking. Stronger variability in block model locking depending on a priori constraints is obtained by McCaffrey et al. (2013), that is, their pn1d and pn2d models have ∼50 per cent different average amplitude even though they fit the data practically equally; similar variation is obtained in the models of Schmalzle et al. (2014). Such variability generally arises because of the poor resolution of shallow locking rate in all considered models when using only land-based interseismic velocities. The viscoelastic-cycle model involves complex deformation around the Gorda plate, with variable slip rate on the Mendocino transform fault (Fig. 11) and reduced half-spreading rate on the Pacific-Gorda Ridge (14 mm yr−1) compared with 27.5 mm yr−1 on the Pacific-Juan de Fuca Ridge to its north (e.g. Wilson 1989). This complexity results from the internal deformation of the Gorda plate, which involves a combination of north–south contraction and east–west extension (Wilson 1989; Gulick et al. 2001; Chayton et al. 2004). We have followed section 4 of Pollitz et al. (2010) in assigning 10 mm yr−1 north–south shortening and 10 mm yr−1 east–west extension across this internal deformation zone. Although these details of the Gorda plate-area deformation are not included in the viscoelastic block model, resulting estimates of Cascadia locking in the southern part of the megathrust are similar (i.e. compare Figs 4 and 8). 7.2 Interplate locking rate uncertainty It is well known that the estimation of locking along a descending slab boundary is uncertain using land-based geodetic data alone, for example, around Japan (e.g. Nishimura et al. 2000; Suwa et al. 2006) and Cascadia (e.g. McCaffrey et al. 2007, 2013; Evans et al. 2015). For the Cascadia megathrust we quantify the uncertainty in locking rate estimates in two ways: (1) using the formal standard deviation in interplate locking resulting from the TVR block model, and (2) using the empirical uncertainty derived from the variance of several available models. In the latter case we use the time-independent models obtained from TVR viscoelastic block model inversion (Fig. 13a), viscoelastic cycle model (Fig. 13b), and models pn1d and pn2d of McCaffrey et al. (2013) (Figs 13c and d). Figure 13. Open in new tabDownload slide Distribution of locking rate on the Cascadia megathrust, obtained from various time-independent models: (a) TVR viscoelastic block model, (b) viscoelastic-cycle model, (c) model pn1d of McCaffrey et al. (2013) and (d) model pn2d of McCaffrey et al. (2013). Superimposed are major plate boundaries and the slab contours of McCrory et al. (2006) in intervals of 10 km. Figure 13. Open in new tabDownload slide Distribution of locking rate on the Cascadia megathrust, obtained from various time-independent models: (a) TVR viscoelastic block model, (b) viscoelastic-cycle model, (c) model pn1d of McCaffrey et al. (2013) and (d) model pn2d of McCaffrey et al. (2013). Superimposed are major plate boundaries and the slab contours of McCrory et al. (2006) in intervals of 10 km. The two estimates of uncertainty are shown in Fig. 14. The uncertainties exhibit the same trend in both cases, that is, locking rate estimates along the shallower part of the megathrust are more uncertain than those at greater depth. It is also apparent that the formal errors resulting from the TVR inversion are smaller than the empirically derived estimates by a factor of ∼3. The lower empirical and formal uncertainties at depth ≳ 30 km imply that the relatively low locking rates inferred at greater depth are robust. Figure 14. Open in new tabDownload slide Distribution of standard error in estimated locking rate on the Cascadia megathrust, obtained from the formal errors in least squares inversion in the time-independent viscoelastic block model (a) and from the variance in four models of the locking rate (Fig. 13). Superimposed are major plate boundaries and the slab contours of McCrory et al. (2006) in intervals of 10 km. Figure 14. Open in new tabDownload slide Distribution of standard error in estimated locking rate on the Cascadia megathrust, obtained from the formal errors in least squares inversion in the time-independent viscoelastic block model (a) and from the variance in four models of the locking rate (Fig. 13). Superimposed are major plate boundaries and the slab contours of McCrory et al. (2006) in intervals of 10 km. 7.3 Balance of time-dependent and time-independent motions In the viscoelastic block model, total interseismic deformation is the sum of a time-dependent term (δvcycle(r, t)) and time-independent terms (vblock(r) and vbackslip(r)). In the viscoelastic-cycle model, the time-dependent term for faulting sources is the first term of eq. (14), that is, that proportional to |$\mathbf {\dot{G}}^{(d)} (\mathbf {r},\mathbf { r^{\prime }},t - t_n + j T_n)$|⁠, and the time-independent terms are those proportional to G(d)(r, r΄, ∞) or [G(d)(r, r΄, ∞) − G(d)(r, r΄, 0+)]. Time-dependent deformation varies over a seismic cycle and generally involves large strain rates near the source fault early in its cycle and low strain rates late in its cycle (e.g. Savage & Prescott 1978; Pollitz 2001; Hearn et al. 2013). If observed contemporary strain rates in a given area are small, then either local long-term slip rates are small or local slip rates may be substantial but the local fault(s) is late in its earthquake cycle, as noted for specific cases by Chuang & Johnson (2011) (Garlock fault), Hearn et al. (2013) (Mojave and San Bernardino segments of the San Andreas fault) and Wang et al. (2012) (major subduction zones). Similarly, if contemporary strain rates are small and the time-dependent component of model deformation has large strain rates, then this must be balanced by the time-independent components of model deformation. All source faults of model deformation must be considered in such a balance. For example, observed east–west extension between southeastern Oregon and Wyoming is very small (Figs 5 and 11). The time-dependent model strain in this area from the megathrust earthquake cycle in Fig. 10(a), which involves substantial east–west extension, must be balanced by the time-independent model strain contributed by all other model components. The time-dependent east–west extension is balanced primarily by time-independent east–west contraction contributed by oceanic ridges and transform faults (Fig. 10c); little model strain is contributed by other faults (Fig. 10d) or distributed deformation (Fig. 10e). The summed model deformation (Fig. 11) consequently has little east–west strain in eastern Oregon and Washington. A similar example involves the amount of north–south contraction in eastern Oregon and Washington. We focus on the decomposition of time-dependent and time-independent velocity fields presented in fig. 12 of Pollitz et al. (2010) because this decomposition was discussed in section 3.5 of McCaffrey et al. (2013). (The time-dependent velocity fields depicted in fig. 12(a) of Pollitz et al. (2010) include contributions of Cascadia megathrust and San Andreas fault (and other minor fault) earthquake cycles; this distinguishes it from Fig. 10(a) of the present study, which includes only the contribution from the Cascadia megathrust earthquake cycle.) Little north–south contraction is observed, but the time-dependent motions in fig. 12(a) of Pollitz et al. (2010) reveal substantial north–south extension. This extension is balanced primarily by time-independent north–south contraction contributed by oceanic ridges and transform faults (fig. 12c of Pollitz et al.2010), and no other sources of deformation, for example, local faulting sources as proposed by McCaffrey et al. (2013), need be invoked to achieve this balance. 7.4 Interseismic vertical motions Vertical crustal motions are constrained by levelling surveys in western Oregon calibrated by estimates of sea level rise (Burgette et al. 2009). The resulting interseismic vertical motions involve coastal uplift up to 4 mm yr−1 (Supporting Information Fig. S11a). Although our inversions do not use the levelling data, the resulting model predictions (Supporting Information Fig. S11b) agree with these observations. Burgette et al. (2009) note that the lesser uplift rate between ∼45° and 43°N requires a shallowing of the locking depth, consistent with the locking models of Figs 4 and 8. Observed uplift rate is much higher south of 43°N because of the lesser distance of the coast to the trench. Our models replicate this observation and suggest that an even shallower locus of locking occurs south of 43°N. 8 CONCLUSIONS The horizontal GPS interseismic velocity field constrains the kinematics of crustal deformation in the Western US. Time-dependent models of the kinematics considered here are the viscoelastic block model and viscoelastic cycle model. The former is formally a special case of the latter in which the lithosphere may be divided into blocks that behave rigidly in the long term. Both models embody the dynamics of a layered lithosphere-asthenosphere system driven by viscoelastic relaxation of the ductile mantle asthenosphere and lower crust from repeated fault slip. The resulting time-dependent models constrain long-term fault slip rates and locking rates throughout the Western US. We focus on the implications of earthquake-cycle effects for the distribution of locking on the slab interface of the Cascadia subduction zone. For a range of mantle asthenosphere viscosity ηasth, we find that maximum locking is concentrated in the 10–20 km depth range of the slab interface and that the depth extent of significant locking (e.g. the 30 mm yr−1 contour line of locking rate) is maximal in the Puget Lowland (Fig. 1) section of the subduction zone. Inferred locking distributions are sensitive to ηasth, with shallow (0–22.5 km deep) locking tending to be greater, and deeper (22.5–45 km deep) locking tending to be smaller, the smaller is ηasth. These results, which are common to the viscoelastic block model and viscoelastic cycle model, suggest that time-dependent effects on the megathrust locking pattern are quantifiable and that epistemic uncertainty on the locking pattern could be reduced if the viscosity structure were sufficiently well known. Epistemic uncertainty remains even without earthquake-cycle effects because, in the context of time-independent models, much variability arises from different methodologies, including possible constraints on the shape or smoothness of the locking pattern. This variability arises to a large extent because of the exclusive use of inland geodetic data to inform these models, and it may be reduced by offshore measurements of interseismic velocities in the future. Acknowledgments We thank Rob McCaffrey for providing the pn1d and pn2d models. This paper benefitted from internal reviews by Roland Bürgmann and Jim Savage. We thank Bill Hammond, Wang Rongjiang and the Associate Editor Dr Kosuke Heki for their constructive criticisms. REFERENCES Adams J. , 1990 . Paleoseismicity of the Cascadia subduction zone: Evidence from turbidites off the Oregon–Washington margin , Tectonics , 9 , 569 – 583 . Google Scholar Crossref Search ADS WorldCat Aki K. , Richards P.G., 1980 . Quantitative Seismology , vol. 1 , Freeman & Co . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Atwater B. , Hemphill-Haley E., 1997 . Recurrence intervals for great earthquakes of the past 3500 years at northern Willapa Bay , Washington , Professional Paper 1576 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Burgette R.J. , Weldon R.J., Schmidt D.A., 2009 . Interseismic uplift rates for western Oregon and along-strike variation in locking on the Cascadia subduction zone , J. geophys. Res. , 114 , B01408 , doi:10.1029/2008JB005679 . Google Scholar Crossref Search ADS WorldCat Chayton J.D. , Goldfinger C., Dziak R.P., Fox C.G., 2004 . Active deformation of the Gorda plate: constraining deformation models with new geophysical data , Geology , 32 , 353 – 356 . Google Scholar Crossref Search ADS WorldCat Chuang R.Y. , Johnson K.M., 2011 . Reconciling geologic and geodetic model fault slip-rate discrepancies in Southern California: consideration of nonsteady mantle flow and lower crustal fault creep , Geology , 39 , 627 – 630 . Google Scholar Crossref Search ADS WorldCat Devries P.M.R. , Meade B.J., 2016 . Kinematically consistent models of viscoelastic stress evolution , Geophys. Res. Lett. , 43 , 4205 – 4214 . Google Scholar Crossref Search ADS WorldCat Evans E.L. , Loveless J.P., Meade B.J., 2015 . Total variation regularization of geodetically and geologically constrained block models for the Western United States , Geophys. J. Int. , 202 , 713 – 727 . Google Scholar Crossref Search ADS WorldCat Field E. et al. , 2013 . Uniform California Earthquake Rupture Forecast, Version 3 (UCERF3) - The Time-Independent Model, USGS Open-File Report, 13-1165 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Gulick S. P.S. , Meltzer A.S., Henstock T.J., Levander A., 2001 . Internal deformation of the southern Gorda plate: fragmentation of a weak plate near the Mendocino triple junction , Geology , 29 , 691 – 694 . Google Scholar Crossref Search ADS WorldCat Hammond W.C. , Thatcher E., 2004 . Contemporary tectonic deformation of the Basin and Range province, western United States: 10 years of observation with the Global Positioning System , J. geophys. Res. , 109 , B08403 , doi:10.1029/2003JB002746 . Google Scholar Crossref Search ADS WorldCat Hammond W.C. , Thatcher W., 2005 . Northwest Basin and Range tectonic deformation observed with the Global Positioning System, 1999–2003 , J. geophys. Res. , 110 , B10405 , doi:10.1029/2005JB003678 . Google Scholar Crossref Search ADS WorldCat Hearn E.H. , Pollitz F.F., Thatcher W.R., Onishi C.T., 2013 . How do “ghost transients” from past earthquakes affect GPS slip rate estimates on southern California faults? , Geochem. Geophys. Geosyst. , 14 , 828 – 838 . Google Scholar Crossref Search ADS WorldCat Hetland E.A. , Hager B.H., 2005 . Postseismic and interseismic displacements near a strike-slip fault: a 2D theory for general linear viscoelastic rheologies , J. geophys. Res. , 110 , B10401 , doi:10.1029/2005JB003689 . Google Scholar Crossref Search ADS WorldCat Johnson K.M. , Fukada J., 2010 . New methods for estimating the spatial distribution of locked asperities and stress-driven interseismic creep on faults with application to the San Francisco Bay Area, California , J. geophys. Res. , 115 , B12408 , doi:10.1029/2010JB007703 . Google Scholar Crossref Search ADS WorldCat Lin F.C. , Ritzwoller M.H., Snieder R., 2009 . Eikonal tomography: surface wave tomography by phase front tracking across a regional broad-band array , Geophys. J. Int. , 177 , 1091 – 1110 . Google Scholar Crossref Search ADS WorldCat Loveless J.P. , Meade B.J., 2011 . Stress modulation on the San Andreas fault by interseismic fault system interactions , Geology , 39 , 1035 – 1038 . Google Scholar Crossref Search ADS WorldCat McCaffrey R. et al. , 2007 . Fault locking, block rotation, and crustal deformation in the Pacific Northwest , Geophys. J. Int. , 169 , 1315 – 1340 . Google Scholar Crossref Search ADS WorldCat McCaffrey R. , King R.W., Payne S.J., Lancaster M., 2013 . Active tectonics of northwestern U.S. inferred from GPS-derived surface velocities , J. geophys. Res. , 118 , 709 – 723 . Google Scholar Crossref Search ADS WorldCat McClusky S.C. , Bjornstad S.C., Hager B.H., King R.W., Meade B.J., Miller M.M., Monastero F.C., Souter B.J., 2001 . Present day kinematics of the Eastern California Shear Zone from a geodetically constrained block model , Geophys. Res. Lett. , 28 , 3369 – 3372 . Google Scholar Crossref Search ADS WorldCat McCrory P.A. , Blair J.L., Oppenheimer D.H., Walter S.R., 2006 . Depth to the Juan de Fuca slab beneath the Cascadia subduction margin Ð A 3-D model for sorting earthquakes, USGS Data Series, DS-91, version 1.2 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Meade B.J. , Hager B.H., 2005 . Block models of crustal motion in southern California constrained by GPSÊ measurements , J. geophys. Res. , 110 , doi:10.1029/2004JB003209 . OpenURL Placeholder Text WorldCat Nishimura T. et al. , 2000 . Distribution of seismic coupling on the subducting plate boundary in northeastern Japan inferred from GPS observations , Tectonophysics , 323 , 217 – 238 . Google Scholar Crossref Search ADS WorldCat Petersen M.D. et al. , 2014 . Documentation for the 2014 Update of the United States National Seismic Hazard Maps, USGS Open-File Report, 14-1091 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Pollitz F. , McCrory P., Wilson D., Svarc J., Puskas C., Smith R., 2010 . Viscoelastic-cycle model of interseismic deformation in the northwestern United States , Geophys. J. Int. , 181 , 665 – 696 . OpenURL Placeholder Text WorldCat Pollitz F.F. , 1996 . Coseismic deformation from earthquake faulting on a layered spherical earth , Geophys. J. Int. , 125 , 1 – 14 . Google Scholar Crossref Search ADS WorldCat Pollitz F.F. , 1997 . Gravitational viscoelastic postseismic relaxation on a layered spherical earth , J. geophys. Res. , 102 , 17 921 – 17 941 . Google Scholar Crossref Search ADS WorldCat Pollitz F.F. , 2001 . Viscoelastic shear zone model of a strike-slip earthquake cycle , J. geophys. Res. , 106 , 26 541 – 26 560 . Google Scholar Crossref Search ADS WorldCat Pollitz F.F. , 2003 . Postseismic relaxation theory on a laterally heterogeneous viscoelastic model , Geophys. J. Int. , 155 , 57 – 78 . Google Scholar Crossref Search ADS WorldCat Pollitz F.F. , 2015 . Post-earthquake relaxation evidence for laterally variable viscoelastic structure and water content in the southern California mantle , J. geophys. Res. , 120 , 2672 – 2696 . Google Scholar Crossref Search ADS WorldCat Pollitz F.F. , Snoke J.A., 2010 . Rayleigh-wave phase-velocity maps and three-dimensional shear-velocity structure of the western US from local non-plane surface-wave tomography , Geophys. J. Int. , 180 , 1153 – 1169 . Google Scholar Crossref Search ADS WorldCat Pollitz F.F. , Vergnolle M., 2006 . Mechanical deformation model of the western United States instantaneous strain-rate field , Geophys. J. Int. , 167 , 421 – 444 . Google Scholar Crossref Search ADS WorldCat Pollitz F.F. , McCrory P., Svarc J., Murray J.R., 2008 . Dislocation models of interseismic deformation in the western United States , J. geophys. Res. , 113 , B04413 , doi:10.1029/2007JB005174 . OpenURL Placeholder Text WorldCat Prescott W.H. , Savage J.C., Svarc J.L., Manaker D., 2001 . Deformation across the Pacific-North America plate boundary near San Francisco, California , J. geophys. Res. , 106 , 6673 – 6682 . Google Scholar Crossref Search ADS WorldCat Puskas C.M. , Smith R.B., 2009 . Intraplate deformation and microplate tectonics of the Yellowstone Hotspot and surrounding western U.S. interior , J. geophys. Res. , 114 , B04410 , doi:10.1029/2008JB005940 . Google Scholar Crossref Search ADS WorldCat Savage J.C. , 1983 . A dislocation model of strain accumulation and release at a subduction zone , J. geophys. Res. , 88 , 4984 – 4996 . Google Scholar Crossref Search ADS WorldCat Savage J.C. , Prescott W.H., 1978 . Asthenospheric readjustment and the earthquake cycle , J. geophys. Res. , 83 , 3369 – 3376 . Google Scholar Crossref Search ADS WorldCat Savage J.C. , Lisowski M., Prescott W.H., 1981 . Geodetic strain measurements in Washington , J. geophys. Res. , 86 , 4929 – 4940 . Google Scholar Crossref Search ADS WorldCat Savage J.C. , Lisowski M., Prescott W.H., 1991 . Strain accumulation in western Washington , J. geophys. Res. , 96 , 14 493 – 14 507 . Google Scholar Crossref Search ADS WorldCat Savage J.C. , Simpson R.W., Murray M.H., 1998 . Strain accumulation rates in the San Francisco Bay area, 1972–1998 , J. geophys. Res. , 103 , 18 039 – 18 051 . Google Scholar Crossref Search ADS WorldCat Savage J.C. , Svarc J.L., Prescott W.H., 1999a . Strain accumulation at Yucca Mountain, Nevada, 1983-1998 , J. geophys. Res. , 104 , 17 627 – 17 631 . Google Scholar Crossref Search ADS WorldCat Savage J.C. , Svarc J.L., Prescott W.H., 1999b . Geodetic estimates of fault slip rates in the San Francisco Bay area , J. geophys. Res. , 104 , 4995 – 5002 . Google Scholar Crossref Search ADS WorldCat Savage J.C. , Gan W., Svarc J.L., 2001a . Strain accumulation and rotation in the Eastern California Shear Zone , J. geophys. Res. , 106 , 21 995 – 22 008 . Google Scholar Crossref Search ADS WorldCat Savage J.C. , Svarc J.L., Prescott W.H., 2001b . Strain accumulation near Yucca Mountain, Nevada , J. geophys. Res. , 106 , 16 483 – 16 488 . Google Scholar Crossref Search ADS WorldCat Savage J.C. , Gan W., Prescott W.H., Svarc J.L., 2004 . Strain accumulation across the Coast Ranges at the latitude of San Francisco , J. geophys. Res. , 109 , doi:10.1029/2003JB002612 . OpenURL Placeholder Text WorldCat Schmalzle G. , McCaffrey R., Creager K., 2014 . Central Cascadia subduction zone creep , Geochem. Geophys. Geosyst. , 15 , 1515 – 1532 . Google Scholar Crossref Search ADS WorldCat Shen Z.-K. , Jackson D., Ge B., 1996 . Crustal deformation across and beyond the Los Angeles Basin from geodetic measurements , J. geophys. Res. , 101 , 27 957 – 27 980 . Google Scholar Crossref Search ADS WorldCat Shen Z.K. et al. , 2003 . The SCEC crustal motion map, version 3.0 . Available at: http://epicenter.usc.edu/cmm3 . Suwa Y. , Miura S., Hasegawa A., Sato T., Tachibana K., 2006 . Interplate coupling beneath NE Japan inferred from three-dimensional displacement field , J. geophys. Res. , 111 , B04402 , doi:10.1029/2004JB003203 . OpenURL Placeholder Text WorldCat Svarc J.L. , Savage J.C., Prescott W.H., 2002a . Strain accumulation and rotation in western Nevada, 1993–2000 , J. geophys. Res. , 107 , 10.1029/2001JB000579 . OpenURL Placeholder Text WorldCat Svarc J.L. , Savage J.C., Prescott W.H., Murray M.H., 2002b . Strain accumulation and rotation in western Oregon and southwestern Washington , J. geophys. Res. , 107 , doi:10.1029/2001JB000625 . OpenURL Placeholder Text WorldCat Thatcher W. , Foulger G.R., Julian B.R., Svarc J., Quilty E., Bawden G.W., 1999 . Present day deformation across the Basin and Range Province, western United States , Science , 282 , 1714 – 1718 . Google Scholar Crossref Search ADS WorldCat Wang K. , Hu Y., He J., 2012 . Deformation cycles of subduction earthquakes in a viscoelastic Earth , Nature , 484 , 327 – 332 . Google Scholar Crossref Search ADS PubMed WorldCat Wech A.G. , Creager K.C., Melbourne T.I., 2009 . Seismic and geodetic constraints on Cascadia slow-slip , J. geophys. Res. , 114 , B10316 , doi:10.1029/2008JB006090 . Google Scholar Crossref Search ADS WorldCat Williams T.B. , Kelsey H.M., Freymueller J.T., 2006 . GPS-derived strain in northwestern California: termination of the San Andreas fault system and convergence of the Sierra Nevada Ð Great Valley block contribute to southern Cascadia forearc contraction , Tectonophysics , 413 , 171 – 184 . Google Scholar Crossref Search ADS WorldCat Wilson D.S. , 1989 . Deformation of the so-called Gorda plate , J. geophys. Res. , 94 , 3065 – 3075 . Google Scholar Crossref Search ADS WorldCat SUPPORTING INFORMATION Supplementary data are available at GJIRAS online. Figure S1. GPS data used in this study with 95 per cent confidence ellipses. There are altogether 1896 velocity vectors. Superimposed are the Juan de Fuca plate boundaries and slab depth contours of McCrory et al. (2006) in intervals of 10 km. Figure S2. Fit of viscoelastic block model to observed horizontal GPS velocity field in the San Francisco Bay area using the reference viscoelastic model (Fig. 3 of the main text). Figure S3. Fit of viscoelastic block model to observed horizontal GPS velocity field in Southern California using the reference viscoelastic model (Fig. 3 of the main text). Figure S4. Fit of viscoelastic block model to observed horizontal GPS velocity field in the Northern Basin and Range using the reference viscoelastic model (Fig. 3 of the main text). Figure S5. Corrected velocity fields vobs(r)–δvcycle(r) at GPS sites obtained using final distribution of Cascadia locking rate from block model inversions for three values of ηasth: (a) 4.4 × 1018 Pa s, (b) 4.0 × 1019 Pa s, and (c) 3.6 × 1020 Pa s. Velocity fields have been subsampled for visual clarity. Figure S6. Distribution of shear modulus perturbation averaged vertically over the elastic plate thickness obtained from the viscoelastic-cycle model optimization for three values of ηasth: (a) 4.4 × 1018 Pa s, (b) 4.0 × 1019 Pa s, and (c) 3.6 × 1020 Pa s. Figure S7. Fit of viscoelastic cycle model to observed horizontal GPS velocity field in the San Francisco Bay area using the reference viscoelastic model (Fig. 3 of the main text). Figure S8. Fit of viscoelastic cycle model to observed horizontal GPS velocity field in Southern California using the reference viscoelastic model (Fig. 3 of the main text). Figure S9. Fit of viscoelastic cycle model to observed horizontal GPS velocity field in the Northern Basin and Range using the reference viscoelastic model (Fig. 3 of the main text). Figure S10. Fit of viscoelastic cycle model to observed horizontal GPS velocity field in the Western US using the reference viscoelastic model (Fig. 3 of the main text). The GPS velocity fields are subsampled for visual clarity. Superimposed are the model faults shaded according to their locking rates. Figure S11. (a) Observed uplift rates along leveling lines in western Oregon (Burgette et al.2009). (b) Model uplift rates obtained with ηasth = 4 × 1019 Pa s. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Implications of the earthquake cycle for inferring fault locking on the Cascadia megathrust
JF - Geophysical Journal International
DO - 10.1093/gji/ggx009
DA - 2017-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/implications-of-the-earthquake-cycle-for-inferring-fault-locking-on-Yfgm6fZbnr
SP - 167
VL - 209
IS - 1
DP - DeepDyve
ER -