Introduction

After a large earthquake the whole Earth would “ring” for up to days, when it undergoes free oscillation with discrete normal modes much like a musical instrument. The elastic normal modes of a spherically symmetric Earth are of two types: spheroidal modes have both radial (or vertical) and transverse (or horizontal) components, and toroidal modes have only horizontal components. These modes exist in (2l + 1)-fold degenerate multiplets, where l is the spherical-harmonic degree. Symmetry-breaking perturbations of rotation and ellipticity as well as general heterogeneity cause frequency splitting to the degenerate multiplets into split singlets with close-by frequencies dependent on the spherical-harmonic order m.

The Earth’s elastic normal modes have periods just short of an hour; ones with periods on the order of minutes and longer have been clearly observed since the 1960s in the spectra of seismic records worldwide (e.g., Gilbert and Dziewonski 1975; Bullen and Bolt 1985). The observations are now made routinely by broadband seismometers and gravimeters (Agnew et al. 1986), and to a lesser extent other geophysical instruments such as strainmeters, tiltmeters (e.g., Park et al. 2008), and even groundwater level (Yan et al. 2016).

In seismological applications, GPS (Global Positioning System) measurements of co-seismic static surface displacement have been commonly used to help determine the slip distributions and faulting mechanisms of earthquakes. On the other hand, when sampled rapid enough down to the range of seconds, GPS would record the actual 3-D ground displacement during an earthquake (e.g., Larson et al. 2003), in parallel to seismometers for the ground velocity or acceleration. A GPS receiver is free from any instrument response and not band-limited in frequency, so in principle it can record normal modes of all periods to the extent supported by its sensitivity and sampling rate.

Large earthquakes can excite free oscillations with displacement amplitudes up to millimeter (mm) scale. For example, the spheroidal modes excited by the 1960 Chilean Earthquake were on the order of 0.1–1 mm (Kanamori and Anderson 1975); the “breathing” mode of 0S0 excited by the 2004 Mw 9.3 Sumatra earthquake was about 0.05 mm (Park et al. 2005). That is at least an order of magnitude smaller than typical precision level of modern GPS observations at their best, which is about 3 mm for horizontal and 6 mm for vertical positions. However, that does not mean that the seismic normal-mode displacements are undetectable in GPS time records. For one thing, the normal modes are harmonic in time and hence enjoy the processing gain in Fourier analysis—the amplitude signal-to-noise ratio (SNR) of a harmonic signal in the frequency domain is proportional to the square root of the number of data points. For example, using 17 years of GPS data, Yuan and Chao (2012) demonstrated the achievement of 0.1 mm horizontal and 0.3 mm vertical precisions for tidal displacements, even for single station record.

Furthermore, the fact that all records contain signals of the same normal-mode frequencies allows further boost of signal detection by proper multiple-station data-stacking procedures. Mitsui and Heki (2012) reported detection of the normal modes excited by the 2011 Tohoku, Japan, earthquake using simple frequency-domain power spectral stacking of GPS records from the Japanese GEONET network. One can take further advantage of the spatial coherence in the displacement field from (hopefully worldwide) network stations, because a normal mode is a (generalized) standing wave having a specific surface spherical-harmonic pattern on Earth surface. A number of effective time-domain stacking schemes have been developed for that purpose, most notably the spherical-harmonic stacking (SHS) (Buland et al. 1979; Chao and Ding 2014) and the optimum sequence estimation (Ding and Shen 2013); see Ding and Chao (2015a) for a review. These stacking techniques not only boost the SNR by ~\(\sqrt{N}\) (N is the number of stations or records used) of target modes, but also suppress the non-target signals. Further gain in SNR can be achieved by simultaneous processing of multiple earthquakes, as they excite the same set of normal modes, differing only in excitation amplitudes.

With the paradigm above, we shall conduct a comprehensive study and present the results of the data-stacking techniques on GPS detection of Earth’s normal modes excited by the 2011 Tohoku earthquake.

Data

We use the 3-D displacement data from three continuous GPS networks—two dense regional and one global networks, as shown in Fig. 1: (i) the Japan Geospatial Information Authority GEONET (GNSS Earth Observation Network) network; we use around 1000 (out of 1230) stations in different cases. (ii) The PBO (Plate Boundary Observatory) Array in the Western USA; we use around 600 (out of 1548) stations. (iii) The IGS (International GNSS Service) global network with around 140 (out of 337) stations. Unused stations are those deemed relatively noisy in their records.

Fig. 1
figure 1

Stations of the three continuous GPS networks used in this paper: a the Japan Geospatial Information Authority GEONET (GNSS Earth Observation Network) network; b the PBO (Plate Boundary Observatory) Array in the Western USA; c the IGS (International GNSS Service) network (with a and b also shown in inset boxes)

Full size image

All these GPS data are solved from GPS 30-s interval raw data by us using the kinematic precise point positioning algorithm (Zumberge et al. 1997) of the GIPSY/OASIS-II (Version 6.1) software. JPL’s reprocessed fiducial satellite orbits, 30-s high-rate clocks, and Earth orientation products are held fixed, including single station ambiguity fixing method (Bertiger et al. 2010). To minimize daybreak effects, the data are processed in 30-h sessions centered daily on 12:00 (GPS time), with the 30-s station displacement estimates corresponding to each unique GPS day (24 h) extracted and concatenated to form displacement time series in three local components in the geographical coordinates (East, North and Up), referenced to the IGS08 reference frame. The other parameter settings follow that of Yuan and Chao (2012). The sampling interval is set at 30 s, hence the Nyquist period is 1 min, corresponding to attainable signal frequency up to ~ 17 millihertz (mHz).

We let the time span (T) of our analysis start at 15 min after the nominal origin time of the Tohoku earthquake—05:46:18 UTC March 11, 2011, by which time all the fault rupture actions had ceased. For study of seismically excited normal modes, typically one lets T to be several days and sufficiently long before the free oscillations decay away to allow desired spectral resolutions. Here, however, the GPS solutions are based on IGS daily updated orbits, effectively introducing an artificial discontinuity into the time series. These daily discontinuities, albeit fine, are found to give rise to artificial repeated spectral peaks in the spectra, seriously corrupting the spectra for practical detection of small signals.

So T needs to be terminated at the end of the first UT calendar day, giving a finite 18 h of timespan. In practice, 21 h are actually available because for ambiguity resolution IGS typically issues 30 h of span for each “daily” solution with two overlapping 3-h segments appended on either end of the 24-h segment. The signal levels at the end of this “day” have dropped significantly according to the normal mode’s Q-values typically of a few hundred. Nevertheless, this T limits the spectral resolution at 1/(21 h) or 0.013 MHz, insufficient to resolve the mode splittings from multiplet to singlets—for instance the frequency splitting of the longest-period mode 0S2 is no more than ~ 0.005 MHz. Therefore, the eigenfrequencies we shall refer to below are the (unresolved) degenerate eigenfrequencies of the multiplets.

Data-stacking techniques

Data stacking with respect to multiple records aims to accentuate, for the benefit of detection, a common signal that is present (but “buried” in random noise) in a temporally and/or spatially coherent manner. A seismic normal mode does constitute such a signal—a harmonic standing wave phase-synchronous in time, having a spatial pattern of a surface spherical harmonic function on Earth surface. Taking advantage of these favorable features of the observed data, we shall now try various data-stacking techniques.

Frequency-domain stacking

The excitation amplitude of a normal mode by an earthquake depends on the magnitude and source mechanism of the earthquake. Proportional to that, each actual amplitude that gets recorded by a given seismometer differs from station to station depending on the location of the station relative to the source, but all seismograms contain the same set of normal modes each manifesting as a harmonic signals of its characteristic eigenfrequency.

On a Fourier power spectrum |F(ω)|2, or periodogram, where F(ω) is the complex Fourier transform of the time record obtained, for instance, by the fast Fourier transform (FFT), a harmonic signal manifests as a spectral peak at the specific frequency. Suppose we add up all the N periodograms from a station network: Σj = 1,N |Fj(ω)|2/N, which is equivalent to forming the arithmetic mean of the N periodogram, we obtain more “stabilized” spectral peaks. While it does not boost the overall SNR, this sum-periodogram procedure “flattens” the noise spectrum by reducing it by the factor of N (see Ding and Chao 2015b). This facilitates the detection of the harmonic signals by enabling them to stand out more prominently against noise level. That is the scheme followed by Mitsui and Heki (2012) in their GPS detections of normal modes, except that they took the amplitude, rather than power, spectra |F(ω)|.

In this paper, following Smylie et al. (1993), we shall adopt the product-periodogram, which is equivalent to forming the geometric mean of multiple records. We form the product of N periodograms: [πj = 1,N |Fi(ω)|2]1/N, which upon taking the logarithm forms the conventional power dB, 20/N × Σj = 1,N log|Fj(ω)|. The product-periodogram has an advantage over the sum-periodogram in that any scaling difference among individual spectra (for example, those owing to different instrument types or calibration) only introduces static offsets (which do not matter) and requires no data normalization.

Computing F(ω) via FFT, we zero-pad our time series beforehand from 21 h (see above) to 11 days, thus allowing interpolation of spectrum at much finer frequency spacing than the elementary Fourier frequency bin. For harmonic signals this practice under certain conditions can be beneficial in spectral resolution, as demonstrated by Buland and Gilbert (1978), in the spirit of the super-resolution raised by Munk and Hasselmann (1964).

Time-domain stacking

Time-domain stacking is to enhance signal detectability using multiple time records by taking advantage of prior knowledge about the spatial pattern of the target signal. An earlier success was reported by Gilbert and Dziewonski (1975) for normal mode detection, where the basic idea is to add up (or stack) multiple time series, but not before they are properly normalized. Thus, while the frequency-domain stacking gives equal weight to each time series, the time-domain stacking imposes different weight (with polarity) on each time series in certain designated way. There are a variety of designated normalizations in the literature that have been developed over the decades (see Ding and Chao 2015a for a review).

Suppose the target signal is identical in all records, such as the “breathing” singlet mode 0S0 (or for that matter all nS0 radial overtones). Then a simple summation of all records suffices to boost the amplitude SNR by \(\sqrt{N}\). That is actually a special trivial case of the SHS stacking below.

Now suppose we know that the signal has a spatial pattern of P(x,y) on Earth surface whose value is a function of x- and y-coordinates (e.g., latitude and longitude) with positive or negative polarities. A simple sum would no longer work as it tends to have the signals canceled in random-walk fashion rather than enhanced straightly. So it makes sense to divide out the P(x,y) value from, i.e., to “normalize”, each time record prior to summation. That constitutes the so-called multi-station experiment scheme proposed by Courtier et al. (2000). However, the division is often a poor numerical operation that ought to be avoided, as it can disproportionately magnify noises especially when near P(x,y)’s nodes where the value approaches zero.

Suppose the pattern P(x,y) is a surface spherical harmonic function Ylm(x,y) of degree l and order m, as is the case for many geophysical phenomena including the seismic normal modes under present study, as well as tides, centrifugal potential, mass loading, and multipoles of the gravitational and geomagnetic fields. Note that P(x,y) is not necessarily real-valued; a complex-valued P(x,y) can signify propagating pattern across the Earth surface. Then, instead of division, one can multiply each time record by Ylm*(x,y), the complex conjugate of Ylm(x,y), before summation or stacking. If the stations being stacked are reasonably globally distributed, then the sum would tend to cancel out (or at least suppress) all other Ylm’s of degree and order different than the target one because of the mathematical orthogonality among the spherical harmonics. The end result is a composite time series, the stack, wherein the target signal is well accentuated. That is the principle of the spherical-harmonic stacking (SHS), first devised by Buland et al. (1979) for enhancing the vertical component of spheroidal singlets, and later extended to horizontal components of both spheroidal and toroidal singlets (making use of the vector spherical harmonics) by Chao and Ding (2014).

In this paper, we mainly report the results of the optimal sequence estimation (OSE), another multiple-station stacking method in parallel. OSE was devised by Ding and Shen (2013) and extended by Ding and Chao (2015a) for, although in principle not limited to, spatially spherical-harmonic signal extraction. Suppose the time record from each station is a superposition of many modes of spherical-harmonics Ylms. In OSE processing, one builds the multiple-station data at any given moment in time in terms of superposition of a selected set of Ylm modes, and then solve for the coefficients of these modes simultaneously by least-squares regression. The number of selected modes need not be large—it typically suffices to include the 2L + 1 Ylms of different orders m belonging to the target degree L as “companions” to absorb “stray” signals. As all harmonics are orthogonal, synthetic experiments have demonstrated that this procedure is quite robust, even when the spatial coverage of the network stations is not quite globally distributed. Note that OSE when formulated in complex form applies not only to standing waves, but also propagating waves in general (e.g., Ding and Chao 2017).

Specifically, OSE consists of the following: for the vertical Up component and for any given point in time t, one construct the matrix equation U(t) = YA(t), where \(U\left( t \right) = \left[ {\begin{array}{*{20}c} {u_{1} \left( t \right)} & {u_{2} \left( t \right)} & \ldots & {u_{N} \left( t \right)} \\ \end{array} } \right]^{{\text{T}}}\) is the N-vector of network data u, \(A(t) = \left[ {\begin{array}{*{20}c} {A_{L( - L)} (t)} & {A_{L( - L + 1)} (t)} & \ldots & {A_{LL} (t)} \\ \end{array} } \right]^{{\text{T}}}\) is the (2L + 1)-vector of the coefficients to be inverted for the target degree L, and

$$Y = \left[ {\begin{array}{*{20}c} {Y_{L( - L)} (\Omega_{1} )} & {Y_{L( - L + 1)} (\Omega_{1} )} & \cdots & {Y_{LL} (\Omega_{1} )} \\ {Y_{L( - L)} (\Omega_{2} )} & {Y_{L( - L + 1)} (\Omega_{2} )} & \cdots & {Y_{LL} (\Omega_{2} )} \\ \vdots & \vdots & \ddots & \vdots \\ {Y_{L( - L)} (\Omega_{N} )} & {Y_{L( - L + 1)} (\Omega_{N} )} & \cdots & {Y_{LL} (\Omega_{N} )} \\ \end{array} } \right]$$
(1)

is the N × (2L + 1) kernel matrix. For number of station N ≥ 2L + 1, one can thus obtain the least-squares solution for A. A composite time series is subsequently constructed along time by repeating the procedure in this optimal manner, simultaneously for the included modes. Similar formulation is developed for the two horizontal East/North components that are combined in terms of complex vector spherical harmonics (i.e., vector harmonics replacing scalar harmonics in matrix Y), and corresponding data stack can be constructed for least-squares regression; for details see Ding and Chao (2015a).

Results

Our results presented below consist of Fourier power spectra (i.e., periodograms) of data stacks, stacked either in the frequency domain or in the time domain as described above. The time-domain stacks are the major results; the frequency-domain ones are mainly for reference and comparison. The spectra consist of discrete peaks representing the normal modes excited by the Tohoku earthquake. Overlain on each spectrum, we plot vertical solid/dash lines to indicate the degenerate eigenfrequencies of the fundamental spheroidal/toroidal multiplet modes 0Sl/0Tl, determined according to the Preliminary Reference Earth Model (PREM; Dziewonski and Anderson 1981; Masters and Widmer 1995), against which the observed spectral peaks can be identified. We shall present the frequency range up to 5 MHz (period ~ 3.3 min), corresponding to l ~ 43 fundamental modes; higher frequencies see either little or non-resolvable normal mode energy.

Frequency-domain stacking

First we duplicate the sum-periodograms of Mitsui and Heki (2012; their Fig. 2), but for power spectra (instead of amplitude). Considering the fact that the selection of stations can make significant differences in the stacked spectrum (Gilbert and Dziewonski 1975), we took from the Japan GEONET network as many stations as feasible—1019 stations (see Fig. 1), compared to just over 300 stations by Mitsui and Heki (2012)). The result (in logarithmic dB) for the three components East/North/Up is given in Fig. 2a. Figure 2b is the same as Fig. 2a but for our corresponding product-periodograms. The E-component has the higher SNR, as expected since the observed co-seismic displacements by the Tohoku earthquake are predominantly eastward of up to 5.3 m (Ozawa et al. 2011).

Fig. 2
figure 2

Frequency-domain power spectral stacks (in dB) to show the normal modes, from 1019 station 3-D (E/N/U-components) records of the Japan GEONET network for the 2011 Tohoku earthquake. a Sum-periodogram (corresponding to Fig. 2 of Mitsui and Heki 2012; no relative vertical shift); b product-periodogram (the low-frequency portion of the middle panel shifted downward to show the upper details). The vertical solid/dash lines indicate the degenerate eigenfrequencies of the fundamental spheroidal/toroidal multiplets according to PREM model

Full size image

In Fig. 2, we see similar sets of spectral peaks in the product-periodogram as the sum-periodogram, but somewhat more clearly given its more effective performance as stated above. Yan et al. (2016) have obtained a similar product-periodogram from 43 groundwater array station records for the same earthquake. Against the PREM-modeled eigenfrequencies (vertical lines in the figures), the nearly uniform-spaced large peaks are identified, not surprisingly, to be the well-excited fundamental modes that correspond to surface waves. The spheroidal modes show up in both horizontal and vertical components, whereas the toroidal modes are absent in the vertical, as noted above. One confirms, as did Mitsui and Heki (2012) that the GEONET GPS displacement has recorded unambiguously the Earth’s normal modes excited by the Tohoku earthquake, most clearly for the fundamental modes in the frequency band of ~ 1.5–3.5 MHz. It should be noted that the excited modes in general have anti-nodes with peak amplitude at the earthquake source, thereby favoring their detection by the near-field GEONET stations not far from the source.

Figure 3 presents the product-periodograms just like Fig. 2b but for the western USA PBO stations—we have selected 685/685/489 stations used for East/North/Up components out of the entire 1548 stations (see Fig. 1). The generally weaker SNR of the peaks relative to the GEONET records is attributed to the fact that the PBO network, being far-field, is generally off from the antinode locations (such as in the near-field case of GEONET) for the normal modes and hence of weaker amplitude.

Fig. 3
figure 3

GPS product-periodogram, same as Fig. 2b, but for the Western US PBO network with 685/685/489 stations used for E/N/U-components

Full size image

We have also obtained the product-periodograms for the global IGS network with 260 selected stations. The spectra (not shown) have similar characteristics as Fig. 3, but further weaker.

Time-domain stacking

We have first tried SHS processing of the three sets of GPS network data. The results (not shown) are somewhat more noisy and less clean than those obtained via the OSE-processing scheme (which we report below). That is not surprising because the performance of SHS relies more on the global distribution of the stations than does OSE (Ding and Chao 2015a), while GEONET and PBO are non-global regional networks.

The OSE scheme yields one single composite time series from multiple-station time series upon optimal weighting with respect to designated spatial function, in our case the scalar or vector surface spherical harmonics as the case may be, as described above. We apply OSE to the near-field network GEONET, and obtain the power spectra of the composite time series targeted for designated Ylm, for the horizontal and vertical components. Figure 4a gives the power spectra for l = 1 and m =  − 1, 0, 1, for spheroidal modes in the U-component. Ideally, if the OSE operation is perfectly effective, we expect to see the (m =  − 1, 0, 1) singlets of the nth radial overtones nS1 that are excited, with all others suppressed to the extent dependence upon the global-ness of the station distribution.

Fig. 4
figure 4

The Fourier power spectra of OSE stack for the spheroidal modes from the GEONET network with 968/910 stations, for: a U-component for l = 1 (m = − 1, 0, 1); b U-component for l = 2 (m =  − 2,  − 1, 0, 1, 2); c N/E-component for l = 2 (m =  −  2,  − 1, 0, 1, 2); d N/E-component for l = 3 (m =  − 3, − 2, −  1, 0, 1, 2, 3). The vertical lines indicate the PREM eigenfrequencies of the spheroidal fundamental modes

Full size image

Figures 4, 5, 6 and 7 present selected examples of Fourier power spectra of the OSE stacks, for spheroidal vs. toroidal modes, for the low-degree l = 1, 2, 3 with respective order m’s, and for the vertical (U) vs. horizontal (N/E) components, from different GPS networks of GEONET, PBO or IGS. The most striking difference of these time-domain stacking results from the power spectral stacking (Figs. 2, 3) is the enhanced SNR for the normal modes and hence the heightened spectral peaks, which match remarkably well in general against the vertical solid/dashed lines that indicate the PREM eigenfrequencies of the spheroidal/toroidal fundamental modes.

Fig. 5
figure 5

The same as Fig. 4, but for toroidal modes: a N/E-component for l = 2 (m =  −  2,  − 1,0,1,2); b N/E-component for l = 3 (m = − 3, − 2, − 1, 0, 1, 2, 3). The vertical solid/dashed lines indicate the PREM eigenfrequencies of the spheroidal/toroidal fundamental modes

Full size image
Fig. 6
figure 6

The same as Fig. 4, but for the western USA PBO network of 685/585 stations for l = 3 (m = − 3, − 2, − 1, 0, 1, 2, 3): a U-component of spheroidal modes; b N/E-component of spheroidal modes; c N/E-component of toroidal modes. The vertical solid/dashed lines indicate the PREM eigenfrequencies of the spheroidal/toroidal fundamental modes

Full size image
Fig. 7
figure 7

The same as Fig. 4, but for the global IGS network of 146/143 stations for l = 1 (m = − 1, 0, 1): a U-component of spheroidal modes; b N/E-component of spheroidal modes; c N/E-component of toroidal modes. The vertical solid/dashed lines indicate the PREM eigenfrequencies of the spheroidal/toroidal fundamental modes

Full size image

Figure 4 targets the spheroidal modes from the GEONET network with 968/910 stations—3(a) for U-component for l = 1 (m =  −  1, 0, 1); 3(b) U-component for l = 2 (m =  − 2,  − 1, 0, 1, 2); 3(c) N/E-component for l = 2 (m =  −  2,  − 1, 0, 1, 2); 3(d) N/E-component for l = 3 (m =  −  3,  −  2,  −  1, 0, 1, 2, 3). One sees that essentially all the spheroidal fundamental modes between 1.5 and 5 MHz, i.e., 0S9–0S43, stand out prominently. Lower-frequency modes are less excited by this earthquake (unlike the 2004 Sumatra earthquakes, see Park et al. 2005), where only 0S4 and 0S5 can be seen relatively well. The overtone modes are barely excited to a detectable level at the GPS sensitivity. The only overtone that is identifiable consistently is the cluster 1S3/3S1 at ~ 0.94 MHz.

The results in Fig. 4 further demonstrate an inadequacy in the sensitivity for discriminating the modes by OSE (or any spatial stacking pending global distribution)—that is, the GEONET network being only regional and far from global in either latitude (pertaining to degree l) or longitude (pertaining to order m), all low-frequency modes stay near constant in phase in the limited region because of their long wavelengths. In such a case the OSE procedure is effective in enhancing the SNR of the spatially coherent signals but cannot effectively discriminate and suppress the non-target modes. That is evident in Fig. 4, where all the well-excited modes, namely the fundamental modes (especially spheroidal) whether or not targeted, would show up in the OSE stacks, although weighted differently dependent on the l and m.

The same things happen with respect to the OSE stacks for the toroidal modes. Figure 5 shows the results in the horizontal N/E-components. Here one sees clear toroidal fundamental modes 0T2–0T12 on the left portion of the spectra, as well as the series of spheroidal fundamental modes from 0S10 up to 0S42 for l = 2 in 4(a), and up to 0S31 for l = 3 in 4(b).

Figure 6 gives the OSE stacks for the western USA PBO network using 685/585 GPS stations, targeted for l = 3 (m =  −  3,  − 2,  − 1, 0, 1, 2, 3) and (a) U-component of spheroidal modes; (b) N/E-component of spheroidal modes; (c) N/E-component of toroidal modes. Again, the results show clear detection of the fundamental modes, albeit not as strong as in the GEONET case. The latter is not surprising as stated above, because the PBO network is no longer near the anti-nodes of each of the excited normal modes. In spectrum 5(a) for the U-component, one sees spheroidal 0S8–0S15 (in frequency band 1.4–2.4 MHz), plus a few toroidal modes of 0T3, 0T5, 0T8, and 0T9. The appearance of toroidal modes in the U-component results from their coupling with nearby spheroidal modes (e.g., Masters et al. 1983). The large peak at ~ 0.4 MHz coincides with that of the overtone 2S1, but cannot be further identified in our analysis. In Fig. 6b, c for the N/E-components, most of the toroidal 0T4–0T22 and spheroidal 0S6–0S22 (0.7–3.3 MHz) show up in the spectrum.

Figure 7 gives some examples of the corresponding results for the global IGS network (see Fig. 1) using 146/143 stations for l = 1 (m =  −  1, 0, 1): (a) U-component of spheroidal modes; (b) N/E-component of spheroidal modes; (c) N/E-component of toroidal modes. The results fall short of what are expected despite the relatively worldwide distribution of the IGS network stations. Significant spectral peaks are much fewer and weaker compared to those from the GEONET and PBO networks. Whereas a few peaks in the lowest-frequency end do not correspond to any PREM-modeled eigenfrequencies, only a handful peaks are identifiable against the PREM eigenfrequencies of the fundamental modes (e.g., 0T2–0T4 and 0S3). This presumably can be attributed to the sparsity and small number of available IGS stations and the fact that the suppression of non-target modes is relatively more effective given the global-ness of the network.

Discussion and conclusions

We conduct a host of numerical experiments in a comprehensive effort of detecting and analyzing the Earth’s normal modes of free oscillation in GPS records collected during the 21 h following the 2011 Tohoku, Japan, earthquake. Three continuous GPS networks are taken: two regional network of the Japan GEONET (~ 1000 stations) and the western USA PBO (~ 600 stations), and the global IGS (~ 140 stations). Overall, we find clear signals for fundamental modes both in GEONET and PBO.

The 3-D GPS time records at 30-s sampling rate are solved by us from raw observational data. Using them, our experiments consist of various multiple-record stacking techniques. The simple frequency-domain power spectrum stacking has the benefit of reducing the variance of the noise spectrum to facilitate the identification of the true signals belonging to the normal modes. We use the sum-periodogram (arithmetic mean) scheme for the GEONET records (to duplicate the previous study by Mitsui and Heki (2012), but using three times as many stations), as well as the product-periodogram (geometric mean) scheme which proves to be more effective. The results are presented in Fig. 2 (for GEONET) and 3 (PBO).

Our main results come from the time-domain stackings that convey two benefits for normal mode detection—boosting the amplitude SNR of the target modes by \(\sqrt{N}\), and suppressing non-target modes via the orthogonality of the modes. We have tried two separate schemes: the SHS and OSE that yield comparable results; here we report those of OSE which is found to be more tolerant to non-globalness of the network and yields somewhat more prominent spectral peaks in the present applications.

Representative examples of the OSE results are presented in Figs. 4–7. They as expected show clearly much higher sensitivity and detectability of the normal modes than the frequency-domain stacking. In the case of GEONET in the region where all excited modes have anti-nodes, all the spheroidal fundamental modes 0S9–0S43 below 5 MHz and some of the lower-degree overtones show up as clean spectral peaks against the PREM model eigenfrequencies in various OSE stack spectra with respect to the target degree and order. So do most of the low-degree toroidal fundamental modes. The regional PBO network sees less strong, but still clearly identifiable spectral peaks of the above-mentioned fundamental modes. The global IGS network data on the other hand detect barely a handful of these modes, presumably because of its sparsity and smaller number of stations.

On individual basis under its given measurement sensitivity, a single GPS record is hardly about to reveal the tiny normal mode signals generated even by a greatest earthquake. It is thus gratifying to be able to detect such signals clearly by means of numerical stacking techniques that take advantage of the multiple records in a GPS network. Conversely, this ability demonstrates that the GPS receivers do actually record the seismic signals loyally, notwithstanding their small amplitudes. Currently as things stand in terms of sensitivity and spectral resolution the GPS technique cannot compete with the conventional broadband seismology, or, say, the dozen-or-so superconducting gravimeters that only record the vertical component, in detecting and analyzing the normal modes. Nevertheless, GPS records the “true” 3-D displacement, particularly important in low-frequencies down to static changes and in complete dynamic range unachievable from other instruments. That has potential implications with respect to the GPS seismology (Larson et al. 2003).

We have focused on the eigenfrequencies of the detected normal modes in this paper, and less on their amplitudes and phases (or complex amplitudes). The information on the latter, or more specifically the “true” excited complex amplitudes of the modes, not those in each station record, actually resides in the OSE formulation for Eq. 1. Hence they in principle can be recovered from the OSE stack (cf. Ding and Chao 2017). (The same is true with respect to the SHS stacking for that matter.) That would provide useful information about the source mechanism of the earthquake (cf. Bogiatzis and Ishii 2014; Lentas et al. 2014). In practice, however, that requires a reasonably good global distribution of the network data, which to a large extent is lacking in our present datasets as we have demonstrated. Another interesting prospect that may become possible if, in addition, the GPS time series are truly continuous (without daily breaks) to enable higher frequency resolution is the following. For an earthquake modeled as a double-couple point source, the spherical-harmonic components of the induced displacement field in the epicentral (as opposed to geographical) coordinates for order greater than 2 would vanish identically (Phinney and Burridge 1973). Achievable through a coordinate transformation of GPS data this would “concentrate” the seismic signal in greatly fewer numbers of components (cf. Chao and Liau 2019), hence presumably rendering the OSE procedure more effective in terms of singlet resolution. These await future developments of the GPS or GNSS system with respect to seismic applications.