Mechanical deformation model of the western United States instantaneous strain-rate field

Summary

1 Introduction

Deformation in continental regions is commonly interpreted in terms of two end-member models (King 1994; Thatcher 2003). The first (‘block model’) views the lithosphere as composed of a number of microplates/blocks that behave rigidly over sufficiently long time intervals, the different blocks being separated by faults. The rigid behaviour of individual blocks is realized over a timescale that is much longer than the earthquake cycle associated with a typical fault. This view, originally conceived to explain geologic and palaeomagnetic data in many regions, has the flexibility to accommodate elastic strain accumulation effects over the interseismic period (e.g. Matsu'ura 1986). The second end-member model, known as the ‘thin sheet model’ (England & McKenzie 1982) accommodates the view that lithospheric deformation over length scales longer than the lithospheric thickness is essentially continuous and that over long time periods the lithosphere behaves as a viscous fluid. This model is generally applied to the thermally defined lithosphere, to which an effective viscosity can be derived that depends on the variation of temperature with depth and an assumed rheology of the lithosphere.

Although the relative merits of each end-member model are ardently debated (e.g. Tapponnier 2001), we believe that the complexity of crustal deformation phenomena over the totality of spatial and temporal scales of relevance demands a compromise between the two. The need for a more unified approach is highlighted by the inherent difference between short- and long-term deformation rates. Constraints on long-term deformation rates through fault slip rates and palaeomagnetic measurements of block rotations are often found incompatible with constraints on short-term deformation rates through GPS measurements and principal stress directions as inferred from seismic focal mechanisms. The existence of a ‘GPS–geologic’ discrepancy is documented in many cases in which the GPS velocity field around a major fault is not in accord with the corresponding geologic slip rate. Examples include the Altyn Tagh fault (Mériaux 2004; Wallace 2004) (GPS inferred rate of ∼9 mm yr⁻¹, geologic slip rate of ∼25 mm yr⁻¹), the Owens Valley fault (Dixon 2000) (GPS rate ∼7 mm yr⁻¹, geologic rate ∼2 mm yr⁻¹), the Garlock fault (Peltzer, 2001 and references therein) (GPS rate <∼2 mm yr⁻¹, geologic rate ∼7 mm yr⁻¹), the Agua Blanca and San Miguel-Vallecitos faults (Hirabayashi 1996; Dixon 2002) (GPS rate ∼2–3 mm yr⁻¹, geologic rate ∼6 mm yr⁻¹ for Agua Blanca fault; GPS rate ∼3–4 mm yr⁻¹, geologic rate ∼1 mm yr⁻¹ for San Miguel-Vallecitos fault), and the Wasatch fault (Friedrich 2003) (GPS rate ∼2.7 mm yr⁻¹, geologic rate 0.2–0.3 mm yr⁻¹).

It has been proposed in several of the above studies that apart from uncertainties in GPS measurements and fault slip rates, these discrepancies are to a large extent explained by the behaviour of a fault system during a viscoelastic deformation cycle. In several parts of the western United States (US) the characterization of active continental crust and mantle in terms of a relatively thin (∼15–30 km) mechanical lithosphere underlain by a ductile, relaxing ‘asthenosphere’ is supported by numerous studies of post-seismic relaxation (e.g. Pollitz 2000, 2001; Nishimura & Thatcher 2003) and crustal response to removal of lacustrine loads (Bills 1994) or lake filling (Kaufmann & Amelung 2000). In the case of a 2-D strike-slip regime, if the fault occupies an elastic upper layer underlain by a viscoelastic substrate, then analytic solutions are available to describe the evolution, accounting for the effects of periodic fault slip and subsequent viscoelastic relaxation compounded over many cycles (Savage & Prescott 1978; Pollitz 2001; Smith & Sandwell 2004). These solutions indicate that during the early part of a (periodic) cycle, crustal velocity around the fault is elevated above the average slip rate, while late in the cycle, crustal velocity is less than the average slip rate (Fig. 1). Many of the above-quoted discrepancies have been rationalized in terms of viscoelastic cycle behaviour, for example, Owens Valley fault: last earthquake in 1872, early in cycle (Dixon 2000, 2003); Wasatch fault: last earthquake ∼1200–1300 yr BP, late in cycle (Malservesi 2003). To these examples could be added many others for which viscoelastic relaxation effects early in the cycle are likely dominant, for example, elevated GPS velocities around the 1992 Landers and 1999 Hector Mine ruptures (Deng 1998; Pollitz 2000; Pollitz 2003a) or rapid uplift around the 1959 Hebgen Lake rupture (Nishimura & Thatcher 2003).

Figure 1

(a) A repeating earthquake occurs at a given location with periodicity T. The expected interseismic velocity field is the average velocity during an interseismic period (i.e. time 0⁺ to T⁻). The long-term velocity is the average over a long time interval, which includes several earthquakes. (b) Viscoelastic deformation cycle for a point located outside of the fault. Accumulated displacements over the interseismic period result from the superposition of the secular displacements and transient post-seismic displacements due to viscoelastic relaxation following periodic faulting events. The slope of the short dashed line represents the long-term velocity, the slope of the long dashed line the average interseismic velocity for high viscosity and the slope of the continuous line the interseismic velocity for low viscosity.

Open in new tabDownload slide

Pollitz (2003b) constructed a viscoelastic deformation cycle model and obtained a simple relationship between short-term deformation rates and average long-term fault slip rates. It was obtained by considering the average (or expected) behaviour of the lithosphere over the viscoelastic coupling cycle of a single deformation component and summing over all dislocation sources. The purpose of this paper is twofold:

1. To generalize the Pollitz (2003b) treatment further to account explicitly for additional types of sources, including viscoelastic relaxation effects on faults for which the slip history is sufficiently known and, in the thin-plate framework, steady slip along creeping portions of faults and

2. To apply the new treatment to a GPS data set in the western US.

The proposed model is a departure from the well-known block model, in which GPS strain gradients are primarily driven by slip beneath numerous defined fault zones and viscoelastic effects from past earthquakes produce little time-dependent behaviour (e.g. McCaffrey 2005; Meade & Hager 2005; d'Alessio 2005) (see also Supplementary Appendix B). Our model is essentially an extension of numerical models in which active continental deformation is produced by secularly increasing, horizontally transmitted tectonic forces either in isolation (e.g. Williams & McCaffrey 2001) or in superposition with transient viscoelastic effects (e.g. Roy & Royden 2000; Lynch & Richards 2001). Smith & Sandwell (2006) also apply a superposition of interseismic, post-seismic (viscoelastic), and steady-creep effects to describe the instantaneous velocity field along the San Andreas fault (SAF). Their model is more comprehensive than ours in accounting for relaxation of numerous moderate M∼ 7 events that may affect the present velocity field. However, their methodology differs in detail, particularly in how the plate boundary zone is loaded over time.

2 Quantification of Continental Deformation Rates from Dislocation Sources and Forces

2.1 Instantaneous velocity field

We assume that the continental lithosphere and underlying substrate may be divided into elastic and viscoelastic portions. The lithosphere deforms through the combined effect of elastic dislocations and forces applied to the elastic portions. The model constructed by Pollitz (2003b) describes dislocation sources in terms of the moment tensor density rate and the force density rate , which is considered to arise from an externally applied force acting on an isolated portion of lithosphere (it could be associated with convergence of two plates along a well-defined boundary). Let r′ be a variable denoting the location of a dislocation source and t denote time. The moment density m(r′, t) and force density f (r′, t) are assumed to be associated with repeating source(s) with definite mean periodicities, that is,

1

and a similar expression for f (r′, t) as a separable function of a space-dependent function f (r′) and a time-dependent function. In eq. (1), H is the Heaviside step function, t₀(r′) is the time of the last dislocation event at point r′, and T(r′) is the corresponding interevent time. In general T(r′) can be a distribution of interevent times that depends on the event index j, but we retain the form of eq. (1) for simplicity. Let V define a volume within which contributing dislocation sources and forces are acting. In the statistical sense, the expected interseismic velocity v(r) is given by (Pollitz 2003b, his eq. 7):

2

Here G^(d)(r, r′, t) and G^(f)(r, r′, t) represent the response of the viscoelastic system to the various components of dislocation sources and forces, respectively, applied at r′ and evaluated at point r and time t. The expressions in brackets represent the difference between the completely relaxed response and the initial elastic response. In a single-fault system, eq. (2) represents the average interseismic velocity between two successive earthquakes on that fault. It is to be distinguished from the long-term velocity field which is the average interseismic velocity plus the average rate of coseismic displacements produced over many earthquake cycles (Fig. 1).

For the forces likely to arise in the context of continental deformation, it is more realistic to work not with periodically increasing but rather continuously increasing forces, that is,

3

where t_initial is a long-past initial time (t−t_initial≫ 0). This force term is the fictitious force that would accumulate if the relative plate motion were not relieved by faulting. The total stress at the interplate boundary where this force is applied (or any other point within the continental interior) is the superposition of that due to the force term in eq. (3) and that from the viscoelastic cycle effects. These two types of contributing stress sources balance each other and produce, in principle, an average stress level that is bounded. With f (r′, t) given by eq. (3), the bracketed term with G^(f) in eq. (2) is replaced with G^(f)(r, r′, −t_initial). In the limit t_initial→−∞ we have

4

Eq. (4) provides a relationship between instantaneous velocity, observed during a relatively short time interval that does not include any dislocation events, and the rate of moment release and tectonic force accumulation applied to elastic portions of the lithosphere. Since it is valid in the statistical sense, it is most appropriate when estimates of past event activity (t₀(r′) and T(r′)) are not available.

However, with knowledge of past event history one can obtain a more general relationship. The space- and time-dependent response to dislocation sources of the form of eq. (1) can be written ( terms dropped for brevity)

5

For a single past dislocation event at a single point (e.g. j= 0 and fixed r′), eq. (5) represents the velocity field associated with relaxation of the viscoelastic earth (observed at point r and time t−t₀(r′) after the dislocation event). As a check upon eq. (4), the expected interseismic velocity field can be derived from eq. (5) by averaging the response over the period (t₀(r′) + 0⁺, t₀(r′) +T(r′)) separately for each dislocation source r′:

6

Carrying out the time integration and j−summation of the bracketed term and equating the quantity m(r′)/T(r′) with , the expression in eq. (6) results in the dislocation component of v(r) in eq. (4).

Note that eq. (5) depends on the viscoelastic structure because of the time dependence of , but eq. (4) does not depend on the entire viscoelastic structure but rather on the decomposition of the earth model into elastic and viscoelastic parts. If the viscoelastic structure is correctly specified, then eq. (5) is a more accurate representation of the velocity field. Depending on our knowledge of past slip history, one may choose to employ eq. (5) when such history is available and eq. (4) when it is unavailable. In general, eq. (4) is applicable to either continuous distributions of and or discontinuous distributions localized on a fault plane. By the nature of intended ‘slip events’, eq. (5) is applicable to discontinuous distributions of m(r′) that are associated with slip on one or more fault planes. It is thus useful to decompose V into a set of fault surfaces Γ_fault and the remaining volume V−Γ_fault (Fig. 2). One may further decompose Γ_fault into those fault surfaces Γ_n for which slip history is known and those fault surfaces Γ_m for which slip history is unknown. Thus Γ_fault={Γ_n}∪{Γ_m}. In the first case we define the triplet (Γ_n, t_n, T_n) to be the nth fault surface and the associated time of the last event and the interevent time. In this framework we form an estimate of interseismic velocity v^inst that accounts for known slip history using eq. (5) and otherwise falls back on the statistical average eq. (4). We may further account for the effect of steady fault creep and/or distributed steady slip in the bulk material. In this case we define a set of discrete creeping surfaces Γ_cr on which moment release rate is temporally constant. The remaining volume V−Γ_cr may accommodate steadily slipping dislocations at a rate that may be laterally variable. Putting together the contributions of all deformation sources yields

7

The first, second, and fifth terms in eq. (7) are to be integrated over a volume surrounding the fault surface Γ_n, Γ_m, or Γ_cr of vanishingly small thickness.

Figure 2

A volume V is assumed to contain both distributed and discrete dislocation sources that contribute to the instantaneous velocity field. Discrete dislocation sources are divided into two classes: (1) Γ_n on which the past slip rate , time of previous earthquake t_n, and interevent time t_n are considered known (or estimated) from geologic information, (2) Γ_m on which the long-term slip rate is possibly known but slip history is poorly constrained.

Open in new tabDownload slide

Eq. (7) may be thought of as a combined deterministic and statistical estimate of instantaneous velocity due to dislocation sources and forces applied to the lithosphere. The first term accounts exactly for past fault movements using space- and time-dependent viscoelastic response functions; m (r′) for r′∈Γ_n may be considered associated with a uniform slip s_n, and the corresponding average slip rate is . In the second term, past fault movements are accounted for with a statistical average; for r′∈Γ_m may be considered associated with a uniform slip rate . The third term represents the contribution of relaxation from distributed faulting with a statistic average, and the fourth term the contribution of externally applied forces. Finally, the fifth and sixth terms account for the effects of steady fault creep on discrete creeping surfaces and steady slip distributed over the remaining volume, respectively.

2.2 Long-term velocity field

At a timescale that is much larger than the interevent time of a typical fault in the system, one may obtain the average long-term velocity as the time-averaged instantaneous velocity field plus the average velocity field produced by the compounded coseismic displacement fields of all dislocation sources. The latter is given by

8

The time-averaged instantaneous velocity field is given by eq. (7) but with the fault elements Γ_n and Γ_m grouped together into the second term [because of the time-averaging process as carried out in eq. (6)]. This results in a long-term velocity field:

9

In contrast with the instantaneous velocity field (eq. 7), the long-term velocity field in eq. (9) depends only on the response functions calculated in the completely relaxed state.

3 Western us Velocity and Strain-Rate Fields

The instantaneous surface velocity field of the western US with respect to fixed North America (NA) is shown in Supplementary Fig. A1. It is a composite of the GPS velocity fields determined in nine separate USGS GPS surveys plus the WUSC velocity field determined by Bennett (1999) (version 002 of the WUSC velocity field, ) using continuous and campaign GPS data and VLBI data. The USGS campaign measurements are extracted from online sources ( and ) and are described in numerous prior publications (Savage 1998, 1999a, b; Thatcher 1999; Prescott 2001; Savage 2001a, b; Svarc 2002a, b; Hammond & Thatcher 2004; Savage 2004). The campaign measurements are generally conducted at intervals of 3 to 4 yr, and the associated velocity field is a composite of such measurements conducted between 1993 and 2003. The velocity field for the San Francisco Bay region is based upon not only USGS campaign measurements but also continuous GPS time-series from the CORS (Continuously Operating Reference Sites) and the BARD (Bay Area Regional Deformation) networks (Prescott 2001).

The WUSC velocity field is a composite of continuous and campaign GPS measurements conducted collectively between 1986 and 2000. Additional VLBI data used in the solution span the period 1979 to 1998. Data from the WUSC velocity field have been corrected by its authors for coseismic offsets of significant earthquakes. No correction for short-term post-seismic deformation has been applied to these data.

Each of the nine USGS campaigns data sets were processed at the USGS using the GIPSY/OASIS II software (Zumberge 1997). Velocities are provided in a fixed North American reference frame based on ITRF2000 (Altamimi 2002). Similarly, the WUSC velocity field is referenced to fixed North America. There are a total of 486 GPS velocity vectors contributed by the USGS campaign data and 260 velocity vectors contributed by the WUSC velocity field. The two data sources (USGS campaign; WUSC) have 84 common sites, and we determined a rotation between the two associated velocity fields that aligns the two velocity fields to within the measurement errors (generally ∼1 mm yr⁻¹ standard deviation in both east and north components for the USGS campaign measurements; ∼0.5 mm yr⁻¹ for the WUSC continuous measurements). The RMS of the difference at the 84 common sites are 1.7 mm yr⁻¹ and 1.1 mm yr⁻¹ for the east and north component. The velocity shift between the two data sets is practically a uniform translation of (−0.6 mm yr⁻¹ east, −1.4 mm yr⁻¹ north).

After correction of the velocity for the estimated effects of steady creep on the SAF system (Supplementary Fig. A2), the associated strain-rate field is depicted in Supplementary Figs A3–A6. We choose to model strain rate rather than the original GPS velocity field for several reasons. First, modelling strain rate avoids the issue of absolute reference frame of the GPS measurements. Second, our quantitative framework is not a kinematic model but rather a dynamic model in which instantaneous deformation rates are related to a series of underlying physical processes. This involves serious challenges when attempting to fit the GPS velocity field in detail. For example, the viscoelastic process that is of first order importance is a diffusive process with long-range effects. The potential inadequacy of the assumption of laterally homogeneous viscoelastic structure can lead to substantial errors in predicted velocity at long wavelength (i.e. far from the earthquake source region). Other unmodelled effects may include basal drag and more spatially variable forcing rates. Omission of such key ingredients can compromise a direct fit of the velocity field, and we judge that neglected processes are associated with long-wavelength effects that carry less weight when working with strain rate.

Because the GPS measurements represent a discrete sample of the total velocity field, it is necessary to smooth the underlying strain field in some manner. The procedure for deriving a strain field from the surface velocity field at a given spatial scale is a slightly modified version of the method of Shen (1996) and is described in Supplementary Appendix A. Supplementary Figs A3 and A4 show the velocity gradient field at spatial scales of 40 and 30 km, respectively, while Supplementary Fig. A5 shows a composite map derived from the previous two figures. Fig. 3 (also presented as Supplementary Fig. A6) represents the velocity gradient field of Supplementary Fig. A5 in terms of the magnitude and directions of the principal horizontal strain rate axes plus the rotational strain rate , defined as

10

where u and v are east and north velocity and γ and β measure distance in the due east and north directions, respectively.

Figure 3

Representation of western US strain-rate field in terms of 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 rotational strain rate (indicated by colour shading). There is no restriction on the standard deviations of the velocity gradient values. Grey lines indicate outlines of planes upon which Juan de Fuca–North America and Pacific–North America forces are imposed.

Open in new tabDownload slide

Fig. 3 reveals that, as one would expect, deformation rates are generally greatest near the North American plate boundary zones adjacent to the JdF and Pacific plates. The JdF–NA plate boundary zone deformation is characterized by primarily ENE–WSW shortening combined with strong clockwise rotation in the Pacific Northwest at ∼50–100 nanostrain yr⁻¹. The Pacific–North American plate boundary zone deformation is characterized by right-lateral shear strain parallel to the strike of the SAF zone and the Eastern California Shear Zone (ECSZ) at rates of ∼180 and ∼50–100 nanostrain yr⁻¹, respectively, combined with strong clockwise rotation, resulting essentially in deformation under simple shear. These characteristics have been remarked in many earlier studies (e.g. Savage 1999b; McCaffrey 2000; Savage 2001a, b; Svarc 2002a). Further in the plate interior in the Basin and Range Province, strain rates are characterized by WNW–ESE extension at relatively small rates ∼25–50 nanostrain yr⁻¹ (Hammond & Thatcher 2004).

4 Model Considerations

4.1 Rheology of the western US continental lithosphere

Most information about depth-dependent rheology in the western US is provided by studies of post-seismic relaxation (Pollitz 2001; Nishimura & Thatcher 2003; Pollitz 2003a; Freed & Bürgmann 2004), removal of lacustrine loads (Bills 1994) or glacial loads (James 2001), or lake filling (Kaufmann & Amelung 2000). Among these studies, those of Bills (1994), Nishimura & Thatcher (2003), and Kaufmann & Amelung (2000), all of which pertain to the Basin and Range province, prefer a rheology involving an elastic upper crust and relatively strong lower crust underlain by relatively weak mantle, implying that the uppermost portion of the crust (or the entire crust) is the strength-carrying portion of the lithosphere. Essentially the same conclusion is reached for the Mojave desert region (Pollitz 2001; Pollitz 2003a; Freed & Bürgmann 2004) and northwestern Washington (James 2001). In northwestern Nevada, the relative strength of the lower crust and upper mantle has not been conclusively resolved, and a broader range of rheologies may be consistent with available post-seismic deformation data (e.g. Hetland & Hager 2003), admitting the possibility of a relatively strong mantle lithosphere.

In Fig. 4 we present four candidate rheology models based on the results of the above studies. They are intended to be representative of the western US as a whole. All four models are characterized by an elastic plate thickness of 20 km, which could correspond to the upper crust plus a boundary layer occupying the upper portion of the lower crust. The remaining lower crust to a depth of 30 km and the underlying mantle are assumed Maxwell viscoelastic with viscosity η_c and η_m, respectively. In all cases, we assume a relatively weak mantle with low viscosity (5 × 10¹⁸ to 2 × 10¹⁹ Pa s) and a range of lower crust strength with low to moderate viscosity (2.5 × 10¹⁹ to 1 × 10²⁰ Pa s). Although lateral variations in rheology are almost certainly present in this large area based on variations in mantle seismic velocity (Humphreys & Dueker 1994) and other physical properties (Lowry 2000; Provost & Chéry 2006), we believe that these simple rheological models provide useful layered starting models that are representative of much of the western US. In addition, by considering how well the observed strain-rate field is matched in specific subregions of the western US, we may discriminate among these rheology models on a region-by-region basis.

Figure 4

Rheological stratification of four candidate models considered in this study, each characterized by an elastic upper layer that includes the upper crust and part of the lower crust underlain by the remaining viscoelastic lower crust and viscoelastic mantle. A Maxwell rheology is assumed for the viscoelastic layers with indicated viscosity values.

Open in new tabDownload slide

4.2 Characterization of forces driving plate boundary zone deformation

The western US exhibits active deformation over a wide range of spatial scales and tectonic regimes (Smith 1978; Zoback & Zoback 1989; Hammond & Thatcher 2004) (Fig. 3). The strike-slip regime of northern California is dominated by the SAF to the west, quasi-rigid block motion of the Sierra Nevada block to its east, and a combination of right-lateral strike-slip strain, with maximal shear trending ∼N35°W, and normal faulting with minimum principal stress axis ∼N70°W, in the northern continuation of the ECSZ. The deformation style becomes increasingly dominated by normal faulting further east in the Basin and Range province, the transition between the strike-slip and normal faulting regimes occurring roughly between the Central Nevada Seismic Zone (CNSZ) in western Nevada and the northern Walker Lane in northeastern California. These patterns yield to ENE–WSW horizontal compression in the Pacific Northwest related to JdF–NA subduction. These tectonic patterns reflect the roles of several tectonic driving forces and their interaction in different regions of the western US crust.

In order to quantify the influence of the bounding Pacific and JdF plates on the deformation of the North American plate, we note that an oceanic plate is characterized by relatively strong lithosphere compared with the continental lithosphere with which it is in contact. The oceanic mechanical plate thickness, assuming that it is associated with the ∼700°C isotherm, is only about 15–20 km for the relatively young oceanic lithosphere adjacent to the Pacific Northwest or California. This is based on the age-heat flow relationships provided by Stein & Stein (1992). This is even slightly less than the ∼15–30 km thickness of the continental mechanical lithosphere estimated in numerous localities in the western US (Section 4.1). However, we note that the rigidity of olivine is about 2.5 times as great as that of crustal materials. This great strength contrast ensures that the oceanic lithosphere will be highly resistant to internal deformation compared with adjacent continental lithosphere, even if the mechanical plate thickness is comparable. We therefore assume that, as a first approximation, the western US continental lithosphere possesses a thin mechanical lithosphere and responds passively to forces exerted by the bounding oceanic plates. An implication of this assumption is that the process of oceanic plate to continental plate interaction along their common interplate boundary may be well described as the response (of the continental lithosphere) to secularly increasing forces exerted on that boundary (Section 2.1, eq. 3).

Fig. 5 shows the two main interplate boundaries that affect western North American plate deformation. The first interplate boundary forms the western boundary of the SAF system which divides the Pacific plate from the NA plate boundary zone. The Pacific moves laterally with respect to the NA plate at a rate of about 48 mm yr⁻¹ (DeMets & Dixon 1999) parallel to the interplate boundary. The force of this interaction is opposed during the interseismic period by strain accumulation among the faults distributed throughout the broad western NA plate boundary zone. Although the role of deep dislocations in driving strain accumulation is often debated (e.g. Savage 1999b), we shall assume that this interaction is described purely through a distribution of horizontal forces along the interplate boundary. This loading mechanism is also implicit in those models, which describe the Pacific–NA interaction as the loading of a shear zone driven from the sides (e.g. Roy & Royden 2000; Lynch & Richards 2001).

Figure 5

Distribution of surfaces on which major oceanic plates (Pacific, Juan de Fuca) exert a force on the North American plate. The exerted force per unit area is assumed uniform along the given surface for each of the interplate boundaries. The corresponding forcing rate parameter (, etc.) is indicated for each force plane. The rectangular outlines indicate the 6 subregions discussed in the text (Section 6) and in Fig. 16.

Open in new tabDownload slide

The second interplate boundary is along the Cascadia subduction zone, which divides the JdF plate from the NA plate. In a similar manner to the SAF system, we assume that loading of the continental lithosphere along the subduction zone is described through a distribution of forces directed parallel to the JdF–NA relative motion resolved onto the slab interface, that is, directed down the slab. This defines a mechanism of interseismic strain accumulation at subduction zones that is similar to, but not identical with, the backslip model of Savage (1983). In the latter model, interseismic deformation within the continental lithosphere is driven by shear dislocations distributed along the interplate boundary, with magnitude equal to the negative of the coseismic slip divided by the recurrence interval (i.e. the negative of the long-term slip rate). Williams & McCaffrey (2001) implemented the framework of distributed forces (which they termed a ‘finite plate model’) to describe JdF to Cascadian forearc interactions. They found important differences between the finite plate model and the conventional backslip model when both are calibrated to fit geodetic data, such that stressing rates transmitted at the JdF–NA interplate boundary are smaller in the finite plate model.

The tectonic loading process is here assumed steady state. Cast in terms of eq. (7), it follows that we may approximate the steady-state loading process for both the Pacific and JdF interactions in terms of a distribution of forcing rates along the idealized continental shelf (SAF system) or the subduction interface (Cascadia system).

Note that the contribution of secularly increasing forces to the stress field is balanced in the long term by those contributions from coseismic and post-seismic stress changes associated with earthquakes. That is, the fourth term of eq. (7) and the remaining terms of eq. (7), each become arbitrarily large with increasing time, but the sum of all terms, in principle, remains finite. (In practice, however, imperfect specification of forcing rates and fault geometries or histories would lead to large stress fields at sufficiently large times.) In addition, horizontal forces arising from lateral gradients in gravitational potential energy likely play a role in driving western US active deformation (Flesch 2000). However, such forces do not play an explicit role in shaping the active deformation in our framework because, according to eq. (7), a constant force (i.e. ) does not contribute to the instantaneous velocity or strain-rate fields. Indirectly, these forces generate an absolute stress field, which promotes Basin and Range normal faulting and fault-perpendicular shortening around the SAF system. The moment release associated with these dislocation sources contribute to the instantaneous velocity field.

4.3 Model parametrization

Instantaneous deformation of the western US continental lithosphere is modelled as a superposition of the effects embodied in the various terms of eq. (7). Here we describe how these terms are parametrized.

Table 1 lists the parameters of force interaction between the oceanic Pacific and JdF plates and the continental NA plate. Note that the force interaction is meant to represent that between the oceanic plate and the portion of western US lithosphere with which it is in contact. In the case of the ‘JdF–NA’ boundary, the force vector is chosen appropriate for relative motion between the JdF plate and Cascadia forearc which is migrating northwards at about 10 mm yr⁻¹ with respect to fixed North America (Wells & Simpson 2001; McCaffrey 2000; Svarc 2002b). As noted by McCaffrey (2000) and Williams & McCaffrey (2001) and seen in Fig. 3 this leads to an horizontal strain rate in the Cascadia region dominated by ENE–WNW shortening. For that reason we have chosen a JdF–Cascadia forearc force interaction that is directed N65°E. In order to account for possible segmentation of the Cascadia margin (e.g. Trehu 1994), we subdivide the JdF–NA plate interface into four subplanes. The corresponding forcing rates, which are assumed uniformly distributed on each respective subplane, are denoted , and , as labelled in Fig. 4. The P-NA boundary is taken to coincide with the North American continental shelf, which is generally the western limit of significant faulting. This boundary of total length 1200 km, upon which P-NA forcing rates are applied, is located ∼100 to 200 km west of the SAF system. It consists of a total of four vertical planes which are grouped into a northern set, which is approximately locally parallel to the Pacific–Sierra Nevada/Great Valley (SNGV) relative motion direction (Argus & Gordon 2001), and a southern set, which is approximately locally parallel to the P-NA relative plate motion direction. The corresponding forcing rates are denoted , and , as labelled in Fig. 4.

Table 1

Parametrization of Tectonic Forces.

Open in new tabDownload slide

Table 2 lists the geometry and slip history associated with selected major and minor faults in the western US, including the date of last major rupture. It includes only those faults whose combined magnitude and slip history as such as to be deemed capable of contributing substantial viscoelastic relaxation signals to the present-day strain-rate field. The contributions of countless other faults will be accommodated in other ways, as described below and in Section 5.4. (It would be inappropriate to use the information in Table 2 to construct any budget of long-term slip. It is commonly the case that faults, which are important contributors to the long-term velocity field are practically opaque with respect to the interseismic velocity field. We refer the reader to Supplementary Appendix B for further elaboration.) Note that the 1999 Hector Mine event is omitted because the southern California GPS observations originate from Release 2 of the SCEC velocity field (), which was released in 1998.

Table 2

Parametrization of Major (Class 1) and Minor (Class 2) Faults.

Open in new tabDownload slide

For purpose of classification we refer to the major (M≥∼8.0) and minor (M≥ 6.8) ruptures as Class 1 and Class 2 faults, respectively. The three Class 1 and nine Class 2 faults considered here are indicated by red and purple line segments, respectively, in Fig. 6, and the corresponding slip parameters of the events are denoted {s_j |j= 1, … , 12} (Table 2). In most cases the recurrence time T is poorly constrained. However, for Class 1 faults such as the northern SAF (accommodating the 1906 earthquake) or the Cascadia subduction zone (accommodating the 1700 earthquake), T is well constrained. Although viscoelastic relaxation effects tend to be dominated by relaxation from the last event, the contributions from preceding events are generally important at great distance from the fault zone. The concept of cyclicity is dubious for many fault zones where earthquakes occur in clusters and that this may affect the expected viscoelastic response (Meade & Hager 2004). However, we believe that our framework accounts well to first order the viscoelastic effects generated by past sequences of earthquakes on at least the major fault zones.

Figure 6

Red and purple lines indicate major faults (Class 1 faults) and minor faults (Class 2 faults) associated with historic earthquakes with known slip history (Table 2). For each fault we calculate the contributions of viscoelastic cycle effects on the instantaneous strain-rate field. Green lines indicate active faults with less certain slip history (Class 3 faults) (Table 3) and corresponding a priori slip rates, for which we calculate averaged interseismic effects on the instantaneous strain-rate field. Blue lines indicate creeping faults slipping at the indicated rates. Superimposed are the outlines of the force planes from Fig. 5.

Open in new tabDownload slide

There are important fault strands parallel to the SAF, indicated by green lines in Fig. 6. There is known large moment release in the past 100 yr on the San Jacinto, Salton trough, and Imperial faults, and in the past 150–300 yr on the Macama and Rodgers Creek–Hayward faults. Several M6.5–7.5 earthquakes are involved on each fault strand, and some of these faults (i.e. Hayward fault; Imperial fault) are thought to be low-friction faults that may accommodate episodic creep. Therefore, we choose to incorporate the contributions of these faults to interseismic deformation using the assumption of uniformly distributed moment release during the seismic cycle on each respective fault segment. For the averaged interseismic velocity to be accommodated along these faults, it is then appropriate to utilize the second term of eq. (7). The moment rate density distribution for one of these fault segments is prescribed by the fault geometry and long-term slip rate. We fix the fault geometry and respective long-term slip rates with parameters as assigned in Table 3, where they are termed ‘Class 3 faults’. For the six listed faults these slip rates are denoted .

Table 3

Parametrization of Class 3 Faults.

Open in new tabDownload slide

Distributed faulting or steady slip within the western US lithosphere is evaluated here using a vertical average over a prescribed depth range and smooth functions to describe the horizontal dependence. Following Pollitz (2003b), one tensor component of moment release rate , which would be associated with a particular dislocation geometry over the volume V−Γ_fault (for faulting) or V−Γ_cr (for steady slip), is assumed laterally variable but uniform in depth from an upper depth d₁ to lower depth d₂. We define to be the vertically integrated moment release rate:

11

It is parametrized in terms of Hermite–Gauss functions. Letting :

12

where l+m≤l_max for fixed l_max= 20, the h_m are normalized Hermite polynomials such that

13

and L₁ and L₂ are proportional to the dimensions of the rectangular grid which, in our application, covers a 1112 × 1051 km² area. We choose values such that 1112 km/L₁ = 1051 km/L₂ equals the last local maximum of the HG function of degree l_max. Most of the HG functions so defined taper off smoothly at the edges of the rectangular area, and only the smaller-wavelength functions contain some signal near the edges. With l_max= 16 the above expansion involves 152 parameters.

4.4 Other sources shaping the active deformation field

Predicted deformation fields are further shaped by fault creep and known lateral variations in viscoelastic structure (i.e. the relatively thick SNGV lithosphere), each of which produces a first order effect on the predicted lithospheric response.

The effect of fault creep is specified on portions of the central San Andreas, Hayward, and Calaveras faults at rates ranging from 12 to 30 mm yr⁻¹. It is based on measured surface creep on the respective fault traces (refer to Pollitz & Nyst 2004, for a more complete description). In the present study creep rates of the central SAF are fixed at rates based on known surface creep rates and assumed constant from the surface down to specified depths. One complication is that when steady fault creep penetrates to the surface, as for the central SAF, the true velocity field is discontinuous, and the method of estimating the continuous velocity gradient field presented in Supplementary Appendix A breaks down. This is an issue when many GPS measurements are present on both sides of a creeping fault, as is the case for much of the San Francisco Bay area. Since the a priori creep considered for the central San Andreas, Hayward, and Calaveras faults all penetrate to the surface, the most practical way to account for a priori creep is to correct the GPS velocity field for this effect prior to estimating the continuous velocity gradient. The velocity field associated with steady creep (Supplementary Fig. A2) has been subtracted from the observed velocity field (Supplementary Fig. A1) prior to estimating the velocity gradient field.

The mantle lithosphere beneath the SNGV is thicker than that of surrounding western US lithosphere based on seismic tomography (e.g. Benz 1993; Humphreys & Dueker 1994). For simplicity, we assume that the mantle lithosphere beneath the SNGV block extends to 40 km depth, which constitutes a large contrast with respect to the surrounding lithosphere of assumed thickness 20 km. The effect of this lateral heterogeneity on the response to applied forces and post-seismic relaxation of the system can be estimated from first order perturbation theory (Pollitz 2003c). The lateral heterogeneity is represented as a contrast in depth-dependent shear modulus in the Laplace transform domain. The perturbation in the deformation field is then prescribed in three steps:

1. 1

   The response to tectonic forces and relaxation following earthquakes on the laterally homogeneous model is evaluated in the volume where the lateral heterogeneity is present,

2. 2

   By converting these ‘incident’ deformation fields into virtual sources of deformation within the laterally heterogeneous volume and

3. 3

   Evaluating the consequent effect on the deformation fields.

We find that this approach allows us to capture the relative rigidity of the SNGV block (e.g. Fig. 3) with a reasonably realistic model of the relatively thick SNGV lithosphere.

5 Interpretation of Western us Strain-Rate Field

We aim to characterize the tectonic forces which act upon the North American lithosphere, the rheology of the western US lithosphere-asthenosphere system as a whole, and the distribution of moment release. A key question is: How well does a laterally homogeneous rheological model capture the principal deformation characteristics of the instantaneous velocity field? To a large extent, evaluation of candidate rheologies is bound to the assumed fault geometry and slip associated with the earthquakes through their corresponding post-seismic relaxation signals. In order to systematically model the velocity field (or, equivalently, the velocity gradient field) with this complexity, we perform a succession of parameter estimations based on least-squares inversion of the observed velocity gradient field. In subsequent sections we describe the inversion procedure, results, and implications for characterizing the sources of western US active deformation.

5.1 Inversion for model parameters

Let α_k denote the collection of model parameters, including: forcing rates , slip values s₁, s₂, … , s₁₂, slip rates , and distributed-moment HG expansion coefficients a_ijlm. Using the notation of Supplementary Appendix A, suppose that we have a velocity gradient field {Ψ₁₁(r_i), Ψ₁₂(r_i), Ψ₂₁(r_i), Ψ₂₂(r_i) –i= 1, …I}. Let Ψ be a vector containing the collection of velocity gradient components at all I points, and let C be the a priori covariance matrix among these observables. In the inverse problem we minimize a functional of the form

14

where

15

In eq. (15), G_ik represents the Greens function response of the system at observable i to model parameter α_k. In eq. (14) the first term represents the data misfit, and the second term represents the integrated roughness of the lateral gradients in vertically integrated moment release rate, weighted by S; the integration in the roughness term is over the regions of distributed faulting V−Γ_fault and/or distributed steady slip V−Γ_cr.

Minimization of eq. (14) with respect to the model parameters leads to the normal equations

16

where q and k span the set of model parameter indices, and the subscript (i) means that the derivatives are evaluated using initial values of α_k, which we assume to be zero. Inversion of eq. (16) yields estimates of the model parameters and associated marginal covariances among them.

5.2 Estimation of tectonic forces

Initial estimates of the tectonic forcing rates may be obtained by jointly inverting the strain-rate field for the force parameters and all three Class 1 {s_j} (Tables 1 and 2). The resulting forcing rates are shown in Fig. 7 for each of the four rheological models. It is noteworthy that estimated forcing rates do not depend greatly on the assumed rheology. The value of is very small for all rheology models. This suggests that the coupling of the JdF and NA plates is very small in central Oregon. The corresponding stressing rates on the JdF–NA subduction interface or the Pacific-NA transcurrent interface are given by the forcing rate divided by the area of the segment, i.e,

17

where L_i and W_i for segment i are given in Table 1. On Model 4, for example, this yields , and . For the JdF–NA interaction {i= 1, … , 4}, these values are only about 10 to 50 per cent of the stressing rates of ∼2 kPa yr⁻¹ estimated by Williams & McCaffrey (2001) for the Cascadia megathrust. The disparity between estimated stressing rates is compensated by the inclusion of post-1700 relaxation in our model. This is seen by considering all of the contributions to the model strain-rate field on Model 4 (Figs 8a and c). At the present stage of the Cascadia seismic cycle, the contractile strain rate perpendicular to the Cascadia coastline is contributed primarily by post-1700 relaxation, and similarly for the stressing rate. This tendency is obtained for all four rheological models considered.

Figure 7

Tectonic forcing rates obtained through inversion of the strain-rate field for the forces parameters and all three Class 1 fault parameters {s_j}. Rheological models are given in Fig. 4.

Open in new tabDownload slide

Figure 8

For the inversion described in Section 5.3, the components of the strain-rate field on Model 4 contributed by: (a) relaxation following earthquakes on Class 1 faults using the ‘revised’ slip amplitudes; (b) relaxation following earthquakes on Class 2+3 faults using estimated slip amplitudes and a priori amplitude, respectively; (c) tectonic forces and (d) sum of all contributions.

Open in new tabDownload slide

5.3 Estimation of repeating slip and average slip rates

Unlike the forcing rates, repeating slip values s_j and average slip rate values are highly dependent on the rheological model. The influence of post-earthquake relaxation is demonstrated by successively adding one or more major faults to the set of deformation sources. In a first test, we hold fixed according to the values found previously (Fig. 7, Section 5.2) and invert jointly for and one Class 1 {s_j}, or for and all three Class 1 {s_j}. The overall fits to the strain-rate field for these cases are shown in Fig. 9. By examining the improvement in fit by the addition of relaxation from a single earthquake sequence, the figure shows that relaxation from repeating 1906 events has the greatest impact among the major faults tested. Fig. 10(a) shows the corresponding slip values obtained in the case where and all three Class 1 {s_j} are inverted simultaneously with held fixed. For Models 1, 2, and 3, estimated slip amplitudes for the 1857 and 1906 sources are generally very large, about 300–600 and 150–350 per cent of a priori slip values, respectively (Fig. 10a, Table 1). For these models, the estimated slip value for the 1700 source is consistently less than the a priori value. For Model 4, however, estimated 1700, 1857 and 1906 slip amplitudes are much closer to a priori values. These tendencies of the inverted slip values lead us to construct a set of ‘revised’ slip amplitudes for the Class 1 faults. The a prioris_j are scaled up by +10 per cent for the 1906 source, +25 per cent for the 1857 source, and −20 per cent for the 1700 source. Inversion for forcing rates with all three Class 1 {s_j} held at the revised values [‘revised’ case in Fig. 9 (left)] improves the fit relative to the case where the slip rates are held fixed at a priori values for all rheological models. Given the tenuous constraints on slip and magnitude of the last major rupture on the major source faults, the revised slip values are a plausible alternative.

Figure 9

(left) Fits to the strain-rate field of models derived by joint inversion for and one Class 1 s_j with held fixed. Two variations are considered: event magnitude is fixed (‘F’) to an a priori value (Table 2) or estimated (‘E’) with the inversion. Also shown is the result of joint inversion for and all three Class 1 faults with held fixed, including a case in which the slip magnitudes {s_j} are held fixed at ‘revised’ values (see text). (right) Fits derived by joint inversion for and {s_j} for Class 2 faults, with held fixed, {s_j} of Class 1 faults held fixed at ‘revised’ values and of Class 3 faults held fixed at a priori values (Table 3). Class 2 faults are included either by fixing (‘F’) the slip magnitude at a priori values (Table 2) or estimating the slip (‘E’). All fits are plotted relative to that fit obtained by inversion of the data set for alone.

Open in new tabDownload slide

Figure 10

(a) Ratio between inverted and a priori slip value of Class 1 faults obtained through inversion of the strain-rate field for and all three Class 1 s_j with held fixed. (b) Ratio between inverted and a priori slip values of Class 2 faults obtained through inversion of the strain-rate field for and Class 2 {s_j|j= 4, … , 12} with slip values of Class 1 faults {s_j|j= 1, 2, 3} held fixed at revised values and slip rate values of Class 3 faults held fixed at the a priori values given in Table 3. Rheological models are given in Fig. 4. (c) Estimated Class 1 & Class 2 slip ratio and standard deviation for inversion described above for Model 4.

Open in new tabDownload slide

The fit of the strain-rate field is further improved by inclusion of post-earthquake relaxation effects from Class 2 and Class 3 faults [Fig. 9 (right)]. For simplicity, the Class 3 slip rates are fixed at their respective a priori values (Table 3) in these tests. More precisely, in these tests and all nine Class 2 {s_j} are inverted simultaneously with fixed at the values determined in the previous section (Table 1), Class 1 slip values fixed at ‘revised’ values, and Class 3 slip rates fixed at a priori values. The resulting predicted strain-rate field is shown in Fig. 8(d). Most of the improvement is from the inclusion of Class 3 faults, with slight additional improvement when Class 2 slip magnitudes are estimated rather than fixed [Fig. 9 (right)]. Estimated Class 2 slip magnitudes (Fig. 10c) are not as well constrained as estimated Class 1 slip magnitude. The standard deviation in a Class 2 s_j estimate can be of the same order as s_j itself, in particular with Model 3 rheology and for Pyramid Lake and Olinghouse faults (s₁₁+s₁₂). It is difficult to discriminate the rheology with the faults located in the CNSZ {s_i |i= 4, 8, 9, 10}. However, regardless of the rheology, estimated slip amplitudes are close to (or exceed) a priori slip amplitudes for strike-slip fault in this area (Cedar Mountain fault, s₉) whereas they are reduced for oblique or normal faulting (Pleasant Valley, Dixie Valley and Fairview Peak faults; s₄, s₁₀, s₈, respectively) (Fig. 10b). These results are in good agreement with recent GPS results showing dextral shear motion with no dilatation in the westernmost Basin and Range (Hammond 2004). The estimated slip amplitude for the three faults in the south ({s_i |i= 5, 6, 7}, Owens Valley, Landers and Kern County faults, respectively) depend more strongly on the rheology, and estimated slip amplitudes with Model 4 show somewhat better agreement with the a priori values (Fig. 10b).

According to these results, estimated slip amplitudes for Class 1 and Class 2 faults are generally closer to the a priori values with Model 4 rheology than with Model 1, 2 and 3 rheologies. In addition, for Model 4 the inverted Class 1 slip values are much greater than the associated standard deviations, and inverted Class 2 slip values are statistically greater than zero to within one or two standard deviations (Fig. 10c). Thus Model 4 is the most consistent with the geological and seismological slip estimates for the twelve events associated with these faults. The Model 4 rheology also yields better fits to the data set than Models 1 and 2, although the differences are not statistically significant. The variance reduction obtained for Model 4 with these slip values, relative to a model with forcing rates alone, is 39.5 per cent. The good overall agreement between estimated and a priori Class 1 and Class 2 slip values for Model 4 leads us to prefer this rheology over the other rheology models. It is convenient to define ‘revised’ slip values for Class 2 faults using the Model 4 inversion results given in Fig. 9 as a guide. The revised slip values {s_i |i= 4 − 12} are assigned 50, 70 or 100 per cent of the a priori values, as indicated in Table 2.

5.4 Unaccounted deformation sources

When only background tectonic forces and relaxation from slip events on specified faults contribute to model strain rates, there are significant misfits to the strain rates. For all rheology models considered, predicted tensor strain rate around the SAF system is much lower than observed, and in the plate interior east of about 119°W it is too large. We admit two main possibilities for explaining the strain that is unaccounted for by the viscoelastic cycle model. The first is that steady deep slip in the lower part of the elastic upper lithosphere may be localized beneath some or all faults. The second is that (periodic) faulting events at locations distributed over the broad areas between the eighteen identified source faults may contribute additional viscoelastic relaxation signals.

The budget of moment accumulation and release along the Pacific–North American plate boundary system also demonstrates the need to complement the set of eighteen considered faults (Table 2) with additional deformation sources. Slip accumulation at a rate of 50 mm yr⁻¹ along a 1000-km long, 20-km wide plate boundary with average shear modulus of 37 GPa yields a moment accumulation rate of 3.7 × 10¹⁹ N m yr⁻¹. This is about 2.5 times larger than the moment release rate of 1.47 × 10¹⁹ N m yr⁻¹ arising from earthquake-cycle deformation based on slip values in Table 2, excluding the Cascadia megathrust. The remainder must be made up with distributed moment release on additional faults, steady deep slip, or a combination of the two. In Section 5.5 we consider the former, in Section 5.6 the latter, and in Section 5.7 a combination of the two. It is worth noting that there is theoretically no requirement in crustal dynamics that there be any deep slip beneath faults. For example, the viscoelastic coupling model of Savage & Prescott (1978) provides a comprehensive explanation of the seismic cycle of an idealized fault system without any deep dislocations. Our introduction of steady slip in the lower elastic lithosphere reflects the inadequacy of viscoelastic-cycle deformation to account for all observed strain localization.

5.5 Relaxation from distributed faulting

The space of possible moment release includes right-lateral strike slip on NW–SE-trending faults and normal slip on NNE–SSW-trending faults in the California-Nevada area. The geometry of involves pure normal slip on a 45°-dipping, N20°E-striking plane; the geometry of involves right-lateral strike-slip motion on a vertical N40°W-striking plane. We use the parametrizations for and given by eqs (11) and (12) in Section 4.3. We assume uniformity of moment release rate from the surface to the base of the elastic lithosphere (20 km) and hence d₁= 0 km and d₂= 20 km. These moment release rates are related to the velocity field, and hence strain rate, field via the third term in eq. (7). We invert the instantaneous strain-rate field for the distributions of and jointly with assuming the Model 4 rheology. Because of the tradeoffs with Class 1 slip magnitudes and Class 3 slip rates, these parameters are held fixed at the ‘revised’ slip values (Table 2) and a priori slip rates (Table 3), respectively. The Class 2 slip amplitudes are held fixed at the revised values specified in Table 2. The roughness weight S in eq. (15) is chosen at a value that results in a modest contribution of distributed faulting to the overall moment release budget. The obtained vertically integrated moment rate distributions are shown in Fig. 11(a). The estimated pattern is not consistent with the tectonic environment around the SAF, but the signal associated with it is small. The much larger is positive over most of the western US and reaches a maximum of ∼0.6 × 10¹⁴ N m (km² yr)⁻¹ around the SAF. The standard error in a point estimate of or is about 0.1 × 10¹⁴ Nm (km² yr)⁻¹. Distributed strike-slip faulting is thus formally well resolved above the noise level, while distributed normal faulting is not well resolved. Figs 11(b) and (c) show the contributions of relaxation from distributed faulting and the Sierra Nevada perturbation, respectively. The latter acts to account for the quasi-rigidity of the Sierra Nevada block lithosphere (Section 4.4).

Figure 11

(a) Distributed moment release rates and after inversion of the strain-rate field for distributed faulting and with Class 1, 2 and 3 faults held fixed to ‘revised’ or a priori values (see text, Sections 5.3 and 5.5). Isolines are each 0.2 × 10¹⁴ N m (km² yr)⁻¹. and are related to the strain-rate field via the third term of eq. (7) and are associated with relaxation following hypothetical slip events distributed throughout the 20 km-thick elastic lithosphere. The contributions of relaxation from distributed faulting (B) and the Sierra Nevada block perturbation (c) to strain-rate field are shown in terms of the amplitudes and directions of the principal strain-rates axes and rotational strain rate.

Open in new tabDownload slide

The variance reduction of this model is 66.5 per cent, an improvement of 48 per cent over the model obtained with forcing rate alone, an improvement of 28 per cent over the modelled obtained with the forcing rate and the revised Class 1 parameters (Section 5.2), and an improvement of 14 per cent over the model obtained with revised deformation parameters without distributed moment release (Section 5.3). The total calculated strain-rate field (Fig. 12) is correspondingly an improvement over that obtained previously, for example, with a priori/revised deformation parameters and without distributed moment release (Fig. 8d). Strain rate magnitudes over the SAF system, the Pacific Northwest, and the Basin and Range Province, while still smaller than observed (Fig. 3), are of higher amplitude with distributed moment release than without it. Similarly, strain rate magnitudes east of 119°W are of smaller amplitude with distributed moment release than without it, in better accord with observation.

Figure 12

Total calculated strain-rate field on Model 4, described in Section 5.5. It takes account for the relaxation following earthquakes on Class 1, 2 and 3 faults, the tectonic forces, and the relaxation from the distributed faulting and Sierra Nevada block perturbation (Figs 11b and c).

Open in new tabDownload slide

5.6 Distributed steady deep slip

In the case of steady deep slip we define moment release rates over the depth range 15–20 km (any depth range concentrated near the base of the elastic crust would suffice). For simplicity we focus on the single distribution ; 15 km, 20 km). This quantity is parametrized with eqs (11) and (12) and is related to the velocity field via the sixth term of eq. (7). We invert simultaneously for and all forcing rates assuming the Model 4 rheology. As in the previous section, Class 1 slip magnitudes and Class 3 slip rates are held fixed at the ‘revised’ slip values (Table 2) and a priori slip rates (Table 3), respectively, and the Class 2 slip amplitudes are held fixed at the revised values specified in Table 2.

The resulting pattern of deep slip in Fig. 13(a) is concentrated primarily around the SAF system. The contributions of distributed deep slip and the Sierra Nevada perturbation are shown in Figs 13(b) and (c). The contribution to the strain-rate field in Fig. 13(b) not only helps localize strain beneath the SAF system but also removes excess strain from the Basin and Range province. In this model, the integrated moment release rate from aseismic slip in the deep lithosphere is ∼1.3 × 10¹⁹ Nm yr⁻¹, roughly equal to the moment release rate from earthquakes (excluding the Cascadia event). This is equivalent to the moment release rate that would be generated by 50 mm yr⁻¹ slip on a 1000-km-long, 5-km-wide vertical dislocation at the base of the elastic layer. This would imply either aseismic slip in the lower 5 km of the lithosphere at the Pacific–North American relative plate motion rate or, invoking perturbation theory (Pollitz 2003c), that the lower 5 km of the lithosphere along an equivalent 1000 km length is absent in the western part of the plate boundary zone. The total predicted strain field in Fig. 14 is similar to that generated with relaxation from distributed faulting (Fig. 12).

Figure 13

(a) Distributed moment release rates after inversion of the strain-rate field for distributed steady deep slip and with Class 1, 2 and 3 faults held fixed to ‘revised’ or a priori values (see text, Sections 5.3 and 5.6). Isolines are each 0.2 × 10¹⁴ N m (km² yr)⁻¹. is related to the strain-rate field via the sixth term of eq. (7) and is associated with steady slip in the lower elastic lithosphere. The contributions of the distributed deep slip (b) and Sierra Nevada block perturbation (c) to strain-rate field are shown in terms of the amplitudes and directions of the principal strain-rates axes and rotational strain rate.

Open in new tabDownload slide

Figure 14

Total calculated strain-rate field on Model 4, described in Section 5.6. It takes account for the relaxation following earthquakes on Class 1 faults (Fig. 8a), the relaxation following earthquakes on Class 2+3 faults (Fig. 8b), the tectonic forces, and the steady slip in the lower elastic lithosphere and Sierra Nevada block perturbation(Figs 13b and c).

Open in new tabDownload slide

5.7 Combined distributed faulting and steady deep slip

We perform a joint inversion for distributed faulting (which contribute to the deformation through relaxation via the third term in eq. 7), distributed steady deep slip (the sixth term of eq. 7), and all forcing rates assuming the Model 4 rheology. As before, Class 1 slip magnitudes and Class 3 slip rates are held fixed at the ‘revised’ slip values (Table 2) and a priori slip rates (Table 3), respectively, and the Class 2 slip amplitudes are held fixed at the revised values specified in Table 2. The roughness weight S in eq. (15) is chosen such that the summed moment release rate of all fault and steady-slip sources equals the moment accumulation rate of 3.7 × 10¹⁹ Nm yr⁻¹.

The distributions of moment release from faulting-related and steady deep-slip are shown in Fig. 15(a), and the corresponding contribution to the strain-rate field is shown in Fig. 15(b). Combined with the contributions of relaxation following discrete faulting events (Figs 15d and e) and tectonic forcing (Fig. 15f), the additional dislocations sum to a total predicted strain-rate field (Fig. 15g) that agrees well with the observed strain field (Fig. 3) in both high- and low-strain regions. The variance reduction achieved by this model is 70.9 per cent. We note that steady slip in the lower elastic layer beneath the SAF in southern California was also inferred by Pollitz (2001) from a GPS profile using a viscoelastic coupling model with repeating 1857 Fort Tejon-type events. Although marginally resolved here as a localized feature, the pattern of deep slip between the Mojave Desert and Owens Valley is consistent with observed strain localization in the southern ECSZ (Peltzer 2001).

Figure 15

(a) Distributed moment release rates and after inversion of the strain-rate field for distributed faulting, distributed steady deep slip, and with Class 1, 2 and 3 faults held fixed to ‘revised’ or a priori values (see text, Sections 5.3 and 5.6). Isolines are each 0.2 × 10¹⁴ N m (km² yr)⁻¹. The combined contributions of distributed faulting and distributed deep slip (b) and Sierra Nevada block perturbation (c) to strain-rate field are shown in terms of the amplitudes and directions of the principal strain-rates axes and rotational strain rate. For the inversion described in Section 5.7, the components of the strain-rate field on Model 4 contributed by: (d) relaxation following earthquakes on Class 1 faults using the ‘revised’ slip amplitudes; (e) relaxation following earthquakes on Class 2+3 faults using estimated slip amplitudes and a priori amplitude, respectively and (f) tectonic forces. Parts D and E are identical to Figs 7(a) and (b), respectively. (g) sum of all contributions shown in parts (c) to (f). Thick and thin green line segments denote contractile and tensile principal strain axes, respectively, at selected points for visual clarity.

Open in new tabDownload slide

The existence of additional faulting sources as depicted in Fig. 15(a) is supported by the occurrence of other historical earthquakes. In southern California additional sources could include repeating 1812-type sources, which ruptured at least the Pallett Creek and Wrightwood localities on the southern SAF (Fumal 2002b). Similarly, the slip budget of southern California south of Wrightwood demands a total of about 50 mm yr⁻¹ long-term slip across the fault system, but the San Jacinto and Landers source faults included here represent only a fraction of the expected total. One could append to these two faults others which accommodate a substantial fraction of the long-term slip rate (and which would complete the budget of expected long-term slip), but the potential viscoelastic relaxation signals from these faults (i.e. their contribution to the present interseismic velocity field) are judged to be very small. Notably, the SAF between Wrightwood and the Coachella Valley, which last ruptured in the late 17th century with modest slip (Sieh & Williams 1990; Fumal 2002a; McGill 2002), is expected to contribute only small viscoelastic relaxation signals to the present deformation field, and, therefore, much of the episodic slip contributed by the southern SAF may be difficult to detect (see Supplementary Appendix B). However, Fig. 15(a) suggests that steady deep slip beneath the southern SAF accounts for part of the slip budget, and it is manifested as a substantial localized strain rate (Fig. 15b).

The slip budget in the San Francisco Bay area also demands the additional sources. Long-term slip rates are (Table 2) 22 mm yr⁻¹ for the SAF and 9 mm yr⁻¹ for the Hayward-Rodgers Creek-Macama fault system; the sum of 31 mm yr⁻¹ is only about 80 per cent of the Pacific-Sierra Nevada relative motion rate (Argus & Gordon 2001), demanding additional deformation sources on the north-central SAF.

Finally, thrust faulting accommodating northward convergence between the Pacific and North American plates in the Los Angeles region is not included in any of our models. Account for numerous M∼ 7 events such as the 1971 San Fernando and 1994 Northridge events would qualitatively increase the predicted north–south contractile strain and better match the observed strain rate in that region.

6 Discussion

In the tests for which Class 1 slip values are variable (Fig. 10a, Section 5.2), the estimated slip amplitude on the 1857 rupture is ∼500 per cent and ∼300 per cent greater than the a priori slip with Model 1 and 2 rheologies, respectively. There is a possibility that 1857 slip was much larger than given by palaeoseismic estimates (Sieh 1978) as proposed by Runnerstrom (2002) for the Cholame segment. Average slip near 14 m rather than 6 m would be consistent with the strain build-up expected between the previous event ∼1480 and 1857, given a long-term geologic slip rate of ∼3.4 cm yr⁻¹ (Sieh 1984). If true, then Models 1 and 2 would be more plausible in southern California. However, there is no direct evidence for slip values near 14 m except possibly on the Cholame segment, and therefore, we prefer the results derived assuming an average 7.5 m slip, which is consistent with palaeoseismic evidence (Grant & Sieh 1993) and which favours Model 4 in southern California.

In addition to their impact on the various slip estimates and global data fits, the candidate rheology models may be evaluated by the fits to specific subregions. In the following we refer to the results of joint inversion for and {s_j} for Class 2 faults, with {s_j} of Class 1 faults held fixed at ‘revised’ values (Table 2) and of Class 3 faults held fixed at a priori values (Table 3). Fig. 16 shows the fits to the strain-rate field in each of six subregions and for the entire data set. The Model 4 rheology (η_m= 2 × 10¹⁹ Pa s and η_c= 1 × 10²⁰ Pa s) fits the entire data set slightly better than the Models 1 and 2 rheologies, and it better fits most of the subregional data sets, particularly western Nevada and southern California around the Landers, Kern County, Owens Valley rupture zones. On the other hand, the Oregon subregion and the San Francisco Bay subregion strain-rate fields are better fit with the Model 1 rheology that has a weaker upper mantle and lower crust (η_m= 1 × 10¹⁹ Pa s and η_c= 2.5 × 10¹⁹ Pa s). This suggests, first, that different regions are governed by different viscosity structures, pointing to lateral heterogeneity in viscosity structure in the western US. Second, the average 1-D long-term viscosity structure in the western US found in this study shows higher upper mantle and lower crust viscosity than in smaller regions in this area over shorter timescale (e.g. Dixon 2004). This suggests possibly multiple material timescales in the relaxing portions of the Earth. Pollitz (2003a) suggests two material time constants ∼0.07 and ∼2 yr for the mantle rheology beneath the Mojave Desert, based on modelling 2.5 yr of post-seismic time series after the 1999 Hector Mine earthquake with a Burghers body rheology. At longer times, however, it is conceivable that a broader spectrum of relaxation times, including times of order 20 yr as implied by Model 4, are needed to adequately describe the viscous component of the complete rheology.

Figure 16

Fits to specific subregions (Fig. 5) derived by joint inversion for and {s_j} for Class 2 faults, with held fixed, {s_j} of Class 1 faults held fixed at ‘revised’ values and of Class 3 faults held fixed at a priori values (Table 3). The fit to the entire data set is shown on the extreme left.

Open in new tabDownload slide

Fig. 17 shows the predicted GPS displacement field based on the combination model described in Section 5.7. Qualitatively it reproduces the gross features of the observed velocity field (Supplementary Fig. A1), particularly the large right-lateral shear strains along the SAF system and the east–west contraction and rotation of the Cascadia region. Many details of the predicted velocity pattern, however, do not conform to the observed velocity field. This includes the observed rapid decrease in velocity as one moves inland from the Pacific Northwest and the observed west–southwest azimuth of motion south of about 36°N. To produce these features in the model velocity field would require, in the absence of additional deformation sources, west-directed basal drag in the Pacific Northwest and southwest-directed basal drag in the southwest US. (An alternative remedy for southern California is that steady slip at depth follows the trend of the Big Bend rather than the local Pacific–American relative motion vector (Lisowski 1991).) Mantle flow fields consistent with such basal drag have been proposed for these respective regions by Williams & McCaffrey (2001) and Liu & Bird (2002). This raises the question as to what extent imposed forces must be balanced in this type of model, how the compensating forces are distributed in the sublithosphere, and what mechanisms maintain them. To properly answer these questions will require consideration of not only the surface velocity field but also the sublithospheric flow field induced by tectonic interactions and mantle convection.

Figure 17

Predicted GPS velocity field based on the combination model of Section 5.7, shifted to produce a negligible velocity field around Great Salt Lake. The effect of background fault creep (e.g. Supplementary Fig. A2) has been added, rendering the resulting velocity field directly comparable with the observed velocity field of Supplementary Fig. A1.

Open in new tabDownload slide

In the Pacific Northwest north of the Mendocino triple junction, the tectonic forcing pattern implies a substantial component N–S contractile strain (Fig. 8c). This is in good agreement with principal stress directions in east-central Washington and northern Oregon (Smith 1978; Zoback & Zoback 1989; Wang 2000), E–W-trending thrust faults and folds in Washington (Wells & Simpson 2001) and the expected long-term northward displacement of the Oregon-Washington forearc ‘block’ with respect with North America. Closer to the subduction front, the relaxation from the repeating 1700 events generate E–W compression, perpendicular to the Cascadia coastline (Fig. 8a). The superposition of these two processes shapes the net strain-rate field (Fig. 8d), which exhibits a N70E compression near coastal Washington and Oregon (Fig. 3). Moreover, the tectonic forcing rates obtained along the strike of the Cascadia subduction zone show variable coupling of the JdF and NA plates, with apparently no coupling in central Oregon. (The forcing rate is , which is statistically indistinguishable from zero.) This is qualitatively consistent with reduced E–W contractile strain perpendicular to the coast from 43.5°N to 46.0°N that is seen in the observed strain pattern (Fig. 3) as well as relatively small GPS velocity vectors between these latitudes after correction for the Cascadia forearc rotation (Fig. 3c of Wang 2003). However, it conflicts with independent inferences of strain accumulation along the entire Cascadia subduction zone (McCaffrey 2000; Wang 2000; Svarc 2002b; Wang 2003). The chosen plane samples portions of the coast with both relatively high E–W contractile strain (south of about 43.5°N) and relatively small E–W strain north of 43.5°N, each with different magnitudes in the landward strain gradient. This suggests that a more detailed analysis is warranted for this region, for example, allowing for additional along-strike and downdip variations in forcing rates. The assumption of uniform slip of the 1700 earthquake also affects inferred forcing rates. Too much coast-perpendicular convergence may be generated with post-1700 relaxation that is artificially large around Oregon. Reduction of prescribed 1700-event slip in this area would be expected to translate into larger and likely positive comparable in magnitude with the other . Nevertheless, low tectonic forcing rates in central Oregon are correlated both with the presence of the thickest part of the Siletz terrane (Trehu 1994) and the presence of the lowest forearc seismicity (Trehu 1994). Our results are qualitatively consistent with the recent results of Verdonck (2004), who suggests a strong coupling north and south of central Oregon and a low coupling in central Oregon on the basis of observed vertical deformation rates across the boundary zone.

In California and Nevada, tectonic forcing and post-seismic relaxation following repeating 1906 and 1857 events constructively interfere around the SAF system but destructively interfere further east in the plate boundary zone (Figs 8a and c). This helps explain the large strain accumulation localized around the SAF fault and much smaller and variable strain accumulation further inland. Contributions from the Class 2 and 3 faults are necessary to localize the deformation in secondary deformation zones such as the ECSZ. Relaxation effects from the modelled events in western Nevada (i.e. CNSZ) are small but tangible (Fig. 8b). Estimated slip amplitudes of the 1954 Fairview Peak and 1932 Cedar Mountain earthquakes are similar to a priori values and statistically well above zero (i.e. inverted s₈ and s₉ in Fig. 9). This agrees qualitatively with Hetland & Hager (2003), who isolated the post-seismic relaxation signal of the central Nevada earthquakes in the GPS velocity field directly. Nevertheless, our methodology may be limited in its ability to detect a strong relaxation signal in the CNSZ because the discretization of the strain-rate field with ∼30 × 30 km² cells may smooth out a very localized post-seismic relaxation signal from the CNSZ earthquakes.

The component of steady slip in the lower elastic plate reinforces the constructive interference of tectonic forcing and post-seismic relaxation along the SAF system. We believe that a revised approach to modelling the steady-deep-slip component would be more effective in localizing predicted strain gradients in areas where they are observed. The most obvious improvement would come from constraining hypothetical deep slip to occur only along the deeper extensions of discrete faults. The excellent fit of the interseismic velocity field to those obtained by block models (e.g. Meade & Hager 2005; McCaffrey 2005; d'Alessio 2005) supports the idea that much of the steady slip that occurs in the lower elastic plate should occur on or near the deeper extension of locked portions of faults. In the context of our model, however, such additional steady deep slip is required only to the extent needed to produce local velocity gradients that cannot be explained by other means (e.g. viscoelastic-cycle effects).

7 Conclusions

The application of a viscoelastic cycle model to a comprehensive GPS data set in western US shows that the instantaneous velocity gradient field is well explained through the physical behaviour of an elastic-viscoelastic coupled system. The system is driven by the forces imparted on the NA lithosphere by the oceanic JdF and Pacific plates, viscoelastic relaxation from well-constrained past fault ruptures, cycle-averaged viscoelastic relaxation from less-constrained faults (including unrecognized faults), and steady creep/deep slip.

Post-earthquake relaxation (combined with steady deep slip) and tectonic forcing constructively interfere near the western margin of the plate boundary zone, producing, locally large strain accumulation along the SAF system. However, they destructively interfere further into the plate interior, resulting in smaller and more variable strain accumulation patterns in the eastern part of the plate boundary zone. The best rheological model that applies to western US as a whole (Model 4) exhibits a higher mantle viscosity than in smaller regions in this area over shorter timescale. This may hint at the influence of multiple material relaxation times in the complete mantle rheology, for example, a transient rheology (e.g. Yuen & Peltier 1982; Pollitz 2003a), as well as lateral variations in rheology as suggested by numerous indicators of lithospheric mechanical properties (Lowry 2000).

A combination of several mechanical processes generates a strain-rate field that agrees in pattern with the observed strain-rate field in the western US: (1) Pacific to North America and JdF to North America tectonic forcing, (2) viscoelastic relaxation cycles over active faults and (3) steady deep slip. When restricted to processes (1) and (2) using 18 identified fault zones, the predicted amplitude around the SAF system is too small, while the predicted amplitude around the Basin and Range province and elsewhere east of about 119°W is too large. Consideration of the budget of moment accumulation and release along the ∼1000 km-long plate boundary system also shows that past events on recognized fault zones represent only a fraction of the needed moment release in the system. Matching of observed strain rates is improved by introducing additional sources of moment release distributed over the broad areas between the eighteen identified source faults. This includes viscoelastic cycle effects on otherwise unaccounted faults (likely repeating M∼ 7 sources) treated in a cycle-averaged sense and steady slip in the lower lithosphere (nominally from 15 km to the base of the elastic plate at 20 km depth). We view the requirement of some degree of steady deep slip in addition to viscoelastic cycle effects as the most important conclusion of this study.

We have assumed throughout that viscoelastic cycle effects dominate the non-steady state component of strain accumulation in the western US, with slip in repeating events generally constrained by palaeoseismology. While this may be a reasonable approach for well-understood large fault systems (e.g. northern and southern SAF), it is unclear how applicable it is to fault zones in the Basin and Range province (Wallace 1987). Dixon (2003) find congruity between viscoelastic cycle effects and geologic slip rates for many faults in the ECSZ. However, Chang & Smith (2002), Hetland & Hager (2003) and Wernicke (2004) point out that non-cyclic viscoelastic effects and/or steady deep slip at rates exceeding recent geologic slip rates appear necessary in other localities.

Numerous issues remain to be clarified in future studies. What are the relative importance of viscoelastic cycle effects and steady deep slip in the context of a thin plate model? To what extent do relaxation and steady slip related to SAF-perpendicular shortening shape the instantaneous strain-rate field? Is a laterally heterogeneous viscoelastic structure and/or a more complicated (i.e. non-Maxwellian) rheology required to explain the instantaneous crustal strain-rate field? Are deformation measurements spatially distributed enough to permit determination of spatially variable forcing rates along the western margin of North America? Do deeper creep rates vary significantly with time, presumably in step with the seismic cycles of the major fault zones? Is basal shear a key component of a dynamic model of instantaneous crustal deformation? The framework adopted here is flexible to permit further investigation of these issues in more detailed studies.

Acknowledgments

This research has made use of Smithsonian Astrophysical Observatory WUSC velocity solution version 002 obtained from. The authors are grateful to Jerry Svarc for generating the USGS campaign GPS results presented here. MV is particularly grateful to the French Ministry of Foreign Affairs, whose support through the Lavoisier grant allowed her to carry out one year of postdoctoral research in the Earthquake Hazard Team, USGS, Menlo Park. This paper was improved by constructive criticisms from Brendan Meade, Bridget Smith and the associate editor. The authors thank Ruth Harris, Jim Savage, and Wayne Thatcher for their comments on a preliminary draft.

References

Author notes