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

Seismological detection methods are traditionally based on picking techniques. These meth-ods 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 conﬁrms the occurrence of this event. We compare the seismic observations with models of the spectral density of the pressure ﬁeld 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.


I N T RO D U C T I O N
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 Global detection methods that use arrays of receivers have previously 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, making 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) identified body waves generated by the hurricane Katrina using beamforming analysis (Rost & Thomas 2002).Similarly, Gerstoft et al. (2008) used a beamforming method based on P-waves backprojection, 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 associated 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;Hasselmann 1963).The analysis of regional data also shows Rayleigh waves traveling from the same location.

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 variation between the stations would reduce significantly the coherence across the array and thus the array-analysis quality.Moreover, using 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).

D E T E C T I O N A N D L O C AT I O N M E T H O D
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 computing 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 sliding 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 normalization, 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 geometry 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 −1 (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 timeslowness 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 −1 in slowness.It was L. Retailleau et al. 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 timeslowness energy is summed within the filter rectangles, to get the final response of the network Ē(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 Ē(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 Ē(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 lowfrequency volcanic tremors.At the end of the process, we obtain the function Ē(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 Ē R (t, f ) as a function of time and frequency after obtaining the maximum values over the different locations R, following eq.( 1): The result is smoothed over 1-hr time windows, as follows: where N is the number of points in time and T = 1 hr is the smoothing window length.
The function Ē R 1h (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.

C H A R A C T E R I Z AT I O N O F T H E 2 0 1 0 F E B RUA RY 1 , E V E N T U S I N G A N T E N N A A N A LY S I S
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 Ē(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 bodywave 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).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.

Time-frequency dependence
More precisely, Fig. 3(c) shows Ē R 1h (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.

Event location
We obtain a function that varies only according to the location.We first average Ē(t, f, r ) over the frequencies 0.09-0.11Hz, the  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).
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: where T is the length of the event, F = 0.09-0.11Hz is the frequency range, and N f is the number of frequencies stacked.The resulting function Ē 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.

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 Ē 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 −1 .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.

C H A R A C T E R I Z AT I O N O F T H E 2 0 1 0 F E B RUA RY 1 E V E N T U S I N G A R E G I O N A L DATA S E T
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: With, 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, β 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: 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  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: 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.

D I S C U S S I O N
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 narrowband 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- 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 USArray 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 hindcasts 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 USArray 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 USArray 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 USArray data (Fig. 4).The seismic amplitudes generated might be too small on January 31 to be resolved by our array-analysis method.In  At the end of February 1, the power spectral density of the pressure field due to the wave-wave interaction weakens and disappear 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.

C O N C L U S I O N S
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.

A C K N O W L E D G E M E N T S
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     Figure S7.Combination of the power spectral density of the pressure (Pa 2 m 2 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) 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.

Figure 1 .
Figure 1.Map of the USArray stations available on 2010 February 1 (black).The 200 stations used in this study are in red.

Figure 2 .
Figure 2. (a) The source locations R tested are represented by the black points.(b) Function Ē R 1h (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).

Fig. 3
Fig. 3(b)  shows the time and frequency variations of the function Ē R 1h (t, f ) obtained from Ē(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 Ē R 1h (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

Figure 3 .
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 Ē R 1h (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 Ē R 1h (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).

Figure 4 .
Figure 4. Location results Ē 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.

Figure 5 .
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.

Figure 6 .
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.
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 8 .
Figure 8. Combination, in logarithmic scale, of the power spectral density of the pressure (Pa 2 m 2 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.

Figure S1 .
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).FigureS4.Examples of the amplitude obtained in the synthetics compared to the distance for several phases represented in Fig.S3.P represents P and P diff .FigureS5.Location results ĒF,T (r) of the antenna analysis without a phase weighting.

Figure S2 .
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).FigureS4.Examples of the amplitude obtained in the synthetics compared to the distance for several phases represented in Fig.S3.P represents P and P diff .FigureS5.Location results ĒF,T (r) of the antenna analysis without a phase weighting.

Figure S4 .
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).FigureS4.Examples of the amplitude obtained in the synthetics compared to the distance for several phases represented in Fig.S3.P represents P and P diff .FigureS5.Location results ĒF,T (r) of the antenna analysis without a phase weighting.

Figure S5 .
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).FigureS4.Examples of the amplitude obtained in the synthetics compared to the distance for several phases represented in Fig.S3.P represents P and P diff .FigureS5.Location results ĒF,T (r) of the antenna analysis without a phase weighting.

Figure S6 .
Figure S6.(a-h) Evolution of the geopotential height at 1000 hPa with time (colour scale in [m]).FigureS7.Combination of the power spectral density of the pressure (Pa 2 m 2 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) 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.