A physical model for strain accumulation in the San Francisco Bay Region

Summary

1 Introduction

The San Francisco Bay region (SFBR) is an active zone of strain accumulation accommodating approximately 37–40 mm yr⁻¹ relative plate motion (Savage et al. 1998; Argus & Gordon 2001; Murray & Segall 2001; Prescott et al. 2001). Located within the Pacific—Sierra Nevada/Great Valley (SNGV) plate boundary zone, it is traversed by several major fault zones (Fig. 1) accommodating long-term slip rates ranging from a few mm yr⁻¹ for faults in the East Bay to as much as 25 mm yr⁻¹ along the San Andreas fault (SAF; Working Group on California Earthquake Probabilities 1999). Although historical observations of seismicity do not span even one SAF recurrence time, seismicity patterns indicate that the rate of moment release along the major faults is roughly in accord with the long-term slip rates (Bakun 1999). Strain accumulation as measured by geodetic measurements is dominated by [plate-boundary(PB)-parallel] right-lateral shear on roughly N33°W trending strike-slip faults and an integrated PB-perpendicular relative motion approximately 2–3 mm yr⁻¹ accommodated within the Pacific–SNGV plate boundary zone. Outstanding questions concerning the active deformation of this region are as follows.

Figure 1

Map of San Francisco Bay Region (SFBR) indicating major faults.

Open in new tabDownload slide

• (i)

  Is strain accumulation in the PB-parallel and PB-perpendicular senses uniform throughout the area or laterally variable?

• (ii)

  Is strain accumulation temporally uniform or are transient strain processes contributing to present or past deformation?

• (iii)

  What are the physical mechanisms that control strain accumulation?

We shall address these questions by constructing a physical model of strain accumulation for the SFBR that is consistent with the present-day velocity field and allows a simple interpretation in terms of the geometry of the driving Pacific–SNGV relative tectonic motion. Unlike earlier interpretations of SFBR geodetic data (Savage et al. 1998; Murray & Segall 2001), the proposed model does not depend explicitly on either the long-term slip rates or locking depths associated with the faults. The new approach is motivated by a few observations. First, the strain accumulation pattern may be explained to first order by simple shear across an ∼135-km-wide shear zone (Savage et al. 1998). Secondly, present strain rate within 20–30 km of the SAF accounts for approximately 60 per cent of the net PB strain (Savage et al. 1998; Murray & Segall 2001; Prescott et al. 2001) and strain rates are similarly elevated near the southern Santa Cruz mountains (Segall et al. 2000). Previous authors have documented that post-seismic relaxation following the M7.8 1906 San Francisco earthquake and M6.9 Loma Prieta earthquakes may have led to localized elevated strain accumulation in the respective source regions for years or decades following the events (Thatcher 1983; Pollitz et al. 1998; Kenner & Segall 1999; Parsons 2002). If transient strain from each of these events persists today, then it may provide an explanation for the aforementioned observations of elevated strain rates.

In the following sections, we examine the deformation pattern in greater detail by deriving the regional strain rate field, describe the elements of our strain accumulation model and discuss how the overall deformation pattern is shaped by the various contributing physical processes.

2 Regional Deformation

2.1 Velocity field

The GPS velocity field for the period 1994 to 2000 is shown in Fig. 2. It contains velocity vectors from 102 sites measured by the US Geological Survey (USGS), Bay Area Regional Deformation Network (BARD), Continuously Operating Reference Stations (CORS) and the International GPS Service (IGS). Details of the data processing may be found in Prescott et al. (2001) and Savage et al. (2004). The velocity field is an ensemble of six velocity profiles, which traverse the SFBR. Plots of velocity along each individual profile are presented in Fig. 3. In these plots PB-parallel and PB-perpendicular velocities are resolved along directions N33.85°W and N56.15°E, respectively for profiles 1, 2 and 3, and N40°W and N50°E, respectively for profiles 4, 5 and 6. These are averages of the local direction of Pacific—SNGV motion at the respective latitudes. The observed velocity field exhibits predominantly simple shear within an ∼135-km-wide plate boundary zone, offsets across creeping segments of the SAF system and, where it exists, distributed PB-perpendicular contraction. The effects of fault creep are most evident along profiles 5 and 6. Profiles 5 and 6 each show large offsets along the Calaveras fault and Profile 6 shows a large offset across the SAF as well.

Figure 2

GPS velocity field with respect to fixed North America (NA) from 1994 to 2000 with 95 per cent confidence regions (Savage et al. 1998). Boundaries between the 135-km-thick plate boundary zone and the Pacific and SNGV plates indicated by yellow lines, each of which is a small circle about a common Euler pole.

Open in new tabDownload slide

Figure 3

(a) Each of the six SW—NE profiles that compose the GPS velocity field are depicted with black line segments and labelled with boxed numerals. (b) Representation of the velocity field in terms of PB-parallel and PB-perpendicular components. These directions are specified by N33.85°W and N56.15°E, respectively for profiles 1–3, and N40°W and N50°E, respectively for profiles 4–6. San Andreas fault (SAF); Rodgers Creek and Hayward faults (RCF); Concord–Green Valley and Calaveras faults (GVF). X = 0 corresponds with the intersection of the particular transect with the San Andreas fault (SAF).

Open in new tabDownload slide

2.2 Velocity gradient field

The velocity field is useful for examining background tectonic motions, but suspected processes such as post-seismic relaxation and regional contraction are more subtle and smaller scale features that are better exhibited in the strain rate field. In order to extract information from the velocity field that is not directly evident in either the velocity vectors or the profiles, we construct an image of the regional strain and rotation rate fields. This is done by fitting the velocity field to a velocity gradient field subject to smoothing constraints in a damped least-squares inversion procedure (Spakman & Nyst 2002). In this method, one converts the velocity field into a spatially continuous velocity gradient field in coherent crustal blocks and spatially discontinuous motion on boundaries between the blocks. The continuous velocity gradient field may be interpreted in terms of tectonic loading and post-seismic relaxation. The discontinuous velocity gradient field is attributable to either coseismic offsets during the period of observation or fault creep.

We parametrize the study region by constructing a triangulation grid of uniform density. The size of the Delaunay triangles is chosen small enough to adequately interpolate between the GPS stations and represent possible small-scale variations in the strain rate field. We compute the four components of the horizontal velocity gradient tensor (∇v) at each node and adopt a linear dependence between the vertices of each triangle. To ensure that creep is not mapped into the continuous velocity gradient field, we parametrize the principal creeping sections (Figs 2 and 11e) as discontinuities in the triangular grid and assume that the velocity gradient field is piecewise continuous. We impose a priori the creep profile that is described in more detail in Section 4.

Figure 11

Components of the average 1994–2000 velocity field calculated from the best-fitting model. (a) Field v_ps, (b) v_shear, (c) field v_rot, (d) v_compr, (e) v_creep. Each velocity field is calculated at the 102 GPS sites as well as 20 points on the plate boundary. Each of the Pacific and SNGV plate boundaries are sampled at 10 uniformly spaced locations.

Open in new tabDownload slide

In addition to the GPS data, we impose an extra constraint that requires ∇×∇v=0 within each triangle, based on the zero-curl property of a continuous gradient vector field (Spakman & Nyst 2002). The inversion procedure selects a solution that fits the data in a least-squares sense and at the same time minimizes the model norm:

1

where m is the model parameter vector, d is the data vector of GPS velocities, A contains the components that link velocity to the velocity gradient and C_d is the data covariance matrix with 1-σ uncertainties of the horizontal GPS velocity components and their correlation coefficients. The trade-off between good data fit and minimal model norm can be regulated by tuning three parameters: α_d and α_s control the influence of amplitude damping (I) and smoothing (D), respectively, and γ_r determines the weight on the ∇×∇v= 0 constraints (hereafter extra data constraints). Because the extra data constraints are treated as data, we weigh them by tuning their error (∼1/γ_r) in the data covariance matrix C_d. A small error (equivalent to a large γ_r) increases the influence of the extra data constraints with respect to the GPS velocity data. In the inversion procedure, the full model covariance (C) and resolution (R) matrices are computed.

For the selection of the final solution, we consider fits of GPS velocity data and extra data constraints and the quality of model covariance and model resolution. For the model covariance, we compute the size of the unit model covariance matrix as a measure for the amount of error amplification mapped from data to solution (Menke 1989, pp. 67–68):

2

with M the number of model parameters. For the model resolution, we compute the size of the resolution spread function to provide some average measure for the independence of the model parameters (Michelini & McEvilly 1991; Eberhart-Phillips & Reyners 1997):

3

R_j stands for the jth component of the diagonal of R. R_kj represents all elements of the corresponding jth row of R, weighted by the distance D_jk between the nodes of the kth and jth model parameters. If S_j is relatively large, then model parameter j is poorly resolved. This may be the result of either a small diagonal resolution value R_j or a strong dependence on other model parameters k, amplified with distance, or a combination of both.

Figs 4(a) and (b) display the well-known trade-off between resolution and model variance for different combinations of α_s and α_d (for constant γ_r). Figs 4(c) and (d) represent the normalized root-mean-square misfit function of the data and the extra data constraints, respectively. The misfit function is defined as

4

N is the number of velocity measurements; d_i and p_i represent the ith velocity measurement and predicted velocity vector, respectively (in Fig. 4d), and the ith extra data constraint and its prediction, respectively (in Fig. 4c); σ_i is the standard deviation of d_i. Two additional trading-off relations exist between resolution (Fig. 4b) and data misfit (Fig. 4c) and between model variance (Fig. 4a) and data fit (Fig. 4c), both as functions of γ_r (for constant α_s and α_d). An increase of the influence of the extra data constraints for constant α_s and α_d reduces the size of the spread function and increases data misfit. The coupling of the standard deviation of the extra data constraints to the inverse of γ_r causes the misfit of the extra constraints to grow with increasing γ_r. The final solution used for further interpretation has size[S]= 1245, size[C]= 3.3 × 10⁻⁵ yr⁻¹, root-mean-square misfits of 1.5 for the GPS data and of 1.6 for the extra data constraints and is indicated by the white dot in Fig. 4.

Figure 4

Sensitivity analysis of inversions of GPS velocity data (Fig. 2) for the velocity gradient field. Results for size[C] (eq. 2, Figs 1a, 2a and 3a), size[S] (eq. 3, Figs 1b, 2b and 3b), the normalized root-mean-square misfit (eq. 4) of the ∇×∇v= 0 constraints (Figs 1c, 2c and 3c) and the GPS data equations (Figs 1d, 2d and 3d) for varying weights on damping (α_d, along y axes), smoothing (α_s, along x axes) and ∇×∇v= 0 constraints (γ_r, subplots 1, 2 and 3). Every contour plot represents 25 inversions.

Open in new tabDownload slide

Reasonably well resolved model parameters have a spread function value between 2 and 6 (Fig. 5c) with a diagonal resolution between 0.2 and 0.4 (or between 20 and 40 per cent, Figs 5a and b). For poorly resolved model parameters, the spread function value is greater than 6. The interior of the study area, despite some localized areas with zero diagonal resolution, is relatively well resolved. As may be expected, the areas with few or no stations are poorly resolved, i.e. north and east of San Pablo Bay and south and east of Monterey Bay (Fig. 5a). In general, the final solution has relatively small 1-σ model errors (Fig. 5d).

Figure 5

Quality in terms of resolution and model variance of the velocity gradient field computed from the GPS velocity data (Fig. 2): diagonal resolution of the model parameters of partial derivatives of the velocity field in (a) longitudinal and (b) latitudinal direction. Contours are plotted at 10 per cent intervals (every 0.1); (c) the resolution spread function with contours at intervals of 2; (d) model standard deviations (i.e. 1-σ errors) with contours at intervals of 1 × 10⁻⁸ yr⁻¹. The model error pattern for derivatives in the longitudinal direction is almost similar to the pattern for derivatives in the latitudinal direction.

Open in new tabDownload slide

The piecewise continuous velocity gradient field derived from the SFBR velocity field is shown in Figs 6(a), (b) and (c). The pattern of the second invariant of the strain rate tensor (i.e. pure shear strain rate, Fig. 6a) combined with rotation rate (Fig. 6c) confirms that the regional deformation is dominated by right-lateral simple shear strain rate. It further reveals that somewhat greater strain accumulation is localized near the SAF than around the East Bay faults (Hayward and Rodgers Creek faults; Greenville fault). The pattern of the first invariant of the strain rate tensor (i.e. areal strain rate, Fig. 6b) confirms an overall small regional contraction as noted by previous authors (Argus & Gordon 2001; Murray & Segall 2001; Prescott et al. 2001). Care should be taken while interpreting the areal strain rate pattern. For example, the absence of GPS sites along the northeastern boundary of the study area may make the detection of the velocity gradient in the east—west direction difficult in the area east and northeast of San Pablo Bay. The relatively good resolution of this part of the solution (Fig. 5a) can be explained by the unambiguous influence of the damping pushing the model parameters towards zero. The strain rate regime in this area shows almost uniaxial north—south contraction, causing a local minimum in the areal strain rates. However, if the real signal in this region is east—west extension, its detection and modelling would, at least partly, neutralize the negative area change found in our results. Similar data distribution problems may be responsible for (at least part of) the dilatational and contractional signals found at Point Reyes and southeast of Monterey Bay, respectively. It is noteworthy that a 50-km-long section of the Central Bay/San Pablo Bay shows no resolvable areal strain rate. This was previously noted by Argus & Gordon (2001), who further noted the low topography of this part of the SFBR and explained the relative lack of a contractional signal as a result of the local geometry of the faults accommodating the long-term strain release. Finally, localization in both shear and areal strain rate may be noted in the epicentral region of the 1989 Loma Prieta earthquake.

Figure 6

Components of the velocity gradient field derived from the GPS velocity field (A, B and C) and the best-fitting model velocity field (D, E and F). Plots A and D represent the second invariant of the strain rate tensor (i.e. pure shear) and indicate the axes of maximum rates of contraction and extension with black and white line segments, respectively, with magnitudes proportional to the lengths of these segments. Plots B and E represent the first invariant of the strain rate tensor (i.e. areal strain rate) with C and E indicative of rates of contraction and extension, respectively, plots C and F represent the rotation rate tensor (i.e. the antisymmetric part of the velocity gradient tensor) with CW and CCW indicative of clockwise and counter-clockwise. Superimposed dots indicate GPS sites. Contours are plotted at intervals of 1 × 10⁻⁷ yr⁻¹, with blue for positive, red for negative and green for zero values.

Open in new tabDownload slide

3 Viscoelastic Stratification

The basis for forward modelling of SFBR deformation here is the response to imposed dislocation sources of a gravitational viscoelastic coupled medium. The viscoelastic structure used here is shown in Fig. 7. It consists of an elastic upper crust of thickness 15 km, a Maxwell viscoelastic lower crust of thickness 15 km and viscosity η_c, and a Maxwell viscoelastic upper mantle of viscosity η_m. This structure is one of two alternative structures derived by Pollitz et al. (1998) on the basis of post-seismic geodetic observations following the 1989 Loma Prieta earthquake. Their model B is characterized by a relatively strong crust and weak mantle with η_c/η_m= 3.3 and η_m= 2 × 10¹⁹ Pa s. We adopt the crust-to-mantle viscosity ratio of that model because it is not well constrained by only horizontal post-seismic data, but we allow η_m to vary in order to attain a better estimation of viscosity structure using the present GPS velocity field.

Figure 7

1-D viscoelastic stratification of the SFBR assumed in this study, following model B of Pollitz et al. (1998).

Open in new tabDownload slide

Source parameters must be specified for the two earthquakes most likely associated with significant relaxation effects: the M7.8 1906 and M6.9 1989 Loma Prieta earthquakes. For this purpose, we adopt the coseismic fault slip and geometry models of Thatcher et al. (1997) and Marshall et al. (1991).

4 Evolution of Deformation

The regional deformation is envisaged to be composed of three principal physical processes: horizontal simple shear and PB-perpendicular uniaxial compression driven by Pacific—SNGV relative motion; post-seismic relaxation of the viscoelastic Earth following major earthquakes; and creep along parts of the San Andreas, Calaveras and Hayward faults. These elements provide a description of the processes of tectonic loading, time-dependent response resulting from earthquakes, and fault creep. In spherical r, θ, φ coordinates, the momentum equation for a 3-D distribution of viscoelastic properties in the Laplace transform domain is (eq. B1 of Pollitz 2003)

5

where ρ₀ and g₀ are density and gravitational acceleration on a spherically symmetric reference model, m_j(r′, s) represents a distribution of dislocation sources (earthquakes or fault creep) on the jth fault surface with area element dS_j, f represents a distribution of forces associated with tectonic loading, T is the stress tensor:

6

7

where u(r, s) is the displacement field. Eq. (5) accounts for the first-order coupling of viscoelastic deformation with the gravitational acceleration of the Earth; the second-order effect of coupling with changes in gravitational potential is neglected. This approximation is sometimes referred to in the seismological literature as Cowling's approximation (Dahlen & Tromp 1998). We assume relaxation of a Maxwell viscoelastic solid:

8

where κ₀, μ₀ are the static elastic constants and η is the viscosity. All quantities dependent on s are evaluated in the Laplace transform domain. Eqs (5)–(8) are to be solved subject to boundary conditions at the surface of the Earth (vanishing traction at the surface of the Earth) and an appropriate interior boundary condition (for example, vanishing displacement at an arbitrarily specified interior interface).

In principle, eqs (5)–(8) should be solved with a 3-D distribution of isotropic elastic parameters κ and μ. However, the main variations of κ and μ are with depth and we assume that the remaining lateral variations are controlled by the presence of thick lithospheric blocks bounding the plate boundary zone, i.e. the Pacific Plate and SNGV lithosphere. [This produces a large contrast in u(r, s) at all depths, including those beneath the elastic plate thickness assigned to the plate boundary zone.] Our approach to the solution is:

• (i)

  define a reference 1-D model with depth-dependent material properties corresponding to those of the shear zone;

• (ii)

  characterize the thick Pacific and SNGV plates as 3-D perturbations with respect to this reference 1-D model; and

• (iii)

  solve eqs (5)–(8) using perturbation theory.

As described in appendix B of Pollitz (2003), the solution for u and T with a 3-D distribution of viscoelastic properties is then the solution to the following equation with a 1-D distribution of viscoelastic properties:

9

where f_eq represents a distribution of equivalent forces that would be applied within the volume of material perturbations (i.e. in the volume outside of the plate boundary zone). For points r located on or within the plate boundary zone, the solution of u and T in eq. (5) on the 3-D model is very nearly the same as the corresponding solution of eq. (9) on the 1-D model.

It is not necessary to know f or f_eq explicitly; for brevity we refer to the sum f+f_eq as simply f. We construct a set of special solutions on the equivalent 1-D model involving specific dislocation sources m_j and physically plausible forces f. These special solutions and the associated velocity fields v=∂u/∂t are described in the time domain as follows.

• (A)

  Velocity field v_ps(r, t). Viscoelastic relaxation following specified earthquakes, m_j≠ 0, f= 0 globally.

• (B)

  Velocity field v_shear(r, t). Horizontal simple shear along vertical planes locally tangent to the curvilinear plate boundary (Fig. 1), m_j= 0, f= 0 within the plate boundary zone, f≠ 0 outside the plate boundary zone.

• (C)

  Velocity field v_rot(r, t). Rigid rotation, m_j= 0, f= 0 globally. This rotation is described by an Euler vector Ω.

• (D)

  Velocity field v_compr(r). Uniaxial horizontal contraction directed perpendicular to vertical planes locally tangent to the curvilinear plate boundary, m_j= 0, f= 0 within the plate boundary zone, f≠ 0 outside the plate boundary zone.

• (E)

  Velocity field v_creep(r). The velocity field associated with steady fault creep, m_j≠ 0, f= 0 globally.

Two points should be noted. First, solutions (B) and (D) are associated with non-trivial f outside the plate boundary zone. We postulate that f is distributed in such a manner that the complete solution of eq. (9) satisfies, to a high degree of accuracy, a kinematic boundary condition of constant Pacific—SNGV relative velocity applied to the edges of the shear zone (Fig. 2). This is critical as it allows us to construct the necessary special solutions using GPS constraints from within the plate boundary zone (where f= 0) and Pacific–SNGV boundary conditions, which are essentially known. Secondly, while both v_shear(r, t) and v_rot(r, t) may have arbitrary time dependence, v_compr(r) is assumed to be independent of time. Specifically, in a local (east, north) Cartesian coordinate system, these two velocity fields take the form

10

11

where δ is the distance measured positive of a point on or within the plate boundary from the SNGV plate boundary, W= 135 km is the width of the plate boundary and φ is the strike of the vertical plane locally parallel to the plate boundary. This strike depends upon the geometry of the plate boundary zone (yellow lines in Fig. 2). The boundaries are meant to approximate the physical boundaries of the plate boundary zone, i.e. the eastern edge of the Coast Ranges on the SNGV side and the offshore faults on the Pacific side. The given boundaries were determined by a process of trial and error. They are specified as small circles about a pole located at angular distances 17.82° and 19.04° from the pole (Fig. 8). The local azimuth of a small circle about this pole through a given point defines φ. We can then write

12

for points located in the plate boundary zone. Similarly, we define a pole to be 90° from the plate boundary zone along an azimuth tangent to it (Fig. 8). We then have

13

Figure 8

The SNGV–Pacific relative plate motion is decomposed into a plate-boundary(PB)-parallel component parallel to small circles about (48.0°N, −100.0°E) and a PB-perpendicular component parallel to small circles about (−41.1°N, −74.1°E). The SNGV–Pacific plate boundary zone is delineated by two small circles about at angular distances of 17.82° and 19.04°. After estimation of the regional PB-perpendicular velocity v₂, the SNGV–Pacific angular velocity vector Ω_SNGV–Pac is obtained from eq. (16): (44.64°N, −99.34°E, 1.081° Myr⁻¹) (triangle). Plotted ellipse indicates the SNGV–Pacific angular velocity vector and 95 per cent confidence ellipse obtained by Argus & Gordon (2001). The upper plot shows a closer view of the boxed region in the lower plot.

Open in new tabDownload slide

Note that the solution for v_shear in eq. (12) is a valid solution of eq. (9) even for time-dependent v₁(t), whereas the corresponding solution for v_compr in eq. (13) is valid only for constant v₂. Time-dependent v₂ would produce transient shear strains, which relax with time and would modify the solution in eq. (13). In order to account for the additional relaxation would require explicit dependence upon the history of v₂. In the subsequent process of matching boundary conditions, it is convenient to have a direct relationship between v_shear or v_comp on the plate boundaries and the corresponding velocity fields within the plate boundary. The choice of time-independent v₂ is thus somewhat limiting but allows this approach to be carried out straightforwardly and greatly simplifies the analysis.

The principal creeping fault segments associated with solution (E) are shown in Figs 2 and 11(e). We prescribe a priori the behaviour on these creeping faults by specifying the depth range and rate of slip as follows: Hayward fault, 0–5 km, 5 mm yr⁻¹ based on Simpson et al. (2001); central SAF, 0–15 km, variable slip rate 12–30 km yr⁻¹ (Rymer et al. 1984); NW creeping segment, 0–15 km, 12 mm yr⁻¹; S Calaveras fault, 0–15 km, 12 mm yr⁻¹ (Oppenheimer et al. 1990). The velocity field produced by steady creep prescribed by these dislocations is evaluated in the fluid limit of the viscoelastic model in a spherical geometry using the method of Pollitz (1996).

The model velocity field at point r in a fixed SNGV reference frame may be written in the time domain as follows:

14

In order to compare the model velocity field with the GPS velocity field, it is necessary to apply an additional SNGV—North America (NA) rotation:

15

5 Model of SFBR Deformation

5.1 Estimation of model parameters

The model velocity field with respect to fixed SNGV in eq. (14) depends upon v₁, v₂, η_m and Ω (three components). The model velocity field with respect to fixed NA further depends upon Ω_SNGV–NA. We require that v(r, t|SNGV) be consistent with constant Pacific–SNGV boundary conditions. To make this more precise, it is convenient to define the SNGV–Pacific angular velocity vector as the composite of PB-parallel and PB-perpendicular motions:

16

where v_pl is the magnitude of Pacific—SNGV relative motion parallel to the plate boundary. We require that

17

for points r_i located on the Pacific Plate boundary, and

18

for points r_i located on the SNGV plate boundary.

The strategy for determining the various parameters is as follows.

• (i)

  Estimate v₁ and rotation Ω for a given η_m by least-squares inversion such that v(r_i, t|SNGV) satisfies the above boundary conditions with optimally small error. Note that v₂ is indeterminate in this step because of its common appearance in eqs (14) and (17).

• (ii)
  Determine additionally v₂ and the rotation Ω_SNGV–NA, which minimizes the reduced χ²

  

  19

  where v is a composite velocity vector consisting of all model or observed velocity components at the N= 102 employed GPS sites and C is the data covariance matrix. M= 6 is the number of independent parameters.
• (iii)

  Repeat the above two steps with a grid search over η_m in order to find the optimum simultaneous fit of both the boundary conditions and the GPS data.

For purposes of illustrating sensitivity to the model parameters, one may remove one of the parameters (say, v₂) from the inversion process in step (i) and include it together with η_m as a grid search parameter in step (iii).

5.2 Results

The above procedure has been carried out using a range of possible relative plate velocities v_pl. Results are very similar for 37 mm yr⁻¹ < v_pl < 40 mm yr⁻¹. We choose the value v_pl= 38 mm yr⁻¹ for further consideration.

The best-fitting model is specified by η_m= 1.2 × 10¹⁹ Pa s, v₁= 32 mm yr⁻¹, v₂= 3 mm yr⁻¹, Ω= (39.326°N, −118.509°E, 0.594° Myr⁻¹) (clockwise rotation) and Ω_SNGV–NA= (48.575°N, −115.128°E, 0.490° Myr⁻¹) (Fig. 9, clockwise rotation). From eq. (16) this yields Ω_SNGV-Pac= (44.64°N, −99.34°E, 1.081° Myr⁻¹) (Fig. 8). The sensitivity of model fit to η_m and v₂ is shown in Fig. 10, which indicates a clear minimum in these parameters.

Figure 9

Estimated SNGV—North America (NA) angular velocity vector and 95 per cent confidence region obtained in this study. Alternative estimates are by Murray & Segall (2001; circle), Argus & Gordon (2001; star) and Dixon (2000; square).

Open in new tabDownload slide

Figure 10

Results of grid search for η_m and v₂ to minimize reduced χ² (eq. 19). The best-fitting model is obtained at η_m= 1.2 × 10¹⁹ Pa s and v₂= 2.9 mm yr⁻¹.

Open in new tabDownload slide

Inference of these parameters is dependent upon the choice of plate boundaries. The Pacific and SNGV plate boundaries in Fig. 2 are small circles about in Fig. 8. Because the v_shear component rotates if is changed, a change in will introduce a trade-off between inferred v_shear and v_compr. Consequently, the value of best-fitting v₂ is very sensitive to the initial choice of . However, the inferred SNGV–Pacific angular velocity vector Ω_SNGV–Pac and, consequently, predicted model deformation are practically insensitive to this choice. For example, if were chosen 1° farther north, inferred v₂ would increase to 4 mm yr⁻¹, but the compounded effect of these changes on Ω_SNGV–Pac when propagated through eq. (16) is negligible. Therefore, the pole location of Ω_SNGV–Pac shown in Fig. 8 is largely insensitive to initial assumptions. It is located just south of the 95 per cent confidence region for this motion obtained by Argus & Gordon (2001). Argus & Gordon (2001) derived plate motions using very different geodetic data, which includes sites on the southern SNGV plate. We attribute the difference between the two pole locations to the limited aperture of the GPS data used in the present study. Considering that the present data set is confined to the northern boundary of the SNGV plate, we note that the formal errors associated with our model Ω_SNGV–NA (Fig. 9) are too small because , the dominant contributor to Ω_SNGV–Pac, (Fig. 8) was chosen subjectively in order to represent the plate boundary. Our formal errors in Ω_SNGV–Pac are thus smaller than those of Argus & Gordon (2001) and Murray & Segall (2001). Consequently, the errors that have been propagated through to errors in Ω_SNGV–NA are too small. An additional factor may be that there are subtle differences between SNGV motion in the SFBR and its motion further south. The second possibility is also suggested by the fact that our Ω_SNGV–NA pole and that of Murray & Segall (2001) both lie to the northeast of the plate boundary zone (clockwise motion) whereas other estimates of Ω_SNGV–NA lie to the southwest of the plate boundary zone (counter-clockwise motion; Fig. 9). Ours and Murray and Segall's estimate are based on data from the SFBR and northern SNGV plate, whereas the other estimates are based on data more broadly distributed over the SNGV plate. In any case, it is clear that geodetic data confined to the SFBR are insufficient to place solid constraints on overall SNGV—Pacific motion. This is further seen by the fact that our model leads to a NA—Pacific Euler pole of (46.1°N, −104.0°E, 1.564° Myr⁻¹), which is substantially different from the well-constrained Argus & Gordon (2001) value of (50.1°N, −75.9°E, 0.778° Myr⁻¹). The local Pacific–NA motion in the SFBR predicted by the two Euler poles are, however, very similar. Predicted rates are within 2 mm yr⁻¹ and azimuths differ by 0° to 3°, the same order as the azimuth misfit between model and data vectors. Murray & Segall (2001) gained greater control on NA–Pacific motion, in good agreement with Argus & Gordon (2001), by including additional GPS data on the northern SNGV plate and Nevada, well to the east of our study area.

The five components of the model velocity field are shown in Fig. 11. It is clear that most (approximately 85 per cent) of the simple shear being accommodated in the plate boundary zone is represented by v_shear. The remaining 15 per cent is accommodated by v_ps primarily through viscoelastic relaxation following the 1906 earthquake. The total model velocity field with respect to fixed SNGV obtained by summing these five components is shown in Fig. 12. There is excellent agreement between the model velocity field and the imposed relative velocity boundary condition given by eq. (17), which is a combination of 38 mm yr⁻¹ PB-parallel and 3 mm yr⁻¹ PB-perpendicular motion. Note that no single component of the model velocity field exhibited in Fig. 11 is consistent with Pacific–SNGV kinematics. However, only the total model velocity field with respect to SNGV, given by their sum, should be compared. To the extent that the total model velocity field satisfies the imposed Pacific–SNGV velocity boundary conditions, the model is consistent with known Pacific–SNGV kinematics.

Figure 12

The total model velocity field v(r, 1994 –2000 | SNGV) (eq. 14) evaluated at the 102 GPS sites and 20 plate boundary sites are shown by the black vectors. The imposed relative Pacific–SNGV velocity boundary condition prescribed by eqs (17) and (18) are shown by the red vectors. In eq. (14) we take v_pl= 38 mm yr⁻¹ and v_compr given by eq. (11) with v₂= 2.9 mm yr⁻¹. (The boundary condition and model velocity vectors are plotted at 10 sample points on both the Pacific and SNGV boundaries. They are negligible on the SNGV boundary.)

Open in new tabDownload slide

6 Discussion

The total model velocity field with respect to fixed NA is shown in Fig. 13, where it may be compared with the GPS velocity field. The two velocity fields are in very good agreement and it is difficult to identify systematic biases, even near creeping faults, suggesting that the creep model described in Section 4 accounts for the main discontinuities observed in the regional velocity field. We compute the piecewise continuous velocity gradient field for the modelled velocity field by application of the same method under the same conditions for parametrization and regularization as described in Section 2.2. We use a uniform value for the a priori data standard deviation of 1.2 mm yr⁻¹, which is the average standard deviation of the GPS velocity data. The inversion characteristics of the final gradient field of the modelled velocity field are very similar to those of the GPS derived solution: size[S]= 1233, size[C]= 3.3 × 10⁻⁵ yr⁻¹, root-mean-square misfits of 1.0 for the GPS data and 1.5 for the extra data constraints. The velocity gradient fields associated with the modelled and observed velocity fields (Fig. 6) agree very well. The greater concentration of pure shear and rotation rates (approximately equal to simple shear strain rate) around the SAF is exhibited in both modelled and observed strain rate fields. The pattern of areal strain rates are also similar. The small differences that exist are mainly in amplitude rather than pattern. Thus, the large contractile strain rates observed in the NE part of transect 1 and near the Loma Prieta rupture zone, as well as the dilatational strain rates observed near Point Reyes, are also present in the modelled strain pattern but with reduced amplitude. Excluding the area near the Loma Prieta rupture zone, the observed dilatational strain rates (Fig. 6b) are uniform over most of the SFBR where resolution is relatively good (Fig. 5a). A similar uniformity in modelled strain rate pattern characterizes our solution (Fig. 6e). Except for the northeastern part of transect 1, where modelled dilatation rate does not match the observed contraction and where no local post-earthquake effects have been accounted for, the assumption of uniform regional contraction (eq. 13) appears to be a good approximation. If regional contraction is indeed manifested over several earthquake cycles within a localized band, this feature is not exhibited during the present interseismic period.

Figure 13

Total model velocity field v(r, 1994 –2000|NA) (black arrows), calculated from eq. (15), and GPS velocity field (red arrows).

Open in new tabDownload slide

The curvilinear boundaries shown in Fig. 2 are intended to represent the average trend of the Pacific—SNGV plate boundary at any given location. The trend defined by those curves varies from N31.7°W in the north (Lake Berryessa) to N36.8°W in the south (Calaveras–SAF junction). The plate boundary so defined closely follows the average trend of the major faults at any given latitude. For example, Argus & Gordon (2001) found an average strike-slip fault trend of N33.9°W across Lake Berryessa (profile I–I′ in their fig. 3), N33.1°W across the San Francisco Peninsula (profile G–G′) and N37.5°W near the Calaveras–SAF junction (profile E–E′). Comparison of such average fault trends with the local relative Pacific–SNGV motion vector determined by those authors leads to PB-perpendicular convergence rates varying from 1.8 to 3.3 mm yr⁻¹ (except in the San Pablo Bay zone, which is undergoing a small amount of PB-perpendicular extension), close to the value of 2.9 mm yr⁻¹ determined in our study. In an independent study, Savage et al. (2004) have determined an average regional PB-perpendicular contraction rate of 0.9 ± 1.1 mm yr⁻¹, i.e. a regional extension that is statistically indistinguishable from zero. This estimate is based on the fact that the best-fitting uniform horizontal strain field for the region is characterized by principal strain rates of ε₁₁= 164.7 ± 7.2 nanostrain yr⁻¹ and ε₂₂=−157.9 ± 6.9 oriented N74.0°W and N16.0°E, respectively. This is equivalent to a combination of pure shear of 161.4 ± 5.0 nanostrain yr⁻¹ oriented N28.4°W and 6.9 ± 7.9 nanostrain yr⁻¹ of uniaxial extension perpendicular to this direction. Savage et al. (2004) suggests that 0.9 ± 1.1 mm yr⁻¹ should be representative of the regional PB-perpendicular compression, which would imply regional extension at a statistically insignificant level. However, we note that the average fault trend is ∼N33.5°W, based on SFBR fault trends summarized in Table 2 of Argus & Gordon (2001). Use of that trend in the presence of the derived uniform strain rate field would result in ∼3 mm yr⁻¹ PB-perpendicular compression, similar to the other estimates.

The agreement between observed and modelled strain rate fields is attributed in part to the effect of post-seismic relaxation. The strain rate field associated with post-1906 and post-1989 relaxation is shown in Fig. 14. It suggests that the greater concentration of shear strain rate near the SAF is a result in large part of post-1906 relaxation, which persists more than 90 yr after the earthquake. Such behaviour is also predicted in fully numerical finite element models that include viscoelastic elements within parts of the lower crust and mantle (e.g. fig. 2b of Kenner & Segall 1999). Similarly, parts of the observed shear and areal strain rate maxima near the Loma Prieta rupture zone are attributed to post-1989 relaxation.

Figure 14

Strain rate pattern as a result of viscoelastic relaxation following the 1906 and 1989 earthquakes, derived from the corresponding velocity field of Fig. 11(e). (a) Pure shear strain rates; (b) areal strain rates; (c) rotation rates. Contours are plotted at 5 × 10⁻⁸ yr⁻¹ intervals, with blue for positive, red for negative and green for zero values.

Open in new tabDownload slide

7 Conclusions

Strain accumulation in the SFBR can be described to first order as the product of simple shear within an ∼135-km-wide plate boundary zone combined with minor PB-perpendicular compression, lingering effects of post-seismic relaxation following the 1906 San Francisco and 1989 Loma Prieta earthquakes, and the effects of steady fault creep. This physical model has been calibrated to optimally fit a set of 102 GPS velocity vectors by a grid search for both the degree of PB-perpendicular compression and the viscoelastic parameters that govern the behaviour of the plate boundary zone. With a mantle viscosity of 1.2 × 10¹⁹ Pa s and a regional PB-perpendicular compression of ∼3 mm yr⁻¹, there is excellent agreement between the observed and modelled velocity field. Slightly greater strain rates centred around the SAF proper (approximately 20 nanostrain greater than the average 140 nanostrain for the entire plate boundary) at present is attributed to post-seismic relaxation following the 1906 earthquake. The driving forces in this model are horizontal forces transmitted by the relatively thick Pacific and SNGV blocks to the plate boundary zone along its edges. This physical model is not unique but, if applicable, it suggests that previously considered loading mechanisms involving continual slip beneath a locking depth are not necessary in order to explain the main features of the strain accumulation pattern.

Acknowledgments

We are grateful to Weijun Gan and Will Prescott for making available to us the processed GPS data set used in this study. We thank Jim Savage and other members of the Crustal Strain Group at USGS for discussions. This paper was improved by constructive comments of two anonymous reviewers and the associate editor.

References