# Detection and analysis of a transient energy burst with beamforming of multiple teleseismic phases

Lise Retailleau, Matthieu Landès, Lucia Gualtieri, Nikolai M. Shapiro, Michel

Campillo, Philippe Roux, Jocelyn Guilbert

► To cite this version:

Lise Retailleau, Matthieu Landès, Lucia Gualtieri, Nikolai M. Shapiro, Michel Campillo, et al. Detection and analysis of a transient energy burst with beamforming of multiple teleseismic phases. Geophysical Journal International, 2018, 212, pp.14-24. 10.1093/gji/ggx410. insu-03589387

HAL Id: insu-03589387

https://insu.hal.science/insu-03589387

Submitted on 25 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Detection and analysis of a transient energy burst with beamforming of multiple teleseismic phases

Lise Retailleau,<fnref n="1" />,∗ Matthieu Landès,<fnref n="2" /> Lucia Gualtieri,<fnref n="3" /> Nikolai M. Shapiro,<fnref n="4" />,<fnref n="5" /> Michel Campillo,<fnref n="1" /> Philippe Roux<fnref n="1" /> and Jocelyn Guilbert<fnref n="6" />

Accepted 2017 September 22. Received 2017 September 20; in original form 2016 November 16

## SUMMARY

Seismological detection methods are traditionally based on picking techniques. These methods cannot be used to analyse emergent signals where the arrivals cannot be picked. Here, we detect and locate seismic events by applying a beamforming method that combines multiple body-wave phases to USArray data. This method explores the consistency and characteristic behaviour of teleseismic body waves that are recorded by a large-scale, still dense, seismic network. We perform time-slowness analysis of the signals and correlate this with the time-slowness equivalent of the different body-wave phases predicted by a global traveltime calculator, to determine the occurrence of an event with no a priori information about it. We apply this method continuously to one year of data to analyse the different events that generate signals reaching the USArray network. In particular, we analyse in detail a low-frequency secondary microseismic event that occurred on 2010 February 1. This event, that lasted 1 d, has a narrow frequency band around 0.1 Hz, and it occurred at a distance of 150° to the USArray network, South of Australia. We show that the most energetic phase observed is the PKPab phase. Direct amplitude analysis of regional seismograms confirms the occurrence of this event. We compare the seismic observations with models of the spectral density of the pressure field generated by the interferences between oceanic waves. We attribute the observed signals to a storm-generated microseismic event that occurred along the South East Indian Ridge.

**Key words:** Body waves; Seismic noise; Wave propagation.

## 1 INTRODUCTION

Dense broad-band seismic networks have recently been deployed at the continental scale, such as the USArray (United States), F-net (Japan) and IberArray (Spain) networks. They provide the opportunity to record continuous signals on a large scale. The installation of networks containing numerous stations closely located began with the need to increase the signal-to-noise ratio for detection of nuclear tests. The first networks that were built with this idea were the Large Aperture Seismic Array (Green et al. 1965) and the Norwegian Seismic Array (Bungum et al. 1971). These

networks also provided the chance to discover new kinds of seismic events. Indeed, some natural seismic sources have an emergent character and thus do not present any clear arrival. This can make them difficult to be detected by traditional methods, as they use the arrival times of the first body-wave phases at individual stations, before they are grouped and associated with a potential event. Networks permitted to discover, for example ‘non-volcanic tremors’ (Obara & Ito 2005) and glacial earthquakes (Ekström et al. 2003).

New methods are based on coherence of short-period seismic radiation backpropagated to the most likely source location in space and time. Beamforming algorithms have been used with Hi-net data to study the rupture of large earthquakes, such as the 2004 Sumatra–Andaman event (e.g. Ishii et al. 2005) and the 2011 Tohoku earthquake (e.g. Satriano et al. 2014).

<footnote n="1">¹Institut des Sciences de la Terre, CNRS, Université Grenoble Alpes, Grenoble, France. E-mail: retaille@stanford.edu</footnote>

<footnote n="2">²Euro-Med Seismological Centre, DAM, DIF F-91297 Arpajon, France</footnote>

<footnote n="3">³Lamont-Doherty Earth Observatory of Columbia University, Palisades, NY 10964, USA</footnote>

<footnote n="4">⁴Institut de Physique du Globe de Paris, Sorbonne Paris Cité, CNRS (UMS 7154), F-75005 Paris, France</footnote>

<footnote n="5">⁵Institute of Volcanology and Seismology FEB RAS, 9 Piip Boulevard, Petropavlovsk-Kamchatsky, 683000, Russia</footnote>

<footnote n="6">⁶CEA/DAM/DIF, F-91297 Arpajon, France</footnote>

<footnote>*Now at: Department of Geophysics, Stanford University, Stanford, California, USA.</footnote>
Global detection methods that use arrays of receivers have previ-
ously been shown to be efficient event detectors (e.g. Shearer 1994;
Ekström 2006). These studies analysed surface waves at long period
using global networks [IDA in Shearer (1994) and GSN in Ekström
(2006)]. In this paper, we analyse the body-wave phases recorded
on a regional network (USArray), following the method detailed in
Retailleau et al. (2015). Notably, Retailleau et al. (2015) introduced
a beamforming method that uses multiple body-wave phases, mak-
ing it possible to analyse events without prior information about
their location, and at all epicentral distances. This means that this
method can locate events automatically.

Beamforming methods are also used to study seismic signals
generated by meteorological events. Gerstoft et al. (2006) identi-
fied body waves generated by the hurricane Katrina using beam-
forming analysis (Rost & Thomas 2002). Similarly, Gerstoft et al.
(2008) used a beamforming method based on P-waves backpro-
jection, and showed that body waves are generated in the ocean by
storms, and that they propagate afterwards within the Earth as P, PP
and PKP seismic phases. Gerstoft et al. (2008) also observed some
variability in the frequency content depending on the storm. More
recently, Landès et al. (2010) studied microseismic sources asso-
ciated with P, PP and PKP seismic phases at the global scale, and
observed seasonal variations using a dynamic ocean wave model
(Kedar et al. 2008; Hillers et al. 2012). Nishida & Takagi (2016)
also recently suggested that SV and SH waves generated by distant
storms can be extracted though their amplitudes are much lower
than the P waves.

The beamforming method described by Retailleau et al. (2015)
uses: (1) the multiplicity of the seismic phases that are generated
during an event and (2) the density of sensors of networks. It does
not need any *a priori* information about the source, which makes it
possible to study continuous data. In this paper, we add a parameter
to the method to consider the relative variation of amplitude of
the different phases with distance. We then perform an analysis of
12 months of data recorded by the transportable component of the
USArray in 2010, to study the consistency of the seismic waves that
reach the array, in terms of space, time and frequency.

We particularly highlight a high-amplitude event that occurred
on 2010 February 1, and we demonstrate that it originated along the
South East Indian Ridge. We compare our results with models of
the power spectral densities of the pressure field generated near the
ocean surface (Tolman 2009; Ardhuin et al. 2011) that is associated
with ocean wave-wave interactions (Longuet-Higgins 1950; Has-
selmann 1963). The analysis of regional data also shows Rayleigh
waves traveling from the same location.

## 2 DATA

We use the 2010 vertical-component records from the Transportable
Array component of the USArray network, which in 2010 included
roughly 400 receivers (black points in Fig. 1). The raw data were
pre-processed on each individual day by correcting them from their
instrumental response and suppressing data that contained possible
glitches and gaps. We kept the data of 200 stations shown in Fig. 1
with the red points, to shorten the array aperture, and thus the ray
parameter variations, across the network. A ray parameter varia-
tion between the stations would reduce significantly the coherence
across the array and thus the array-analysis quality. Moreover, us-
ing only the 200 stations in Fig. 1 (red points) allows us to reduce
the azimuthal variation of the array aperture (Fig. S1, Supporting
Information S1).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Map of the USArray stations available on 2010 February 1 (black). The 200 stations used in this study are in red.</figcaption></figure>

## 3 DETECTION AND LOCATION METHOD

We perform array processing to detect and characterize seismic
events, combining the information offered by a large network and
the different phases that are generated. This is done by comparing
the result of the Vespagram analysis to the different body-wave
phases that are expected for a specific event.

The location method is described by Retailleau et al. (2015).
It consists in a grid search on a set of potential source locations
to determine which of them matches the observed data. For each
potential location of the grid, we perform a two-step process on the
data.

As a first step, we construct a time-distance gather by comput-
ing the great-circle distances between the potential location and
all the network stations. We then compute a Vespagram (Rost &
Thomas 2002) to represent the time-distance gather as a function
E(t, s, r) of the traveltime t, apparent slowness s and potential source
location r. The Vespagram is smoothed in time with a 30-s slid-
ing window, to suppress high-frequency variations. The Vespagram
amplitude is normalized to emphasize the first arriving body-wave
phases, rather than surface waves, whose consistency in time and
space arrivals is much less precise. Retailleau et al. (2014) also
showed that the late-arriving phases (e.g. PPS, SSS and others) are
more influenced by the heterogeneities in their travel paths than the
earlier phases (e.g. PKP, PP and others). To perform the normal-
ization, we compute the average Vespagram over all the slownesses
and smooth this function in time within a 600-s sliding window. We
then divide the Vespagram by the resulting function. The slowness
resolution expected of each phase depends on the network geome-
try and can be characterized using the network response (Fig. S1,
Supporting Information S1). Around 10 s of period, we expect a
resolution of about 0.0065 s km⁻¹ (more details can be found in the
Supporting Information S1).

To further focus on body waves, as a second step, we design a
time-slowness filter of the different body-wave phases. The purpose
behind this filter is to determine if the arrivals observed in the
Vespagram correspond to expected arrivals in a one dimension Earth
model. The time-slowness pairs are generated using Crotwell et al.
(1999) and Buland & Chapman (1983) methods in the IASP91
Earth model (Kennett & Engdahl 1991). The filter (in the time-
slowness domain) consists of a set of rectangular areas centered at
the arrival times and the slowness of the selected phases. The size of
these rectangles is 20 s in time and 0.01 s km⁻¹ in slowness. It was
determined empirically for frequency bands of analysis between 0.01 and 0.3 Hz. For each source location, the filter is continuously shifted over the Vespagram along the time coordinate. The time-slowness energy is summed within the filter rectangles, to get the final response of the network $\tilde{E}(t, r)$ as a function of time $t$ and the source location $r$.

Along a given azimuth, the Vespagrams are the same for different distances. This could potentially lead to matching the filter to a wrong phase arrival on the Vespagram and thus to a mislocation of the source. To reduce the ambiguity of the source location, we add a new parameter to the method detailed in Retailleau et al. (2015) that accounts for the seismic waves amplitude (details can be found in the Supporting Information S2). While the phase filter is shifted over the Vespagram, a weight is applied to each phase, and thus to each rectangle of the filter. The weights represent each phase amplitude as a function of epicentral distance (Supporting Information S2). To calculate these weights, we compute synthetic seismograms with the spectral-element method AXISEM (Nissen-Meyer et al. 2014), using the isotropic version of the Preliminary Reference Earth Model (Dziewonski & Anderson 1981), using a single-point force as a source (details in Supporting Information S2). In future applications on different types of sources, the source mechanism for building the synthetics should be modified.

Performing a grid search on the Earth surface, we test different source positions to locate potential events in time and space. When the test position corresponds to a real source of location $r_0$, the amplitude of $\tilde{E}(t, r = r_0)$ shows a maximum at the arrival time of the first phase arrival. Several examples of earthquake detection and location are described in Retailleau et al. (2015).

In this study, we choose to perform the analysis after filtering the pre-processed data on different frequency bands. This leads to explore the frequency content of the events under examination. In order to analyse all possible locations and several frequency bands, it is necessary to compute a Vespagram for every combination of location and frequency, making the process heavily time consuming. Moreover, the analysis is performed on a large amount of data (one year). We thus employ parallel computation to perform our analysis. Additionally, in the one-year continuous analysis presented in this paper, we constrain the depth at zero, which results in biased estimations for deep earthquakes. To analyse the depth variations, it is necessary to add the depth as a supplementary parameter, increasing furthermore the computational time.

The process is performed here continuously over one year of data to obtain a function $\tilde{E}(t, r, f)$ of time of the first arriving phase, location of the event and frequency. In this analysis of the 2010 USArray data, we test 3173 potential sources (ensemble of locations $R$), as represented in Fig. 2(a) (with a horizontal spacing between sources of 250 km). The time-series are filtered into 71 frequency bands between 0.0205 and 0.225 Hz (distributed logarithmically in frequency), to observe the events variations of behaviour as opposed to the time, space and frequency. We expect to observe different types of events, such as earthquakes, microseisms and low-frequency volcanic tremors. At the end of the process, we obtain the function $\tilde{E}(t, r, f)$ as a function time $t$, keeping a sampling of 1 point per 10 s, location $r$ and frequency $f$.

We plot the resulting $\tilde{E}^R(t, f)$ as a function of time and frequency after obtaining the maximum values over the different locations $R$, following eq. (1):

$$ \tilde{E}^R(t, f) = \max_{r \in R}[\tilde{E}(t, r, f)] \quad (1) $$

The result is smoothed over 1-hr time windows, as follows:

$$ \tilde{E}_T^R(t, f) = \frac{1}{N} \int_{t-\frac{T}{2}}^{t+\frac{T}{2}} \tilde{E}^R(t, f) dt \quad (2) $$

where $N$ is the number of points in time and $T=1$ hr is the smoothing window length.

The function $\tilde{E}_{1h}^R(t, f)$ is represented in Fig. 2(b). The amplitude tends to increase with decreasing frequency because of the constant rectangle size in the body-wave filter (described earlier in this section). The phase arrivals have shorter duration at higher frequencies compared to the filter box size. As a consequence, the time averaging weakens the energy at high frequencies. There might also be some residual contamination from surface waves.

Fig. 2(b) highlights different types of events. The vertical lines show several earthquakes. These events generate energy at all frequencies, and have short durations. Other events also appear with narrower frequency ranges and longer source lengths. These events may be associated with the interactions between the ocean and the solid Earth or with low-frequency tectonic or volcanic events. In this paper, we do not analyse all different events observed in Fig. 2(b). However, a short comparison can be made with the list of events computed by Obrebski et al. (2013), who applied a beamforming technique to the Southern California Seismic Network data also recorded in 2010. Some events listed by Obrebski et al. (2013) have been also observed in this study by using the USArray data [two examples are highlighted in Fig. 2(b) with continuous pink circles], some others have not been observed [one example is highlighted in Fig. 2(b) with a dashed pink circle]. On the other hand, some events were detected here using the USArray analysis, but are not in the list of Obrebski et al. (2013), because they only list events within 30°–90° of epicentral distance from the Southern California Seismic Network.

An advantage of the method is the detection of events with no a priori information about the parameters of their source, which allows detection of possible non-catalogued events. This is the case with the energetic event that occurred on 2010 February 1. It appears as the most energetic coherent signal recorded by USArray during the period of analysis for periods about 10 s (Fig. 2b). In this paper, in order to exploit the capacity of the new method to detect emergent signals, we analyse in detail this long-duration high-amplitude event.

## 4 CHARACTERIZATION OF THE 2010 FEBRUARY 1, EVENT USING ANTENNA ANALYSIS

The analysis of the USArray data highlights several long-lasting events (of about 1 d) that have not been catalogued, and in particular an event that occurred on 2010 February 1. Our first global analysis suggests that the source is located in the southern Indian ocean. To describe this event more precisely, we analyse the results for 2 d of data centred on 2010 February 1. We perform a narrower grid search in an area around the South East Indian Ridge (ensemble of locations $R$ in Fig. 3a), and with frequencies $f$ between 0.011 and 0.315 Hz. This leads to a new function $\tilde{E}(t, r, f)$ spatially focused to analyse the selected event. This study focuses on $P$ waves as it has been demonstrated by polarization analysis that body-wave microseismic noise is dominated by $P$ waves on the vertical component (e.g. Landès et al. 2010; Zhang et al. 2010; Gualtieri et al. 2014).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) The source locations *R* tested are represented by the black points. (b) Function $\tilde{E}_{1h}^R(t, f)$ is obtained for 2010 as a function of time and frequency. The arrows show notable earthquakes and the rectangles show other notable increases in energy. The continuous pink circles represent events detected by Obrebski *et al.* (2013) and the dashed pink circle an example of event not detected by Obrebski *et al.* (2013).</figcaption></figure>

## 4.1 Time-frequency dependence

Fig. 3(b) shows the time and frequency variations of the function $\tilde{E}_{1h}^R(t, f)$ obtained from $\tilde{E}(t, r, f)$ using eqs (1) and (2).

We observe an increase in energy around the beginning of February 1 that lasts for about a day. This strong energy is concentrated around 0.1 Hz and is consistent in time.

More precisely, Fig. 3(c) shows $\tilde{E}_{1h}^R(t, f)$ over two frequency bands. The black curve shows the mean over all of the frequencies ([0.011–0.315] Hz) and the grey curve over a selection of

frequencies ([0.09–0.11] Hz). The event strongly emerges in the [0.09–0.11] Hz frequency band. Considering the broad frequency band (black curve), we only observe short peaks, which correspond to five earthquakes of magnitude over 4.5.

## 4.2 Event location

We obtain a function that varies only according to the location. We first average $\tilde{E}(t, f, r)$ over the frequencies 0.09–0.11 Hz, the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) The ensemble of source locations *R* tested for the thinner grid search are represented by the white dots. (b) Time-frequency energy $\tilde{E}_{1h}^R(t, f)$ as a function of time and frequency after obtaining the maximum of the energy over all of the locations tested. (c) Results shown in (b) after averaging $\tilde{E}_{1h}^R(t, f)$ over the two frequency bands of [0.011–0.315] Hz (black) and [0.09–0.11] Hz (grey) (delimited by the white rectangle in (b). In (b) and (c), the dashed lines show the earthquakes of magnitudes >4.5 that occurred during the selected time (magnitude 4.5, yellow; 5, orange; 6, purple). These events were obtained from the CMT catalogue (Dziewonski & Anderson 1981; Ekström *et al.* 2012).</figcaption></figure>

range of frequencies highlighted in Figs 3(b) and (c). Then, for each location, we obtain the value of the highest amplitude over the length of the event:

$$ \tilde{E}^{F,T}(r) = \max_{t \in T} \left[ \frac{1}{N_f} \sum_F \tilde{E}(r, t, f) \right] \quad (3) $$

where *T* is the length of the event, *F* = 0.09–0.11 Hz is the frequency range, and *N*<sub>*f*</sub> is the number of frequencies stacked. The resulting function $\tilde{E}^{F,T}(r)$ is represented in Fig. 4. We locate the event along the South East Indian Ridge, South of Australia, at about 150° dis-

tance and 236° azimuth with respect to the USArray. This location remains consistent for the whole duration of the event. The location spot, which is spread on over 1000 km, is wider than the resolution expected with the array response, which is of about 250 km, suggesting that the source might be spread.

## 4.3 Slowness observations

The Vespagram *E*(*t*, *s*, *r*) is computed as explained in Section 3, and shown in Fig. 5. The location *r* is the location obtained with
<figure><img src="image_6.png" /><figcaption>Figure 4. Location results $\tilde{E}^{F,T}(r)$ of the antenna analysis. Ridges are shown with the thin grey lines, coasts with thick grey lines. The maximum of amplitude is represented with the black circle.</figcaption></figure>

$\tilde{E}^{F,T}(r)$, at an azimuth of 236° and a distance of 150° with respect to the USArray. The signals were again filtered between 0.09 and 0.11 Hz, the frequencies that contain the event (Fig. 3). Emerging energy is observed during 1 d, as seen in Fig. 3. The earthquakes are not visible in the Vespagram because the length of their occurrence is short compared to the 2 d studied.

The Vespagram shows that the emergent signal has a stable range of slowness over the whole day (Fig. 5). This range is relatively wide, and goes from 0.01 to 0.03 s km⁻¹. The time-slowness couples predicted are represented by the seismic phases PKPab, PP and PPP as white dots in Fig. 5. These suggest that the waves traveled to the network mainly as PKPab. Lower increase in energy is observed at larger slownesses, which can be associated with the PP or PPP phases. The wide range of slownesses observed might indicate perturbations of the medium along the travel path.

In summary, the antenna analysis shows a high-amplitude event that occurred on 2010 February 1, with a dominant frequency of 0.1 Hz along the South East Indian Ridge.

## 5 CHARACTERIZATION OF THE 2010 FEBRUARY 1 EVENT USING A REGIONAL DATA SET

There is no seismic event catalogued at the time and location described in the previous section. To confirm its occurrence and location, and obtain further information about the event, we analysed the vertical-component records by 20 stations located around the source (Fig. 6a). Fig. 6(b) shows the vertical-component seismograms filtered between 0.09 and 0.11 Hz, for 2010 February 1, as a function of the distance to the source location (estimated in the previous section). On average, we see that the signal amplitude decreases when the stations are further away from the source.

The amplitude of a seismogram can be approximated as a simple function of the distance to the source:

$$A_i = A_0 \frac{e^{-Br_i}}{r_i^n} \quad (4)$$

With,

$$B = \frac{\pi f}{Q\beta} \quad (5)$$

where $A_0$ is the source-radiated seismic amplitude, $A_i$ is the seismic amplitude recorded at receiver i, $r_i$ is the distance between receiver i and the source, $f=10$ s is the frequency, $\beta$ is the shear wave velocity and $Q=250$ is the attenuation quality factor. The exponent $n$ is $n=1$ for body waves and $n=0.5$ for surface waves. In considered distance and frequency ranges, we assume that the signals are mainly composed of surface waves (and thus $n=0.5$).

The amplitude ratio between the different receivers is:

$$\frac{A_i}{A_j} = \left(\frac{r_j}{r_i}\right)^n \exp(-B[r_i - r_j]) \quad (6)$$

which is independent of the source seismic intensity $A_0$, and where i and j denote a pair of receivers. The amplitude difference is a direct function of the distances of the stations to the source. The analysis of the amplitude ratios between the different receivers then leads to an approximate location of the event, with regional data (Battaglia & Aki 2003; Taisne et al. 2011).

The aim is to minimize the difference between the amplitude ratio observed and the amplitude ratio predicted with the distance

<figure><img src="image_7.png" /><figcaption>Figure 5. Vespagram $E(t, f = [0.09–0.11]Hz, r_0)$ obtained after filtering the data between 0.09 and 0.11 Hz for a source $r_0$ located at 236° azimuth and 150° epicentral distance with respect to USArray. Theoretical time-slowness arrivals at this distance are showed with white dots for the seismic phases PKPab, PP and PPP.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6. (a) Map of the stations used for the regional data analysis. (b) Signals of the different stations on 2010 February 1, as a function of their distance to the averaged location obtained. (c) Mean location obtained over all of the events. (d) Location obtained at the beginning of the event (at 03:00 hr). (e) Location obtained at the end of the event (at 18:00 UTC). The magenta contour shows the 80 per cent maximum energy location of the Vespagram analysis method.</figcaption></figure>

ratio (eq. 6). Hence, we perform a grid search to compute the misfit function of the difference between the two members in eq. (8) for each pair of stations i and j, as a function of the possible source location for different time windows:

$$M_{ij}(x, y) = \exp \left( - \left| \frac{A_i^{\text{obs}}}{A_j^{\text{obs}}} - \frac{A_i^{\text{pred}}}{A_j^{\text{pred}}} \right| \right) \quad (7)$$

$$= \exp \left( - \left| \frac{A_i^{\text{obs}}}{A_j^{\text{obs}}} - \left( \frac{r_j}{r_i} \right)^n \exp(-B[r_i - r_j]) \right| \right) \quad (8)$$

where $A_i$ and $A_j$ are the signal amplitudes at receivers i and j, respectively, at the time of maximum energy.

All $M_{ij}$ obtained from different station couples are stacked to obtain the event location. The results are showed in Fig. 6.

Figs 6(c)–(e) show the location obtained from eq. (8) as intensity. Both the results averaged over the whole event (Fig. 6c) and over the end of the event only (Fig. 6e) show a high amplitude close to the source location obtained using body waves recorded by the US-Array, although it is significantly shifted to the east. The amplitude observed at the beginning of the event (Fig. 6d) is located further north. These results confirm the occurrence of a microseismic event in that region, but they do not give a precise indication of the source location. The stations used for the surface-wave-based location are very spread (Fig. 6a), and the seismic waves reaching them propagate through different media that can be only oceanic crust or partly continental crust. Therefore, the uncertainties of the results are very large and the obtained location can only be used as an approximate.

## 6 DISCUSSION

The analysed event is interesting because it generates continuously strong seismic signal during a period of about 20 hr. These signals have an average amplitude over 10 s segments that is of the same order as the one produced by the magnitude 6 earthquake that occurred in the same area at the end of the day in the Solomon Islands (Fig. 6). Moreover, Fig. 3(c) shows that in the frequency band [0.011–0.315] Hz, the waves traveling from the earthquake and this event appear as coherent on the USArray. Although it is difficult to compare the broad-band radiation of earthquakes with the narrow-band radiation of this observed microseism, Figs 3(b) and (c) show that the February 1 event produced a large-energy radiation, which appears as coherent as earthquakes of magnitude up to 6.

Different hypotheses can be defined concerning the origin of the phenomenon observed in this study. The location presented in Fig. 4 might suggest that the origin of the event is linked to the South East Indian Ridge activity, such as underwater volcanism. It has indeed been shown that volcanism could generate low-frequency events (around 1 Hz) and even very-long-period events that have periods ranging from 3 to 100 s or longer (see McNutt 2005, for a review on this subject). Xia et al. (2013) thus argued that the 26 and 28 s microseisms originate from volcanism. Another possibility is that the observed phenomenon is linked to microseisms generated by ocean wave-wave interactions (Longuet-Higgins 1950).

Storm winds generate gravity waves at the ocean surface. The non-linear interaction among ocean waves that come from nearly opposite directions and have similar frequency generates a pressure that is transmitted at the ocean bottom and converted to elastic waves (e.g. Kedar 2011; Ardhuin & Herbers 2013). This gen-
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Ocean wave height computed using ocean wave model WAVEWATCH III(R) (Tolman 2009; Ardhuin et al. 2011): (a) January 31 at 12:00 UTC. (b) February 1 at 00:00 UTC. (c) February 1 at 12:00 UTC. (d) February 2 at 00:00 UTC. The black circle represents the location obtained with the USAarray.</figcaption></figure>

erated signal is twice the frequency of the interacting ocean waves (Longuet-Higgins 1950, 1952; Hasselmann 1963) and it constitutes a secondary microseism with periods between 3 and 10 s.

Body waves have been first identified in the secondary microseism by Backus et al. (1964), Toksöz & Lacoss (1968), Haubrich & McCamy (1969), Vinnik (1973), Koper & de Foy (2008) and have been associated to specific storms (Gerstoft et al. 2008; Farra et al. 2016). To investigate the hypothesis that a storm generated the observed body waves, we compare the results obtained with our analysis of the USAarray data with the output of the numerical ocean wave model WAVEWATCH III(R) (Tolman 2009; Ardhuin et al. 2011). The first quantitative modeling of surface waves from secondary microseisms using ocean wave model hincasts was performed by Kedar et al. (2008).

Fig. 7 shows the ocean wave height from January 31 to February 2 (Ardhuin et al. 2011); available at ftp://ftp.ifremer.fr/ifremer/cersat/products/gridded/wavewatch3/pub/HINDCAST/NC/2010). We observe the presence of a wave height increase along the South East Indian Ridge, close to the locations obtained from the seismic signals using the USAarray and the method in Section 3. These wave heights do not explain the seismic signals, however they suggest the occurrence of a storm in the region. The wave heights approximately correspond to the geopotential height reanalysis during this period time (Supporting Information S3), which confirms that our event should be linked to a relatively circular and close-shaped low atmospheric pressure in the region, that is, to a mid-latitude cyclone.

We observe that the location of the maximum ocean wave height moves from west to east, whereas our source remains relatively stable in time.

The generation of the secondary microseisms is not directly linked to the ocean wave height but to the ocean wave-wave interaction. Fig. 8 shows the power spectral density of the pressure field linked to the wave interactions for ocean waves at 0.055 Hz (i.e. seismic waves at 0.11 Hz, eq. 2 in Ardhuin et al. 2011) generated with the numerical ocean wave model WAVEWATCH III(R)

(Tolman 2009; Ardhuin et al. 2011), with (Figs 8a-c) and without (Figs 8d-f) reflection at the coast. We corrected this result by applying the coefficient that accounts for the effect of the ocean layer on body waves generated by a source at the ocean surface (called ‘ocean source site effect’ in Gualtieri et al. (2014). This ocean source site effect was computed at 0.1 Hz for a take-off angle related to PKPab waves at these epicentral distances (eqs 4–8 in Gualtieri et al. 2014). We analyse the power spectral density of the pressure field linked to the interaction of ocean waves corrected by the source site effect related to the arrival of the PKPab phase. Fig. 8 shows the power spectral density of the pressure field with and without reflection at the coast on January 31 at 12:00, February 1 at 00:00 and February 1 at 12:00. The power spectral density of the pressure field calculated every 3 hr with reflection at the coast can be found in Fig. S7 in the Supporting Information S4.

The maximum wave height moves eastward along the ridge (Fig. 7). Contrarily, the amplitude of the power spectral density of the pressure field is low on January 31 and increases at the beginning of February 1, when it is located South of Australia (Fig. 8 and Fig. S7, Supporting Information). This observation is consistent with the results obtained with the USAarray data (Fig. 4 and previous section). Moreover, the amplitude is much higher when the reflection at the coast is added to the model (Figs 8d-f). This implies that the signals might have been generated by the interactions between ocean waves reflected at the coast and ocean waves traveling to it. This could explain why no seismic signal is observed west from the main source, although an increase in wave height is observed.

The small eastward movement of the pressure field (Fig. 8) on January 31 is not observed using the USAarray data (Fig. 4). The seismic amplitudes generated might be too small on January 31 to be resolved by our array-analysis method. In Figs 4(a)–(c), we do not consider coastal reflection, while in Figs 4(d)–(f), we consider a large amount of coastal reflection (10 per cent). This value is larger than values found so far (e.g. Stutzmann et al. 2012) and therefore can well represent an upper bound for the reflection coefficient.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Combination, in logarithmic scale, of the power spectral density of the pressure (Pa²m⁻²s) at the ocean surface related to the ocean wave-wave interaction for waves at 0.055 Hz (i.e. elastic waves at 0.11 Hz) (a)–(c) with the ocean source site effect acting on PKP waves without reflection and (d)–(f) with 10 per cent reflection at the coast. (a) and (d) January 31 at 12:00 UTC. (b) and (e) February 1 at 00:00 UTC. (c) and (f) February 1 at 12:00 UTC. The grey shadow represents inland regions and distances at which the PKPab phase is not seen by the USArray.</figcaption></figure>

At the end of February 1, the power spectral density of the pressure field due to the wave-wave interaction weakens and disappears when moving eastward. This is because the PKP phase does not exist at epicentral distances smaller than 145° (grey shadow in Fig. 8). When we correct the pressure field due to wave-wave interactions with the source site effect we use the take-off angle of the PKP phase, which is the strongest phase generated by the event. This implies that we have no model observation towards the east (highlighted by the grey shadow in Fig. 8 and Fig. S7, Supporting Information). This matches well the location obtained with the USArray. On the other hand, with the regional data set, surface waves are still observed to the east, which might explains the shift of the high-amplitude peak to the east observed with the regional stations compared to the observation with USArray.

These observations suggest that the event analysed from the USArray network is linked to the generation of seismic waves by a storm. However, some questioning persists concerning the dominant frequency of the event. Indeed, the results obtained with the US-Array data show signals at a quite narrow frequency band (around 0.1 Hz). The power spectral density of the pressure field corrected by the source site effect for the PKPab phase gives a location in agreement with our results, considering ocean waves at 0.055 Hz (i.e. seismic waves at 0.11 Hz), but it also shows higher amplitudes of the pressure field at higher frequencies, which are not observed with the USArray. This effect might be linked to attenuation along the wave path. Furthermore, the 10 per cent coefficient of reflection at the coast is probably too high and this is why some of the discrepancies persist between the model and observations. For example, Stutzmann et al. (2012) observed at station (CAN), a maximum reflection coefficient of 6 per cent. In this paper, we compare two extreme cases, a case without reflection and a case with a high re-

flection coefficient (10 per cent) to observe that reflection is needed to explain better the observations.

## 7 CONCLUSIONS

The pressure field associated with the ocean wave-wave interaction shows a good agreement with the results obtained with the array analysis and only some minor discrepancies. For this reason, we conclude that the generation of the observed signals by a storm is the most probable hypothesis.

The simultaneous antenna analysis of multiple body-wave phases efficiently detects and locates various sources. It can therefore be used to scan continuous records from large networks. This leads to a catalogue of previously undetected events to be created. Most studies on particular events have focused on the most well-known climate events (e.g. Gerstoft et al. 2006), thus a catalogue would permit to characterize the signal generation from meteorological events in more various locations on Earth, and to understand more globally their behaviour.

## ACKNOWLEDGEMENTS

While Jocelyn Guilbert contributed to the development of the approaches presented in this paper, the authors regret that he was not able to read the last versions of the manuscript. We thank Alexandre Fournier and Tarje Nissen-Meyer for their help on the computing of synthetics with Axisem. Data from the Transportable Array network are made freely available as part of the EarthScope USArray facility, operated by Incorporated Research Institutions for Seismology. Numerical computations were performed on the S-CAPAD
platform, IPGP, France. This work was supported by the Commisariat à l'Energie Atomique (CEA, France) and by the European Research Council under contract FP7 ERC Advanced grant 227507 (WHISPER), and by the project DataScale funded by the French BigData program. The work of NMS was supported by the Russian Science Foundation (grant no. 14-47-00002). ISTERe is part of Labex OSUG@2020. LG acknowledges support from a Lamont-Doherty Earth Observatory Postdoctoral Fellowship. This is the Lamont-Doherty Contribution Number 8148. We thank the editor and two anonymous reviewers for their useful comments.

## REFERENCES

Ardhuin, F. & Herbers, T.H., 2013. Noise generation in the solid earth, oceans and atmosphere, from nonlinear interacting surface gravity waves in finite depth, *J. Fluid Mech.*, **716**, 316–348.

Ardhuin, F., Stutzmann, E., Schimmel, M. & Mangeney, A., 2011. Ocean wave sources of seismic noise, *J. geophys. Res.*, **116**, C09004, doi:10.1029/2011JC006952.

Backus, M., Burg, J., Baldwin, D. & Bryan, E., 1964. Wide-band extraction of mantle p waves from ambient noise, *Geophysics*, **29**(5), 672–692.

Battaglia, J. & Aki, K., 2003. Location of seismic events and eruptive fissures on the Piton de la Fournaise volcano using seismic amplitudes, *J. geophys. Res.*, **108**, 2364, doi:10.1029/2002JB002193.

Buland, R. & Chapman, C., 1983. The computation of seismic travel times, *Bull. seism. Soc. Am.*, **73**(5), 1271–1302.

Bungum, H., Husebye, E.S. & Ringdal, F., 1971. The norsar array and preliminary results of data analysis, *Geophys. J. Int.*, **25**(1–3), 115–126.

Crotwell, H. & Owens, T.J.R., 1999. The TauP Toolkit: flexible seismic travel-time and ray-path utilities, *Seismol. Res. Lett.*, **70**(April), 154–160.

Dziewonski, A. & Anderson, D., 1981. Preliminary reference Earth model, *Phys. Earth planet. Inter.*, **25**, 297–356.

Ekström, G., 2006. Global detection and location of seismic sources by using surface waves, *Bull. seism. Soc. Am.*, **96**(4A), 1201–1212.

Ekström, G., Nettles, M. & Abers, G.A., 2003. Glacial earthquakes, *Science*, **302**(5645), 622–624.

Ekström, G., Nettles, M. & Dziewoski, A., 2012. The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes, *Phys. Earth planet. Inter.*, **200–201**, 1–9.

Farra, V., Stutzmann, E., Gualtieri, L., Schimmel, M. & Ardhuin, F., 2016. Ray-theoretical modeling of secondary microseism p waves, *Geophys. J. Int.*, **206**(3), 1730–1739.

Gerstoft, P., Fehler, M.C. & Sabra, K.G., 2006. When Katrina hit California, *Geophys. Res. Lett.*, **33**, L17308, doi:10.1029/2006GL027270.

Gerstoft, P., Shearer, P.M., Harmon, N. & Zhang, J., 2008. Global P, PP, and PKP wave microseisms observed from distant storms, *Geophys. Res. Lett.*, **35**, L23306, doi:10.1029/2008GL036111.

Green, P., Frosch, R. & Romney, C., 1965. Principles of an experimental large aperture seismic array (LASA), *Proc. IEEE*, **53**(12), 1821–1833.

Gualtieri, L., Stutzmann, E., Farra, V., Capdeville, Y., Schimmel, M., Ardhuin, F. & Morelli, A., 2014. Modelling the ocean site effect on seismic noise body waves, *Geophys. J. Int.*, **197**(2), 1096–1106.

Hasselmann, K., 1963. A statistical analysis of the generation of microseisms, *Rev. Geophys.*, **1**(2), 177–210.

Haubrich, R.A. & McCamy, K., 1969. Microseisms: coastal and pelagic sources, *Rev. Geophys.*, **7**(3), 539–571.

Hillers, G., Graham, N., Campillo, M., Kedar, S., Landès, M. & Shapiro, N.M., 2012. Global oceanic microseism sources as seen by seismic arrays and predicted by wave action models, *Geochem. Geophys. Geosyst.*, **13**, Q01021, doi:10.1029/2011GC003875.

Ishii, M., Shearer, P.M., Houston, H. & Vidale, J.E., 2005. Extent, duration and speed of the 2004 Sumatra-Andaman earthquake imaged by the Hi-Net array, *Nature*, **435**(7044), 933–936.

Kedar, S., 2011. Source distribution of ocean microseisms and implications for time-dependent noise tomography, *C. R. Geosci.*, **343**(8–9), 548–557.

Kedar, S., Longuet-Higgins, M., Webb, F., Graham, N., Clayton, R. & Jones, C., 2008. The origin of deep ocean microseisms in the North Atlantic Ocean, *Proc. R. Soc. A: Math. Phys. Eng. Sci.*, **464**(2091), 777–793.

Kennett, B. & Engdahl, E., 1991. Traveltimes for global earthquake location and phase identification, *Geophys. J. Int.*, **105**(2), 429–465.

Koper, K.D. & de Foy, B., 2008. Seasonal anisotropy in short-period seismic noise recorded in South Asia, *Bull. seism. Soc. Am.*, **98**(6), 3033–3045.

Landès, M., Hubans, F., Shapiro, N.M., Paul, A. & Campillo, M., 2010. Origin of deep ocean microseisms by using teleseismic body waves, *J. geophys. Res.*, **115**, B05302, doi:10.1029/2009JB006918.

Longuet-Higgins, M.S., 1950. A theory of the origin of microseisms, *Phil. Trans. R. Soc. Lond.*, **243**(857), 1–35.

Longuet-Higgins, M.S., 1952. On the statistical distributions of sea waves, *J. Mar. Res.*, **11**(3), 245–265.

McNutt, S.R., 2005. Volcanic seismology, *Annu. Rev. Earth Planet. Sci.*, **32**, 461–491.

Nishida, K. & Takagi, R., 2016. Teleseismic s wave microseisms, *Science*, **353**(6302), 919–921.

Nissen-Meyer, T., van Driel, M., Stähler, S.C., Hosseini, K., Hempel, S., Auer, L., Colombi, a. & Fournier, a., 2014. AxiSEM: broadband 3-D seismic wavefields in axisymmetric media, *Solid Earth*, **5**(1), 425–445.

Obara, K. & Ito, Y., 2005. Very low frequency earthquakes excited by the 2004 off the kii peninsula earthquakes: a dynamic deformation process in the large accretionary prism, *Earth Planets Space*, **57**(4), 321–326.

Obrebski, M., Ardhuin, F., Stutzmann, E. & Schimmel, M., 2013. Detection of microseismic compressional (P) body waves aided by numerical modeling of oceanic noise sources, *J. geophys. Res.*, **118**(8), 4312–4324.

Retailleau, L., Shapiro, N.M., Guilbert, J., Campillo, M. & Roux, P., 2014. Antipodal focusing of seismic waves observed with the USArray, *Geophys. J. Int.*, **199**(2), 1030–1042.

Retailleau, L., Shapiro, N.M., Guilbert, J., Campillo, M. & Roux, P., 2015. Detecting and locating seismic events with using USArray as a large antenna, *Adv. Geosci.*, **40**, 27–30.

Rost, S. & Thomas, C., 2002. Array seismology: methods and applications, *Rev. Geophys.*, **40**(3), 1008, doi:10.1029/2000RG000100.

Satriano, C., Dionicio, V., Miyake, H., Uchida, N., Vilotte, J.-P. & Bernard, P., 2014. Structural and thermal control of seismic activity and megathrust rupture dynamics in subduction zones: lessons from the Mw 9.0, 2011 Tohoku earthquake, *Earth planet. Sci. Lett.*, **403**, 287–298.

Shearer, P.M., 1994. Global seismic event detection using a matched filter on long-period seismograms, *J. geophys. Res.*, **99**(B7), 13 713–13 725.

Stutzmann, E., Ardhuin, F., Schimmel, M., Mangeney, A. & Patau, G., 2012. Modelling long-term seismic noise in various environments, *Geophys. J. Int.*, **191**(2), 707–722.

Taisne, B., Brenguier, F., Shapiro, N.M. & Ferrazzini, V., 2011. Imaging the dynamics of magma propagation using radiated seismic intensity, *Geophys. Res. Lett.*, **38**, L04304, doi:10.1029/2010GL046068.

Toksöz, M.N. & Lacoss, R.T., 1968. Microseisms: mode structure and sources, *Science*, **159**(3817), 872–873.

Tolman, H.L., 2009. User manual and system documentation of WAVE-WATCH III TM version 3.14, Technical note, MMAB Contribution, 276: 220.

Vinnik, L., 1973. Sources of microseismic p waves, *Pure appl. Geophys.*, **103**(1), 282–289.

Xia, Y., Ni, S. & Zeng, X., 2013. Twin enigmatic microseismic sources in the Gulf of Guinea observed on intercontinental seismic stations, *Geophys. J. Int.*, **194**(1), 362–366.

Zhang, J., Gerstoft, P. & Shearer, P.M., 2010. Resolving P-wave travel-time anomalies using seismic array observations of oceanic storms, *Earth planet. Sci. Lett.*, **292**(3–4), 419–427.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.
**Figure S1.** (a) USArray, plotted in cartesian coordinates. (b) Array response as a function of the wavenumber (c) zoom in the white square shown in (b).

**Figure S2.** Beamforming computed for an azimuth to the USArray network of 236 degrees. The circles show the three first arriving phases for two potential distances to the network.

**Figure S3.** Synthetic signals computed for a vertical force at the surface as a function of epicentral distance and time. The different theoretical time arrivals of the P-wave phases are plotted in blue (Crotwell *et al.* 1999).

**Figure S4.** Examples of the amplitude obtained in the synthetics compared to the distance for several phases represented in Fig. S3. *P* represents *P* and *P*<sub>diff.</sub>

**Figure S5.** Location results $\tilde{E}^{ET}(r)$ of the antenna analysis without a phase weighting.

**Figure S6.** (a–h) Evolution of the geopotential height at 1000 hPa with time (colour scale in [m]).

**Figure S7.** Combination of the power spectral density of the pressure ($Pa^2m^2s$) at the ocean surface related to the ocean wave-wave interaction for waves at 0.055 Hz (i.e. elastic waves at 0.11 Hz) with the ocean site effect acting on *PKP* waves and with 10 per cent reflection at the coast in logarithmic scale. The shaded region represent onland region and distance with no *PKPab* phase.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.