Imaging of aseismic fault slip transients recorded by dense geodetic networks

Summary

1 INTRODUCTION

Most subduction zones accommodate more than half of their fault-slip aseismically (Pacheco et al. 1993). This large component of plate-motion has been generally inferred to occur as steady-state fault creep. However, as dense, continuous geodetic networks have started to record data on the overriding plates of a number of subduction zones it has become clear that much of the aseismic deformation occurs during relatively brief periods. There are now large-scale permanent GPS arrays above subduction zones in Japan (GEONET) (Miyazaki et al. 1997; Mazzotti et al. 2000), the Pacific Northwest (PANGA)(Miller et al. 2001), and elsewhere. The episodes of rapid aseismic deformation observed by these arrays are sometimes associated with earthquake afterslip (Kawasaki et al. 1995; Burgmann et al. 2001), but increasingly aseismic slip events are being identified that are not associated with any seismic earthquake making them ‘silent earthquakes’ (Beroza & Jordan 1990). Silent earthquakes have been observed by continuous GPS networks along numerous subduction zones including southwest Japan (Hirose et al. 1999; Ozawa et al. 2001), Cascadia (Dragert et al. 2001), the middle America Trench (Lowry et al. 2001), and the Aleutions (Freymueller et al. 2001). Prior to the advent of large scale GPS arrays, aseismic slip was documented on the San Andreas Fault with strain and creep meters (Linde et al. 1996). Aseismic deformation transients resulting from both fault-slip (Cervelli et al. 2002) and dyke intrusion (Aoki et al. 1999) have also been observed in volcanic areas. Thus, with the advent of continuous geodetic networks a new class of deformation events has become available for detailed study.

There are numerous important reasons for engaging in detailed studies of these events which go beyond the initial documentation stage. Aseismic slip can have a strong influence on the moment budget of faults which needs to be quantified for accurate seismic hazard estimates. Aseismic slip also redistributes stress in the crust affecting the locations of future earthquakes. Transient deformations proceed some volcanic eruptions (e.g. Linde et al. 1993; Owen et al. 2000), and deep aseismic slip, similar to the recently observed transients, preceded the largest subduction zone thrust earthquake ever, the 1960 Chile Mw 9 earthquake (Kanamori & Cipar 1974; Cifuentes & Silver 1989; Linde & Silver 1989). Finally, detailed imaging of aseismic fault-slip transients may help constrain fault zone constitutive laws.

There are major challenges in developing inverse methods that can determine detailed information about the space–time variations in the aseismic slip-rate on faults from geodetic data. First, we do not know a priori the nature of the time variations that we are trying to detect. Secondly, space geodetic measurements contain spatially and temporally correlated errors owing to atmospheric path delays, and in the case of GPS multipath and random benchmark motions (Wyatt 1989). This is a fundamentally different situation than the earthquake seismology inverse problem where the only significant contribution to the ground-displacement time-series is that resulting from fault-slip. Thirdly, owing to the imperfect resolution of the data, any inversion procedure will involve some spatial and temporal smoothing. It is very desirable that the amount of smoothing be determined by some rigorous procedure, not simply the analysts prejudice about what the slip-distribution should look like.

Segall & Matthews (1997) developed a recursive Kalman filter algorithm for determining the space–time variation in fault slip-rates from a set of geodetic time-series. This method, referred to as the Network Inversion Filter (NIF), incorporates several features, such as a stochastic description of local benchmark motion, that are advantageous when inverting geodetic data. The principal advantages of the NIF over other methods include: a non-parametric description of the slip rate as a function of time, and the fact that the filter operates on raw position time-series rather than derived quantities such as average displacement rate. Subsequent improvements in the technique have included separating the spatial and temporal smoothing of the slip-rate distribution (Segall et al. 2000) and estimation and removal of reference frame errors in the GPS data (Miyazaki et al. 2003). These studies have utilized linear descriptions of the underlying physical system and measurement process. The linear formulation has lead to two primary limitations of the filter in its application to real data sets. First, non-linear constraints on the model, in particular enforcing slip-rate non-negativity, have not been incorporated often leading to the need to over smooth the solutions to avoid unphysical solutions (i.e. left-lateral slip on the San Andreas Fault). Secondly, the hyper-parameters which control the extent of spatial and temporal smoothing have had to be estimated separately using a maximum likelihood procedure. Using a prediction error decomposition, Segall & Matthews (1997) showed that the likelihood could be computed using a recursive filter that inverted the state-covariance matrix at each epoch. While relatively fast, compared to inverting the full data covariance matrix, this approach proved to be computationally burdensome—the full data set needed to be run through a forward filter for every choice of smoothing (hyper) parameters. Moreover, the maximum-likelihood estimates of the smoothing hyper parameters were often disregarded in favour of a significantly oversmoothed solution as a means of enforcing slip-rate positivity. In this paper we seek to correct these two related problems by employing a non-linear system and measurement model. Our new formulation allows the incorporation of non-negativity constraints through the use of pseudo-data, and incorporates the hyper-parameters directly into the state vector so that they are estimated by the filter simultaneously with the slip-rate. The new system and measurement models are solved using an extended Kalman Filter. This non-linear approach provides significantly improved spatial and temporal resolution, lower computational cost, and is equivalent to the maximum-likelihood estimate in high signal-to-noise ratio (SNR) cases which include the larger aseismic transients observed in real data.

2 METHOD

We model GPS station positions as a function of time t, as

1

where the first term on right hand side represents coseismic offsets X^cos(k) at times t_eq(k), and H(t) is the Heavyside function. The second term on the right hand side represents deformation due to transient aseismic slip. We relate the time dependent site motions to slip on faults in an elastic medium as a function of position ξ and time, s_p(ξ, t−t₀) through Green's functions G^r_pq(x, ξ). To date, Green's functions relating slip to surface displacement have been computed using analytical expressions for dislocations in uniform half-spaces, although there is no impediment to computing Green's functions in heterogeneous or layered media. In (1)p, q, r= 1, 2, 3, and summation on repeated indexes is implied, n_q(ξ) is the unit normal to the fault surface A(ξ).

In (1) the fault slip includes both stead-state and transient components. In some cases where modelling the steady-state deformation field is not possible (at least to sufficient accuracy), we have found it useful to separate the steady station velocities by adding a term V_sec·(t−t₀) to (1). In this case s_p in (1) represents only the transient component of fault slip. Given sufficient data it is possible to estimate the secular velocity V_sec along with the transient fault slip. In other cases we have estimated V_sec separately and subtracted it from the position time-series prior to filtering.

The remaining terms are related to measurement and reference frame errors. The term Ff(t) represent reference frame errors, where F is a linearized Helmert transformation and f(t) is a vector of rigid body translations, rotations, and a scale factor (Miyazaki et al. 2003). The fourth term on the right hand side of (1), L(x, t−t₀), represents random benchmark wobble, which we model as a Brownian random walk with scale parameter τ (units length time^−1/2) (Wyatt 1989). While some studies have advocated a flicker noise model for benchmark motion (Mao et al. 1999; Zhang et al. 1997), Langbein & Johnson (1997) using over 10 yr of data show that a random walk fits the data quite well. Without a long time-series it is very difficult to distinguish between flicker noise (with a 1/f spectral decay) and random walk (with a 1/f² spectral decay). A key feature of the Network Inversion Filter is that it can distinguish spatially correlated transient signal from site specific coloured noise. The reason for this is that elastic deformation causes a spatially coherent signal, whereas the local benchmark motions are, by definition spatially incoherent. The fact that the local benchmark motions are spatially incoherent may be more important than the precise form of the spectral decay. The final term, ε, represents observation error, which we take to be normally distributed with zero mean and covariance σ²Σ^GPS, where Σ^GPS is the covariance matrix of the GPS positions, and σ² is a scale factor to account for unmodelled errors such as mulitpath, or azimuthally varying path delays.

To estimate the time-dependant quantities s_p(ξ, t−t₀), f(t), and L(x, t−t₀) we incorporate (1) into a state–space model of the underlying system and measurement processes. The state–space system is most generally specified by a non-linear measurement equation:

2

and a non-linear state-transition equation:

3

where the k subscripts indicate observation epoch, d_k is the data vector at epoch k, x_k is the state vector at epoch k which contains all the quantities necessary to describe the system. The state transition equation predicts the state of the system at time k+ 1 given an estimate of the state at time k. The measurement and process noise are described by the covariance matrices R and Q respectively (see Appendix A).

Following Segall & Matthews (1997) we model the local benchmark motion at the jth station as a random walk with scale parameter τ

4

where dw is formal white noise with unit variance. In this paper, τ is assumed to have the same value for all sites, usually around 1–2 mm yr^−1/2 (Langbein & Johnson 1997). Similarly the slip-rate of the ith subfault is modelled as a random walk with scale parameter α

5

Thus, the cumulative slip on any subfault is a integrated random walk and the slip accelerations between observation epochs are modelled as a white noise process. The hyper-parameter α controls the temporal roughness of the slip-rate function (Segall & Matthews 1997).

2 .1Pseudo-data

There are two types of constraints we would like to enforce on the solution. The first is that the slip-rate distribution be one-sided (i.e. no left-lateral slip on the San Andreas fault). The second is that the slip-rate distribution be spatially smoothed. Following Segall et al. (2000) we incorporate these constraints into the state estimation process through the use of pseudo-data

6

7

Eq. (6) has been used successfully in previous implementations of the NIF to enforce spatial smoothing (Segall et al. 2000; Burgmann et al. 2002; Miyazaki et al. 2003). Eq. (7) can be reformulated as an equality constraint by introducing the dummy variable λⁱ_k which is estimated simultaneously with (see Hel-Or & Werman 1996)

8

The hyper-parameters γ and ρ control the extent to which the two constraints are enforced in this stochastic formulation. Since one-sidedness is a physical constraint we typically fix ρ at some very small value (∼1/1000 of the expected slip-rates) to strictly enforce it. However, in cases where the secular velocity has been subtracted from the data time-series we are actually solving for deviations in slip-rate relative to ‘steady state’. In these cases, values as low as −1 times the plate motion rate are physically allowable as they would be interpreted as a region that is normally fully uncoupled becoming completely coupled. The hyper-parameter γ which controls the extent of spatial smoothing applied at any epoch has been estimated in previous NIF studies by maximum-likelihood using a prediction error decomposition (Segall & Matthews 1997). This method is ideal in a low signal-noise-ratio case, but it is computationally expensive. In this paper, γ is estimated directly by the filter algorithm as discussed below.

2 .Direct estimation of hyper-parameters

In previous linear formulations of the NIF the hyper-parameters which control spatial and temporal smoothing of the slip-rate distribution, γ and α, have appeared in the covariance matrices of the pseudo-data noise and the process noise respectively. Since we will need to use a non-linear state–space model to incorporate the positivity constraint in (8) there is also the possibility of directly estimating γ and α in the filter. Hyper-parameters which control the process noise, such as α and τ can be incorporated into the measurement equation by defining the associated random-walk process as:

9

with the slip and slip-rate amplitude given by s=αW and . The contribution of slip on the ith subfault to displacement at the jth station is G_ijs_i (see appendix A). , and Wⁱ are all included in the state vector and estimated by the filter.

Hyper-parameters which specify the covariance matrices of the data and pseudo-data, such as γ, ρ, and σ, can be directly incorporated in the observation equation by simply dividing equations like (6) by the appropriate standard deviation

10

Using these direct relationships between the hyper-parameters and the data/pseudo-data the hyper-parameters can then be considered part of the state-vector (see Appendix A). Augmenting the state-vector with the unknown parameters of the system model and solving them using a linearized (extended) filter is a common process in the data assimilation literature referred to as adaptive estimation (see Wunsch 1996, section 6.6.2).

2 .3Extended Kalman filter algorithm

The non-linear state–space system described by eqs (2) and (3) (and detailed in appendix A) can be solved for estimates of the state vector at each of the observation epochs using an algorithm known as the extended Kalman filter (see Gelb 1974). The filter is initialized with an a priori estimate of the state vector, x_0|0 and its covariance Σ_0|0 (where the subscripts x_k|j indicate the estimate of the state vector at epoch k given data through epoch j). Typical values of these priors will be discussed in the synthetic and data examples in subsequent sections. The a priori estimates are then predicted forward one time step (there is no requirement that subsequent time-steps be of the same size). The predicted state is then compared with the data from the first epoch and any residual is used to update the estimate of the state at epoch 1. The prediction and update process are then iterated through the entire data set.

2 .3.1Prediction step

11

12

2 .3.Update step

13

14

15

where

16

After reaching the last observation epoch we have estimates of the state vector at each epoch given data up to that epoch, x_k|k. To obtain estimates of the state vector at each epoch given all the data, x_k|N (where N is the number of observation epochs), it is necessary to run a similar recursive procedure backward in time. This algorithm is known as smoothing (see Rauch et al. 1965; Segall & Matthews 1997).

3 ȀSIMULATION

In this section we present the results of inverting a data set for a simulated aseismic slip transient on a dipping thrust fault recorded by an array of GPS receivers on the overriding plate (see Fig. 1). The transient begins near the southern downdip corner of the fault and propagates both updip and along strike with a rupture velocity of 3 km d⁻¹ and a rise-time at any one point of 17 d. The total duration of the event is about 45 days. Data time-series are computed for the East and North components for the 42 sites shown in Fig. 1. Each data time-series is the sum of the signal resulting from the fault-slip, the contribution from local benchmark motion, and the contribution from measurement error (assumed to be a white noise process). While the absolute scale of the slip and resultant displacements (mm, cm, etc.) is arbitrary, we have chosen the relative size of the fault-slip, benchmark motion, and measurement error to be representative of the best recorded transients (e.g. Miyazaki et al. 2003). The measurement errors are normally distributed with a scale of 1, the benchmark motion is a random walk with a scale parameter of 5 yr^−1/2, and the maximum surface displacement from fault slip is 29. An example of a data time-series and its contributions from fault-slip and benchmark motion is shown in Fig. 3. The following sections discuss the inversion of this data set with positivity constraints using both the original maximum likelihood algorithm of Segall & Matthews (1997) and the new algorithm which estimates the hyper-parameters directly by including them in the state vector.

Figure 1

Map of the station distribution (asterisk) and fault grid (coloured rectangles) for our synthetic tests. The shading indicates the final slip distribution. The fault strikes North-South and dips 15° to the East. The upper (Western) edge is at 10 km depth and the lower (Eastern) edge is at 30 km depth.

Open in new tabDownload slide

Figure 3

The top panel shows the East component data (blue asterisks) for the station located at (−5 km, 15 km) in Fig. 1 which is the sum of the fault-slip contribution (middle panel), the benchmark wobble contribution (bottom panel), and white noise measurement error. The black lines in each panel show the fit from running the extended filter on the entire data set.

Open in new tabDownload slide

3 .1ȀMaximum likelihood estimation of hyper-parameters

In previous implementations of the Network Inversion Filter, the hyper parameters which control the spatial and temporal smoothing of the slip-rate distribution, γ and α, have been estimated before running the filter using maximum likelihood (Segall & Matthews 1997). The best estimates of the hyper parameters are found by maximizing the log likelihood L over the set of hyperparameters Θ=[σ, α, γ]. Maximizing L for the simulated data set results in an estimate of α= 1.9e4 and γ= 4.8 (see Fig. 4).

Figure 4

Likelihood surface (−2L) as a function of γ and σ. Asterisk marks the minimum, corresponding to the maximum likelihood estimate.

Open in new tabDownload slide

Running the synthetic data through a version of the filter which used hyper-parameters fixed to the maximum likelihood values and incorporated positivity constraints, produced the slip-rate estimates shown in Fig. 5. In general the estimated slip and slip-rate distributions are similar to, but somewhat more smoothed than, the true distribution. We measured the accuracy of the inversion result in terms of the summed squared residual between the true slip and slip-rate distributions and the recovered distributions:

17

18

where s and are the true distributions and and are the estimated ones. Thus the quantities r₁ and r₂ are analagous to a variance reduction. For the maximum likelihood inversion r₁= .24 and r₂= .04 (a value of 0 indicates perfect recovery of the true distribution).

Figure 5

Maps of the estimated slip-rate using maximum-likelihood estimates of the hyper-parameters, similar to Fig. 2.

Open in new tabDownload slide

Figure 2

Maps of the true slip-rate distribution every 4 d from the start of the transient to its end. The event nucleates near the southern downdip edge of the fault and propagates updip and to the North.

Open in new tabDownload slide

3 .2ȀDirect estimation of hyper-parameters

As described in Section 2.2 and Appendix A, an alternate approach to [maximum-likelihood] estimation of γ and α is to directly incorporate them into the state vector so that they are determined as part of the estimation process. This is useful because the filter need only be run once rather than once for each combination of hyperparameters, Θ_i, to be tested. The key additional piece of information that is necessary to include the hyper-parameters in the state vector is an a priori estimate, x_0|0 and it's covariance, Σ_0|0. We assign apriori values of α and γ that would result in a very smooth slip-rate distribution, but also assign very large variances so that the a priori information will contribute little to the final estimate. We also assume that the hyper parameters do not vary as a function of time during the observation period, i.e. there is no process noise for these state vector elements and that there is no a priori covariance between them. Thus the best estimate of the hyperparameters comes from the state vector at the final epoch, x_N|N. A similar approach is taken with the prior information for the slip-rate which is set to an a priori value of zero with a large variance. The hyper-parameters which control the benchmark motion, τ, and the scale parameter for the GPS covariance matrix, σ, are assumed known (say based on an analysis of a quiet portion of the time-series), and fixed to their correct values in these simulations as was also done in the maximum-likelihood simulation above. In principle there may be sufficient information in a long time-series to estimate the values of these parameters as well, but we have not experimented with this.

For the synthetic data set we first ran the filter with an apriori values of α= 1 ± 1000 and γ= 1 ± 1000, this produced estimates of α= 2.2e3 and γ= 116 indicating that our prior variance on α was too small. So the filter was rerun with an a priori values for α and γ of 1e3 and 1e2, and a priori variances of 1e6 and 1e6 respectively. This resulted in estimates of α= 1.3e4 ± 402 and γ= 780 ± 2e3. (runs with higher a priori values of α and γ did not result in significantly different final estimates). These results are similar to the maximum likelihood estimate in that α estimates are similar and that γ is poorly constrained (see Fig. 4). The estimates of the slip-rate distribution from the extended filter are shown in Fig. 6. The slip-rate estimates from this filter provide a better recovery of the true slip-rate values (r₁= .14, r₂= .07) than the maximum-likelihood filter as evidenced by the higher peak slip-rates in Fig. 6 than in Fig. 5. The datafits for one station are shown in Fig. 3. Although the input signal was purely thurst slip, we inverted for both the thrust and strike-slip components of motion on each subfault. The estimated erroneous amount of strike-slip motion estimated by the filter corresponds to only about 3 per cent of the total slip indicating that by using both horizontal components of GPS displacement, we are able to effectively resolve the azimuth of slip. In general, the filter can separate the spatial coherent signal due to fault-slip from the spatially incoherent signal resulting from benchmark-motion. Thus, the extended filter provides nearly optimal estimates of the slip-rate distribution in the high signal-to-noise ratio case without the computational cost of the maximum likelihood filter.

Figure 6

Maps of the estimated slip-rate using estimates of the hyper-parameters determined by the filter, similar to Fig. 2.

Open in new tabDownload slide

3 .ȀLow signal-to-noise ratio example

The previous simulation incorporated a transient with a high signal to noise ratio which is representative of the largest aseismic transients. In this section we address the problem of detecting a transient in a noisy data set where the fault-slip signal is not unambiguously present in the GPS time-series by eye. We applied the filter to a case where the other time-dependent contributions to the GPS time-series are larger than the fault-slip contribution. The amplitudes of the fault-slip and benchmark wobble terms were chosen to have approximately equal maximum amplitude contributions to the GPS time-series. The duration of the slip-transient was also increased in this episode to increase its similarity to the bench-mark wobble signals. The resulting data is shown in Fig. 7. The transient is contained between days 50 and 100 and is visible in the trench perpendicular component of some of the most sensitive stations, but numerous stations exist where the net displacement over the 150 d time period is of the opposite sign than what would be expected from fault slip owing to the contribution from the benchmark wobble term. In this low signal to noise ratio case, the filter is more unstable than the previous case in that using very large a priori uncertainties on α and γ can result in (unphysical) negative values of these hyper-parameters. Thus it was necessary to start with a smooth prior for α and γ with low uncertainty. This resulted in an estimate of α that was larger than 2σ above the prior, so the filter was rerun a few times gradually increasing the apriori estimate of α and it's uncertainty. The final run, when the estimate of α was within 1σ of the a priori values, was initiated with prior values for α and γ of 4000 and 100 with associated uncertainties of 1e3. The final estimates of the hyper parameters were α= 4466 ± 320 and γ= 226 ± 549.

Figure 7

Data time-series for the low signal-to-noise ratio simulation.

Open in new tabDownload slide

The true and estimated slip-rate distributions resulting from the low SNR test are shown in Figs 8 and 9. The estimated slip-rate distribution is qualitatively quite similar to the true distribution although the peak slip-rates are not quite recovered (r₁= .71, r₂= .69). The estimated moment-rate distribution captures the duration and even some of the internal structure of the transient (Fig. 10), but a small amount of the benchmark wobble and measurement noise is mapped into the fault-slip components of the state vector when no transient is present. Again, the input signal was purely thrust slip erroneous amount of strike-slip motion estimated by the filter remained small (about 10 per cent of the total slip). Thus even in the case with a SNR less than 1 for the overall data set, the filter is capable of pulling out an accurate and useful description of the transient's spatial extent, duration, sense of slip, and propagation.

Figure 8

Input slip-rate distribution for the low SNR simulation.

Open in new tabDownload slide

Figure 9

Recovered slip-rate distribution for the low SNR simulation.

Open in new tabDownload slide

Figure 10

Input (dashed) and recovered (solid) source-time functions for the low SNR simulation.

Open in new tabDownload slide

4 APPLICATION TO THE 1999 CASCADIA SLOW EARTHQUAKE

In 1999 a cluster of continuous GPS stations run by the Gelogical Survey of Canada and the Pacific Northwest Geodetic Array (PANGA, (Miller et al. 2001)) network recorded a 40 day period of anomalous displacements, on the order of a few mm, resulting from slow-slip on the deep portion of the Cascadia thrust interface (Dragert et al. 2001). The stations which showed the anomalous displacements are contained within a ∼200 km by ∼200 km region extending from Seattle to the central portion of Vancouver Island (see Figs 11 and 12). Dragert et al. (2001) matched the anomalous displacements with a forward model consisting of a uniform zone of 2.1 cm of slip between about the 30 and 40 km depth contours of the thrust interface with slip tapering to zero by about ∼25 km depth. This model did not attempt to fit the time dependence of the signal beyond breaking the slip-event into three time-periods which showed general northward propagation consistent with the observed delay in the arrival of the event at the northern stations (Dragert et al. 2001). Here we present the results of applying the version of the NWIF that includes both positivity constraints and direct estimation of the hyper-parameters to the Cascadia data set.

Figure 11

Map showing the station distribution (black triangles and 4 letter labels) and fault-geometry (rectangles) for the 1999 Cascadia slow earthquake.

Open in new tabDownload slide

Figure 12

East-component time-series (asterisks) for the time-period used in the filter estimation. The contributions from secular velocity, reference frame error, and benchmark wobble have been removed. The solid-lines show the fit to the data from the estimated slip distribution.

Open in new tabDownload slide

We model the daily station positions (see Fig. 12) determined by (Dragert et al. 2001) as the sum of the contributions from transient fault slip, reference frame error, secular velocity, benchmark wobble, and measurement error. Secular velocities at these stations vary on the time-scales of years (Dragert et al. 2001) so to avoid mapping any of this variability into our slip-rate estimates we inverted data from a relatively short time period centered on the transient, 1999.5 to 1999.85, in which the secular-rate could be considered constant. The secular velocities were determined ahead of time using a Kalman filter that estimated the velocity, reference frame error, and benchmark wobble components of the time-series in the period between 1999.2 and 1999.6. The contribution from secular velocity only was then extrapolated and subtracted from the data time-series.

We used a 3-D fault geometry derived from the model of the Cascadia thrust interface published by Flueck et al. (1997) (see Fig. 11). Rectangular fault segments were fit to a set of points specifying the thrust interface at approximately 5 km depth intervals between 20 and 50 km depth. The location and dip of the interface are specified in a manner that best approximates an elastic half-space by accounting for topography and preserves the dip-angle of the thrust interface (Flueck et al. 1997). There were 72 subfaults with average lengths of about 25 km and widths of about 15 km. Green's functions were then calculated for the subfaults assuming a homogeneous isotropic half-space (Okada 1985). To account for station DRAO being held fixed in the GPS processing (Dragert et al. 2001), the Green's function for this station was subtracted from those of other stations. Despite station DRAO being held fixed in the GPS processing, short term random translations of the network are still present in the time-series (Dragert et al. 2001). To account for these errors, we estimated the translations of the GPS reference frame, but owing to the relatively small scale of the network, we did not solve for rotations or scale. Additional process noise was introduced in the benchmark wobble terms at a number of epochs to account for known antenna/hardware changes (Dragert, personal communication). The temporal smoothing parameter, α was initialized at a value of 50 ± 50 and had a final estimate of 63 ± 30. Larger values of the a priori variance of α led to unstable hyper parameters, so the a priori estimate and its uncertainty were gradually increased until the final estimate was within 1σ of the a priori value, similar to the low SNR simulation. The spatial smoothing parameter was initialized with a value of .1 ± 1000 and had a final estimate of 0.8 ± 72. The other hyper-parameters, τ, σ, ρ, were fixed at 2 mm yr^−1/2, 10, and 1e-4 (m yr)⁻² respectively.

The results of the filter are shown in Figs (12 and 13). The slip transients begins at depths of 30–35 km just south of station ALBH. It then expanded to the south before propagating both updip and along-strike to the north. The peak slip was about 8 cm in the region just south and updip of ALBH, while the peak slip-rates were about 1 m yr⁻¹ (though this value trades off with grid-size). We estimated both the thrust and strike-slip component of slip on each subfault, the magnitude and direction of the final slip is shown in Fig. 14. On most faults the direction of slip is resloved to be roughly parallel to plate motion direction, but for the deepest row of subfaults, resolution of the azimuth (and amplitude) of slip is poorer. The moment-release distribution was not smooth in time. There are two peak periods of high slip-rates centered around day 240 and 270 separated by a roughly a 15–20 d period between about julian day 250 and 265 where slip-rates temporarily decreased below the 0.5 m yr⁻¹ level on all subfaults (roughly 1999.71 in Fig. 13). The inferred late moment release results from the late changes at stations UCLU and NANO (Fig. 12). The total moment-release is equivalent to about an M_w 6.9 event and the average stress-drop was approximately .01 MPa. The duration of the period of high slip-rate on any one subfault ranged from about 5–30 d.

Figure 13

Maps of the estimated slip-rate as a function of time and station distribution (black triangles).

Open in new tabDownload slide

Figure 14

Final slip distribution. Arrows denote magnitude and direction of slip on each fault segment for which the slip was larger than 1.0 cm. Faults with less slip than this are considered unresolved. In general the slip vectors are closely aligned with the plate-motion vector except for the deepest row of subfaults, which are poorly resolved. Grey triangles show station locations.

Open in new tabDownload slide

5 DISCUSSION

The algorithm described in the previous sections has several advantages for interpreting continuous geodetic data and may be suited for real-time fault monitoring. In particular, the inclusion of a positivity constraint on slip-rate increases the spatial resolution of slip-rate anomalies over previous versions of the NIF. This feature can be accommodated with either direct or maximum-likelihood estimation of the filter hyper-parameters. The incorporation of the hyper-parameters into the state-vector may become quite useful in real-time applications with high sample-rate data when the maximum-likelihood algorithm is too expensive to implement. Including the hyper-parameters may also suggest an event detection algorithm, as they are only forced away from the smooth a priori values during periods of anomalous slip.

The slip-rate anomalies shown in Fig. 13 show several characteristics that are analogous to slip in ordinary earthquakes, but with longer time-scales. The slip-rate on a portion of the fault takes a finite amount of time, from a few to ∼20 d, to reach its peak. This property is much more similar to unstable earthquake slip than the exponential decay that is usually observed in stable earthquake afterslip. The 1999 Cascadia slow-event also shows something like slip-pulse propagation where the southern portion of the fault had finished slipping before the northern portion began. The key differences between the Cascadia event and a normal earthquake are: (1) the longer time-scales, both for the total event and local rise-time and associated the lower propagation and slip velocities. (2) The much lower average stress drop resulting from the relatively small amount of slip over a large area. However, these two characteristics are very similar to typical afterslip distributions. Thus if earthquakes represent fault-slip governed by unstable velocity-weakening friction and afterslip is an example of fault-slip governed by stable velocity-strengthening friction, then silent earthquakes appear to have behaviour somewhat in-between, sharing important characteristics with both.

The aseismic slip transients that have been observed recently at numerous subduction zones Miyazaki et al. (2003), Dragert et al. (2001), Lowry et al. (2001), Freymueller et al. (2001) have several features in common. Primarily that they all rupture the portion of a subduction zone thrust interface immediately below the ‘locked’ zone, and that while the fault slip-rate clearly accelerates over a period of time, there is a limit to the acceleration and the slip does not becoming unstable. This portion of the thrust interface is predicted to be ‘conditionally stable’ in the rate-state friction formulation owing to the temperature approaching the onset of plasticity in feldspar at about 450°C (Scholz 1998). The conditionally stable regime of rate-state friction has been well characterized analytically and in single degree of freedom numerical simulations (Rice & Gu 1983). Under steady-state loading, faults in the conditionally-stable regime slip aseismically unless there is a sudden perturbation to the system. If the perturbation is large, the fault can be pushed over a stability boundary causing it to become unstable and slip seismically. However, if the perturbation is such that the fault approaches but does not cross the stability boundary, a period of stable sliding with increased velocity will follow as the fault evolves back to its steady-state behaviour. This type of behaviour represents one potential mechanism for explaining slow earthquakes. However the types of triggering events that could lead to a perturbation in the state of the fault, primarily the shear and normal stress, are hard to constrain for slow events. There were no nearby large earthquakes to trigger either the Cascadia slow event, and there is no clear connection between the 1996, M_w 6.7, Hyuganada earthquakes and the 1997–1998 Bungo-Channel slow earthquake (Miyazaki et al. 2003). However, stress-state changes that initiate the slow event could result from a number of other mechanisms that are difficult to observe from the surface, such as changes in pore fluid pressure. There is also a portion of the conditionally stable rate-state parameter space where oscillatory behaviour between stable sliding and locked behaviour is observed in 1-D simulations. This has been proposed as a potential explanation for the creeping section of the San Andreas fault (Scholz 1998) where repeated creep events have been observed. If the deep slow-slip transients on subduction zones are roughly periodic events as has been suggested for Cascadia (Miller et al. 2002), and the 1-D simulations can be extrapolated to real faults this may provide a mechanism that doesn't require a triggering event. Additionally, recent simulations of rate-state friction on a 2-D subduction interface have generated slow earthquakes confined to the conditionally stable region below the seismogenic zone (Shibazaki & Iio 2002).

ACKNOWLEDGMENTS

This work was supported by NSF-EAR0074181 and NASA SENH NAG5-6167. We thank Herb Draggert for giving us access to his processed station coordinates which were derived from GPS data collected by the Geological Survey of Canada and the PANGA network. We thanks Stephane Mazzotti and Takeshi Sagiya for helpfull reviews.

References

Appendix

Appendix A: STATE VECTOR SYSTEM

This appendix details the non-linear state–space system used in the extended Kalman filter inversion algorithm for the case where the secular velocity term has been subtracted from the data prior to running the filter. The system is described by the observation equation,

(A1)

and the state-transition equation:

(A2)

The state vector, x_k contains all the parameters necessary to describe the state of the system, i.e. the terms on the right-hand side of eq. (1). The state-vector has elements corresponding to four types of model-parameters, fault-slip , benchmark motion (βⁱ_k), reference frame (f), and hyper-parameters (τ, σ, α, γ, ρ),

(A3)

There are three types of ‘data’ in the vector d_k=[d^GPS_k, d^∇_k, d⁺_k], the GPS observations at epoch k, the spatial smoothing pseudo-data, and the positivity pseudo-data. The observation equation for each type of data can be written as:

(A4)

(A5)

(A6)

where i and j indicate the station and subfault index respectively as in Section 2. Thus, the covariance matrix of the concatenated data vector is,

(A7)

The extended Kalman filter algorithm uses the full non-linear h(x_k), i.e. eqs (A4)–(A6), to calculate the residual at each epoch. The update step uses the partial derivatives of eqs (A4)–(A6) with respect to the elements of the state-vector to form the matrix H.

The state transition equation is specified by the stochastic models assumed for each process, which are discussed in the main text. In our case, it is actually linear, x_k+1= Tx_k+δ_k+1, where

(A8)

(A9)

The upper-left submatricies result from the random-walk and integrated random-walk models assumed for and S respectively. The values for the elements of Q corresponding to the positivity dummy variable λ were determined from numerical simulations of random-walk time-series. The hyper-parameters are assumed to be time-invariant, and Σ_f is the time-independent covariance assumed for the reference frame errors (Miyazaki et al. 2003).