Detecting and locating seismic events with using USArray as a large antenna

We design an earthquake detection and location algorithm that explores coherence and characteristic behavior of teleseismic waves recorded by a large-scale seismic network. The procedure consists of three steps. First, for every tested source location we construct a time-distance gather by computing great-circle distances to all stations of the network and aligning the signals respectively. Second, we use the constructed gather to compute a Tau-P transform. For waves emitted by teleseismic sources, the amplitude of this transform has a very characteristic behavior with maxima corresponding to different seismic phases. Relative location of these maxima on the time-slowness plane strongly depends on the distance to the earthquake. To explore this dependence, in a third step, we convolve the Tau-P amplitude with a time-slowness filter whose maxima are computed based on prediction of a global travel-time calculator. As a result of this three-step procedure, we obtain a function that characterizes a likelihood of occurrence of a seismic event at a given position in space and time. We test the developed algorithm by applying it to vertical-component records of USArray to locate a set of earthquakes distributed around the Globe with magnitudes between 6.1 and 7.2.


Introduction
Traditional earthquake detection and location algorithms use the arrival times of short period body wave phases such as P or S. The waves are identified at every individual station when the wavefront generated by the source passes through it.The measured arrival times are then grouped and inverted to locate the event.An alternative approach is to study seis-mic events with using a network of seismic receivers as an antenna.The records by large regional-scale seismic networks installed during recent decades such as USArray in the United States, HI-Net Array in Japan or VEBSN in Europe have been recently used to study sources of large earthquakes such as the 2004 Sumatra-Andaman event (Ishii et al., 2005), the 2011 Tohoku earthquake (e.g., Satriano et al., 2014) and the 2012 Wharton Basin earthquake (Satriano et al., 2012) based on beam forming algorithms.More recently, Retailleau et al. ( 2014) used a network-based analysis of the US-Array data to study the antipodal focusing of seismic waves emitted by an earthquake in the Indian Ocean.In these studies, the authors focused on characterization of details of the source processes of large earthquakes with a priory known locations.In the present paper, we explore another approach when records from large networks are used to detect and to locate epicenters of seismic events without a priory knowledge about their occurrence.We first test the method with a set of synthetic seismograms computed in the spherically symmetric Earth model PREM (Dziewonski and Anderson, 1981) and then apply it to a few real earthquakes.

Data
We use vertical-component records from the Transportable Array (TA) component of the USArray, which contained roughly 400 receivers during 2010-2011.We selected 200 closely located stations forming a compact group that can be used as an antenna (Fig. 1).A quality check was performed to suppress data with high amplitude glitches.We selected five earthquakes, described in Table 1, to illustrate the method.

L. Retailleau et al.: Detecting and locating seismic events
Table 1.Description of the 4 earthquakes on which the method was tested.These characteristics were obtained from the CMT catalog (Dziewonski and Anderson, 1981;Ekström et al., 2012) For every earthquake we used records of 4 h starting at the source time of the event.
We also created a synthetic dataset to test the method.For this purpose, we used the focal mechanism from the Global CMT catalog for event as our source.Its location was set to longitude 140 and latitude 0 as shown in Fig. 2a.The synthetic seismograms were computed with the method AXISEM, a 2-D spectral-element solver for 3-D elastodynamics in global, spherically symmetric background models (Nissen-Meyer et al., 2007), using the isotropic version of the PREM model (Dziewonski and Anderson, 1981) and with including the viscoelastic attenuation.Both the real signals and the synthetic seismograms were bandpassed in the period band between 20 and 50 s.

Method
The basic idea of our method is to use simultaneously distinctive characteristics of different seismic body-wave phases recorded by an array of receivers.Every seismic phase is characterized by its slowness and back-azimuth that can be measured with an array-based analysis.For body waves emitted by the same earthquake, the back-azimuths measured at a given time are similar, because the seismic structure of the Earth is well approximated by a spherically symmet- ric model.The slowness of every seismic phase is varying with the epicentral distance.When using N seismic phases, the combination of their slownesses is becoming a N values vector whose dependence on the epicentral distance is much more selective comparing to a single phase.
Main steps of our method are illustrated in Fig. 2 with using the synthetic dataset.First we select a tested source position (Fig. 2a) and construct a time-distance gather by computing great-circle distances to all stations of the network and aligning the signals respectively (Fig. 2b and c).Then, we compute a Tau-P transform to present the gather energy as function of travel time and apparent slowness (Fig. 2d and  e).
The Tau-P transform is smoothed in time with a 30 s sliding window to suppress high frequency variations.Then, we decided to emphasize the energy of the first arriving body waves that is lower than the energy of the latter arrivals (especially the surface waves).The early body-wave phases have narrow maxima on the time-slowness plane and are, therefore, more useful for determining the distance with our method.To increase their relative contribution, we apply a normalization, consisting in dividing the Tau-P transform by its average over all slownesses smoothed in time with a 600 s sliding window.Finally, we design a time-slowness filter whose maxima are computed based on predictions of the global travel-time calculator of Crotwell et al. (1999).
This filter consists of a set of rectangles centered at arrival times and slowness of selected phases as illustrated by yellow lines in Fig. 2d and e.The filter is continuously shifted along the time coordinate and the time-slowness energy within the boxes is summed.When the tested source position coincides with the epicenter and the time coincides with the arrival time of all phases, the summed energy E sum shows a clear maximum (Fig. 2f).If the tested source position is far from the epicenter, the values of E sum are relatively weak (Fig. 2g).We perform a grid search and test different source positions located on a geographical grid.The results of this grid search for the synthetic test are shown in Fig. 3 and show a clear maximum located very close to the true source position.

Discussion
After testing the method with the synthetic seismograms, we applied it to records of real earthquakes shown in Table 1.The results shown in Fig. 4 demonstrate the example of three events of magnitude 7 and one of magnitude 6 which were well located.These events occurred at different locations on Earth with different depths and magnitudes providing a first view of the abilities of the method.
Further assessment of the method performance requires its testing with more events and with signals filtered in different frequency bands.Indeed these different frequency bands should highlight different events in the signals and a higher frequency would locate the events more precisely.Another important issue that remains out of scope of this short paper is the possible depth sensitivity of the proposed methods.In the presented tests, the source depth used to compute travel times and to design the time-slowness-filters was fixed to their depth (obtained from the Global CMT catalog, Table 1).We expect that for shallow events (above 50 km) the difference in depth should not have much consequences in the shape of the time slowness-filters and the resulting E sum values.However, the depth might become important for deep earthquakes (such as event 2 with its depth of 586 km) when the pP and sP phases appear.Including these phases into the analysis may help to discriminate the depth.
Methods to study the slowness and back azimuth of a source in the time domain have been developed at Norwegian Seismic Array (NORSAR) long ago (King et al., 1976).Their method generates a function of energy of the waves as opposed to time, azimuth, and slowness.Nonetheless our method also takes into account the distance of the source to the network, and thus the location of the source of signal.Indeed our mask study obtains the decays of arrivals between the different phases, which give a good idea of the location of the event.We are able to get more information from the Tau-P transform.
Moreover, we do not only focus on a precise location as Ishii et al. (2005) did to study the Sumatra rupture.Adding several phases increases the precision and enables us to observe events spread worldwide, in all the range of distances.
Basically our method does not need any a priori information such as the time and location of the event we want to detect.It is possible to detect two events that occurred at the same time in different places.The phases also do not need to Since this approach does not need any information about the events it should lead to the detection of events we do not expect.For this purpose, we intend to apply the method to a few years of continuous USArray records to evaluate its capacity for automatic detection and location of different types of seismic events.

Figure 2 .
Figure 2. (a) Map showing the network used and two points showing the correct location (in red) and a wrong one (in blue).Details of the method for two potential locations, the correct (b, d, f) and a wrong one (c, e, g).(b) and (c) show the alignment of the seismograms with respect to their respective potential source.(d) and (e) Tau-P transforms and theoretical masks (in orange) generated for each case.(f) and (g) The energy functions E sum generated.The theoretical arrival of the first phase is shown with a dashed line.

Figure 3 .
Figure 3. Application of the described network-based location method to the synthetic dataset.Colors show the maximum of E sum at every location.The black point shows the correct location of the source.

Figure 4 .
Figure 4. (a) Results of the location method for four earthquakes described in Table 1, which locations are represented in (a).The dashed black line shows the 90 % energy around the maximum of E sum .(b) Vanuatu 27 May 2010 magnitude 7.2 earthquake, (c) Ryukyu islands (Japan) 26 February 2010 magnitude 7.0 earthquake, (d) Santiago del Estero (Argentina) 1 January 2011 magnitude 7.0 earthquake, and (e) Turkey 8 March 2010 magnitude 6.1 earthquake.Table 1 describes the events.
et al.: Detecting and locating seismic events be selected and identified, nor their time arrivals measured a priori. .
Figure 1.Map of the United States where the red circles show locations of the USArray Transportable seismic stations used in this study.
Table 1 describes the events.