Locating an earthquake from P and S arrival times

An earthquake has four unknowns: its origin time and the three coordinates of its hypocentre, the point where the rupture began. Each station adds one equation, its arrival time. The equations are not linear in the unknowns, so the standard method guesses a location, linearizes the travel times around the guess, solves a small least-squares problem for a correction and repeats. That method is Geiger’s, from 1912, and it still sits inside the USGS location programs.

Sandbox Physics · 2026-10-30

The S-minus-P time measures distance

The P wave and the S wave leave the hypocentre at the same moment. Since S is slower, the gap between the two arrivals grows with distance, and it does not depend on when the quake happened. In rock where the S speed is the P speed divided by the square root of three, a 6 km/s P speed gives about 8.2 km of distance for every second of gap.

One station with both picks therefore puts the hypocentre on a sphere around it. Three spheres meet at a point, at least in principle. Real picks have errors and the speed model is only approximate, so a better method uses every pick at once.

t_{\rm S}-t_{\rm P}=r\left(\frac{1}{\beta}-\frac{1}{\alpha}\right),\qquad r=\frac{t_{\rm S}-t_{\rm P}}{1/\beta-1/\alpha}

Linearize the travel time around a guess

Lee and Lahr summarize Geiger’s method in Section 4.1 of the HYPO71 report. Pick a trial origin time and hypocentre. For each station compute the predicted arrival time and the residual, the observed time minus the predicted one (their Eq. 1).

If the guess is close, each residual is nearly a linear combination of four small corrections: the shift in origin time plus each spatial shift times the partial derivative of travel time along that direction (their Eq. 2). For a straight ray in uniform rock, the derivative along any direction is minus the component of the unit vector from source to station, divided by the speed.

Minimizing the sum of squared misfits leads to four normal equations (their Eqs. 3 to 5). Their solution moves the trial solution (Eq. 6), and the step is repeated until the corrections are small. If S-minus-P times are used instead of absolute times, the origin time drops out (Eq. 7).

R_i=T_i-t_i\approx\delta t_0+\frac{\partial t_i}{\partial x}\delta x+\frac{\partial t_i}{\partial y}\delta y+\frac{\partial t_i}{\partial z}\delta z,\qquad \mathbf G^{\mathsf T}\mathbf G\,\delta\mathbf m=\mathbf G^{\mathsf T}\mathbf R

Depth and origin time pull against each other

Stations far away from a shallow quake all see rays leaving at a low angle. Moving the source down then delays every arrival by about the same small amount, which looks almost the same as starting the quake slightly later. The columns of the matrix for depth and origin time are nearly parallel, and the normal equations become ill-conditioned.

Lee and Lahr name this problem and describe HYPO71’s answer: a step-wise regression that only adjusts the variables the data can determine at each iteration. The location Lab uses a different common remedy, a Levenberg–Marquardt step that adds a small multiple of the matrix diagonal and grows it only if a full step would raise the misfit.

S picks help more than extra P picks. An S wave from the same point arrives with a different slope against distance, so it separates depth from origin time. A station close to the epicentre helps too, because its ray leaves almost vertically.

\delta\mathbf m=\left(\mathbf G^{\mathsf T}\mathbf G+\lambda\operatorname{diag}\mathbf G^{\mathsf T}\mathbf G\right)^{-1}\mathbf G^{\mathsf T}\mathbf R

Worked example

One distance and one row of the matrix

Use the Lab’s rock: a P speed of 6 km/s and an S speed of 6 km/s divided by the square root of three, 3.46 km/s. The gap per kilometre is 0.289 minus 0.167 seconds, 0.122 s. Its inverse is 8.20 km/s.

A station records the S wave 3.00 s after the P wave. The hypocentral distance is 3.00 times 8.20, 24.6 km. If the quake is 10 km deep, the epicentre lies 22.5 km from the station, on a circle around it.

Now take a trial hypocentre at 10 km depth directly below the origin of the map and a station 30 km east of it. The distance is 31.6 km and the P travel time 5.27 s. Moving the source 1 km east shortens the ray, so the derivative along east is minus 30 over 31.6 times 6, minus 0.158 s/km. Moving it 1 km down lengthens the ray by the depth over the distance, so the depth derivative is 10 over 31.6 times 6, 0.0527 s/km. With the origin-time column of ones, that station’s row of the matrix is 1, minus 0.158, 0, 0.0527.

\frac{1}{1/\beta-1/\alpha}=8.20\ \mathrm{km/s},\quad r=24.6\ \mathrm{km};\qquad \mathbf G_i=\big(1,\ -0.158,\ 0,\ 0.0527\big)\ \mathrm{s/km}

An experiment to think through

Predict, change, explain

In the location Lab, start from Ring network: a tight ellipse, then switch S picks to P only. Predict which grows more: the epicentre ellipse or the depth uncertainty.
Choose P picks only: depth goes loose and compare the depth section with the ring network preset.
Explain the change with the depth and origin-time columns of the matrix.

Locate quakes and test their error ellipses

The Labs use flat, isotropic, elastic or acoustic rock with made-up but realistic speeds; none of them reproduces a real survey. Each Lab states its model and limits in its notes and can export the data behind its readouts.

Try explaining it yourself

1. A station records P at 10.0 s and S at 14.0 s on its clock. With the Lab’s speeds, find the hypocentral distance and the origin time.

Read an explained answer

The gap is 4.0 s, so the distance is 4.0 times 8.20, 32.8 km. The P travel time is the gap divided by the speed ratio minus one, 4.0 over 0.732, which is 5.46 s. The quake therefore started at 10.0 minus 5.46, which is 4.54 s.

r=4.0\times 8.20=32.8\ \mathrm{km},\qquad t_P=\frac{4.0\ \mathrm s}{\sqrt3-1}=5.46\ \mathrm s,\qquad t_0=10.0-5.46=4.54\ \mathrm s

2. Why does a residual that is small at every station not prove the location is right?

Read an explained answer

The residuals are computed with the assumed speed model. A wrong model can fit the picks well while putting the source in the wrong place. Lee and Lahr warn that small residuals and standard errors are not sufficient to guarantee an accurate hypocentre.

3. All stations lie east of a quake. Which direction of the epicentre is worst constrained, and why?

Read an explained answer

East–west. Every ray leaves towards the east, so moving the source east or west changes all travel times by nearly the same amount, which the origin time can absorb. North–south moves change the times differently from station to station.

Sources to read alongside

MIT OpenCourseWare, 12.201 Essentials of Geophysics (Fall 2004), lecture notes, Chapter 4: Seismology
Section 4.6 separates the elastic wave equation into P and S waves and gives their speeds (Eqs. 4.34 and 4.37). Sections 4.14 to 4.16 derive Snell’s law from Fermat’s principle (Eqs. 4.91 to 4.96) and show that the ray parameter is the slope of the travel-time curve (Eq. 4.98). Section 4.19 notes that in travel-time tomography the matrix entries are the ray lengths in each block.
W. H. K. Lee and J. C. Lahr, HYPO71 (revised): a computer program for determining hypocenter, magnitude, and first motion pattern of local earthquakes, U.S. Geological Survey Open-File Report 75-311 (1975)
Section 4.1 sets out Geiger’s method: expand the arrival-time residuals to first order in the origin time and the three coordinates (Eq. 2), solve the least-squares normal equations (Eqs. 3 to 5), update the trial hypocentre (Eq. 6) and repeat. The introduction warns that small residuals and standard errors do not guarantee an accurate hypocentre.
F. W. Klein, User’s guide to HYPOINVERSE-2000, a Fortran program to solve for earthquake locations and magnitudes, U.S. Geological Survey Open-File Report 02-171 (2002)
Appendix 4 builds the 4 by 4 covariance matrix from the eigenvalues and eigenvectors of the location problem, scaled by the variance of the arrival times, and takes the error ellipsoid from its 3 by 3 spatial block. It notes that the ellipse scales with the assumed timing error, and that a true error estimate would also need the uncertainty of the crustal model, which is nearly impossible to know.