# Bulletin of the Seismological Society of America

- Copyright © 2006 Bulletin of the Seismological Society of America

## Abstract

I prove three theorems on earthquake location in the laterally homogeneous, radially concentric layered earth showing that (1) epicenters of earthquakes are locatable without using numerical travel-time information, (2) focal-depth residuals change sign at some epicentral distance that is a function of the focal depth but not of the distance residuals, and (3) the distance Δ between two random points on the surface of a sphere is distributed as sinΔ. Theorem (1) implies that the epicenter can be located without simultaneously solving for hypocentral depth or origin time. Theorem (2) suggests that the focal depth should never be determined by joint regression on all four hypocentral parameters. Theorem (3) suggests a likely reason why late readings tend to overpower accurate readings when a standard least-squares technique is used. Together these theorems suggest that the problem of earthquake location is ill posed.

Joint least-square regression will lead to the true hypocenter only if well-posedness is restored by restricting the class of admissible solutions using *a priori* knowledge, such as a travel-time table. But if the travel-time table is derived from the hypocentral solutions there is a feedback between location errors and errors in the travel-time tables, and local minima cannot be eliminated. Gutenberg instructively attempted to correct for bias due to azimuthal clustering, by assuming that residuals were nonnegative vectors.

## Introduction

Large mislocations of earthquake hypocenters are not uncommon. The published confidence bounds of earthquake location errors are widely believed to be unrealistic. However, an objective comparison between different location techniques, or a calibration of errors arising from various approaches, is hampered by the scarcity of standard measurements from controlled sources.

Jeffreys (1970, p. 490) pointed out that “the standard error is a minimum estimate of uncertainty.” Mislocations are common even in large, well-recorded events. For example, the 1960 Chile earthquake (*M*_{w} 9.5) was found to be mislocated by about 75 km landward after its relocation by seismic-imaging techniques (Krawczyk *et al.*, 2003). Menke and Schaff (2004) tested routine location procedures on simulated data and concluded that the errors were as large as or larger than those found from relative locations for pairs of earthquakes.

Current procedures of earthquake location may be self- consistent without being necessarily accurate, because they are not checked against data from controlled sources. In this article we re-examine some of the fundamental assumptions that underlie location procedures. We show that the problem of locating a point inside a radially symmetrical sphere from observations on the surface is an ill-posed problem in four dimensions (three spatial dimensions plus origin time). Linear regression in the time domain may be useful as an approximation, especially when the flat-earth approximation may be assumed to be valid. But in the general case the errors in the hypocentral depth may contaminate the epicentral solution.

Many seismologists are aware of these problems, yet they continue to use the present location methods because they yield self-consistent solutions. This consistency is due to bootstrap procedures, as the hypocentral solutions are fed back into the travel-time tables. Calibration of networks by large controlled explosions is still rarely performed.

One may wonder why we tolerate such admittedly unsatisfactory procedures. The answer, as Sir Harold Jeffreys once explained to me, is that the location problem has been seen as a subsidiary problem to the construction of travel- time tables. Errors in the location of individual events were considered of negligible importance, as “it will all come out in the wash.” The underlying—but unproved—assumption is that the errors are normally distributed and unbiased.

The statistical theory of earthquake observations is not highly developed. For example, the fact that the focal-depth component of the travel-time residuals changes sign in the near field seems to have been overlooked (see theorem 2). Generations of station seismologists faced the thankless chore of trying to pick an emerging signal in the presence of noise, yet the fact that this operation is inherently biased was ignored. When I once attempted to point this out, Jeffreys was upset. He claimed that no such bias existed because seismogram readers would attempt to overcompensate for late readings by picking some spurious early noise for the true signal. I have never found a seismologist who actually did that.

Without a reliable set of controlled sources observed over a wide range of distances, bias is hard to prove. The problem has not been widely studied, but it seems reasonable to expect that a *P* arrival will be easier to pick if the frequency of the signal exceeds the frequency of the noise. As distance of propagation selectively damps out the higher frequencies, distant arrivals are normally correlated with lower frequencies and larger positive errors.

In Jeffreys (1970), the heading “Determination of epicenters” is immediately followed by the statement “This is done by the method of least squares” (p. 71). Thus it has been taken for granted (starting with Geiger, 1910) that all four hypocentral parameters (the latitude *φ*, the longitude *λ*, the focal depth *h*, and the origin time *t*) are independent and normally distributed random variables. Otherwise the method of least squares will not work. But Jeffreys (1970) cautiously added that “it is usually desirable to restrict the solution to the best stations” (p. 71) or to use weighted least squares.

## Theorem (1), on Locating Seismic Events from Surface Observations on a Sphere

Let *E* be the epicentral area of some event *O* at depth *h* below the surface of a radially layered sphere. Let *t*_{i}, *i* = 1, …, *n* be the measured arrival times at *n* stations located on the surface of the sphere. Then *E* is locatable to any desired accuracy without using additional travel-time information. However, the focal depth *h* remains indeterminate.

*Proof*

(see Anderson, 1978; Lomnitz, 1982). Suppose, without loss of generality, that the *n* arrivals are ordered by increasing lateness, so that *t*_{1} is the earliest arrival time and *t*_{i} < *t*_{i+1} for all *i.* Let *M*_{i} be the great circle that bisects the arc [*P*_{i}, *P*_{i+1}] between any sequential pair of stations (Fig. 1). Let *H*_{i} be the left hemisphere bounded by *M*_{i}, that is, the hemisphere that contains all stations *P*_{1}, *P*_{2}, …, *P*_{i} and therefore all surface points with arrival times earlier than *t _{i}*.

*H*must therefore contain the epicentral area

_{i}*E*, which is the area of minimum travel times on the surface of the earth. Therefore

*E ⊂ H*.

_{i}

We now repeat this operation for all *n* − *1* sequential pairs of stations. In every case we obtain a left hemisphere *H* that contains the epicenter. Note that the actual observations *t* do not come into the procedure, as only the order of arrivals is used. Now consider the intersection *Q* = *H*_{1} ∩ *H*_{2} ∩ … ∩ *H*_{n−1} of all left hemispheres. This polygonal area on the surface of the sphere has the property that 1 using the basic theorem of set theory, which says that if *A ⊂ B*, then *A* ∩ *B* = *A*. Thus the intersection of all left hemispheres contains the epicenter. The size of the convex polygon *E* depends on the number and configuration of stations. The larger the number of observations, the smaller is the probable size of *E*. Thus, by increasing the density of stations we may locate the epicenter to any desired accuracy without using the numerical values {*t _{i}*} of the arrival-time observations. Only the rank sequence is used.

*Corollary 1.*

Note, however, that the focal depth *h* cannot be located in this manner. For the intersection *Q* between meridian planes *M _{i}* is a pyramidal volume bounded by the epicentral polygon

*E*, with its vertex at the center of the sphere. The focus must be somewhere within this volume, but its depth

*h*cannot be ascertained, because all

*M*go through the center of the sphere. This indeterminacy may be attributed to the fact that seismic observations are restrained to the surface. Thus no matter how small the area of the epicentral polygon, the focal depth remains indeterminate.

_{i}The efficiency of this location procedure may be slightly improved if instead of bisecting the arcs [*P _{i}*,

*P*

_{i+1}] we partition them in proportion to the differences in the travel times

*T*−

_{i}*T*

_{i+1}. But even this simple idea involves serious problems. The travel time is a positive increasing function of the epicentral distance Δ, but we don't know how to estimate the intercept at the origin.

*Corollary 2.*

In principle, one might estimate a provisional travel-time curve *T*(Δ) by fitting the arrival-time differences to the epicentral distance Δ. Next we may try to find the velocity structure of the interior by inverting *T*(Δ). In this way, one should be able to obtain increasingly accurate estimates of the focal depth *h*. But such estimates will be affected by errors in measuring the arrival times and in estimating the travel-time curve—particularly, the zero-distance intercept—in the presence of lateral heterogeneity. The estimation of the epicentral parameters (*φ, λ*) might remain reasonably free of these errors as shown by theorem (1), as long as we keep it separate from the procedure of focal-depth estimation.

## Theorem (2), on Estimating the Focal Depth

Let *E* be the epicenter of an event *O*, determined as in theorem (1). Consider the focal depth *h* below *E* in the radially layered sphere. Let Δ_{i}, *i* = 1, … *n* be the measured epicentral distances at *n* stations located on the surface of the sphere. Now consider the travel-time error *τ*(Δ, *δh*) associated with an error *δh* in focal depth. Notice that this error changes sign at some critical epicentral distance Δ_{c} (Fig. 2). The area within Δ < Δ_{c} may be called the epicentral cap (or the near field) for a focal depth *h*. When observations exist both inside and outside the epicentral cap the error in focal depth will contaminate the epicentral location if joint- regression procedures are used.

*Proof.*

Assume that the epicenter *E* is known exactly, but that the focal depth *h* is unknown. From symmetry the travel- time error *τ* associated with an error *δh* in focal depth depends only on the epicentral distance Δ (Fig. 2). We define the critical distance Δ_{c} as the epicentral distance where a ray through the equator of the focal sphere intersects the surface of the earth. Since the error *δh* is normal to this ray the path to Δ_{c} changes by an amount that is second-order in *δh*. In other words, observations made at an epicentral distance Δ_{c} have a negligible effect on the focal-depth determination.

In the near field (Δ < Δ_{c}) all arrival times have the property that *τ* > 0, where *τ* is the travel-time error corresponding to an error *δh* in focal depth, because all rays exit the focal sphere through the upper hemisphere. Conversely, in the far field all the rays will exit the focal sphere through the lower hemisphere and we have *τ* < 0 (Fig. 2). This means that a given travel-time residual *τ* should be interpreted in opposite ways, depending on whether the observation is made inside or outside the epicentral cap.

But the critical distance Δ_{c} depends on the focal depth *h*, which is unknown. Suppose that *h* = 30 km; then Δ_{c} might be somewhere around 75 km. But if *h* = 100 km, Δ_{c} could be around 300 km. If the focal depth is unknown all stations at distances between 75 km and 300 km would be ambiguous in terms of their influence on focal-depth determination. During joint regression, the sign of their contribution would depend on the result of the location procedure. Yet these stations could well be crucial for locating the epicenter. Two stations with identical residuals located on either side of Δ_{c} will cancel out in a joint-regression procedure.

Suppose that the location of the epicenter *E* is known from a separate procedure such as described in theorem (1). Let *n* and *m* be the numbers of observations in the near field and the far field, respectively. Without loss of generality assume *n* < *m*. As any pair of residuals *ρ* on either side of Δ_{c} will have opposite effects on the estimation of the focal depth, their joint effect may be reduced to a remainder *ρ*_{inside} − *ρ*_{outside} at the near-field station.

Consider now the problem of joint hypocenter determination, where the epicenter *E*, the focal depth *h*, and the origin time *t* are located jointly by multiple regression on all four independent variables. There are *m*!/(*m* − *n*)! possible configurations of pairs of stations, which means that the problem is nonunique since any remainder will contaminate the epicentral solution and every configuration will produce a different remainder. This may be explained as follows.

Suppose that the location of the epicenter *E* is unknown and we attempt to determine the epicenter and the focal depth by joint regression in four-dimensional parameter space. Thus the travel-time residuals *ρ* will be partitioned between four normal equations. For any given epicentral solution the pairs of focal depth residuals *δh* that straddle the boundary between the near field and the far field tend to cancel out, and the leftover residuals will be split between an epicentral correction and a focal-depth correction. But there is no unique way of carrying out such a partition.

For every station configuration there will be a functional dependence of the epicentral parameters on the focal depth. In general there will be several local minima, and the solution will converge to a local minimum in the neighborhood of the trial epicenter selected at the beginning of the iterative procedure. A combination of positive residuals in the far field with few negative residuals in the near field may even pull the focus upward into the air (i.e., lead to a negative estimate of *h*). On the other hand, where there are many near- field stations, as in California, the focal depths may be systematically overestimated.

## Theorem (3), on the Distribution of Random Observations on a Sphere

The distance Δ between random points on the surface of a sphere is distributed as sinΔ.

*Proof*

(Fig. 3). Unlike the preceding theorems, this theorem does not require an assumption of radial symmetry, or any other special property of the internal structure of the Earth. Consider two arbitrary points *O*, *P* on the surface of the sphere. Let Δ be the angular distance *OP*. The shaded annular area of width *d*Δ at a distance Δ equals 2*πR*sinΔ*d*Δ. Then the probability of a random point *P* falling within the shaded area is proportional to this area, normalized over the area of the sphere (Fig. 3). Hence the probability of a point *P* to be found at a distance between Δ and Δ + *d*Δ from point *O* is 2 and 3 where *F* is the cumulative probability distribution of the angular distance *x*. By differentiation, the probability density is found to be *f*(*x*) = 0.5 sin*x*.

*Corollary 1.*

The epicentral location of earthquakes ought not to depend on the location of seismic observations (though a case might be made for the converse to be true!). Yet the actual sample of arrival-time measurements is always likely to be dominated by distant observations.

For any given earthquake of epicenter *O* the probability of a station *P* being located at some distance Δ looks like Figure 4. This may be easily confirmed by plotting the actual distribution of epicentral distances in a seismic bulletin. Thus we should expect the observations to increase in number with increasing distance, up to a distance of 90 degrees, and then decay back to zero at 180 degrees. Actually small earthquakes are much more frequent than are large ones, and small shocks are not recorded at large distances. Thus the actual distribution begins to fall off well before 90 degrees: it looks like the dotted line in Figure 4. The peak in the distribution depends on the magnitude of the earthquake.

In conclusion, the probability of an observation is not uniformly distributed in distance. Therefore the sample of observations always tends to be dominated by distant observations. But the amplitude and the dominant frequency of the signal decrease with distance while the noise remains independent of distance, as it depends only on local conditions. Hence the signal/noise ratio tends to decrease with distance. If distant observations also tend to be more numerous, we have a serious bias. The location procedure is likely to be dominated by low-accuracy observations at large distances.

*Corollary 2.*

The joint-regression procedure may tend to introduce a systematic error into the epicentral solution when distant readings are more numerous, and less accurate, than are readings in the near field.

This corollary may seem unpersuasive at first sight, especially for a reader who is unfamiliar with observational seismology. The problem is inherent to location procedures. When distant observations tend to be late, near-field observations must tend to be early as a result of regression. The bias may not seem to be obvious, but it is a result of plate tectonics. Most earthquakes occur at plate boundaries, and most plate boundaries are subduction margins. Most subduction margins are also ocean–continent transitions. For obvious reasons, many more seismic stations are located on continents than in ocean basins. All these factors gang up against the valiant seismologist in his or her vain effort to produce accurate, unbiased locations.

For an epicenter near Acapulco, Mexico, for example, most Mexican and American stations are in the northwest quadrant while European stations are in the northeast quadrant. The nearest stations in the southwest quadrant are somewhere in Samoa or Australia. There are a few stations in the southeast quadrant, but they are an order of magnitude less numerous and powerful than are those in the north. The bias is systematic for this location, yet it is less tragically so than it would be, say, in South America or Indonesia.

## An Ill-Posed Problem

A problem is well posed when a solution exists, is unique, and depends continuously on the initial data. It is ill posed when it fails to satisfy at least one of these criteria (Hadamard, 1923). A well-known example of an ill-posed problem is scene analysis, where one attempts to extract 3D information from a 2D image.

The above theorems, and others like them, suggest that the earthquake-location problem may be ill posed. For example, a least-square location procedure minimizes Σ(*ρ*^{2}), where *ρ* is the travel-time residual. This implies that Prob{+*ρ*} = Prob{−*ρ*}, as assumed by C. F. Gauss (1823) in his derivation of the bell-shaped curve named after him. But when a seismic signal arrives at a station it gradually emerges from the local noise. The arrival time is picked with some delay +*ρ*, which reflects inversely the signal/noise ratio. On the other hand, picking a signal in advance of its actual arrival time is a purely mental phenomenon, since there is nothing there to be picked. Thus a negative residual is not the same kind of object as a positive residual. In one case there is a signal, in the other the signal does not yet exist. The bias is a result of causality.

Somewhat similarly, in joint-regression procedures an error in focal depth does not affect the solution in the same way as an error in latitude or in longitude. The bias is due to spherical geometry. All our observations are on the surface of the sphere, and negative depths do not exist.

Well-posedness may be restored by restricting the class of admissible solutions using *a priori* knowledge, such as a travel-time table. In this case the solution space is assumed to be convex, but local minima cannot be eliminated and there is a feedback between location errors and the errors in the travel-time tables.

The problem is compounded by the fact that we locate an object in space by using measurements in the time domain. By reducing the problem to the time domain we blind ourselves to the fact that the residuals are actually vectors. A 1-sec residual is not the same thing (i.e., should not have the same effect on the solution) as a residual of 1 sec to the north.

## Discussion

I realize that I have not provided a solid theory of earthquake location, or a remedial approach supported by real data. The reason is clear. The location problem is ill posed and there are not enough observations of controlled sources.

Controlled explosions are rarely large enough to calibrate a regional network. The early results using nuclear tests for amending travel-time tables have been disappointing, precisely because of serious discrepancies with tables in common use. In the end, seismologists were tempted to merge explosion data with earthquake data. Travel times from earthquakes have even been used to correct explosion data!

Statisticians are apt to suggest that travel-time residuals should be transformed into distance vectors on the surface of the sphere. But statistics on the sphere is an arcane subject (Mardia and Gadsden, 1977). In theory there is no obvious way of translating a least-square algorithm in the time domain to the space domain. Should we minimize the vectorial sum of the distance residuals? This is not the correct answer either.

Beno Gutenberg used a graphical location method that was routinely utilized in Pasadena for years, particularly by Charles F. Richter in the classical joint research work *Seismicity of the Earth* (Gutenberg and Richter, 1954). Many solutions are on file at the Caltech Seismological Laboratory where Gutenberg's famed yellow pads are still kept. While the criteria used by Gutenberg to weigh stations were never set down in writing, they could be inferred from the solutions, as the residuals were normally labeled.

Gutenberg realized, even before plate tectonics, that seismic observations were not scattered at random about the epicenter. They were clustered in certain azimuths and distance ranges. Stations in the near field or in the southern continents were rare and might easily be overwhelmed by distant stations, for example, in the United States or in Europe, during regression. More importantly, Gutenberg never attempted to include the focal depth in his location procedure. Focal depth was estimated separately from reflected phases such as *pP*, *sP*, and *PcP*—never by joint regression. This was because Gutenberg was keenly aware of the fact that the sign of the focal-depth residual changes over small distances in the near field.

Consider an earthquake epicenter *E(λ,φ)*, where *λ* is the latitude and *φ* is the longitude. Let *E*′ be a trial epicenter in the neighborhood of *E*. For any set of *n* stations consider the travel times *t _{i}*,

*i*= 1, …,

*n*. Gutenberg proposed to graph the azimuth

*α*,

_{i}*i*= 1, …,

*n*of the observations against the corresponding travel-time residuals

*ρ*,

_{i}*i*= 1, …,

*n*; both referred to the trial epicenter

*E*′ (Fig. 5).

Now, Gutenberg fitted a sine-curve envelope to the early residuals as shown in Figure 5. Gutenberg and Richter always performed this crucial step in such a manner as to weigh the observations in terms of their position relative to the epicenter, and in terms of their personal opinion about the reliability of the station. Often an important station was used to “anchor” the sine fit in such a way that other arrivals would turn out having been late. Early observations were a problem. When they belonged to an “unreliable” station they were simply discarded.

Clearly, such a procedure would be difficult—though not impossible—to program on a computer. On the other hand, the number and quality of seismic observations have improved considerably since Gutenberg's time. Fitting a sine envelope to a set of points offers no major challenge.

Gutenberg's procedure was iterative. In order to obtain the next trial epicenter *E*″ he would simply add the amplitude and phase of the fitted sine curve to the latitude and longitude of *E*′. Then the residuals were recomputed for *E*″ and the procedure would be repeated. The iteration procedure would stop when the sine curve became a straight horizontal line. The method might be adapted to interactive graphics.

## Conclusion

Theorem (1) shows that the epicentral parameters *λ* and *φ* may be determined on the sphere without using travel times or travel-time tables. This suggests a location method that assumes only radial homogeneity and is free of any bias introduced by travel-time measurements on the sphere. The location can be made as accurate as desired. It could be adopted as the gold standard for earthquake locations.

The focal depth remains indeterminate and must be determined separately. This is because all the measurements are on the surface. In joint regression the errors in the focal depth will contaminate the epicentral solution and vice versa.

Theorem (2) shows that the estimation of focal depth from travel-time residuals is not unique, because the rays that exit the focal sphere at the equator divide the surface of the sphere in two regions. Residuals will have opposite effects on focal-depth estimation depending on which of these two regions contains the observation. Since the areas of the two regions are very unequal there is a strong bias in focal- depth estimation.

Theorem (3) is intended to show that epicentral location is biased too, but in a different way. In a random field of stations the number of observations increases as the sine of the distance. This means that most samples are dominated by distant observations. But at distant stations the frequency content of the initial *P* wave becomes increasingly similar to the frequency of the background noise. This means that detection of early arrivals becomes increasingly difficult, and residuals are increasingly positive. As stations tend to cluster in certain areas, such as California or Europe, epicentral solutions tend to be pulled in the direction of these areas.

These results are intended to illustrate the fact that the location of earthquakes from observations on the Earth's surface is an ill-posed problem. We do not know how to locate earthquakes. Therefore we don't know the structure of the earth's interior either, as it is derived from the inversion of travel-time tables obtained from earthquake locations. Some scientists have been aware of this situation for years, and have bravely attempted to restore well-posedness by inventing new ways of locating earthquakes, or by restricting the class of admissible solutions using *a priori* information.

Unfortunately, the location problem is the fundamental problem of seismology. If we do not know how to locate earthquakes, we cannot be certain of our information on the Earth's interior—not so much because it matters whether some earthquake is mislocated a few miles seaward or landward but because locations have been used since Zoeppritz (1907) to estimate travel-time tables. Location errors creep into travel-time tables and from there into models of Earth structure. There does not seem to be any good way of preventing this, since we are dealing with the effects of bias. Earthquakes are not homogeneously distributed in space and neither are stations. Plate tectonics make sure that observations will be skewed by the position of the major continents. This bias is systematically built into the travel-time tables.

- Manuscript received 4 March 2005.