# New approach to detect seismic surface waves in 1Hz-sampled GPS time series

N. Houlié, G. Occhipinti, T. Blanchard, N. Shapiro, P. Lognonné, M.

Murakami

► To cite this version:

N. Houlié, G. Occhipinti, T. Blanchard, N. Shapiro, P. Lognonné, et al. New approach to detect seismic surface waves in 1Hz-sampled GPS time series. Scientific Reports, 2011, 1, 10.1038/srep00044. insu-03606462

HAL Id: insu-03606462

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

Submitted on 12 Mar 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>

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

# New approach to detect seismic surface waves in 1Hz-sampled GPS time series

SUBJECT AREAS:

GEOPHYSICS

SEISMOLOGY

GEOLOGY

EARTH SCIENCES

N. Houlié¹, G. Occhipinti², T. Blanchard³, N. Shapiro², P. Lognonné² & M. Murakami⁴

¹UC Berkeley, Seismological Laboratory, Berkeley, California, USA, ²Institut de Physique du Globe de Paris (Université Paris Diderot-CNRS), ³School of Earth and Environment, University of Leeds, LS2 9JT, Leeds, UK, ⁴Geographical Survey Institute, Tsukuba, Japan.

Received
3 February 2011

Accepted
1 July 2011

Published
22 July 2011

Correspondence and
requests for materials
should be addressed to
Nicolas Houlié.
(houlie@seismo.
berkeley.edu)

Recently, co-seismic seismic source characterization based on GPS measurements has been completed in near- and far-field with remarkable results. However, the accuracy of the ground displacement measurement inferred from GPS phase residuals is still depending of the distribution of satellites in the sky. We test here a method, based on the double difference (DD) computations of Line of Sight (LOS), that allows detecting 3D co-seismic ground shaking. The DD method is a quasi-analytically free of most of intrinsic errors affecting GPS measurements. The seismic waves presented in this study produced DD amplitudes 4 and 7 times stronger than the background noise. The method is benchmarked using the GEONET GPS stations recording the Hokkaido Earthquake (2003 September 25th, Mw=8.3).

Our capability to characterize both the seismic wave propagation and the seismic sources depends on the quality (accuracy and variety) of the datasets available. Since the 1970s, seismic networks were deployed to achieve these goals at the global scale¹,². The rise of broadband and very-broadband seismometers allowed tremendous progress. These networks made many discoveries possible, amongst which anisotropy, study of the core-mantle boundary, mapping of subducting plates in the mantle, etc. Nowadays, the rising cost of operations and the difficulty of accessing best locations for new sites, make it arduous to expand these networks much further. GPS receiver, with its cheaper installation and maintenance, offers an alternative to seismic sensors when it comes to the collection data in the range of periods from 10 to 200 seconds. Other advantages of GPS receiver are 1) that it does not have instrument response, 2) that it does not clip during strong motion and 3) that it provides ground motion time-series in displacement. The challenges of integrating the acceleration or velocity records³ are thus avoided. If the GPS waveforms were demonstrated to be as reliable as seismic waveforms, they could be integrated into global Earth tomographies at long-period (T=20 seconds or more).

Today, a few studies using GPS take advantages of GPS waveforms. Mostly, daily processing studies rather focus on the study of coseismic (first days)⁴ to post-seismic (first weeks to months) deformation⁵ in near- and far-field (distance larger than rupture length), and with remarkable results. These successes caused us to forget that the accuracy of the ground deformation measurements rely on the distribution of satellites. For positions computed over hours or days, the distribution of satellites is not an issue because the precise orbits help defining an average daily position using the satellites that are visible during the sidal day. Nevertheless, for shorter periods (minutes to hours), the lack of satellites in the sky might decrease the quality of the solutions. This effect can be seen when computing static displacement across very short sessions, but might be even more critical when GPS waveforms are computed. Yet recent studies showed that GPS receivers with high rate capabilities (e.g., 1Hz sampling rate) are sensitive to ground motion induced by seismic events including the 2002 Denali⁶⁻⁹, the 2003 Tokachi-Oki¹⁰,¹¹, the 2003 San Simeon¹², ¹³, the 2004 Sumatra-Andaman¹⁴, the 2008 Iwate-Miyagi Nairiku¹⁵, the 2010 Maule Earthquake¹⁶ earthquakes. In some cases, GPS observations collected were positively correlated to co-located seismic records.

The events detected by GPS stations are generally large ($M_w$ > 6.0). This is due to the reduced GPS sensitivity to ground displacement compared to seismometers, as well as to most networks being designed to survey continental deformation (with a spacing of 40 km or more). There is no theoretical limitation to the detection of small events ($M_w$ < 6.0). GPS phase aliasing (related to internal receiver filtering) appears to be the main parameter in reducing the sensitivity of GPS to describe ground vibrations¹⁷. To minimize this effect, operators converted their network to high sampling rate (1Hz or more). Today, a growing number of GPS receivers could be converted to 1Hz sampling rate (Figure 1) at minimal cost (change of telemetry), which encourages the support of seismic studies using long-period (Period T> 20s) GPS waveforms.
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 1** | Global high-rate GPS network operating today. Among the Global GPS network 1905 GPS sites providing data on a daily basis to SOPAC archive⁴², we highlight the GPS sites that belong to the NTRIP (red triangles) and IGS (black circles) networks that are recording freely available 1Hz data. For clarity reasons, the GEONET, Plate Boundary Observatory (1000+ sites) and BARD⁴³⁻⁴⁵ networks are not fully plotted on this map.</figcaption></figure>

Demonstrating that the GPS is at least as reliable as seismic sensors requires 1) to calibrate successfully GPS waveforms against seismic records, 2) to address computational issue(s) that could affect the GPS accuracies, and 3) to reliably detect small amplitude ground vibrations at long distance (500 km or more i.e. *distance that is comparable with previous works*⁷).

So far, some intrinsic errors have limited the reliability of records: the propagation of the GPS signal through the ionosphere and troposphere; satellite visibility; the clock delay; local site obstacles such as the multipath produced by reflection and/or diffraction of the GPS signal at the receiver site environment; and lastly, the accuracy of the displacement time-series. Because of these, GPS has failed to be integrated with seismic data into global tomographies¹⁸ yet. These issues, though, have been solved, and time-series of daily solutions produced, already. Most of the effects, related to the propagation GPS signals from the satellite to the site, are now eluded by the use of adequate data analysis based on a GPS network approach¹⁹,²⁰. These methods are based on the regional spatial filtering²¹ or/and the modified sidereal filtering (MSF)²².

In this study, we investigate the impact of a low number of satellites in the sky on the shape of GPS waveforms. We then test a new method to study seismic waves in 1Hz GPS time series based on the well-known computation of GPS double difference (DD) (see Methodology section). This allows for the detection of 3D ground shaking produced by surface wave propagation. The method is benchmarked using the GEONET²³,²⁴ GPS stations recording the Hokkaido Earthquake (2003 September 25th, Mw=8.3, Figure 2). The source mechanism was favourable for surface waves generation, and local static displacement at the epicentre was estimated at around 5–9 metres¹⁰,²⁵,²⁶.

## Results

Daily processing is efficient at eliminating most of the GPS uncertainties (see Introduction). Hence, we processed the GPS data collected during the hours following the 2003 Hokkaido event (the data processing is described in detail in the Methodology section). We then extracted the residuals (LC ionosphere-free) from the daily solutions for each line-of-sight (LOS), to recover the motion of each site. The phase residuals correspond to all the phase changes that cannot be explained by the motion of satellites above the network.

We use them to extract the short-lived disturbances of the phase received across the network. This methodology has already been applied to the study of volcanic atmospheric plumes in order to constrain the troposphere variations with respect to the standard atmospheric model used during the GPS processing²⁷,²⁸. In that case,

<figure><img src="image_6.png" /><figcaption>**Figure 2** | USGS Location of the main shock (black star). GEONET permanent sites are indicated with squares while International GPS Service (IGS) sites are symbolized using triangles. A black circle indicates the location of the GEOSCOPE station INU.</figcaption></figure>
the troposphere state changes were perturbing the GPS signals. As the GPS receivers were not moving during the explosions, the phase residuals were fully attributed to the presence of the volcanic plume. The problem we address here is similar in that the GPS phase is changed. But not by a volcanic plume; here it is the receiver's change of positions that affect the GPS phase. Each phase change is valid for a LOS (defined by the azimuth and elevation angles of one satellite from a given site).

The number of LOS residual histories equals the numbers of satellites visible at a single site. We use all the LOS residual histories available to compute the East, North and Vertical motion (E,N,U) at a given time:

$$ \begin{pmatrix} u_1 & v_1 & w_1 \\ u_2 & v_2 & w_2 \\ \vdots & \vdots & \vdots \\ u_n & v_n & w_n \end{pmatrix} \begin{pmatrix} E \\ N \\ U \end{pmatrix} = \lambda \begin{pmatrix} \phi_1 \\ \phi_2 \\ \vdots \\ \phi_n \end{pmatrix} $$

where the u, v and w are the components of the unit vector of each satellite LOS and φ the phase change in the direction of LOS (e.g.¹⁷). λ is the wavelength of the residuals carrier.

Following this approach, the quality of the displacement time-series inferred from the inversion of residuals in the line-of-sight (LOS) do rely on the distribution of the satellites in the sky (Figure 3 shows these distributions of satellites for various sites) at the time of the seismic wave arrival at the site. As the GPS system is based on the time arrival of the GPS signal at the receiver, if fewer satellites are located along a North-South azimuth the accuracy of the north-south component will be affected. The quality of the ground displacement solutions will be function of a parameter that describes the geometry of the tetrahedrons defined by 3 satellites and a site. Closer to 90 degrees the angles of tetrahedron will be, better the solutions will be. We build an operator that we name K, which is the average of all the mixed products of LOS unit vectors we can create at a site:

$$ K = \frac{1}{N} \sum_{1}^{N} (\vec{u} \times \vec{v}) \cdot \vec{w} \quad (1) $$

where N is the maximal number of triplets and $\vec{u}, \vec{v}, \vec{w}$ are the three unit vectors of three given satellites LOS. The equation (1) says that in the best conditions (elevation angle = 0; no masks), K=1.

<figure><img src="image_7.png" /><figcaption>Figure 3 | Skyview plot of the satellites above the sites 0051, 0222, 0292, MIZU and TSKB at the time of the seismic waves passage. The satellite configuration is not optimal because the number of satellites is limited along the axis N0/N180. Although, this lack of satellites is not due to a satellite tracking issue, because the position of satellite vary little from site to site (over 300 km distance). All sites will have poor accuracy on the north component of the vibration as there are very few sites in the quadrants N135-N225 and N315-N45</figcaption></figure>
<table><caption>**Table 1 |** Component of a triplet of vectors in our computation (minimum elevation of 15 degrees above the horizon). The maximum value of *K* for this triplet is ~ 0.93.</caption><thead><tr><td>Unit vector</td><td>Azimuth (°)</td><td>Elevation (°)</td></tr></thead><tbody><tr><td>u&#773;<sub>0</sub></td><td>0</td><td>15</td></tr><tr><td>v&#773;<sub>0</sub></td><td>90</td><td>15</td></tr><tr><td>w&#773;<sub>0</sub></td><td>90</td><td>90</td></tr></tbody></table>

However, as the satellites lower than 15 degrees elevation cut-off are usually not considered in the processing, the value of *K* can only be < 1. *α* is the maximum value for *K* in our computing:

$$ \alpha = (\vec{u}_0 \times \vec{v}_0) \cdot \vec{w}_0 \quad (2) $$

The value of *K* ranges from 0 (all the satellites are aligned along a perfect line in the sky) to a best possible value (*K* ~ 0.93 in the case of

<figure><img src="image_8.png" /><figcaption>**Figure 4 |** Potential GPS constellation for three IGS (AIRA, MIZU and TSKB) sites located in Japan during the DOY 268 (25 September 2003). The sites are respectively indicated using gray, red and black lines. We characterize this potential by plotting the number of satellites visible at each site (top panel) and the parameter *K* (computed every 15 min) at the site 0292 (lower panel) during the day. At the time of the earthquake (vertical red line) the number of satellites available over Japan is close to 10 satellites. Still, we stress that the orientation of these satellites is important too. The parameter *K* (see text for an accurate description) at the time of the seismic wave propagation is not close (*K* ~ 0.2) to the theoretical value of *K* (~ 0.93). We indicate the maximal theoretical value with a horizontal red line on the lower panel).</figcaption></figure>

the distribution described in Table 1). At the time of the seismic wave propagation in Japan, the value of *K* does not exceed 0.25, one of the lowest values for that day (Figure 4).

The analysis showing the frequency content of the LOS residuals showed phase residuals could be used to recover ground motion (Figure 5). In Figure 6, we compare the motion recovered for the site 0292 with the seismic records of the GEOSCOPE¹ very-broadband (STS-1[29,30]) seismometer located in INUyama, Japan (Figure 2). The seismic records have been corrected from instruments response (into displacement) and decimated to 1Hz. The distance between the two sites is ~ 10 km. As we consider periods ranging from 30 to 50 seconds, we do not expect any difference in the motion of the two sites.

The agreement between both sets of horizontal records (Figure 6) is not optimal ($RMS_{east} = 9.7310^{-2}$, $RMS_{north} = 5.8710^{-2}$,

$$ RMS = \sqrt{\frac{1}{N} \left( \frac{\text{synth} - \text{data}}{\text{data}} \right)^2} \quad \text{with } N = 400). $$

The horizontal amplitudes are not fully recovered by the inversion of the GPS LOS. Even so, the phase of the east component is in good agreement with the phase of the east component of the seismic record. The phases and amplitudes recovered for the north component are not matching

<figure><img src="image_9.png" /><figcaption>**Figure 5 |** Periodograms of the residual time-series of the satellite PRN 4 at the site 0292. The seismic signal is only visible in the period range 10–80s with a strongest amplitude in the ~ 30–50s (20 to 35 mHz) period range.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6** | Results of the inversion computed for the 0292 site (Figure 2). We compare the time-series of GPS displacement (in red) with the seismograms recorded at the GEOSCOPE site INU (in black). We note that while the agreement between the two sets of waveforms is reasonable for eastern and vertical components it could be more problematic for the northern one. Instrument responses have been provided by GEOSCOPE. Both datasets have been band-pass filtered between 30 and 50 seconds (20 and 35 *mHz*) according to the periodogram presented in Figure 5.</figcaption></figure>

the north component of the seismic record. We note a better agreement (phase match) for east direction which is better constrained by the satellites positions in the sky (Table 2). The vertical GPS component compares well ($RMS_{vert} = 2.3210^{-2}$, N=400) with the vertical component of the seismic sensor. We attribute the slight phase delay between the two time-series to an effect of the band-pass Butterworth filter used to produce the seismic time-series. We compare the frequency content of both seismic and GPS vertical components (Figure 7). While the agreement is good in the frequency range from 0.0125 to 0.1 Hz (10 to 80 seconds), unexpected high amplitude

<table><caption>**Table 2** | Positions of the satellites visible from the site 0292 at the time of the seismic wave arrival.</caption><thead><tr><td>PRN</td><td>Az (deg. N)</td><td>Elevation (deg.)</td></tr></thead><tbody><tr><td>4</td><td>118.95</td><td>49.52</td></tr><tr><td>8</td><td>127.90</td><td>21.93</td></tr><tr><td>10</td><td>290.72</td><td>48.10</td></tr><tr><td>13</td><td>42.77</td><td>22.80</td></tr><tr><td>17</td><td>309.21</td><td>17.49</td></tr><tr><td>24</td><td>357.60</td><td>79.12</td></tr><tr><td>27</td><td>92.48</td><td>28.88</td></tr><tr><td>29</td><td>213.30</td><td>28.70</td></tr></tbody></table>

noise is visible for longer periods. We have not fully explained such discrepancy yet.

We showed that the limitations of the GPS are due at least partially to the distribution of the satellites in the sky. We need to set up a processing strategy that enables us to estimate the ground motion whatever the number of visible satellites and their distributions in the sky are. This methodology must enable the detection of both static motions close to seismic ruptures but also of ground displacements produced by propagation of surface waves, namely Rayleigh and Love waves.

Our approach is based on the computation of the double difference (DD) from a pair of satellites and stations (Figure 8). GPS processing packages commonly use this method for computing GPS data because it allows working with few satellites. The DD allows an instantaneous quasianalytical reduction of the major errors (ionosphere, troposphere and clocks delays) assumingmost of the errors are seen by the other components of the DD (satellite or site) at the same time.

Here, in order to prove the sensitivity of DD to ground motion, we employed 2 GPS stations (sites 0051 and 0222) located in the far field at ~ 6° from the epicenter (Figure 2). We picked these two sites because they are located on a small circle and their DD residuals are expected to be minimal. Success with these two stations will make our experiment meaningful for all the other pairs of stations.

We detected oscillations for all the stations at the expected time of seismic arrival at the sites. We focused on the two sites (0051 and 0222), 128 kilometers apart (Figure 2), along the propagation front. We show all DD residuals between the 2 sites in Figure 9. We have compared the Fourier spectra of GPS DD residuals with and without seismic transients. During the hour including the earthquake occurrence, we observe a strong peak between 20 and 35 *mHz* (T ~ 30 to 50s) (Figure 10).

Our next aim was to determine whether the GPS DD observations are related to the seismic source only, or to some other process (such as site effect, monument oscillation, etc.) as well. We thus performed a test, comparing the DD observations with synthetic DD. In order to obtain synthetic DD, we first computed the synthetic waveforms expected at the GPS sites by using a combined discrete-wave number and reflectivity method³¹⁻³³. This is completed by combining synthetic seismograms computed using a CMT solution³⁴, and a crustal model (Table 3) from a regional average of the global tomography³⁵.

<table><caption>**Table 3** | Crustal structure used in the numerical simulation</caption><thead><tr><td>ThicknessΔ (km)</td><td>α(km/s)</td><td>β(km/s)</td><td>ρ(g/m³)</td></tr></thead><tbody><tr><td>11</td><td>5.65</td><td>3.43</td><td>2.68</td></tr><tr><td>9</td><td>6.55</td><td>3.70</td><td>2.90</td></tr><tr><td>10</td><td>7.16</td><td>3.84</td><td>3.06</td></tr><tr><td>halfspace</td><td>7.98</td><td>4.47</td><td>3.31</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Figure 7 | Comparison of the time-frequency estimates of a) Seismic data and b) GPS Data. Time-frequency amplitudes are calculated using the s-transform<sup>46</sup> (resolution factor=1) and are normalized using the maximum amplitude for each frequency near the event arrival time. Relative arrival times of the peak amplitudes (c) are picked from both seismic and GPS data. Good correlation is found in the 15–50 second (0.015–0.07Hz) period band.</figcaption></figure>

The source parameters used in the computation are based on the seismic moment of ~3 × 10<sup>21</sup> Nm, the source depth of 13 km, and the focal mechanism: strike 234°, Dip 7°, and Rake 103°. The synthetic ground displacements are then converted into synthetic DD residuals (S<sub>i</sub><sup>k</sup>) by projecting the local motion (E, N and U) onto the LOS of the satellite that is defined by the unit vector (u,v,w):

$$S_i^k = E.u + N.v + U.w \quad (3)$$

The synthetic seismic DD, S<sub>ij</sub><sup>kl</sup>, is using four projections created at two sites k and l:

$$S_{ij}^{kl} = \frac{1}{\lambda_{LC}} (S_i^k - S_i^l) - (S_j^k - S_j^l) \quad (4)$$

where λ<sub>LC</sub> is the wave length of the L<sub>C</sub> signal (~48cm per cycle) and S<sub>i</sub><sup>k</sup> the seismic vibration 154 projected along the LOS between satellite x and the satellite y.

We compare the synthetic residuals S<sub>ij</sub><sup>kl</sup> with the observed phase time-series (Figure 11). All seismic and GPS waveforms have been

<figure><img src="image_12.png" /><figcaption>Figure 8 | Schematic illustration of the ground shaking propagating along the surface and perturbing four ray-paths arriving at the two GPS sites i and j. The four residuals of the four signals received at site i and j are computed using GAMIT<sup>36</sup>. Because Rayleigh waves are not arriving at the same time to both sites, the residuals of the Double Difference are not equal to zero.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 9 | Examples of the residuals observed on the DD between the benchmarks 0222 and 0051 and each couple of satellite (as in the sketch of Figure 8). The wavelength of the L~C~ combination is ~ 48 cm per cycle. Time is indicated in seconds. Pseudo-Range numbers (PRN) are indicated for each DD.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 10** | Frequency spectra of the residuals presented in Figure 9. a) Amplitudes of the Fourier spectra of the GPS residuals computed for data recorded before the seismic wave arrival. Peak frequencies are located below 25mHz. b) Amplitudes of the Fourier spectra of the GPS residuals plotted in Figure 9. Peak frequencies are located around 30mHz.</figcaption></figure>

band-pass filtered between 20 and 35 mHz, coherently with the frequency signature of the transient displacement (Figure 10). The good agreement between the observed and the predicted residuals (Table 4) confirms that the DD are detecting seismic surface waves at around 700 km from the epicentre. We list the values of the RMS in Table 4.

## Discussion

We show that GPS inversion of phases in the local reference (East, North and Vertical directions) is limited by the distribution of the GPS constellation. The start of the GALILEO program, by adding more satellites in the sky, brings great hopes in the matter. Additionally, we show that the inversion of phases is meant to be limited to the local reference frame (East, North, Up) but could be expanded to the reference frame defined by the LOS of the visible satellites, using the DD method.

Whatever its capability to recover local motion histories (East, North, Vertical components), GPS can provide constraints by using LOS residuals within DD geometry. We suggest that GPS observations can provide additional constraints on the seismic moment.

This successful comparison between synthetic ground displacement and the DD GPS observations suggests that the use of GPS networks in regional or even global tomographies is possible. Our methodology also offers an opportunity to improve the GPS sensitivity to ground shaking, and constitutes a powerful tool for a better understanding of long-period geophysical processes through the lithosphere.

<table><caption>**Table 4** | RMS of synthetic DD plotted on Figure 11. The number of observations N=400. RMS formula is given in the text.</caption><thead><tr><td>DD</td><td>RMS</td></tr></thead><tbody><tr><td>04-08</td><td>0.0715</td></tr><tr><td>10-04</td><td>0.0526</td></tr><tr><td>10-08</td><td>0.1339</td></tr><tr><td>13-04</td><td>0.1426</td></tr><tr><td>13-08</td><td>0.1700</td></tr><tr><td>13-10</td><td>0.0734</td></tr><tr><td>24-04</td><td>0.0083</td></tr><tr><td>24-08</td><td>0.0543</td></tr><tr><td>24-10</td><td>0.0557</td></tr><tr><td>24-13</td><td>0.1574</td></tr></tbody></table>

In validating the comparison between a seismic 1D model and the GPS observations, we suggest the GPS might now be considered a supplement to broadband seismic networks. Furthermore, GPS networks could provide complementary constraints on the surface displacement by increasing the spatial resolution of the observed surface wave fronts and/or to the kinematic rupture modeling. The agreement between data and synthetics suggests that DD can be integrated in tomographic problems instead of fixed local three components, taking advantage of the higher sensitivity of DDs.

## Methods

**Daily Data processing.** To test the sensitivity of the network to lateral motion, we select GPS receivers located approximately at 700 km radius from the epicentre (indicated by a circle in Figure 2). Working with such sub-network configuration allows for detecting both the arrival time of the surface waves, and the wave polarity and amplitude variations with the site.

The local permanent GPS stations (KSMV, MIZU, TSKB, USUD, YSSK [not on Figure 2, N47.03, E142.72] in Figure 2 belonging to the International GPS Service (IGS) have been integrated in the analysis to stabilize the processing using accurate orbits.

All calculation have been completed using the GAMIT/GLOBK Tool suite³⁶, ³⁷, following a standard procedure for large static networks (the network is not affected by co-seismic offsets) We have extracted residuals that could not be explained by tides, GPS constellation cycle slip, or tropospheric and ionospheric perturbations.

In order to minimize the noise due to the very low elevation LOS, the minimal elevation angle was fixed at 15 degrees above the horizon and the ambiguity resolution was done assuming a baseline length lower than 500 km. The troposphere delay was estimated using the standard meteorological model³⁸ but as the weather is supposed to be constant during the wave propagation over the network, we do not expect oscillation in the computed GPS waveforms. Additionally, we have used the LC (ionosphere free) frequency, which minimizes the effect of the ionosphere on the GPS observables, and assures the results are free of high frequency post-seismic perturbation taking place in the ionosphere³⁹,⁴⁰.

**The double difference.** A Double Difference (DD) is computed using signals from two given satellites received at two selected stations. The residuals on the DD are the difference between the residual phases of four signals received at two different stations. The DD are not oriented in the local reference frame but in a reference frame defined by the four lines of sight of the two satellites at the two GPS receivers. The DD represents the sum of displacements of two sites seen by two different satellites. For each couple of available sites, the residual on the double difference between two sites has been investigated. The final DD residuals $\phi(L1/L2)_{ij}^{kl}$ can be described as follows:

$$ \phi(L1/L2)_{ij}^{kl} = \frac{f(L1/L2)}{c} (\delta\rho_j^l - \delta\rho_j^k + \delta\rho_i^k - \delta\rho_i^l) + \epsilon_{ij}^{kl} \quad (5) $$

where $\delta\rho_j^l$ is the residual distance between the station j and the satellite l, c the speed of light and f the frequencies of the signal for the two carrier frequencies L1 and L2. $\epsilon_{ij}^{kl}$ is a scalar number comprising the unmodeled processes (orbits errors, troposphere, ionosphere, noise on the measurement, etc).
<figure><img src="image_15.png" /><figcaption>**Figure 11** | GPS Observables plotted in Figure 9 compared with synthetics. The data and the synthetic data are filtered by using a bandpass filter 0.020 – 0.035Hz. This bandpass filter is in agreement with the spectra made using unfiltered GPS residuals data (Figure 10). There is a good matching between the residuals GPS filtered data and the filtered seismic time-series.</figcaption></figure>
<figure><img src="image_16.png" /></figure>

The residuals computed with equation 5 are free of all instrumental time stamp errors. Finally, we computed the ionosphere free residual $\Phi(LC)_{ij}^{kl}$ from L1 and L2 carrier residuals⁴¹.

$$ \Phi(LC)_{ij}^{kl} = \Lambda_{L_1} (\phi_{ij}^{kl})_{L_1} + \Lambda_{L_2} (\phi_{ij}^{kl})_{L_2} $$

$\Phi_{ij}^{kl}$ represents the sum of all processes not modeled such as the measurement noise of the GPS receiver.

1. Romanowicz, B., Cara, M., Fels, J., and Rouland, D. GEOSCOPE: A French Initiative in Long Period, Three Component, Global Seismic Networks. *Eos Trans. Am. Geophys. Union* **65**, 753–754 (1984).

2. Roult, G., Lépine, J., Bonaimé, S., and Riveira, L. The GEOSCOPE program: state of the art in 2005. In *AGU 2005 Fall Meet. Suppl.*, volume 86, (2005).

3. Boore, D. M. Analog-to-Digital Conversion as a Source of Drifts in Displacements Derived from Digital Recordings of Ground Acceleration. *Bull. Seismol. Soc. Am.* **93**(5), 2017–2024 (2003).

4. Simons, M., Fialko, Y., and Rivera, L. Coseismic deformation from the 1999 M w 7.1 Hector Mine, California, earthquake as inferred from InSAR and GPS observations. *Bull. Seismol. Soc. Am.* **92**(4), 1390–1402 (1999).

5. Marone, C. J., Scholz, C., and Bilham, R. On the mechanics of earthquake afterslip. *J. Geophys. Res.* **96**, 8441–8452 (1991).

6. Kouba, J. Measuring seismic waves induced by large earthquakes with GPS. *Stud. Geophys. Geod.* **47**, 741–755 (2003).

7. Larson, K., Bodin, P., and Gomsberg, J. Using 1-Hz GPS Data to Measure Deformations Caused by the Denali Fault Earthquake. *Science* **300**, 1421–1424 (2003).

8. Bock, Y., Wrawirodirdjo, L., and Melbourne, T. Detection of arbitrary large dynamic ground motions with a dense high-rate GPS network. *Geophys. Res. Lett.* **31**(L06604) (2004).

9. Bilich, A., Cassidy, J., and Larson, K. GPS Seismology: Application to the 2002 Mw 7.9 Denali Fault Earthquake. *Bull. Seismol. Soc. Am.* **98**, 593–606 (2008).

10. Miyazaki, S. et al Modeling the rupture process of the 2003 September 25 Tokachi-Oki (Hokkaido) earthquake using 1-Hz GPS data. *Geophys. Res. Lett.* **31**, L21603 (2004).

11. Emore, G., Haase, J., Choi, K., Larson, K., and Yamagiwa, A. Recovering absolute seismic displacements through combined use of 1-Hz GPS and strong motion accelerometers. *Bull. Seismol. Soc. Am.* **97**, 357–378 (2007).

12. Ji, C., Larson, K., Tan, Y., Hudnut, K., and Choi, K. Slip history of the 2003 San Simeon earthquake constrained by combining 1-Hz GPS, strong motion, and teleseismic data. *Geophys. Res. Lett.* **31**(17) (2004).

13. Wang, G., Boore, D. M., Tang, G., and Zhou, X. Comparisons of ground motions from collocated and closely spaced one-sample- per-second Global Positioning System and accelerometer recordings of the 2003 M 6.5 San Simeon, California, earthquake in the Parkfield Region. *Bull. Seismol. Soc. Am.* **97**(1B), 76–90 (2007).

14. Ohta, Y., Meiano, I., Sagiya, T., Kimata, F., and Hirahara, K. Large surface wave of the 2004 Sumatra-Andaman earthquake captured by the very long baseline kinematic analysis of 1-Hz GPS data. *Earth Planets Space* **58**, 153–157 (2006).

15. Yokota, Y., Koketsu, K., Hikima, K., and Miyazaki, S. Ability of 1-Hz GPS data to infer the source process of a medium-sized earthquake: The case of the 2008 Iwate-Miyagi Nairiku, Japan, earthquake. *Geophys. Res. Lett.* **36**(L12301) (2009).

16. Delouis, B., Nocquet, J.-M., and Vallée, M. Slip distribution of the February 27, 2010 Mw = 8.8 Maule Earthquake, central Chile, from static and high-rate GPS, InSAR, and broadband teleseismic data. *Geophys. Res. Lett.* **37**(L17305) (2010).

17. Houlié, N. and Allen, R. The Instantaneous Displacement (ID) Method: Application of rapid displacement estimates to earthquake early warning alerts. *Eos Trans. Am. Geophys. Union* **89**(23) (2007).

18. Romanowicz, B. Global mantle tomography: Progress status in the past 10 years. *Annual Review of Earth and Planetary Sciences* **31**(1), 303–328 (2003).

19. Altamimi, Z., Sillard, P., and Boucher, C. ITRF2000: A New Release of the International Terrestrial Reference Frame for Earth Science Applications. *J. Geophys. Res.* **107**(B10), 2114 (2002).

20. Trota, A., Houlié, N., Sigmundsson, F., and Jónsson, S. Deformations of the Furnas and Sete Cidade Volcanoes. Velocities and furt her investigations. *Geophys. J. Int.* **166**(2), 952 (2006).

21. Wdowinski, S., Bock, Y., Zhang, J., Fang, P., and Genrich, J. Southern California Permanent GPS Geodetic Array: Spatial filtering of daily solutions for estimating coseismic and postseismic displacements induced by the 1992 Landers earthquake. *J. Geophys. Res.* **102**(B8), 18057–18070 (1997).

22. Choi, K., Bilich, A., Larson, K., and Axelrad, P. Modified sidereal filtering: implications for high-rate GPS positioning. *Geophys. Res. Lett.* **31**(22) (2004).

23. Miyazaki, S., Hatanaka, Y., Sagiya, T., and Tada, T. The nationwide GPS array as an Earth Observation System. *Bull. Geogr. Surv. Inst.* **44**, 11–22 (1998).

24. Hatanaka, Y. et al Improvement of the Analysis Strategy of GEONET. *Bull. Geol. Surv. Inst.* **49** (2003).

25. Coffin, M. and Hirata, N. Large Earthquake Strikes Hokkaido, Japan. *Eos Trans. Am. Geophys. Union* **84**(42) (2003).

26. Yagi, Y. Source rupture process of the 2003 Tokachi-oki earthquake determined by joint inversion of teleseismic body wave and strong ground motion data. *Earth & Planet. Sc. Lett.* **56**, 311–316 (2004).

27. Houlié, N., Briole, P., Nercessian, A., and Murakami, M. Sounding the plume of the 18 August 2000 eruption of Miyakejima volcano (Japan) using GPS. *Geophys. Res. Lett.* **32**(L05302) (2005).

28. Houlié, N., Briole, P., and Nercessian, A. 2005, March 9th Mount St Helens explosion. *Eos Trans. Am. Geophys. Union* **86**(30), 277–284 (2005).

29. Wielandt, E. and Streckeisen, G. The leaf-spring seismometer; design and performance. *Bull. Seismol. Soc. Am.* **72**(A), 2349–2367 (1982).

30. Wielandt, E. and Steim, J. A digital very broadband seismograph. *Ann. Geophys.* **4b**(227) (1986).

31. Bouchon, M. A simple method to calculate Green's functions for elastic layered media. *Bull. Seismol. Soc. Am.* **71**, 959–971 (1981).

32. Kennett, B. L. N. Seismic Wave Propagation in Stratified Media. Cambridge University Press, Cambridge, England, (1985).

33. Coutant, O. *Program of Numerical Simulation AXITRA. Res. Report LGIT-Grenoble, in French*, (1989).

34. Dziewonski, A., Chou, T.-A., and Woodhouse, J. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. *J. Geophys. Res.* **86**, 2825–2852 (1981).

35. Shapiro, N. and Ritzwoller, M. Monte-carlo inversion for a global shear velocity model of the crust and upper mantle. *Geophys. J. Int.* **151**, 88–105 (2002).

36. King, R. and Bock, Y. *Documentation of the GAMIT software*. MIT/SIO, (2006).

37. Herring, T. GLOBK: Global Kalman Filter VLBI and GPS Analysis Program, version 10.2 (2005).

38. Saastamoinen, J. Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites. *AGU Geophysics Monograph Series* **15**, 247–251 (1972).

39. Ducic, V., Artru, J., and Lognonné, P. Ionospheric remote sensing of the Denali Earthquake Rayleigh surface waves. *Geophys. Res. Lett.* **30**(18), 1951, doi:10.1029/2003GL017812 (2003).

40. Lognonné, P. et al Ground based GPS imaging of ionospheric post-seismic signal. *Planet. Space. Science* **54**, 528–540 (2006).

41. Bender, P. L. and Larden, D. GPS carrier phase ambiguity resolution over long baselines. In *Proceedings of the First International Symposium on Precise Positioning with the Global Positioning System*, Goad, C., editor, volume 1, 357–361, (1985).

42. Jamason, P. et al SOPAC web site (http://sopac.ucsd.edu). *GPS Solutions* **8**(4), 272–277 (2004).

43. King, G. C. P., Stein, R. S., and Lin, J. Static stress changes and the triggering of earthquakes. *Bull. Seismol. Soc. Am.* **84**, 835–953 (1994).

44. Murray, M. H., Dreger, D., Neuhauser, D., Baxter, D., Gee, L., and Romanowicz, B. Real-time geodesy. *Seis. Res. Lett.* **69**(145) (1998).

45. Houlié, N. and Romanowicz, B. Asymmetric deformation across the San Francisco Bay Area faults from GPS observations in Northern California. *Phys. Earth Planet. Inter.* **184**(3–4), 143–153 (2011).

46. Stockwell, R., Mansinha, L., and Lowe, R. Localization of the complex spectrum: the S transform. *IEEE transactions on signal processing* **44**(4), 998–1001 (1996).

## Acknowledgments

We thank Geographical Survey Institute for providing the 1Hz GPS data. This work was partially supported by Berkeley Seismological Laboratory (BSL) and by US Geological Survey Cooperative Agreement to BSL #G10AC00141. It is BSL Contribution #11-05

## Author contributions

N.H. and G.O. initial idea, conception and writing; T.B. calibration with seismic sensor; N.S. computed seismic synthetics presents and helped with the writing; P.L. provided the dataset and helped with the writing; and M.M. provided the dataset.

## Additional information

Competing financial interests: The authors declare no competing financial interests

License: This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivative Works 3.0 Unported License. To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-nd/3.0/

How to cite this article: Houlié, N. et al. New approach to detect seismic surface waves in 1Hz-sampled GPS time series. *Sci. Rep.* 1, 44; DOI:10.1038/srep00044 (2011).