How do “ghost transients” from past earthquakes affect GPS slip rate estimates on southern California faults?

In this study, we investigate the extent to which viscoelastic velocity perturbations (or “ghost transients”) from individual fault segments can affect elastic block model‐based inferences of fault slip rates from GPS velocity fields. We focus on the southern California GPS velocity field, exploring the effects of known, large earthquakes for two end‐member rheological structures. Our approach is to compute, at each GPS site, the velocity perturbation relative to a cycle average for earthquake cycles on particular fault segments. We then correct the SCEC CMM4.0 velocity field for this perturbation and invert the corrected field for fault slip rates. We find that if asthenosphere viscosities are low (3 × 1018 Pa s), the current GPS velocity field is significantly perturbed by viscoelastic earthquake cycle effects associated with the San Andreas Fault segment that last ruptured in 1857 (Mw = 7.9). Correcting the GPS velocity field for this perturbation (or “ghost transient”) adds about 5 mm/a to the SAF slip rate along the Mojave and San Bernardino segments. The GPS velocity perturbations due to large earthquakes on the Garlock Fault (most recently, events in the early 1600s) and the White Wolf Fault (most recently, the Mw = 7.3 1952 Kern County earthquake) are smaller and do not influence block‐model inverted fault slip rates. This suggests that either the large discrepancy between geodetic and geologic slip rates for the Garlock Fault is not due to a ghost transient or that un‐modeled transients from recent Mojave earthquakes may influence the GPS velocity field.


Introduction
[2] Elastic block models are generally used to infer slip rates on fault segments in tectonically complex areas, such as southern California [e.g., McCaffrey, 2005: Meade and. These models implicitly assume steady state deformation. However, owing to viscoelastic coupling, deformation rates and patterns around major faults are expected to vary with time between large earthquakes [e.g., Savage and Prescott, 1978]. Where viscoelasticity has been incorporated into block models, differences in inferred slip rates have resulted [Huang et al., 2010;Johnson and Fukuda, 2010;Pollitz et al., 2010;Chuang and Johnson, 2011]. [3] Here we investigate the extent to which viscoelastic velocity perturbations (or "ghost transients") from individual fault segments can affect elastic block model-based inferences of fault slip rates from GPS velocity fields. By "ghost transients," we mean perturbations in the GPS velocity field relative to the cycle-average velocities (see also Hetland and Hager, 2005). Ghost transients are not the same as post-seismic deformation; they are nonzero (and may be significant) throughout the inter-seismic interval. Our approach is distinct from modeling approaches in which the post-seismic deformation from recent known earthquakes is estimated and subtracted from the velocity field prior to using elastic models to invert for fault slip rates [e.g., Hammond et al., 2011]. Since velocities and strain rates around a fault segment late in the inter-seismic interval may be significantly smaller than their cycle-average values [e.g., Thatcher, 1983, Figures 2-6;Cohen and Kramer, 1984, Figure 4], GPS velocity corrections may be necessary throughout the seismic cycle. [4] We focus on the southern California GPS velocity field (SCEC CMM4.0) [Shen et al., 2011], exploring the effects of several known large earthquakes for end-member rheological structures. The GPS position time series data used to construct the SCEC CMM4.0 velocity field were corrected for early post-seismic deformation following the Landers, Northridge, and Hector Mine earthquakes, using a logarithmic function with a characteristic decay time of 10 days [Shen et al., 2011]. SCEC CMM4.0 is based in part on GPS position time series data from 1986 to 2004. SCEC CMM4.0 velocities differ little from the Plate Boundary Observatory solution, which was derived from GPS data from 2004 to 2008 [Shen et al., 2011].
[5] For selected fault segments, we construct an idealized earthquake history consisting of a sequence of periodic, identical repeating events ending with the most recent large earthquake. We first calculate average velocities and time-dependent perturbations relative to this average at all GPS sites in the neighborhood of the fault segment. (We deal with perturbations because to recover the complete GPS velocity field, we would have to compute and sum cycle-average velocities and perturbations for all fault segments in the region, many of which are unknown or poorly constrained.) Next, we invert two GPS velocity fields for slip rates using a block modeling approach-one field that has been corrected for the perturbation and one which has not-and we compare the resulting slip rates. For this study, the viscoelastic models we use to calculate the perturbations are simple (layered models with linear [Maxwell or Burgers Body] rheologies), and the locking depth is fixed in the block models. Our approach differs from that of Chuang and Johnson [2011] in that we assess the effect of earthquake cycles on individual fault segments, rather than simultaneously modeling the entire ensemble of southern California faults. This enables us to isolate the effects of individual faults and well-constrained earthquakes, and identify segments which contribute to the differences in inferred slip rates between elastic and viscoelastic block models.

Modeling Methods
[6] To the extent possible, we have represented the same lithosphere in both the block models and the viscoelastic earthquake cycle models. We assume a uniform elastic Earth with l = m = 30 GPa. The block models are elastic so the layered viscosity structure is described solely for earthquake-cycle models. We use the code VISCO1D [Pollitz, 1997] to model viscoelastic earthquake-cycle deformation. VISCO1D uses the first two terms in equation (1) of Pollitz et al. [2010] to compute time-dependent and cycle-average velocities, respectively, at points on the Earth's surface. We are interested in the perturbation (i.e., the first term) only. Segall [2010] also presents a velocity solution for a viscoelastic earthquake-cycle model, making use of both a steady state, cycle-average term and a time-dependent, viscoelastic perturbation term (his equations 12.27 and 12.28, respectively).
[7] The models require that we specify a layered rheological structure, with a viscosity and two elastic modulii set for each layer. The middle to upper crust, extending from 0 to 25 km depth, is assumed to be elastic. Viscosities for our low-and high-viscosity Earth models (M1 and M2, respectively) are shown in Table 1 and Figure 1. Model M1 viscosities are taken from the ranges given by Johnson et al. [2007] which are inferred from inter-seismic deformation around the San Andreas Fault Zone (SAFZ). Model M2 viscosities are based on post-seismic deformation following the 1999 Izmit earthquake sequence on the North Anatolian Fault Zone (NAFZ) [Hearn et al., 2002a[Hearn et al., ., 2009, and on late inter-seismic deformation around the central NAFZ [Hearn et al., 2002b[Hearn et al., , 2009. Hearn et al. [2009] found that the late interseismic NAFZ deformation required steady state viscosities of 10 20 Pa s or greater to match the observed pre-1999 era GPS velocities across the fault, though Kokum and Johnson [2011] obtain a lower minimum viscosity (5 Â 10 19 Pa s) when modeling more realistic, finite-length earthquake ruptures. (Hereinafter, we will refer to Johnson et al. [2007] as JEA07, Chuang andJohnson [2011] as CJ11, andHearn et al. [2009] as HEA09.) [8] Though HEA09 assume a layered elastic structure in their models, we do not do so here. Given In the model of Chuang and Johnson [2011], the elastic crust extends from 0 to 20 km depth and the lower crustal layer from 20 to 30 km depth.

Geochemistry Geophysics
Geosystems G 3 G uncertainties in the viscosity estimate from HEA09, and the fact that we wanted an end-member lithosphere model with a long Maxwell time (/m), we deemed the elastic layer over a uniform half-space (UHS) approximation acceptable for our purpose here. When we began this project, JEA07 and HEA09 represented our best guesses at low-and high-viscosity end-members for the region. Since then, CJ11 have inferred a viscosity for the uppermost mantle (between depths of 30 and 50 km) which is lower than that of M1, though below 50 km, their viscosity value is 2 times the model M1 value (Table 1).
[9] This is not meant to be a comprehensive study of all fault segments in southern California capable of producing large earthquakes and perturbing the GPS velocity field. We have chosen to model earthquakecycle deformation for the fault segments that produced the 1857 SAF earthquake, the~1640 Garlock Fault earthquake, and the 1952 Kern County earthquake. These segments were selected because they produce large earthquakes, their earthquake parameters are well defined, and (for the SAF and the Garlock Fault) they are where discrepancies between GPS and Holocene geologic slip rates have been documented (see summary in CJ11).
[10] Earthquake rupture properties for each of the three segments we model are shown on Table 2. For each segment, we assume slip from 0 to 25 km depth; this is meant to include both shallow coseismic slip and rapid afterslip at greater depths in the upper to middle crust. The geometries and lengths for the 1857, 1952, and 1640 ruptures ( Figure 2 and  The Kern County earthquake was modeled as oblique, with 3.9 m of strike-slip and 1.9 m of dip-slip offset, following Bawden [2001]. L is length of the rupture, W is rupture depth (interval is from 0 to W km), D is elastic plate thickness, U is coseismic slip, T c is inter-seismic interval, and t is time since the last large earthquake.

Viscoelastic Perturbation Calculation
[11] The nature of viscoelastic post-seismic surface deformation is to generate a surface velocity field that is more rapid than the average cycle velocities early in the earthquake cycle and slower than average later in the cycle [Savage and Prescott, 1978]. This means that corrections of observed GPS velocities at any epoch for the effects of viscoelastic earthquake-cycle deformation must be referenced to the average inter-seismic cycle, that is, the deformation averaged over the entire earthquake cycle. This is shown schematically in Figure 3 for an idealized M w = 8 earthquake (not unlike the 1857 event) with a fault length of 300 km, average slip of 5 m, and faulting depth of 25 km. The assumed rheology is that of Model 1 (M1), with the lower upper mantle velocities that show larger earthquake cycle effects. Earlier in the cycle (up to~100 years after the major cycle event), transient velocities are higher than the cycle average; later in the cycle, the velocities are lower. The perturbation relative to the cycle average (ghost transient) is positive, that is, of the same sign as the cycle-average velocity, early in the inter-seismic interval. The ghost transient is negative (   block modeling approach where GPS velocities are fit by rotation of spherical caps and specified locking at the cap boundaries that accounts for elastic strain accumulation on cap-bounding faults [McCaffrey, 2005;Meade and Hager, 2005]. We assume a locking depth of 15 km. Our test models based on the uncorrected SCEC CMM4.0 velocity field confirm that slip rate estimates are insensitive to modest changes in the assumed locking depth (i.e., from 15 to 25 km).

Results
[13] Depending on the rheological model, the ghost transient associated with the 1857 SAF rupture segment may or may not be large enough to significantly affect block model-inferred slip rates along the SAF. Inferred slip rates on other faults in our model are not significantly altered regardless of lithosphere rheology. For rheology model M1, inversion of the corrected GPS velocity field yields SAF slip rates that are about 5 mm/a greater than slip rates estimated from the uncorrected field (Table 3 and Figure 4). Velocity perturbations exceed 1 mm/a over a region hundreds of kilometers in dimension, centered on the 1857 rupture ( Figure 4). When the M2 rheological structure is assumed, the ghost transient contribution to the GPS velocity field is similar in appearance but the perturbation amplitudes are small. In this case, the ghost transient increases the inferred SAF slip rates by 1.2-2.5 mm/a (Table 3). We find that for both the M1 and M2 rheologies, ghost transients associated with the White Wolf Fault (i.e., the 1952 Kern County earthquake segment) and the Garlock Fault do not affect inferred slip rates significantly. That is, inversions of the corrected and uncorrected GPS velocity fields in these cases yield similar slip rates (to within less than 0.5 mm/a) on all of the southern California faults represented in our model. [14] The fit of the models to the GPS data is similar for both inversions using SCEC CMM4.0 velocities and for those velocities corrected for transient effects. However, for the Garlock Fault inversions using GPS data corrected for viscoelastic perturbations we find that larger-magnitude transients (i.e., from model M1) increase the misfit residuals without significantly changing the predicted slip rate on the Garlock Fault. For example, if the predicted perturbation for the M1 rheology model is increased by a factor of 4, the predicted Garlock fault slip rate is unchanged while the misfit increases by 41%.
[15] It should be noted that our viscoelastic correction applied to the SCEC CMM4.0 velocities depends on our ab initio assumptions of geologically estimated slip rates on the San Andreas (34 mm/a), Garlock (7 mm/a), and White Wolf (2 mm/a) faults as described above. If these slip rates are smaller than the true (unknown) fault slip rates, our viscoelastic correction will be correspondingly underestimated, and the resulting effects on the SCEC CMM4.0 velocity field will be too small. However, GPS slip rates for the San Andreas and Garlock faults are lower than geologic rates (e.g., CJ11). Therefore, our calculated viscoelastic perturbations (based on the geologic rates) are likely to be upper bounds that maximize the ghost transient effects on the slip rates. GPS and geologic rates for the White Wolf fault agree within uncertainties, and predicted viscoelastic corrections are small, so we are less concerned with the geologic slip rate assumed in our calculations for the effects of the 1952 earthquake on White Wolf fault slip rate.

Slip Rates and Ghost Transients
[16] If the JEA07 rheology is applicable to southern California, then ghost transients from the 1857 SAF rupture segment can perturb inferred slip rates by about 5 mm/a and may partly explain differences between geologic and GPS-based slip rates on the Carrizo and Mojave sections of the SAF. However, if a stiffer rheology (such as that inferred for the lithosphere surrounding the NAFZ in Turkey by HEA09) applies, the ghost transient correction is small. In this case, simple elastic block models should be suitable for inferring southern California fault slip rates from the SCEC CMM4.0 GPS velocity field, and an alternative explanation must be sought for discrepancies between geodetic and geologic fault slip rate estimates. We note that if the CJ11 viscosity structure applies, then the ghost transient velocities are larger than shown on Figure 4. When the SCEC CMM4.0 velocity field is corrected for this transient, the inferred SAF rate increases by 5 to 7 mm/a along the segments shown on Figure 2. [17] Ghost transients from the other two fault segments we investigated (the Garlock and White Wolf faults) had no effect on inferred southern California fault slip rates; the velocity perturbations were too small. We note that for a hypothetical model with smaller blocks (bounded by low-slip rate faults), this conclusion could change. For short rupture segments, the length scale of velocity contributions becomes short ( Figure 5). A similar scaling applies to both the cycle-average velocity contributions and the perturbations. Small-wavelength velocity field perturbations tend to be ignored by models with large blocks, merely adding to the model RMS error but not changing the preferred slip rate solution. If the blocks were small (i.e., of the order of the wavelength of the velocity field perturbations), then the best solution might require altered slip rates.

Comparison with Other Studies
[18] CJ11 explain the Garlock GPS-geologic slip rate discrepancy in terms of viscoelastic earthquake cycle effects. We and CJ11 assume similar earthquake recurrence intervals, rupture geometry, and earthquake size. Our block models are nearly identical. CJ11 assume a 30 km thick lithosphere-that is, their lowest viscosity layer extends downward from 30 km (Table 2). We assume a 50 km thick lithosphere, consistent with JEA07. This tends to reduce the contribution of viscoelastic relaxation to the GPS velocity field, with the greatest differences being in the near field. We re-ran the Garlock viscoelastic earthquake cycle models to calculate the ghost transient using the CJ11 rheology but obtained a small-magnitude velocity perturbation. When we inverted the velocity field corrected for this ghost transient, the Garlock fault slip rate was unaffected. Hence, we cannot duplicate the findings of CJ11. We note that CJ11 use a modeling approach that differs from ours in that they incorporate earthquake cycles on all southern California faults at one time, rather than looking at earthquake cycle effects for individual segments. Our inferred slip rates could be similar to those of CJ11 if we used the same linear viscosity structure and corrected the SCEC CMM4.0 GPS velocity field for earthquake cycles on all of the faults represented in the block model.

Effects of Parameter Choices and Simplifications
[19] Our modeling approach involves several simplifications. Among the most obvious are assuming a Maxwell viscoelastic lithosphere and representing rapid afterslip as part of the (spatially uniform) coseismic slip in the viscoelastic models, and simplifying the southern California fault geometry in the block models. Assuming that the same segments rupture repeatedly in the same manner at regular intervals is the simplest possible representation of the characteristic earthquake hypothesis. Non-periodic earthquake cycles can influence inter-seismic deformation and hence ghost transients [Hetland and Hager, 2005;Meade and Hager, 2004].
[20] If smaller earthquakes occur along the modeled segments (together with the rarer larger events), we suggest that their effects are likely to be unimportant. Smaller earthquakes (on shorter rupture segments) produce small amplitude, short wavelength perturbations. Furthermore, if these earthquakes occur frequently, their Savage parameter values (recurrence time/Maxwell time) are very small and nearly cycleinvariant viscoelastic deformation should result [e.g., Savage and Prescott, 1978]. [21] We assume an elastic plate thickness of 25 km in the viscoelastic earthquake cycle models. Though this is a simplification, we feel that it is justified in part by non-volcanic tremor depths of 20-25 km along the San Andreas Fault in the study area [Shelly, 2010]. Models with thinner elastic plates produce largeramplitude ghost transients that are more localized around the rupture. As noted above, we re-ran the Garlock earthquake-cycle model with a 15 km plate (and a low substrate viscosity, i.e., the CJ11 rheological structure), and the current-day ghost transient velocity vectors differed little from those computed using our model M1. [22] Our test models and other block model results show that locking depth is not well constrained by geodetic data, and that inferred slip rates are insensitive to variations in this parameter [McCaffrey, 2005;Meade and Hager, 2005]. We used a locking depth of 15 km in our block modeling calculations.
[23] The very simple models we present here provide a "worst-case" scenario for possible time dependence of viscoelastic deformation around faults. Welldeveloped plate boundary faults with large cumulative displacement likely have low viscosity ductile shear zones extending through the lithosphere. Seismic studies along the central SAF clearly show that this structure fully penetrates the crust [Shelly, 2010] and geologic evidence for viscous strain localization along exhumed fault zones is common [e.g., Cole et al., 2007]. Models of earthquake cycles along strike-slip faults with viscous shear zones generate steadier deformation throughout the earthquake cycle than layered viscoelastic models that yield similar near-field post-seismic deformation [Vaghri and Hearn, 2012]. Layered models incorporating nonlinear viscoelastic rheologies tend to localize viscous strain within a channel below faults, producing similar behavior [Hearn et al., 2006, Takeuchi andFialko, 2012].
[24] Three-dimensional viscoelastic earthquake cycle models (with finite-length ruptures) display less time variability and more localized strain late in the earthquake cycle than two-dimensional models [Chuang and Johnson, 2011]. This means that lithosphere viscosities that are lower than past estimates may be inferred from inter-seismic deformation. For the North Anatolian Fault, models incorporating 3-D earthquake-cycle effects [Kokum and Johnson, 2011] suggest a broader range of admissible lower crust and upper mantle viscosities than had been inferred from 2-D earthquake cycle models [Hearn et al., 2002b]. The flip side of this important finding is that to explain discrepancies between geologic and geodetic slip rates, very low substrate viscosities are required (e.g., our model M1) [Chuang and Johnson, 2011].

Conclusions
[25] Our models suggest that if asthenosphere viscosities are low (3 Â 10 18 Pa s), the current GPS velocity field is significantly perturbed by the 1857 M 7.9 San Andreas Fault (SAF) earthquake; that is, current strain rates around the SAF are significantly lower than their average values. Correcting the GPS velocity field for this perturbation (ghost transient) adds about 5 mm/a to the inferred SAF slip rate along the Mojave and San Bernardino segments, consistent with the results of JEA2007 and CJ11. GPS velocity perturbations due to the Garlock Fault [M 7.5 1640 rupture segment] and the White Wolf Fault [M 7.3 1952 Kern County earthquake segment] are smaller. For both, inversions of the corrected and uncorrected GPS velocity fields yield near-identical slip rates for all of the faults in our block model. This suggests that either the large discrepancy between geodetic and geologic slip rates for the Garlock Fault is not due to a ghost transient or viscoelastic effects from faults we did not model (e.g., from the 1990s Mojave earthquakes) are significant. If southern California lithosphere viscosities are of the order of 5 Â 10 19 Pa s or greater (at all depths), viscoelastic earthquake-cycle effects cannot explain apparent discrepancies between geodetically and geologically inferred fault slip rates.