# Bulletin of the Seismological Society of America

- Copyright © 2006 Bulletin of the Seismological Society of America

## Abstract

Traditional local-earthquake location using a horizontally layered homogeneous velocity model is limited in its resolution and reliability due to the existence of frequently overlooked 3D complexity of the real Earth. During traditional 3D seismic tomography, simultaneous earthquake relocation using the resultant 3D velocity model has produced reliable earthquake locations; however, only a small subset of events are typically used and thus relocated in the inversion. The rest of the events in a catalog must then be relocated using the 3D models. The repeated calculation of travel times across 3D *V _{P}* and

*V*models is also not efficient and not practical for a routine network earthquake location when the very time-consuming exact 3D raytracing is used. Because high-resolution earthquake data are now available from many modern seismic networks, representative high-resolution 3D

_{S}*V*and

_{P}*V*models for a region can be better determined. By taking advantage of recently available high-speed computer technology and large disk space, we implemented a simple algorithm to efficiently locate every local earthquake using the best available regional 3D

_{S}*V*and

_{P}*V*models. Once the

_{S}*V*and

_{P}*V*information for all cubic cells in a 3D grid model are determined,

_{S}*P*and

*S*travel times from each grid point to all seismic stations can be calculated and stored on disk files for later usage. During the iteration process for earthquake location, travel times from a trial hypocenter to all recording stations can be determined simply by a linear interpolation from those of the adjacent eight grid points available in the previously stored disk files without the need for raytracing. The iterations continue until the hypocenter adjustments at the end of the last iteration are below the given criteria and the travel-time residual, or the difference between the observed and the calculated travel times, is a minimum. Therefore, any local earthquake can be efficiently and reliably located using the available 3D velocity models. This simple location program has been applied to relocate earthquakes in the New Madrid Seismic Zone (nmsz) of the central United States and in the central eastern Taiwan region. Preliminary results in both regions reveal that earthquake hypocenters can be efficiently relocated in spite of the very significant lateral structural variations. Tests with data from Taiwan further demonstrate that the resolution of seismic tomography and the relocated seismicity is sensitive to relative distribution of seismic-network stations and background seismicity. Thus, this single-event location program can be applied to relocate all earthquakes in a seismic-network catalog and, more importantly, to allow routine earthquake location for any seismic network using the available 3D velocity models.

## Introduction

Velocity models with homogeneous horizontal layers are routinely used for travel-time calculation in most popular local-earthquake-location programs, for example, HYPO71 (Lee and Lahr, 1975), hypoinverse (Klein, 1978, 2002), and hypoellipse (Lahr, 1999). However, the real Earth structure is far from being horizontally layered, and thus different degrees of mislocation for earthquakes in most seismic-network catalogs are inevitable. The effect of the 3D complexity of the real Earth on earthquake location has long been overlooked at most seismic networks mainly due to insufficient resolution to produce reliable 3D *V _{P}* and

*V*models from the existing data and network configuration. The application of 3D velocity models for single-earthquake location is also difficult due to the lack of an efficient algorithm for fast calculation of travel times using exact 3D raytracing. Because reliable earthquake locations are essential for a comprehensive study of seismogenic and tectonic processes, and the seismic hazard of a region, it has been a continuous and long-term research focus for seismologists to try to improve seismic-network configuration, derive a better velocity model, and develop a better algorithm for earthquake location. Taking advantage of the modern improvements in seismic instrumentation and seismic-network configuration, the determination of reliable 3D

_{S}*V*and

_{P}*V*velocity models for a region becomes a reachable research target. A single-event location algorithm using 3D velocity models can now be implemented.

_{S}A few relative earthquake location techniques have also been developed in modern times to overcome the discrepancies between the given velocity models and the real Earth. Among them, the joint hypocenter determination (jhd) method has been widely applied to relocate groups of clustered earthquakes using arrival-time information (e.g., Pujol, 1988, 2000). The *P*- and *S*-wave station corrections produced during jhd relocation can be used to quantify the lateral variations of velocities in the region (Pujol, 1992). More recently, the double-difference (dd) method has been developed (Waldauser and Ellsworth, 2000) by cross-correlating waveforms to minimize errors in arrival times and by relocating a group of clustered earthquakes using their relative arrival- time differences. Both the jhd and dd methods have successfully provided reliable relative earthquake locations for many regions. However, the jhd and dd methods require clustered earthquakes as an *a priori* condition. A realistic Earth model is not needed, involved, or resolved in the location process. Thus, neither the traditional earthquake- location programs (e.g., HYPO71, hypoinverse, and hypoellipse) nor the relative earthquake-location techniques (jhd and dd) are practical and efficient for a single- event location when a 3D Earth model is involved.

Recently, Hauksson (2000) relocated all earthquakes from 1981 to 1998 in southern California using representative 3D *V _{P}* and

*V*/

_{P}*V*models derived from a 3D inversion of arrival times from selected local earthquakes (4.3% of total events in the catalog) and from a few surface shots. Lomax

_{S}*et al.*(2000) described a probabilistic earthquake- location methodology and introduced an efficient Metropolis–Gibbs, nonlinear, global sampling algorithm to determine local earthquake location over 3D and layered models. Using probabilistic earthquake-location techniques and a 3D

*V*model with a constant

_{P}*V*/

_{P}*V*ratio, Husen

_{S}*et al.*(2003) relocated selected events recorded by the Switzerland seismic network and showed similar results as those obtained from a tomographic inversion. In a routine application for local earthquake location, the travel-time field from a seismic station across the 3D

*V*model was computed and stored on hard disk (e.g., Husen

_{P}*et al.*, 2003). In this study, we have developed independently a simple single-earthquake- location algorithm with similar ideas as Hauksson (2000), Lomax

*et al.*(2000), and Husen

*et al.*(2003) to take care of travel time of seismic waves across complicated 3D

*V*and

_{P}*V*models.

_{S}In a traditional 3D tomographic inversion of travel-time data, a selected subset of earthquakes are simultaneously relocated during inversion using the resultant 3D *V _{P}* and/or

*V*model. However, it is very common that only a small subset of high-quality earthquakes are selected from the entire earthquake catalog for inversion. Therefore, a large number of smaller events, which are essential for a comprehensive study of spatial and temporal distributions of seismicity, are not selected and must be relocated separately after inversion.

_{S}Modern progress in 3D tomographic-inversion techniques and modern advances in seismic networks have offered an opportunity to determine reliable 3D velocity models and earthquake locations (e.g., Shen, 1999; Kim, 2003; H. Chen *et al.*, unpublished manuscript, 2005). Many seismic networks, even the one in southern California, where they could be using Hauksson's 3D models (Hauksson, 2000), are still using a 1D model in their routine earthquake location; and part (but not all) of what is holding them back is how computationally intensive it would be to do all the raytracing during travel-time calculations. In order to reliably and efficiently relocate all earthquakes from a regional earthquake catalog as well as routinely locate every earthquake for any seismic network, it is essential to develop an efficient and stable algorithm to handle travel-time calculation across a 3D Earth model for a single earthquake location, either for previously archived catalogs or for near real- time reports on earthquakes.

## Single-Earthquake Location Using 3D Velocity Models

The single-earthquake-location technique implemented in this study is an extension of the classical Geiger's method (Geiger, 1912). A 3D grid model established during a 3D tomographic inversion is used to represent the local and regional velocity structure for earthquake location. In this 3D grid model, every group of eight adjacent grid points defines a basic cubic cell or block. Inside each basic cubic cell, *V _{P}* and

*V*are constant so that the associated ray paths inside the cube are straight lines.

_{S}*P*- and

*S*-wave travel times and ray information from all seismic stations to all grid points across the 3D

*V*and

_{P}*V*models can be calculated using any available 3D raytracing methods (e.g., Cerveny, 1987; Um and Thurber, 1987; Prothero

_{S}*et al.*, 1988; Podvin and Lecomte, 1991; Thurber and Kissling, 2000) and stored in computer files. Since a trial hypocenter will be located inside one of the cubic cells, travel times and ray information from the trial hypocenter to all recording stations can be easily determined by a linear interpolation from those at the adjacent eight grid points. Thus, repeated travel-time calculations during earthquake location do not require any extensive and very time-consuming 3D raytracing, and any single earthquake within the region can be quickly located using this simple algorithm and 3D velocity models. A brief description of this newly implemented single-earthquake-location algorithm is discussed in subsequent sections.

### Location Method

Let us assume that an event has been recorded by *n* stations and that the arrival time at the *i*th receiver is . Let *x̄ _{j}*,

*j*= 1, 2, 3, and

*[t with macron]*

_{0}represent the actual unknown hypocentral location and origin time of the event, and let

*x*,

_{j}*j*= 1, 2, 3, and

*t*

_{0}be the corresponding initial estimates. Under the assumption that the differences between the actual and initial values are small, we can write 1 where 2 3 and

*T*is the travel time from the intial estimate of the source location to the

^{i}*i*th station. Let us introduce the difference 4 where

*δt*

^{i}is the arrival-time residual 5 Introducing matrix notation, equation (4) becomes 6 where 7 8 9 Note that

**and**

*δ*t**x**are two column vectors; the superscript

^{T}indicates matrix transposition.

In equation (6) the only unknown is **x**. After solving for **x** we use equations (2) and (3) to compute new initial estimates, which in turn gives rise to a new equation (6). This process is repeated iteratively until the residual is small enough. To solve equation (6) we use the least-squares method. The cost function is 10 where **V** is the covariance matrix of the errors. Here it is assumed that **V** is diagonal with *v*_{jj} = 1/*w*_{j}^{2}, where *w*_{j} is the quality weight for the *j*th arrival-time pick. Letting **V** = **C**^{2}, equation (10) becomes 11 where 12 (see, e.g., Draper and Smith, 1981). Minimizing equation (11) with respect to **x** gives the standard least-squares solution 13

Introducing the singular value decomposition (svd) **H**′ = **UΛV**^{T} (e.g., Aki and Richards, 1980), equation (13) becomes 14

Since the initial values used in equation (1) may be far from the actual values, we used Levenberg's (1944) damped least squares, which has the following cost function: 15 where *θ*^{2} is a scalar that controls the length of the solution vector. The damped least-squares solution is given by 16 where **I** is the identity matrix. In terms of the singular value decomposition of **H′**, equation (16) becomes 17

This formulation is also useful because the condition number of **H′** may be very large (i.e., **H′** is ill conditioned), in which case damping decreases it, thus leading to a more stable solution. Damping, however, has the well-known disadvantage that it decreases the resolution of the solution (e.g., Aki and Richards, 1980).

### Centering and Scaling

Centering and scaling are commonly used in statistical regression analysis (e.g., Draper and Smith, 1981) and has been used by Lee and Lahr (1975) and Lienert *et al.* (1986). Centering of the data is accomplished by removing the mean values from ** δt** and from each of the columns of

**H′**. Using the subscript

*c*to indicate centering, we have 18 where the notation < > indicates mean value and

**1**is a vector with all its elements equal to one. Because

**δt′**is a weighted vector (see equation 12) and the weights have been normalized such that their sum is equal to one, <

**δt′**> is a weighted mean. For the element of

**H′**in the

*i*th row and

*j*th column we have 19 and for

*t*

_{0}we have 20 In view of equation (20), the last column of

**H**

_{c}and the last row of

**δt**can be ignored, which means it is sufficient to solve for the hypocentral coordinates only.

We define the scaling matrix **S** as 21 Then we apply the scaling matrix **S** to : 22 After applying the SVD to , we compute **x**_{cs} as follows: 23 where the subscript *cs* stands for the “centering” and “scaling.” We then recover **x** using **x** = **Sx**_{cs}. The origin time is then obtained from equation (20).

### Calculation of Travel Times

Many different methods have been developed and applied to compute travel times for seismic waves. In general, the selection of a method for travel-time calculation depends on the complexity of the velocity model. In this study, the method of Podvin and Lecomte (1991) using the finite- difference technique is selected for travel-time calculation. The method of Podvin and Lecomte (1991) can be applied to a 3D velocity model with very significant lateral and vertical velocity contrasts. The partial hypocenter derivatives can be shown via a variational argument to satisfy 24 where *u _{s}* is the slowness at the earthquake source,

*dl*is an element along the ray path, and

*dx*/

*dl*,

*dy*/

*dl*, and

*dz*/

*dl*are the components of the unit vector tangent to the ray path from the hypocenter and pointing in the direction of ray propagation (Lee and Stewart, 1981). A posterior raytracing (Podvin and Lecomte, 1991) can be used to derive the ray gradients and the ray paths.

Because the finite-difference calculation of travel times through a 3D velocity model is very time-consuming, *P*- and *S*-wave travel-time information from every grid point in a 3D grid model covering the study region to every surface station (receiver) can be first calculated and then saved on disk files for later usage. Note that the travel time for a seismic wave from a source to a receiver is the same as that from the receiver to the source. The velocity model and receiver locations will remain the same during earthquake location, and so will the travel-time parameters from all receivers to every grid point. Therefore, travel-time information for any ray path, that is, from a trial hypocenter to a recording station, can be determined efficiently by interpolation from those of the adjacent rays hitting the neighboring grid points, whose travel times have previously been calculated and stored on disk. In this study, a linear interpolation of eight grid points of a cube that contains the trial hypocenter is used for the arrival time and partial derivatives estimation. The simple linear interpolation to compute the arrival time and partial derivatives with regard to the unknowns has been widely used with the finite-difference calculation of travel times (e.g., Benz *et al.*, 1996). When the grid size is carefully chosen, numerical experiments show the inversion is very consistent. By doing so, the computation time required in travel-time calculation for a 3D velocity model during the iteration process of single-earthquake location can be significantly reduced.

To estimate the uncertainty of the hypocenter parameters, we follow this approach to minimize root mean square (rms) travel-time residual. The uncertainty in origin time, longitude, latitude, and depth of a hypocenter can be estimated by perturbing independently each of the hypocenter parameters. When a parameter is perturbed (e.g., *t* → *t* + Δ*t*) a new rms travel-time residual is calculated with respect to the perturbed hypocenter parameters (e.g., *t* + Δ*t*, *x*, *y*, *z*). The perturbation of a parameter is determined such that the differences between the new rms travel-time residual and the final rms travel-time residual reaches ±20%. Assuming the uncertainty of the origin time and earthquake location follows a normal distribution, the minimum amount of the positive and negative perturbation is chosen as the representative uncertainty of the associated parameter.

## Test Examples from the nmsz in the Central United States

Local-earthquake data recorded by three-component seismic stations in the nmsz of the central United States are selected to test the newly developed single-event location algorithm using recently available 3D *V _{P}* and

*V*models. These data include those recorded by the panda (Portable Array for Numerical Data Acquisition) experiment between 1989 and 1992 (Chiu

_{S}*et al.*, 1992) and by the regional seismic network between 1995 and 2000 (Chiu

*et al.*, 2001). Figure 1 shows the location of panda and regional seismic- network stations. All earthquakes were located preliminarily using the panda velocity model, which consists of a thin layer of sediments (650 m) with extremely low

*V*and

_{P}*V*overlying five crustal homogeneous layers and a half-space upper mantle (Chiu

_{S}*et al.*, 1992). Figure 2 shows the epicenters of the locations obtained using the panda model and the hypoellipse program (Lahr, 1999).

The above-mentioned nmsz earthquakes with more than five reliable *P* and *S* picks were selected for a 3D tomographic velocity inversion using the program of Benz *et al.* (1996), which was revised by P. Shen (personal comm., 1999) to independently invert for *V _{P}* and

*V*simultaneously. The selected earthquakes are simultaneously relocated after velocity inversion (H. Chen

_{S}*et al.*, unpublished manuscript, 2005) and the relocated epicenters are shown in Figure 3. All earthquakes from the original database are then relocated by the newly developed single-event relocation algorithm discussed in this article using the resultant 3D

*V*and

_{P}*V*models of H. Chen

_{S}*et al.*, (unpublished manuscript, 2005). The relocated epicenters determined with the new location program are shown in Figure 4. While the distribution patterns of epicenters shown in Figures 2–4 seem to be very similar, there are visible differences in depth views as demonstrated in a few along-strike projections of hypocenters (Figs. 5, 6, and 7). Along the most active central segment of the nmsz, only randomly selected hypocenters are shown in the along-strike projection (Fig. 6) to make it possible to visualize the hypocenter shift for each individual event. Basically, the relocated hypocenters from the single earthquake location using 3D velocity models are consistent with those obtained from the simultaneous velocity inversion and earthquake relocation. However, the depths of some relocated earthquakes shift visibly from those of the preliminary hypocenters using the horizontally layered panda model, particularly in Figures 5 and 7.

The uncertainties in the earthquake locations for the nmsz along the *x* axis (longitude), *y* axis (latitude), and *z* axis (depth) are estimated and shown as a function of earthquake number in chronological order in the database (Fig. 8). It is apparent that location uncertainties are, in general, very small, with an average of 0.24 km, 0.30 km, and 0.48 km for erx, ery, and erz, respectively. Location uncertainties are significantly reduced, mostly below the average value, after the time of earthquake number 500 when the expansion, upgrade, and densification of the modern regional seismic network in the nmsz was completed. Thus, the shifts in hypocenters shown in Figures 5, 6, and 7 are most likely the consequence of the application of the 3D velocity models in earthquake relocation.

The 1D panda velocity model was built upon the information of well-constrained layer thickness and layer seismic velocity obtained from a few deep and many shallow seismic reflection/refraction lines and from a 1D *V _{P}* and

*V*velocity inversion using panda data (Chiu

_{S}*et al.*, 1992). Chiu

*et al.*(1997) showed in terms of statistical results (i.e., rms, vertical error erz, and horizontal error erh) that the panda model has provided a very significant improvement in hypocenter determination over other previous models (e.g., Nuttli

*et al.*, 1969; Mooney

*et al.*, 1983; Andrews

*et al.*, 1985; Nicholson

*et al.*, 1984). For the first time, images of active faults in the nmsz could be depicted from the better relocated hypocenters (Chiu

*et al.*, 1992, 1997). The homogeneous-layer velocity model beneath the Upper Mississippi Embayment is an appropriate first-order approximation. However, the layer model fails to consider the lateral variations of thickness and seismic velocities inside the crustal layers, especially in the uppermost sedimentary basin. From a jhd analysis of panda data in the nmsz region, Pujol

*et al.*(1997) demonstrated that

*P*- and

*S*-wave station corrections correlate exceptionally well with the lateral variations of the thickness of the thin sediments beneath the stations, revealing that the thin sedimentary basin has a significant impact on earthquake location.

In the Upper Mississippi Embayment, the thickness of the sediments gradually increases from 0 at the surrounding Paleozoic outcrop boundary to about 1000 m thickness near Memphis (e.g., Chen *et al.*, 1996). Although the thickness of sediments is relatively thin compared to that of all other crustal layers, the travel times of both *P* and *S* waves inside the sediments cannot be overlooked because of its extremely low seismic velocity (*V*_{P} ∼ 1.8 km/sec and *V*_{S} ∼ 0.6 km/ sec) and the large velocity contrast with the underlying Paleozoic basement (*V*_{P} ∼ 6.0 km/sec and *V*_{S} ∼ 3.6 km/sec). Independent of the depth and location of earthquakes, all seismic ray paths will impinge almost vertically at seismic stations in the nmsz due to the extremely low-velocity sediments and the high-velocity contrast across the bottom of the sediments. For example, a misestimation of 100 m of sediment thickness in the velocity model will introduce ∼0.16 sec of travel-time residuals for the *S* wave, which will have a significant effect, particularly on earthquake depth. As an interim solution to accommodate the lateral variation of sediment thickness beneath the seismic stations, Chiu *et al.* (2001) adapted 10 velocity models of different thickness for the uppermost sediment layer to relocate earthquakes in the nmsz region.

The thickness of the sediments varies from 200 to 400 m beneath most seismic stations in the northern Mississippi Embayment. Therefore, the thickness of the low-velocity sediments was overestimated (Δ*h* ≈ 250–450 m) in the panda model for the northern nmsz. The hypocenters beneath the northern nmsz are expected to be located shallower than what they should be when the panda model is used. An along-strike cross-sectional view of hypocenters in the northern nmsz (Fig. 5) shows that the relocated hypocenters are deeper in most cases, as expected. The various horizontal and vertical shifts of the relocated hypocenters between different events are most probably a reflection of the localized discrepancies between the panda model and the more realistic 3D *V*_{P} and *V*_{S} models (H. Chen *et al.*, unpublished manuscript, 2005). Therefore, hypocenters, especially their depths, can be improved either by a simultaneous tomographic inversion and relocation (H. Chen *et al.*, unpublished manuscript, 2005) or by single earthquake relocation using 3D velocity models (this study). The thickness of the sediments is roughly 650 m in the central nmsz, that is, the sediment is properly represented in the panda model. Thus the relocated hypocenters as shown in Figure 6 are very similar to the original locations using the panda model, as expected. On the contrary, the thickness of the sediments ranges from 700 to 1000 m underneath those stations in the southern nmsz, that is, the sediment layer is underestimated (Δ*h* ≈ 50–250 m) in the panda model. The relocated hypocenters shown in Figure 7 for a short section of the southern nmsz are, in general, shifted toward shallower depths, again as expected. Therefore, the relocated hypocenters (Figs. 5, 6, and 7) are shifted properly, with only a few exceptions, from the original panda locations to reflect the deficiencies of the oversimplified homogeneous- layer panda model, which does not represent the sedimentary basin correctly.

The evolution of earthquake location in the nmsz can be briefly summarized in the following five additional transverse cross-sectional views across the northeast, northwest, north-central, south-central, and southwest segments of the seismicity (Fig. 9). The geometry of the active faults in the nmsz can be depicted from the better relocated earthquake locations presented in this article. The active faults in the nmsz consist of the vertical northeast, northwest, and southwest segments and the gently southwest dipping central segment. There are no apparent improvements in the earthquake clusters shown in the AA′ and BB′ sections even with the 3D models (Fig. 9b). This is probably because of very few local earthquakes along the northeast and northwest segments of the nmsz and adjacent areas in the northern margin of the nmsz seismic network. Therefore, there is no sufficient coverage of seismic ray paths to produce high-resolution 3D velocity models beneath the northeast and northwest segments. It is, however, apparent that there are significant lateral variations of fault-zone geometry from north to south in the central segment of the nmsz, where seismic-network coverage is excellent and seismicity is the highest in the entire region. In order to apply the single-event location program presented in this study for routine earthquake location in the nmsz seismic network, additional spatial coverage by seismic stations in the northeast and northwest segments of the nmsz will be needed to improve the resolution of 3D velocity models.

## Test Example from the Hualien Area in Central Eastern Taiwan

The single-earthquake 3D location program has also been applied to the Taiwan region to relocate all earthquakes recorded from 1991 to 2002 after representative 3D *V _{P}* and

*V*models were determined (Kim, 2003). In particular, we focus on the Hualien region in central eastern Taiwan, where data from the island-wide seismic network and from a dense temporary local seismic array are available. This test aims to investigate the significance of seismic-network configurations on the resultant 3D velocity tomography and earthquake relocation. The Hualien area is chosen for the test because it is the most seismically active area of the entire region. However, the spatial correlation between the tectonic structure and seismicity beneath the Hualien area are poorly known due to complicated subsurface structures and lack of reliable earthquake locations.

_{S}Since 1991, a modernized island-wide seismic network of 78 three-component stations operated by the Central Weather Bureau (cwb) is responsible for earthquake monitoring in the Taiwan region. In addition to a few cwb stations distributed in the Hualien area, mostly along the north– south-trending tectonic structural belts, a 30-station panda II array was deployed in the area (Fig. 10) from 1993 to 1995 (Chen, 1995a). Selected data collected by the island- wide seismic network and the temporary seismic arrays were used to determine high-resolution 3D *V _{P}* and

*V*models for the entire Taiwan region (Kim, 2003; Kim

_{S}*et al.*, 2005).

All earthquakes in the cwb catalog from 1991 to 2002 were originally located by the island-wide seismic network using a horizontally layered velocity model simplified from a 3D velocity model of Chen (1995b). Chen's model consists of a 3D *V _{P}* model and a constant

*V*/

_{P}*V*ratio determined using pre-cwb seismic-network data recorded mostly before 1991. These earthquakes were relocated by the single-event location algorithm presented in this study using the 3D

_{S}*V*and

_{P}*V*models of Kim (2003). The original and relocated cwb catalog data are projected into two transverse cross- sectional views in the Hualien area (Fig. 11). Because of data and methodology differences in the 3D tomographic inversion, the spatial resolution of the 3D

_{S}*V*and

_{P}*V*models of Kim (2003) is superior to that of Chen (1995b) in three major aspects. These aspects include the following: (1) the number of seismic stations covering the same area is more than doubled and all data are three-component in the study of Kim (2003) as compared to that used in Y. Chen (1995), which were mostly single-component; (2) both

_{S}*V*and

_{P}*V*were determined simultaneously and independently in Kim (2003) as compared to

_{S}*V*only in Chen (1995b); and (3) local dense seismic-array data were used in the study of Kim (2003), which significantly increases the spatial resolution of the resultant velocity models, particularly beneath the Hualien region. Because only a few cwb stations are located in the Hualien area, the spatial coverage in the Hualien area is basically poor for local-earthquake location.

_{P}

For this test local-earthquake data were selected from the cwb catalog for the period from 1993 to 1995, that is, the same time period when the panda ii array was deployed. The selected earthquakes were relocated using the single-earthquake-location algorithm presented in this study, the 3D *V _{P}* and

*V*models of Kim (2003), and cwb seismic stations. The hypocenters from the original cwb catalog and from the 3D relocation using only cwb stations are projected into two transverse cross-sectional views shown in Figures 11a and 11b, respectively. In spite of the significant improvement of the velocity models and location technique, it is difficult to argue that there is an obvious improvement in hypocenter locations from Figure 11a to Figure 11b when only a few cwb network stations are used. The only differences are that the relocated hypocenters shown in Figure 11b seem to be slightly more clustered and that the clustered seismicity extends to larger depths than in Figure 11a.

_{S}While the panda ii array was deployed from 1993 to 1995 in the Hualien region, it recorded two to three times more earthquakes than the island-wide cwb network during the same time period. Most local earthquakes were very small and either not detected or not locatable by the cwb network alone. Local earthquakes recorded by the panda ii were located by the best available 1D layered model for central Taiwan (Yeh and Tsai, 1981) and are shown in Figure 11c along the two cross-sections used for Figures 11a and 11b. These local earthquakes have been relocated using the algorithm presented in this study and the 3D *V _{P}* and

*V*models of Kim (2003), and using cwb and panda ii stations. Therefore, these local earthquakes have been relocated using the best ever seismic-network configuration and the best possible 3D

_{S}*V*and

_{P}*V*models available for the region. The westward-dipping planar seismicity from the surface to a depth of ∼30 km is slightly more clustered and extends deeper in Figure 11d than in Figure 11c, depicting the geometry of an active fault. This planar seismicity marks the eastern-boundary fault separating the Central Mountain Range of the Eurasian plate from the Coastal Range of the Philippine Sea plate (Kim 2003; Kim

_{S}*et al.*, 2005). From the relocated local seismicity, the width of this boundary fault has further been reduced to a minimum and its geometry has been better defined (Fig. 11d) than before (Figs. 11a and 11b). More than half of the horizontal ground velocity between the converging Eurasian plate and Philippine Sea plate has been accommodated along this westerly dipping fault (Kim

*et al.*, 2005). Thus, the seismic-network configuration, the 3D

*V*and

_{P}*V*models, and an efficient and stable single- event-location algorithm constitute three of the most essential components for efficient and reliable earthquake location.

_{S}## Discussion

Although the Earth is far from being made by homogeneous layers, it is a common practice to use a best- determined 1D layered velocity model to simulate Earth's velocity structure for earthquake location. This 1D model is usually determined based on geological and geophysical data and is capable of providing the best results with minimum statistical uncertainties on travel-time residuals and hypocenter error estimations. However, earthquake hypocenters can be mislocated, particularly in a region of complicated subsurface structure with very significant lateral structural variations. Seismologists are continuously trying to improve the representative velocity models for a region to obtain the best possible earthquake locations. In the test example of the nmsz, we have demonstrated that local earthquake location can be significantly improved when the lateral variations of crustal velocity structure, especially the sedimentary basin, are properly taken into account. However, the 3D velocity structure is not resolved evenly over the coverage area of the seismic network, so the quality of earthquake location is not even. Hypocenters are expected to be determined not as well in the area of poor structural resolution (e.g., near the margin of the network such as the northeast and northwest segments of the nmsz) as in the other areas in the center of the network (Fig. 5).

As a rule of thumb in traditional earthquake location, the quality of earthquake location is poor for earthquakes located outside or near the margin of a seismic network. However, Chiu *et al.* (1997) presented a test case in the nmsz using panda data and concluded that an earthquake outside of a seismic network can still be located reasonably well only if the local and regional crustal velocity models are well determined, even though statistical error estimations for hypocenter location and origin time may indicate otherwise. Thus, the single-earthquake-location technique presented in this article should be applicable not only to local earthquakes but also to regional earthquakes near a seismic network as long as representative 3D *V*_{P} and *V*_{S} models for the region are available.

The resolution of a 3D velocity model for a region can be improved, for example, by increasing the spatial coverage of a seismic network, by using a smaller grid size in velocity model, and by improving algorithms for 3D tomographic inversion and for 3D raytracing. Although 3D velocity models are available for many areas (e.g., Hauksson, 2000; Lomax *et al.*, 2000; Husen *et al.*, 2003), the majority of modern seismic networks in the United States and around the world still depend on a horizontally layered 1D velocity model for routine earthquake location. A layered model is always a good first-order approximation. However, modern improvements in seismic instrumentation and an increasing number of seismic stations in most seismic networks have resulted in the collection of earthquake data with high spatial resolution and wide spatial coverage. Such improvements are essential to provide data adequate for the determination of representative regional 3D velocity models. Any future improvement in regional 3D velocity models or raytracing techniques can be easily adapted and applied to the single- earthquake-location algorithm presented in this article. The test example from the Hualien area of central eastern Taiwan demonstrates further that the spatial resolution of 3D *V*_{P} and *V*_{S} models and earthquake relocation can be significantly improved by the addition of a dense local seismic array/ network.

In addition, slow computer speed and lack of sufficient disk space have previously prohibited an effective application of a 3D raytracing technique on repeated travel-time calculations across 3D velocity models during earthquake location. Modern advances in computer technology have dramatically improved the speed of computing travel times across 3D models. Computation times required for travel- time calculation using 3D raytracing (Podvin and Lecomte, 1991) and using the simple method of searching disk files proposed in this study are estimated for a region of 230 km × 150 km × 17 km. The region was modeled in terms of 3D models with cubic blocks having sides equal to 0.25, 0.5, 0.75, and 1 km to test the time efficiency and dependence on block size of each technique. The computations were carried out using a pc computer running under Linux with a Pentium IV 2.4 GHz cpu and 1 gb of memory. As summarized in Table 1 and presented in Figure 12, the time required for one travel-time calculation using the 3D raytracing technique of Podvin and Lecomte (1991) is inversely proportional to the block size, that is, it varies from ∼2 sec to ∼83 sec for block sizes of 1 km and 0.25 km, respectively. However, the time required for 1000 travel-time calculations using the disk-searching method proposed in this study is almost a constant, ∼0.003 sec, in spite of different block sizes. The computation times required for both methods may vary significantly when different 3D raytracing methods, different block sizes, or different computers are used. Nevertheless, the above tests demonstrate that travel-time calculations using 3D raytracing can be very fast in the modern computational environment, and that the method proposed here is extremely efficient, with an increase of about 10^{6} times in computational speed.

Recent progress in disk-storage size and faster access time to disk files have made the implementation of the proposed single-earthquake-location algorithm possible. Providing accurate 3D *P*- and *S*-wave velocity models for a study region are available, the simple algorithm presented in this article can be efficiently applied for local and regional earthquake location. Most importantly, this single-earthquake-location algorithm using 3D models can be easily adapted by any seismic network for routine earthquake location to produce a high-quality earthquake catalog, which previously seemed to be an impossible goal for any seismic network around the world. With reliable earthquake locations from a routine network operation, it is possible to quickly correlate seismicity with active faults, investigate characteristic features of active faults and their associated seismic hazard, study spatial and temporal variations of seismicity and their implication to precursory studies, and explore regional structural models for tectonics studies.

## Conclusions

An earthquake can be easily mislocated in the range from a few to a few tens of kilometers if the lateral and vertical velocity variations in the Earth are not properly taken into account in the travel-time calculations. Reliable 3D *V _{P}* and

*V*models, an efficient 3D raytracing technique for travel-time calculation, and a stable algorithm for single- earthquake location using 3D models are essential for the improvement of local and regional earthquake location. Modern progress in seismic instrumentation and improvement in seismic-network configuration have allowed the collection of excellent earthquake data that can be used to significantly improve the spatial resolution of 3D

_{S}*V*and

_{P}*V*models for a region. The recent advances in high-speed computer technology and large capacity of storage media permit the development of a simple algorithm for a fast, efficient, and practical travel-time calculation using 3D velocity models. Test examples from the nmsz in the central United States and in the Hualien area in central eastern Taiwan demonstrate that local and regional earthquakes can be reliably and efficiently located by the single-earthquake-location algorithm using 3D

_{S}*V*and

_{P}*V*models presented in this study. In both regions, the relocated hypocenters are more clustered than the original hypocenters and depict the geometry of active faults that is essential for seismic-hazard assessment and regional tectonic studies. Any future improvement in seismic-network configuration or raytracing technique will improve the resolution of 3D

_{S}*V*and

_{P}*V*models as well as the earthquake locations. This simple algorithm for single- earthquake location using 3D velocity models can be easily adapted by any seismic network for routine earthquake location to produce a high-quality earthquake catalog if high- resolution 3D velocity models are available.

_{S}## Acknowledgments

We thank Dr. Harley M. Benz of the U.S. Geological Survey for allowing us to use his 3D tomographic-inversion software in this and related studies. Efforts of Shen (1999) on modifying this 3D tomographic-inversion software to invert simultaneously *P*- and *S*-wave velocity models are highly appreciated. Critical reviews by Dr. Fred Klein of the U.S. Geological Survey and an anonymous reviewer have significantly improved the content of this manuscript. Sincere thanks also go to Dr. Jeanne Hardebeck of the U.S. Geological Survey, whose comments have made this a better article and are highly appreciated. This study was sponsored by U.S. Geological Survey contracts 98HQGR1012 and 99HQGR0053 and by the Center of Excellent program at the Center for Earthquake Research and Information (ceri) of the University of Memphis. All figures in this article were generated using the Generic Mapping Tools software (Wessel and Smith, 1991, 1995), which really makes a big difference in the high-quality presentation of our results. This is ceri contribution number 480.

- Manuscript received 1 June 2004.