Double-Difference Tomography: Method and Application to

the Hayward Fault, California

 

Haijiang Zhang, Clifford Thurber

 

 

Department of Geology and Geophysics

 University of Wisconsin-Madison

 

 

Submitted to BSSA

 

August 22, 2002

 

Abstract

We have developed a double-difference seismic tomography method that makes use of both absolute and relative arrival times. By reducing scatter in event locations using the more accurate relative arrival times, the method produces an improved (sharper) velocity model. Simultaneously, it yields event locations of a quality equivalent to those of the double-difference location method, while avoiding simplifying assumptions of that method. We test this method on a Hayward fault, California, earthquake dataset spanning 1984 to 1998. The earthquakes relocated by this method collapse to a thin line along the fault trace, forming a number of “streaks”, consistent with previous results. The velocity model yields a sharper velocity contrast near the source region than the conventional tomography method with only absolute catalog data.

 

 

Introduction

Local earthquake tomography (LET) has become a relatively routine application for use in seismically active regions covered by a dense seismic network. Common LET algorithms, however, do not take advantage of the many recent developments in earthquake location techniques. Many of these developments are aimed at improving relative and/or absolute location accuracy.

The accuracy of event hypocenters is determined by several factors, including the network geometry, available phases, and arrival time accuracies (Pavlis, 1986). If more phases are available, the event hypocenters will generally be more stable. For example, S-wave phases will be very useful in constraining event depths and providing information that helps to decouple the hypocenters from the structure in the inversion (Roecker, 1982; Gomberg et al., 1990). Due to the presence of noise, the arrival times picked either manually or automatically generally have errors (Douglas et al., 1997). Recent studies have shown substantial improvements in location precision for earthquakes and explosions when waveform cross-correlation (WCC) and event clustering techniques are used to improve arrival time estimates or determine high-precision relative arrival times (VanDecar and Crosson, 1990).  Example applications include Mount St. Helens (Fremont and Malone, 1987), Hawaii (Got et al., 1994; Rubin et al., 1998), California (Poupinet et al., 1984; Nadeau et al., 1994; Shearer, 1997; Rubin et al., 1999; Waldhauser et al., 1999; Waldhauser and Ellsworth, 2000, 2002), the Coso geothermal field, Nevada (Lees, 1998), the Soultz geothermal field, France (Rowe et al., 2002), and explosions at the Balapan test site (Phillips et al., 2001; Thurber et al., 2001). These studies are based on the assumption that waves generated by two similar sources, propagating along similar paths, will generate similar waveforms, and waveform cross-correlation can then be used to determine precise relative arrival times. These and other studies have demonstrated substantial improvement in the definition of seismogenic features and in the accuracy of relative locations of ground-truth events that is possible using multiple-event methods with high-precision absolute or relative arrival-time data.  In addition, they have also provided new insight into tectonic processes, earthquake recurrence, and earthquake interaction.

We can make an important distinction between the two fundamentally different ways WCC data has been used: (1) by directly using relative arrival times to determine relative event locations (e.g., Fremont and Malone, 1987; Got et al., 1994; Waldhauser and Ellsworth, 2000), or (2) by adjusting absolute arrival time picks to minimize discrepancies among relative arrival times (Dodge et al., 1995; Dodge, et al., 1996; Shearer, 1997; Rowe et al., 2002).  The advantage of the former approach is that it incorporates all the available information contained within the multitude of relative arrival time differences with a direct measure of quality (the correlation value) associated explicitly with each datum.  A disadvantage is that some simplifying assumption needs to be made to derive the locations from the arrival time differences. For example, in the methods of Fremont and Malone (1987) and Got et al. (1994), all events in a cluster have precisely the same take-off angle and azimuth to each station. As a result, the derived locations are ultimately relative, not absolute, so that some assumption must be made to end up with useful event coordinates (e.g., final locations are computed relative to a catalog-based cluster centroid). Waldhauser and Ellsworth (2000) proposed a different location algorithm, in which the spatial partial derivatives for a set of events are evaluated at the current location of each event. It is assumed, however, that the path anomalies from velocity heterogeneity are location independent. This assumption is valid for closely spaced events, but is not true for far apart events. As a result, such far apart event locations may be biased due to velocity heterogeneity (Got et al., 1994; Waldhauser and Ellsworth, 2000; Wolfe, 2002).  In contrast, the latter approach uses the relative arrival times to determine a much smaller number of adjusted arrival time picks, but these picks are absolute arrival times and so can be used to determine absolute locations (in an existing velocity model, or using tomography).

We have developed a new method that combines the advantages and avoid the disadvantages of the above approaches.  It is based on the code hypoDD of Waldhauser (2001), and makes use of both absolute and relative arrival time data.  The method determines a three-dimensional (3D) velocity model jointly with the absolute and relative event locations. This approach has the advantage of including relative arrival times with their quality values along with absolute arrival times, thereby not discarding valuable information by only using adjusted picks, and at the same time dispensing with simplifying assumptions about ray path geometries or path anomalies and producing absolute locations, not just relative locations. The velocity model obtained with double-difference tomography should also be superior to that from conventional tomography. With conventional tomography, event locations will be somewhat scattered due to imprecise picks, but in double-difference tomography, the use of the WCC time difference data removes most of this scatter, which will in turn remove some “fuzziness” from the velocity model. To demonstrate the effectiveness of the method, we have applied it to the Hayward fault dataset of Waldhauser and Ellsworth (2002) to resolve and to study the detailed local 3D seismic velocity structure.

 

Double-difference tomography

The body wave arrival time  from an earthquake  to a seismic station  is expressed using ray theory as a path integral

                                                                                        (1)

where is the origin time of event , is the slowness field and is an element of path length. The source coordinates  origin times, ray paths, and the slowness field are the unknowns.  The relationship between the arrival time and the event location is highly nonlinear, so a truncated Taylor series expansion is generally used to linearize Equation (1). This linearly relates the misfit between the observed and predicted arrival times  to the desired perturbations to the hypocenter and velocity structure parameters

                                                                      (2)

A similar expression exists for event , which is also observed at station

                                                                    (3)

Subtracting Equation (3) from Equation (2), we have

               (4)

Assuming that these two events are near each other so that the paths from the events to a common station are almost identical and the velocity structure is known, then Equation (4) can be simplified as

                                    (5)

where  is the so called double-difference (Waldhauser and Ellsworth, 2000).  is the difference between observed and calculated differential arrival times for the two events, and can also be written as

                                                                 (6)

The observed differential arrival times can be calculated from both waveform cross-correlation techniques for similar waveforms and absolute catalog arrival times. Equation (5) is known as the double-difference earthquake location algorithm (Waldhauser and Ellsworth, 2000).

In this approach, earthquake locations may be biased when inter-event distances exceed the scale length of velocity variations. Waldhauser and Ellsworth (2000) apply a distance-weighting factor to reduce or exclude data from event pairs that are far apart. Although the arrival difference data from such events may be excluded, they can still be linked in the inversion via a series of intermediate pairs (Got et al., 1994). For example, the pair  can be linked if the two pairs  and  are included.

To overcome this limitation, we use Equation (4) directly. It is known that there is a coupling effect between the event hypocenters and the velocity structure (Thurber, 1992). Our purpose is to determine not only the relative event locations, but also their absolute locations and the velocity structure. For this reason, we include the absolute arrival times in the inversion. By doing this, we can jointly determine the velocity structure and the relative event locations as well as the absolute event locations accurately. The equations for the double-difference tomography algorithm are:

         (7)

In our simultaneous inversion for velocity structure and event locations, velocity anomalies are constrained by seeking a first-order smooth model. Smoothing regularization should provide a minimum-feature model that contains only as much as structure as can be resolved above the estimated level of noise in the data (Constable et al., 1987; McCaughey and Singh, 1997). We apply the same smoothing weight to the horizontal and vertical directions. The system of linear equations (7) is solved by means of the LSQR algorithm (Paige and Saunders, 1982) for the damped least squares problem, with each equation weighted according to the a priori data uncertainty and misfit during each iteration, with large residuals rejected or downweighted by a biweight function (Waldhauser and Ellsworth, 2000).

 

Application to the Hayward fault

We applied the double-difference tomography method to earthquakes recorded between 1984 and 1998 by the Northern California Seismic Network (NCSN) on the Hayward fault, California, the same as those used by Waldhauser and Ellsworth (2002). The Hayward fault accommodates about 9 mm/yr of right-lateral relative plate motion between Northern American plate and the Pacific plate. Franciscan rocks, consisting mainly of graywacke mélange with smaller amounts of chert, shale, mafic volcanic rock, and limestone, form most of the upper crust to the west of the Hayward fault (Bailey et al., 1964). In the upper crust to the east of the Hayward fault, there exists folded and thrusted shallow marine sandstones and shales of the Great Valley Sequence overlying the Franciscan terrane (Hole et al., 2001). Rocks of the Franciscan terrane are faster than the Great Valley and younger sediments (Hole et al., 2001).

1489 earthquakes with magnitudes from M0.2 to M4.5 and 52 stations are utilized (Figure 1). We use 17,955 P-wave and 9,531 S-wave differential arrival times from waveform cross-correlation, in addition to 767,127 P-wave arrival-time differences directly from catalog data and 20,257 absolute catalog P-wave arrival picks. The catalog arrival time differences are selected so that any event is linked to a maximum of 10 neighboring events by at least 8 pair-wise observations (Waldhauser and Ellsworth, 2002). The cross-correlation data are measured using the cross-spectral method for events with similar waveforms, which are only obtained along the fault.

Figure 2a shows the catalog locations, which are scattered along the fault zone. Figure 2b and 2c show the event relocations by using the double-difference location method (Waldhauser & Ellsworth, 2002) and the double-difference tomography method, respectively. We can see that, after relocation, both double-difference methods provide a sharp picture of the seismicity. However, the absolute event locations are different between the two methods. The cluster centroid for the DD location method is located at latitude 37.721o, longitude –122.064o and depth 7.467 km, while it is located at latitude 37.717o, longitude –122.068o and depth 8.646 km for the DD tomography method. That is, there is a shift of about 632 m in the cluster centroid. The differences in the absolute event locations between these two methods are due to the assumption made in the double-difference location method that the catalog-based cluster centroid is constrained. In the double-difference location method, the RMS residuals for cross-correlation data and catalog data decrease from 112 ms to 5 ms and from 167 ms to 40 ms, respectively. For the double-difference tomography method, the RMS residuals for cross-correlation data and catalog data decrease from 112 ms to 5 ms and from 185 ms to 40 ms, respectively (Note that the DD tomography solution includes the 20,257 absolute P-wave arrival times, which is why the starting RMS value is higher). For comparison, we also relocate the events by using the conventional tomography method with only the absolute catalog data (Figure 2d). We see that the event locations are still scattered, comparable to the catalog locations (Figure 2a).

To compare the event locations more clearly, we zoom in on the region of latitude 37.83o to 37.91o and longitude –122.26o to –122.18o. The coordinate center is located at latitude 37.72 o and longitude –122.06 o. The positive z-axis is downward, and the coordinate system is rotated anticlockwise 35 o so that the y-axis is nearly parallel to the strike of the Hayward fault, with y increasing to the NW and x increasing to the NE. Figures 3, 4, 5, and 6 show map views and cross-sections of event locations for catalog data, the double-difference location method, the double-difference tomography method, and the conventional tomography method, respectively. The seismicity inside the box forms a northwest-striking zone associated with the Hayward fault, and that outside the box a diffuse zone of earthquakes about 2 km northeast of the fault zone.

In map view, both double-difference methods collapse the on-fault seismicity to a thin line (Figure 4a and 5a), with most of the events aligning in depth along linear, horizontal streaks (Figure 4d and 5d). Although most of the relative event locations from the two double-difference methods are quite similar, there are some differences between them (Figure 4c and 5c). The event relocations define a nearly planar, vertical fault zone striking in the direction of the surface trace of the Hayward fault in the double-difference locations (Figure 4c). Double-difference tomography locations, however, show a slight NE dip between 4 and 8 km depth and a SW dip between 8 and 12 km depth. We project the event locations in the box into finer slices (1-km separation in the Y-direction) and find that the bend at the depth of 8 km mainly happens from Y=20 km to 23 km. This is consistent with the velocity complexity in the same area (Figure 10). In contrast, the event locations from conventional tomography method have very little improvement compared to the catalog locations (Figure 3 and 6). We notice, however, that the absolute event locations obtained from the double-difference tomography method are slightly different from those corresponding to the double-difference location method. The shift is about 94 m in the x-direction and 110 m in the y-direction (Figure 4b and 5b). The events also tend to be shallower for the DD tomography method with average depth 8.358 km, while it is 8.646 km for the DD location method (Figure 4d and 5d).

Figure 7 and 8 show horizontal slices and across-strike vertical sections, respectively, through the velocity model resulting from the double-difference tomography method. The velocity grids in the x-direction are located at x=-35, -20, -10, -5, -4 -3, 0, 5, 10, 20, and 35 km and in the y-direction at y=-30, -20, -10, 0, 10, 20, and 30 km. In the z-direction, the grids are located at z=0, 4, 8, 12, and 20 km. We see that the Hayward fault is marked by a strong velocity contrast (Figure 7, 8), and this contrast persists vertically beneath the surface trace of the fault to the maximum depths constrained by the model, similar to the results shown in Hole et al. (2001). Higher velocity rocks are located to the west, consistent with the local geological setting. The strong velocity contrast near the surface is due to the boundary between the Franciscan terrane and the Great Valley Sequence (Hole et al., 2001). Figure 9 and 10 show the horizontal slices and across-strike vertical sections through the velocity model from the conventional tomography method with only absolute catalog data. Comparing Figures 7 and 8 with 9 and 10, we note that the velocity contrast boundary is more consistent with the event locations and the velocity contrast is sharper in the model from the double-difference tomography method. Figure 11 shows the velocity difference between the double-difference tomography and the conventional tomography. It clearly shows that the velocities from the double-difference tomography are faster (about 0.2 km/s) than those from the conventional tomography to the west of the fault and slower (about –0.2 km/s) to the east of the fault. This indicates that the double-difference tomography produces a sharper velocity boundary. It also makes the event locations more linear and concentrated (less “fuzzy”) than those from the conventional tomography method. This removes one component of the error in the tomography model due to the “random” mislocations. Consequently, this makes the velocity contrast sharper and more consistent with the event locations.

 

Conclusion

The double-difference tomography method introduced here is efficient in relocating large numbers of earthquakes accurately as well as characterizing the local velocity structure finely. This method collapses the scattered catalog locations into horizontal lineations of seismicity on the northern Hayward fault, the same as by the double-difference location method (Waldhauser et al., 1999). By including the absolute arrival times in addition to the relative arrival times, our method obtains the absolute event locations without the assumptions made in the double-difference location method.

By using the accurate WCC-derived relative arrival times directly, the double-difference tomography method is able to sharpen the image of the velocity structure, especially near the source region. In particular, we obtain a sharper velocity contrast along the Hayward fault, compared to conventional tomography method. This contrast is also more consistent with the event locations and the surface trace of the Hayward fault.

 

Acknowledgements

We are extremely grateful to Felix Waldhauser for sharing the Hayward data and for valuable comments on early drafts of this manuscript. This research has resulted from a project supported in part by the National Science Foundation under grant EAR-0125164. Additional support was provided by Defense Threat Reduction Agency (contract number DSWA01-98-1-0008), US Department of Defense; the content does not necessarily reflect the position or the policy of the US Government, and no official endorsement should be inferred.

 

 

References

Bailey , E. H., W. D. Irwin, and D. L. Jones, Franciscan and related rocks, and their significance in the geology of western California, Bull. Calif. Div. Mines. Geol., 183, 177 pp., 1964

Constable, S. C., R. L. Parker, and C. G. Constable, Occam’s inversion: A peractical algorithm for generating smooth models from electromagnetic sounding data, Geophysics, 52, 289-300, 1987.

Dodge D. A., G. C. Beroza, and W.L. Ellsworth, Foreshock sequence of the 1992 Landers, California, earthquake and its implications for earthquake nucleation, J. Geophys. Res., 100, 9865-9880, 1995.

Dodge, D. A., G. C. Beroza, and W. L. Ellsworth, Detailed observations of California foreshock sequences: Implications for the earthquake initiation process, J. Geophys. Res., 101, 22,371-22,392, 1996.

Fremont, M.-J., and S. D. Malone, High precision relative locations of earthquakes at Mount St. Helens, Washington, J. Geophys. Res., 92, 10,223-10,236, 1987.

Got, J.-L., J. Frechet, and F. W. Klein, Deep fault plane geometry inferred from multiplet relative relocation beneath the south flank of Kilauea, J. Geophy. Res., 99, 15,375-15,386, 1994.

Hole, J. A., T.M. Brocher, S. L. Klemperer, T. Parsons, H. M. Benz, and K.P. Furlong, Three-dimensional seismic velocity structure of the San Francisco Bay area, J. Geophys. Res., 105, 13,859-13,874, 2001.

Lees, J. M., Multiplet analysis at Coso geothermal, Bull. Seism. Soc. Am., 88, 1127-1143, 1998.

McCaughey, M. and S. C. Singh, Simultaneous velocity and interface tomography of normal-incidence and wide-aperture seismic traveltime data, Geophys. J. Int., 131, 87-99, 1997.

Nadeau, R., M. Antolik, P. A. Johnson, W. Foxall, and T.V. McEvilly, Seismological studies at Parkfield III: microearthquake clusters in the study of fault-zone dynamics, Bull. Seism. Soc. Am., 84, 247-163, 1994.

Paige, C. C., and M. A. Saunders, LSQR: Sparse linear equations and least squares problems, ACM Transactions on Mathematical Software, 8/2, 195-209, 1982.

Pavlis, G. L, Appraising earthquake hypocenter location errors: a complete practical approach for single-event locations, Bull. Seism. Soc. Am., 76, 1699-1717, 1986.

Phillips, W. S., H. E. Hartse, and L. K. Steck, Precise relative location of 25 ton chemical explosions at Balapan using IMS stations, Pure Appl. Geophys., 158, 173-192, 2001.

Poupinet, G., W. L. Ellsworth, and J. Frechet, Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California, J. Geophys. Res., 89, 5719-5713, 1984.

Roecker, S. W., The velocity structure of the Pamir-Hindu Kush region: Possible evidence of subducted crust, J. Geophys. Res., 87, 945-949., 1982.

Rowe, C.A., R.C. Aster, W.S. Phillips, R.H. Jones, B.Borchers and M.C. Fehler, Using automated, high-precision repicking to improve delineation of microseismic structures at the Soultz geothermal reservoir, Pure Appl. Geophys. 159, 563-596, 2002.

Rubin, A. M., D. Gillard, and J. L. Got, A reinterpretation of seismicity associated with the January 1983 dike intrusion at Kilauea Volcano, Hawaii, J. Geophys. Res., 103, 10,003-10,015, 1998.

Rubin, A. M., D. Gillard, and J. L. Got, Streaks of microearthquakes along creeping faults, Nature (London), 400, 635-641, 1999.

Sambridge, M.S., Nonlinear arrival time inversion: constraining velocity anomalies by seeking smooth models in 3-D, Geophys. J. Int., 102, 653-677, 1990.

Shearer, P. M., Improving local earthquake locations using the L1 norm and waveform cross correlation: Application to the Whittier Narrows, California, aftershock sequence, J. Geophys. Res., 102, 8269-8283, 1997.

Thurber, C. H., Hypocenter-velocity structure coupling in local earthquake tomography, Phys. Earth Planet. Int., 75, 55-62, 1992.

Thurber, C.H., C. Trabant, F. Haslinger, and R. Hartog, Nuclear explosion locations at the Balapan, Kazakhstan, nuclear test site: the effects of high-precision arrival times and three-dimensional structure, Phys. Earth Planet. Int., 123, 283-301, 2001.

VanDecar, J.C. and R. S. Crosson, Determination of teleseismic relative phase arrival times using multi-channel cross-correlation and least squares, Bull. Seism. Soc. Am., 80, 1548-1560, 1990.

Waldhauser, F., W. L. Ellsworth, and A. Cole, Slip-parallel seismic lineations on the Northern Hayward Fault, California, Geophys. Res. Lett., 26, 3535-3528, 1999.

Waldhauser, F., and W. L. Ellsworth, A double-difference earthquake location algorithm: Method and application to the Northern Hayward Fault, California, Bull. Seism. Soc. Am., 90, 1353-1368, 2000.

Waldhauser, F., hypoDD: A computer program to compute double-difference hypocenter locations, U.S. Geol. Surv. Open File Rep., 01-113, 25 pp., 2001.

Waldhauser, F., and W. L. Ellsworth, Fault structure and mechanics of the Hayward Fault, California, from double-difference earthquake locations, J. Geophys. Res., 107, ESE 3-1-3-15, 2002.

Wolfe, C. J., On the mathematics of using difference operators to relocate earthquakes, Bull. Seism. Soc. Am., in press, 2002.

 

 

Figure Captions

Figure 1. Stations and Hayward fault seismicity between 1984 and 1998 located by the Northern California Seismic Network (NCSN).

Figure 2. Hayward fault event locations. (a) NCSN Catalog locations. (b) Relocations by the double-difference location method (Waldhauser and Ellsworth, 2002). (c) Relocations by the double-difference tomography method. (d) Relocations by the tomography method with only absolute catalog data.

Figure 3. Catalog locations in the region of latitude 37.83o to 37.91o and longitude  to . (a) Map view of event locations in the latitude and longitude system. (b) Map view of event locations in the coordinate system rotated anticolockwise 35o with the center located at latitude 37.72 o and longitude , positive x-axis directing to the northeast, positive y-axis to the northwest, and positive z-axis downward. (c) A-A’ cross section. (d) B-B’ cross section.

Figure 4. Same as Figure 3 but for the event locations relocated by the double-difference location method (Waldhauser and Ellsworth, 2002).

Figure 5. Same as Figure 3 but for the event locations relocated by the double-difference tomography method.

Figure 6. Same as Figure 3 but for the event locations relocated by the tomography method with only absolute catalog data.

Figure 7. Horizontal slices through the 3-D seismic velocity model from the double-difference tomography method. The depth of the slice is labeled above each plot. Black dots indicate the earthquake hypocenters within half the grid-size of the slice.

Figure 8. Across-strike vertical slices through the 3-D velocity model from the double-difference tomography method at locations shown in the lower left corner in each plot. The final hypocenters within 5 km of each slice are included as black dots.

Figure 9. Same as Figure 7 but for the 3-D velocity model from the conventional tomography method with only absolute catalog data.

Figure 10. Same as Figure 8 but for the 3-D velocity model from the conventional tomography method with only absolute catalog data.

Figure 11. Same as Figure 8 but for the 3-D velocity difference between the double-difference tomography and the conventional tomography. The event locations shown in this plot are from the double-difference tomography method.