Seismic Analysis of the 7 January 1998 Chemical Plant Explosion at Kean Canyon, Nevada

GENE A. ICHINOSE1, KENNETH D. SMITH, and JOHN G. ANDERSON1
University of Nevada, Reno, Seismological Laboratory
UNRSL Logo
Mackay School of Mines, Mail Stop-174, Reno, NV, 89557-0141 phone (775) 784-4265 fax (775) 784-1833 email ichinose@seismo.unr.edu
1Also at the University of Nevada, Reno, Department of Geological Sciences
Bulletin Seismological Society of America, Volume 89, No 4, pgs. 938-945 August 1999

Abstract

An accident at the Sierra Chemical Company Kean Canyon plant, 16 km east of Reno, Nevada, resulted in two explosions 3.52 seconds apart which devastated the facility. An investigation into a possible cause for the accident required the chronological order of explosions. We resolved the high-precision relative locations and chronology of the explosions using a cross-correlation method applied to both seismic and air-waves. The difference in relative arrival times of air-waves between the explosions indicated that the first explosion occurred at the northern site. We then determined two station centroid separations between explosions which average about 73 meters with uncertainties ranging from ± 17 to 41 meters depending on the alignment of station pairs. We estimated a centroid separation of 80 meters using P-waves with a larger uncertainty of ± 340 meters. We performed a grid search for an optimal separation and the azimuth by combining air-wave arrivals from 3 station pairs. The best solution for the relative location of the second explosion is 73.2 meters S35°E from the first explosion. This estimate is well within the uncertainties of the surveyed crater locations by the US Chemical Safety and Hazard Investigation Board (CSB). The CSB reported a separation of approximately 76.2 meters S33°E. The spectral amplitudes of P-waves are 3 to 4 times higher for the second explosion relative to the first explosion but the air-waves have similar spectral amplitudes. We suggest that this difference is due to the partitioning of energy between the ground and air caused by downward directivity at the southern explosion, and upwards directivity at the northern explosion. This is consistent with the absence of a crater for the first explosion and a 1.8 meter deep crater for the second explosion. 
Picture of accident scene

The Sierra Chemical Plant at Kean Canyon
near Mustang, Nevada. Photos courtesy of the
U.S. Chemical Saftey Hazard Review Board.

Introduction

On the morning of January 7, 1998, two explosions occurred 3.52 seconds apart at the Sierra Chemical Company facility, 16 km east of Reno, Nevada. Unfortunately, four people lost their lives and six more individuals were injured in the explosions that devastated the facility. The US Chemical Safety and Hazard Investigation Board (CSB), with the cooperation of several state agencies, initiated an investigation into the cause of the accident. The chronological order of the explosions was an important aspect of the investigation because it may help determine the cause of the accident. The CSB requested that the Seismological Laboratory independently estimate the chronology based on air and seismic wave arrival times from seismograms to compare with other results obtained by the CSB. 

The two explosions were recorded at several stations of the western Great Basin seismic network and by 3-component, telemetered, digital broadband stations (Figure 1). From these records, the estimated origin time of the first explosion is 15:54:03.0 ± 0.3 sec UTC (7:54 AM local time), with an approximate location of 39° 31.8'N, 119° 38.1'W (± 1.5 km). However, based on a survey by the CSB, the location of the northern crater's center is 39° 32.5'N, 119° 38.1'W (John Piatt, personal communication). When treated as an earthquake, the coda magnitude of the largest explosion is estimated to be Md2. 

Picture of accident scene

Booster room 1 was location between the sites
of the two explosions. Photos courtesy of the U.S.
Chemical Saftey Hazard Review Board.

An examination of the seismograms in Figure (2a) show that there are two similar shaped P-wave arrivals (P1 and P2) at station PAH separated in time by about 3.5 seconds. This is followed by two similar shaped N-waves (N1 and N2) also separated in time by about 3.5 seconds and arriving approximately 70 seconds after the P-waves. The N-wave is created by the shock wave traveling in air (Kanamori et al., 1991). The identification of these phases were based on the phase arrivals which do not exhibit a distance moveout expected for reflected air or seismic waves. The coherency between P- and N-waves also indicates two explosions at the same location. 

The (P2) phase in Figure (2a) appears to be an S-wave of the initial explosion but is actually the P-wave arrival from the second explosion which is nearly coincident with the expected S-wave arrival from the initial explosion at both PAH and WCN. The larger P-wave amplitudes for the second explosion suggest that this was the larger of the two explosions, although the S-waves from the initial explosion at both PAH and WCN may contribute some interference. 



Photos courtesy of the U.S. Chemical Saftey Hazard Review Board.

View of crater where PETN storage building used
to be located. This is the site southern explosion.

View of remains of booster room 2. This is the
site of the northern explosion.

Identifying pot locations in booster room 2.

Damaged mixing pot from booster room 2.

See Figure 5 for map of explosion sites. The onsite investigation revealed that the first explosion ocurred at the northern site in boster room 2 where boosters were manufactured. About 3 seconds later (3.54 seconds) falling debris generated a second explosion in an explosives (PETN) storage building about 75 meters to the south at the southern site. Investigators believe the first explosion most likely occurred in boster room 2 when a worker left materials in a mixing pot by mistake overnight. When the mixing pot was turned on at 7:54 am in the morning, its blade may have hit solidified explosives which triggered the electric shock wave that set off the initial explosion.




The CSB surveyed the accident site and reported that two explosions were separated by a horizontal distance of approximately 250 feet (76.2 meters), along a strike of 147° or S33°E (John Piatt, personal communication). Uncertainties in these measurements are due to uncertainty on where the "centers" of the explosions were located. The northern explosion did not leave a crater and occurred in a building that was formerly about 40 feet by 40 feet in dimension, and the southern explosion left a kidney-shaped crater that was about 30 feet across and 50 feet long (John Piatt, personal communication). A circular approximation would have a radius of 40 feet. Based on these dimensions, the separation between the centers of the explosions relative to the centroid could be uncertain by as much as several meters, and the azimuth could be uncertain by a few degrees. According to testimony to the CSB (John Piatt, personal communication), there were about 7500 to 8000 pounds of explosives (TNT or COMP-B) at the northern site, and about 15000 pounds of explosives (PETN) at the southern site. Based on several independent lines of evidence, the CSB has come to the conclusion that the northern explosion occurred first. 

Poupinet et al. (1984), Fremont and Malone (1987), VanDecar and Crosson (1990), Deichmann and Fernandez (1992), Dodge et al., (1996), and Lees (1998) used a cross-correlation method on nearly identical seismograms called multiplets for the high-precision relative hypocenter locations of earthquakes. Fremont and Malone (1987) directly tested the precision of this method with ground truth information from explosions and concluded that explosions within 250 meters could be relocated with an accuracy better than 20 meters. We have also observed numerous cases of multiplets in routine earthquake monitoring activities and would like to determine the resolving power and errors associated with this methodology for our regional seismic network. 

Figure 1. 

Figure 1.Western Nevada digital seismic stations, location of the explosion, and the location of the Reno urban area. 

Table 1. Station locations, and their distances and azimuths to the explosion site.

Code  Name  Latitude  Longitude  Elevation(km)  Distance(km)  Azimuth ° 
WCN Washoe City, NV  39° 18.1'N 119° 45.4'W 1.50 28.6 201.5
VIP Virginia Pk., NV 39° 45.2'N 119° 27.7'W 2.49 27.9 32.2
PAH Pah Rah Range,NV 39° 42.2'N 119° 23.1'W 1.50 28.3 49.4
BEK Bekwourth, CA  39° 52.0'N 120° 21.5'W 1.74 71.8 300.4
WAK Walker, CA  38° 30.3'N 119° 26.2'W 1.89 116.3 171.5
Figure 2a. 
Figure 2. (a) Vertical component seismograms from station PAH and WCN. The P1 and P2 labels indicate P-wave arrivals for explosion A and B. N1 and N2 labels indicate N-wave arrivals for explosion A and B. S1 is the predicted S-wave arrival. (b) Seismograms of vertical component N-waves (air-wave) and 3 component seismic body waves used in this analysis. The arrows point to the arrival of P-wave of explosion A and B. The N-waves are low pass filtered at 20 Hz.

Analysis and Results

We located the initial explosion (explosion A) from the P-wave arrivals recorded at the UNR Keck digital stations and one helicorder record from an analog station (Table 1). The existing regional network of analog stations did not trigger on the explosion. We used a 1D velocity model (Table 2) to estimate the absolute location of the first explosion. Since the location determined from the P-wave arrivals of the first explosion is only 1.2 km from the known mapped location of the explosions, our confidence in the velocity model in Table 2 is increased. 

Finding the relative centroid locations of explosions A and B, requires the precise time difference in P-wave or N-wave arrivals between explosions, tb-ta, where tb is the phase arrival time from explosion B and ta is from explosion A. To determine tb-ta, we performed cross-correlations on windows of the P- and N-wave arrivals in the frequency domain (Fremont and Malone, 1987). The frequency domain technique can establish relative time difference estimates that are below the limit imposed by the sampling interval. This allows relative locations with a precision on the order of several meters even with a system sample interval of 0.01 seconds/sample. We also compute the coherency of the phase arrivals to determine the frequency range over which the slope of the cross-spectrum phase is analyzed. 

Table 2. Velocity model used in location of explosion. 
P-wave(km/s) Depth(km)
3.0 0.0
4.5 1.0
5.5 2.0
6.0 4.0
6.1 7.0
6.2 12.0
6.4 18.0
6.8 28.0
7.8 38.0
 
 
The time difference, DELTA t, between two phase arrivals is proportional to the slope of the phase of the cross spectrum, 
Equation 1.(1) 

where phi(delta f) is the phase of the cross spectrum over a frequency range delta f, and n is the phase unwrapping uncertainty. 

The normalized coherency between two time series measures the similarity of their shapes, ranging between 0, when they are completely dissimilar, and 1 when they are identical. The coherency in the frequency domain, C(f), between the Fourier transforms of the seismograms s1(f) and s2(f) is defined here following Menke et al. (1990), 

Equation 2.(2) 

where f is frequency, < > denotes boxcar averaging over a frequency interval DELTA f centered at f, and s* denotes complex conjugation. 

The seismograms of both explosions were first windowed using 2 second windows around each of the P- and N-wave arrivals and then cosine tapered. Rough estimates of tb-ta were then determined by hand picking a common peak or trough of explosion A and B. The windowed seismograms were then shifted by this tb-ta before computing the coherency and cross phase spectra. We find that the coherency between phase arrivals falls-off at high frequencies, and therefore we only used frequencies between 5 and 20 Hz. A line is fit to the slope of the cross phase spectra in equation (1) by least squares, with an intercept fixed at 0 Hz and an n equal to 0. We finally correct the initial tb-ta value by adding it with DELTA t thus giving the precise time separation of the two phases. The cross-coherence plots and phase spectra from PAH are shown in Figure 3, and the apparent time lags with uncertainties derived from these cross spectra are given in Table 3. 

The results were not affected by different initial tb-ta but estimating it as close as possible to the actual time shift helps to reduce the uncertainty due to phase unwrapping. There is always a 2 pi n uncertainty in unwrapping the phase spectra but since an initial shift was performed, we expect n to be 0 or 1 and the maximum uncertainty in phase to be less than 2 pi. The fall-off in coherency at high frequencies is more likely due to the first explosion superimposed upon the record of the second, or possible differences in the details of the two source time functions rather than being due to slight differences in the travel paths caused by the difference in the source locations. The loss in coherency is reflected in the large uncertainties in estimating the slope of phase of the cross spectra. 

Table 3. Time after the first explosion until the maximum of the cross correlation of the first and second explosion.

Station  Component  P-wave  Uncertainty1  Air-Wave  Uncertainty1 
tb-ta(sec)  (ms)  tb-ta(sec)  (ms) 
WCN Z 3.599 36 -
WCN E - 3.403 27
PAH Z 3.582 38 3.582 27
PAH N - 3.542 36
WAK Z - 3.330 9
WAK N - 3.330 14
1 Standard Deviation of phase spectra converted to apparent time separation between explosions. 
Figure 3a.
Figure 3b. 
Figure 3 (a) Coherency plots of the P- and N-wave windows used in the analysis for 3 stations. (b)The slope of the phase of the cross phase spectra are fit from 5 to 20 Hz. 
Based on the time separations from the cross-correlation method, we estimate L, the distance separation of the second explosion relative to the first as a function of beta, the hypothetical direction from the first source to the second, using: 
Equation 3.(3) 

where c is the air velocity, equal to 343 m/s, DELTA tij is the difference in the arrival times between the explosions at stations i and j, and thetai is the azimuth from the first explosion to the ith station. A grid search for the best solution which approximately satisfies equation (3) is made over values of beta ranging from 0 to 180° for the three station pairs. The optimal values are shown in Figure 4 where the three lines approximately cross at beta approximately 145° (S35°E) and L approximately 73.2 meters. The relative location estimated using N-waves is shown in Figure 5 and unambiguously indicate that the initial explosion was at the northern site, which is consistent with the analysis of the CSB. The origin time of explosions B relative to A, Tb, is determine by, 

Tb = ( tb-ta ) + L cos ( beta - phi )/c (4) 

where phi is the azimuth of the station from explosion A. We fix L equal to 73.2 meters and beta = 145 ° from the optimal values determined above. Using tb-ta values from Table 3 and phi values from Table 1, results in a Tb equal to 3.52 seconds. 

 
Figure 4. 
Figure 4. A plot of L, the source separation versus beta, the azimuth between sources. The arrow points to the best solution of L

Discussion and Uncertainties

The optimal relative location of the second explosion, based on 3 pairs of stations, is within the uncertainties of the CSB investigation. The difference in separation between our estimate (73.2 meters) and the CSB estimate (76.2 meters) is small compared to the source dimensions, and the difference in azimuth between our estimate (145°) and the CSB estimate (147°) is also within the range of angles that is allowed by the source sizes. There is no straight forward way of computing an uncertainty for the optimal values of L and beta, but since the three solution lines in Figure 4 do not exactly intersect, the area created by the intersection illustrates an approximate uncertainty. 

We compute uncertainties for individual station pairs by the standard deviation, sigma, which is the square root of the variance between the best fit slope of equation (1) to the phase of the cross spectra. These standard deviations are converted to tb-ta uncertainties and listed in Table 3. We then propagate tb-ta uncertainties through equation (3) by arbitrarily fixing $beta$ equal to 145°, and using ± 1 sigma in tb-ta to convert to the maximum uncertainty in the source separation, L, for each station pair. The importance of receiver geometry on uncertainty is shown by the difference in separation uncertainties between different station pairs. Station pair PAH and WAK has a path which is almost along strike of the explosions and has an uncertainty of ± 17 meters. The station pair WCN and PAH has a path relatively more perpendicular to the strike of the explosions and has a larger uncertainty of ± 33 meters. The path which is more parallel to the strike of the explosions has a better resolution of the separation. 

 
There is no significant source of error associated with timing in the recorder itself. The digital stations maintain absolute timing by synchronizing with a GPS time signal that is broadcast from the Seismological Laboratory. The GPS signal is broadcast every second and a high precision oscillator in the seismograph unit is phase-locked to UTC by this pulse. A radio frequency delay of 44 ms, which occurs in the electronics and telemetry systems, is accounted for in establishing absolute time of the recorded waveforms. A timing mismatch of 1 msec between the GPS time and the internal clock time results in a clock correction that is reported by the instrument. Timing errors during regular operation rarely exceed several msec. Because the two explosions occurred within 3.52 sec, the absolute timing of the instrumentation can affect the location that we gave for the first explosion, but not the relative locations of the first and second explosion. Only the error in the digitization rate is relevant to uncertainties in the relative locations. The manufacturer reports that errors in the digitization rate for the internal oscillator do not exceed 1 msec for any one sample and are expected to be on the order of 50 microsec. If the oscillator would drift more than 1 msec over any recording period, then a timing correction would be initiated by the instrument and its results would be recorded in the instrument log. 

Atmospheric conditions that affect the speed of sound can shift the estimated separations slightly. We used an air velocity of 343 m/s (Kinsler et al., 1982). For a 10% uncertainty in the assumed air velocity of 343 m/s, which is greater than expected, the relative source separation error would be ± 7.0 meters, which does not impact our conclusion as to the relative source locations or the chronology. A constant wind velocity across the array on the morning of the blasts would not be significant; a 10 MPH wind has only 1.3% of the speed of air. 

Figure 5. 
Figure 5. Schematic map of the explosion site with the two explosions denoted as stars. L and beta is the best estimate separations of the two sources, and theta is the azimuth from explosion A to the station [see equation (3)]. The variable r is the distance between explosion A and the station but drops out in the derivation of equation (3).

Explosion A ocurred at the northern site, boster room 2, where boosters were manufactured. Seconds later falling debris generated a second explosion, "Explosion B", in an explosives storage building at the southern site. Investigators believe Explosion A most likely occurred when a worker left materials in a mixing pot by mistake overnight. When the mixing pot was turned on in the morning, its blade may have hit solidified explosives which triggered the electric shock wave that set off the initial explosion.

One of our objectives was to see if it is feasible to measure the separation between the explosions using the P- or S-wave waves. From Table 4, we see that the best estimate of the separation using the P-waves is 80 m, which is consistent with the estimates using the N-waves. However, the uncertainty in the time separation for the P-waves from WCN and PAH leads to a large uncertainty on the separation. The P-wave at WAK, 116 km epicentral distance, was too weak to provide a reliable separation. These results confirm that high-precision relative locations of closely spaced multiple events can be resolved, given an adequate signal-to-noise ratio, and two or more stations distributed along the strike of the events (e.g. Fremont and Malone, 1987). 

Table 4. Geometry and results of source separation estimation.

Path P-wave (sec) * Air-wave (sec) *  ° Theta sub i° Separation (m) Uncertainty (m) 
WAK-PAH 0.2126 -26.4 95.7 72.3 ± 17
WAK-WCN 0.0752 -26.4 -56.4 73.3 ± 41
WCN-PAH 0.1374 -56.4 95.7 73.2 ± 33
WCN-PAH 0.0144 -56.4 95.7 79.9 ± 340
t sub ij is the difference in tb-ta between station pairs along path. 
 
Figure 6. 
Figure 6. Uncorrected P-wave and N-wave spectra of the two explosions. The spectral curves are the smoothed log average of three components of motion. The arrow points to the peak spectral value of a noise window before explosion A. 

Directivity and Source Dimension

Figure 6 shows the uncorrected spectra of the P- and N-waves from the two explosions at PAH. The spectra of the three components were log averaged and then smoothed. The P-wave spectra at PAH for explosion B is 3 to 4 times larger than explosion A, which is consistent with the report that the site of explosion B contained more explosives however, the amplitude ratio is twice the reported ratio of explosive mass. In comparison, the N-wave spectral amplitudes for B are similar to explosion A, which is inconsistent with the P-wave spectra. The larger P-wave spectral ratio between explosions and a smaller N-wave spectral ratio suggests that explosion B was more coupled to the ground, thus allowing more energy to be partitioned into the ground than into the air. We speculate that explosion B initiated at the top of the stockpile, resulting in downward directivity, partitioning more energy into the ground than into the air, and forming a crater 1.8 meters deep, while explosion A had upwards directivity consistent with the absence of a crater. The CSB hypothesized that the first explosion triggered the second explosion when debris from the first explosion entered the building through the ceiling or skylight at the second site. This mechanism is consistent with the observed directivity. The S-wave from the first explosion is predicted to arrive 4.2 seconds after the P-wave at PAH (Figure 3). The S-wave may contaminate the P-wave spectra of the second explosion, therefore making it higher in amplitude. This could suggest an alternative scenario in which two explosions had the same yields, as suggested by the similar N-wave spectral ratio. However the spectra for the first 0.5 seconds still differ by a factor of 2, making this scenario less likely. 
 
The physical dimensions of the explosions, like the dimensions of earthquakes, should be related to the corner frequency, fc, measured from the Fourier spectrum. To test this hypothesis, we used an equation for the earthquake source radius r from Brune (1970): 
Equation 5.(5) 

where c is either the P-wave velocity or the speed of sound in air. The measured crater of the southern explosion was about 12 meters across and 1.8 meters deep, suggesting that the expected source radius is 6 meters. The northern explosion did not leave a crater but the building was formerly 40 feet by 40 feet, giving an upper limit to the source radius of about 6 meters. 

The Fourier spectra from the N-waves, shown in Figure 6b, are relatively flat from 1 Hz to above 30 Hz. The high frequencies of the N-waves are limited by the anti-aliasing filters in the recorder (at 40 Hz). The spectrum from the northern explosion might suggest the presence of a corner at about 30 Hz, which from equation (4) would give a source radius of 4 meters. Such a result is reasonable considering the independent information about the size of the building. 

The P-wave spectra fall off rapidly above 6 Hz, so we take 6 Hz to represent the corner frequency of these spectra. In equation (4) we use the P-wave velocity of 3000 m/sec from the topmost layer of the velocity model in Table 1 which is reasonable for weathered bedrock. With this combination of parameters, equation (4) gives r equal to 60 meters, which is about factor of of ten larger than the estimate from the N-wave and from ground observations. We therefore suggest that attenuation along the path has played a major role in decreasing the amplitude of high frequency P-waves. The P-waves spectra have decreased to amplitudes comparable to the pre-event noise above 20 Hz, implying that the attenuation eliminates the chance to use P-waves to estimate the source dimension for small events, whether earthquakes or explosions. For earthquakes, the P-waves would pass through the near surface zone of severe attenuation only once, so resolution of a high corner frequency should be somewhat better than in this case. 

 

Conclusions

We have analyzed the relative arrivals of both similar looking P- and N-waves at a number of seismic stations to estimate the spatial separation and orientation of two closely-spaced explosions. Analysis of seismograms indicate that two explosions occurred separated in time by 3.52 seconds. The best solution for the source separation and azimuth is determined by the combination of N-waves from 3 station pairs. We conclude that the initial explosion occurred at the northern site and that the southern explosion occurred about 73.2 meters S35°E. Individual estimates using only a pair of stations have relatively higher uncertainties but still supports our conclusions. The separation and orientation of the two explosions were consistent with, and well within uncertainties of the data provided by the CSB. P-wave arrivals recorded only for one source-receiver path yields a separation of 80 m with an individual uncertainty of ± 340 meters. 

From the relative spectral amplitudes of P- and N-waves, we speculate that explosion B may have had a downward directivity, whereas explosion A may have been more upwardly directed. From the viewpoint of forensic seismology, this experiment was successful, in that the N-waves unambiguously demonstrate that the northern of the two explosions occurred first. We confirm that the relative separation of sources can be determined precisely using only a pair of regional seismic stations with an optimal alignment. 

Acknowledgments

Our condolences goes out to the families of the Sierra Chemical workers who were killed in the accident. David von Seggern estimated the magnitude of the explosion. The reviews provided by D. von Seggern, D. Dodge, and an anonymous reviewer were extremely helpful and led to significant improvements. We thank the Keck Foundation for their generous gift that allowed installation of the digital stations used in this research. This work was made possible through financial support provided by U.S. Geological Survey NEHRP grant 1434-94-G-2479. 

References

  • Brune, J. (1970). Tectonic stress and the spectra of seismic shear waves from earthquakes, J. Geophys. Res 75, 4997-5009. 
  • Deichmann, N., and M. Garcia-Fernandez (1992). Rupture geometry from high-precision relative hypocentre locations of microearthquake clusters, Geophys. J. Int. 110, 501-517. 
  • Dodge, D. A., G. C. Beroza, and W. L. Ellsworth (1996). Detailed observations of California foreshock sequences: Implications for the earthquake initiation process, J. Geophys. Res 101, 22371-22392. 
  • Fremont, M.-J., and S. D. Malone (1987). High precision relative locations of earthquakes at Mount St. Helens, Washington, J. Geophys. Res. 92, 10233-10236. 
  • Lees, J. M. (1998). Multiplet analysis at Coso Geothermal, Bull. Seism. Soc. Am. 88, 1127-1143. 
  • Kanamori, H., J. Mori, D. L. Anderson, and T. H. Heaton (1991). Seismic excitation by the space shuttle Columbia, Nature 349, 781-782. 
  • Kinsler, L. E., A. R. Frey, A. B. Coppens, and J. V. Sanders (1982). Fundamentals of Acoustics, 3rd Ed., John Wiley and Sons, New York. 
  • Menke, W., A. L. Lerner-Lam, B. Dubendorff, and J. F. Pacheco (1990). Polarization and coherence of 5 to 30 Hz seismic wave fields at a hard-rock site and their relevance to velocity heterogeneities in the crust, Bull. Seism. Soc. Am. 80, 430-449. 
  • Poupinet, G., W. L. Ellsworth, and J. Frechet (1984). Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras fault, California, J. Geophys. Res. 89, 5719-5731. 
  • VanDecar, J. C., and R. S. Crosson (1990). Determination of teleseismic relative phase arrival times using multi-channel cross-correlation and least squares, Bull. Seism. Soc. Am. 80, 150-169. 

Links