Detection of Atmospheric–Ionospheric Disturbances in TEC Time Series From Large GNSS Networks Using Wavelet Coherence

Acoustic-gravity waves in the neutral atmosphere induce disturbances in the distribution of electrons in the ionosphere that can be observed in the total electron content (TEC) of global navigation satellite system (GNSS) signals. Large GNSS networks, such as Japan's GEONET, provide spatially dense samples of TEC time series, which can be cross-correlated to detect the coherent structure of these disturbances. Identifying the atmospheric wave-induced perturbations from the short-term background fluctuations can help further our understanding of geophysical processes. In this article, we introduce a method based on wavelet coherence and cross-correlation for detecting and analyzing earthquake/tsunami-induced ionospheric disturbances. We apply this method to detect and isolate the potential disturbances from the data of large GPS networks in Japan, USA, and New Zealand for the 2011 Tohoku-Oki earthquake and the 2015 Illapel earthquake and their resulting tsunamis. We then filter the regions corresponding to each strong coherence structure in time–frequency space to extract the separately identified ionospheric disturbances. The speeds and directions of arrival of these disturbances are found to be compatible with the acoustic-gravity waves from main-shock and after-shock epicenters of each event, offshore tsunami propagation, and seismic Rayleigh waves. Furthermore, we introduce a method to determine the observing heights for the far-field observations.

Abstract-Acoustic-gravity waves in the neutral atmosphere induce disturbances in the distribution of electrons in the ionosphere that can be observed in the total electron content (TEC) of global navigation satellite system (GNSS) signals.Large GNSS networks, such as Japan's GEONET, provide spatially dense samples of TEC time series, which can be cross-correlated to detect the coherent structure of these disturbances.Identifying the atmospheric wave-induced perturbations from the short-term background fluctuations can help further our understanding of geophysical processes.In this article, we introduce a method based on wavelet coherence and cross-correlation for detecting and analyzing earthquake/tsunami-induced ionospheric disturbances.We apply this method to detect and isolate the potential disturbances from the data of large GPS networks in Japan, USA, and New Zealand for the 2011 Tohoku-Oki earthquake and the 2015 Illapel earthquake and their resulting tsunamis.We then filter the regions corresponding to each strong coherence structure in timefrequency space to extract the separately identified ionospheric disturbances.The speeds and directions of arrival of these disturbances are found to be compatible with the acoustic-gravity waves from main-shock and after-shock epicenters of each event, offshore tsunami propagation, and seismic Rayleigh waves.Furthermore, we introduce a method to determine the observing heights for the far-field observations.Index Terms-Global navigation satellite system (GNSS), ionospheric disturbances, tsunami, wavelet analysis.

I. INTRODUCTION
I ONOSPHERIC disturbances associated with tsunamis and other solid Earth events have been studied to understand the ionospheric response to atmospheric waves resulting from solid Earth surface processes.Traveling ionospheric disturbances (TIDs) are manifestations of waves in the neutral ionosphere.
These disturbances may occur spontaneously, or could be induced by atmospheric acoustic-gravity waves (AGWs) or by geomagnetic storms.
Theoretical studies and experimental observations since the early 1960s have improved our understanding of the ionospheric response associated with atmospheric-gravity waves [1].Wave-like disturbances in the ionosphere with periods of approximately 1 min-1 h, and propagation speeds ranging from 50 to a few hundred m/s, have been observed by different techniques, such as ionosondes, Doppler sounding, and global navigation satellite system (GNSS) total electron content (TEC).
More recently, ionospheric disturbances caused by earthquakes [2], mine blasts [3], tsunamis [4], volcanic explosions [5], [6], solar flares [7], chemical explosions [8], underground nuclear tests [9], [10], and rocket launches [11] have been observed using data from dense GPS networks.In contrast to radar sounding, GPS measurements are sensitive to the TEC along the signal propagation path.Methods for detecting ionospheric disturbances using GPS data have included statistical tests on the TEC time series [12], the statistical angle-of-arrival and Doppler method [13], and cross-correlation methods [14], [15].These techniques all approximate ionospheric disturbances as monochromatic or narrowband waves.Bandpass filtering with predefined passbands or polynomial detrending fitting processes can then be applied to extract the perturbation signals from TEC time series.
The aforementioned methods use filtering processes only in frequency space [15] or for detrending polynomial fitting in time domain [14].These methods assume that only a single-source perturbation is observable at a given time, and that the disturbance can be approximated as a quasi-monochromatic wave.However, ionospheric perturbations from multiple sources may be present within the same region simultaneously.Specifically, in the case of natural hazards, such as earthquakes and tsunamis, different ionospheric disturbances can be induced by coseismic AGWs near the epicenter.In previous observations of the 2011 Tohoku tsunami and earthquake events, both surface Rayleigh waves and AGWs from the epicenter and tsunamis were identified in TEC disturbances using various methods [16], [17].TEC depletions as well as large-and medium-scale concentric waves were also observed in TEC after the earthquake event [18].The perturbations in TEC from the Tohoku earthquake and resulting tsunami have been observed in near-field data, such as regions above Japan, and in far-field observations, including those made above the continental United States [19], [20].
For observations that use data from GPS TEC time series, cross-correlation methods are a frequently used technique in studying the natural hazards associated with ionospheric responses (for example, [21] and [22]).However, TEC time series containing different types of ionospheric perturbations not only fail detection with cross-correlation methods due to the assumptions of a single source quasi-monochromatic wave, but the superimposed TEC waveforms result in difficulties of estimating the TIDs propagation speeds [21].The potential to isolate different types of TIDs from a single TEC time series is a key improvement in the detection and identification of atmospheric-induced ionospheric disturbances over crosscorrelation methods.Our previous research [23] demonstrated the improvement in cross-correlation detection of different TIDs by using wavelet coherence analysis to design a filter band to extract perturbations present in different frequency bands.Prior research has applied a similar wavelet method using the Paul wavelet to slant TEC in the case of the Haida Gwaii tsunami in 2012 [24] but did not apply wavelet coherence as a detection method.
Here, we advance the wavelet-based analysis method [23] to reconstruct the TEC perturbation signals directly from wavelet coefficients indicating highly coherent time-space structures.We demonstrate that the wavelet coherence method can successfully detect different TEC disturbances from different perturbation sources.The wavelet coherence method is applied to isolate and detect disturbances in the near-field and far-field networks resulting from both the Japan Tohoku-Oki (11 March, 2011) earthquake/tsunami and the Illapel (16 September, 2015) earthquake/tsunami.Moreover, we introduce a method for estimating the vertical propagation of tsunami-AGW induced TIDs in the far-field regions.
We introduce GNSS TEC, the wavelet transform (WT), wavelet coherence, and the TID detection process in Section II.
We then define a subnetwork-based wavelet coherence analysis for the purpose of isolating potential TIDs within subareas of a large GNSS network in Section II-C.An analysis of the statistical significance of the detection is provided in Section II-D and Section II-E.Then, we apply this method to analyze the ionospheric responses on the earthquake/tsunami events in Section III.Finally, Section IV concludes this article.

A. GNSS Observations of AGW-Induced Ionospheric Disturbances
Perturbations of TEC measurements have been commonly used to study the AGW-induced ionospheric disturbances [25].The TEC of radio signals along the propagation path is defined as the total number of electrons within a 1-m 2 column around the path of radio signal propagation where n e (l) is the electron density at a distance l along this path.TEC can be calculated from GPS phase measurements from a dual-frequency GNSS receiver using the method developed in [26].An early mathematical representation of TEC fluctuations in terms of atmospheric wave and ionospheric parameters was introduced by Georges and Hooke [27].

B. Wavelet Transform
The WT is a commonly used time-frequency analysis technique for examining frequency contents locally in time.The WT can be expressed as inner products of a signal and a set of analyzing functions (or coherence states [28]) to measure their similarity locally in the time-frequency domain.As compared with the shifted and modulated window function used in the short-time Fourier transform, the WT uses a window function ψ, also known as the mother wavelet function, as a basis for the signals.The continuous wavelet transform (CWT) W h (a, b) of a TEC time series I h (t) from receiver h can be expressed as inner products of I h (t) with a set of dilated (by wavelet scale a) and translated (by localized time parameter b) daughter wavelet function ψ as follows: We choose the Morlet wavelet in which the frequency parameter, m, adjusts the time and scale resolution and controls the number of oscillations present.A proper value of m will ensure a good frequency resolution.The correction term e −m 2 /2 can also be assumed to be negligible when m is greater than 5.In addition, the Morlet wavelet function with m >= 6 satisfies admissibility, dω < ∞, which is a requirement for the wavelet reconstruction ( [29], as described in Section II-E).Moreover, the Morlet wavelet (with m = 6) has a special feature in that the inverse of the wavelet scale corresponds to frequency in Hz [30].For these reasons, we choose m = 6 to have a balance in time-frequency resolution for isolating different types of TIDs induced by atmospheric acoustic and gravity waves from tsunamis and earthquakes.The CWT of a discrete TEC time series I h [k] = I h (kT ) with a Morlet wavelet and a sampling time T (in seconds) can be expressed as We evaluate this CWT at discretized scale values given by a j = a 0 2 j/M , with index j = 0, 1, . .., J and a 0 representing the smallest scale.J can be found from the length and the sampling rate of the TEC time series, as given in [30].M = 2 provides adequate sampling in scale with the Morlet wavelet, but we choose M = 4 for finer resolution.It is important to note that the CWT using the shifted and scaled Morlet wavelets does not represent a projection of the original time series onto an orthogonal basis as the scaled and shifted versions of the Morlet mother wavelet are not orthonormal.The reconstruction from Fig. 1.Signal processing algorithm for TID detection and propagation velocity estimation.First, the GNSS network is divided into multiple subareas.Next, the wavelet coherence within each subarea is computed and the detection threshold is applied to detect coherent structures.The detection results are then applied to a wavelet filtering process on the TEC time series.Wave propagation velocity and direction estimation is then performed on the filtered time series using cross correlation methods.
wavelet coefficients using Morlet wavelets can carry information from a future point in time and, therefore, is noncausal.
From the definition of the WT, the area of each timefrequency atom remains constant [29].The time width (aΔt) and frequency bandwidth (Δω/a) of each atom is tuned by the scale parameter a corresponding to the wavelet functions ψ( t−b a ).Thus, by adjusting the sharpness using the scale a, the WT can provide adequate resolution at higher and lower frequencies.
Our algorithm applying the WT to detect disturbances in the TEC time series is shown in Fig. 1.
A large and dense GNSS network is first partitioned into several subareas (see examples in Fig. 2).We approximate the TIDs as planar waves propagating through each small subarea.TEC time series from different receiver stations are first detrended by removing a linear trend and then a fourth-order polynomial fit.The CWT is then used to transform the detrended TEC time series into time-frequency space at discrete steps in a and b, referred to as wavelet coefficients. 1Wavelet coherence is computed between pairs of TEC time series within each subarea and a threshold based on noise properties is set for detection of coherence between pairs of signals.In each small subarea, wavelet coefficients exceeding this threshold indicate the presence of coherent signals in multiple time series.The corresponding disturbance signals are then reconstructed using only these wavelet coefficients.A cross-correlation technique [15] is then applied to the reconstructed TEC signal in the time series to estimate the propagation direction and speed of the corresponding disturbance.A new approach (see Section II-H) is defined for estimating the height and the vertical propagation speed of the disturbance.The complete procedure is described in more detail in the following sections.

C. Subdivision of Network Into Subareas
We first divide the full GNSS network into many small subareas, assuming that each disturbance can be approximated as a planar wave within each subarea.The coherence detection and estimation are performed independently for each subarea.For example, we divided 1235 Japan GEONET stations into 32 subareas (approximately 1 • × 1 • grid), selected to provide at least 50 stations within each subarea, shown in Fig. 2(a).Wavelet coherence analysis is performed between the TEC time series from two small subsets of GPS stations within each subarea.The resulting coefficients are used to design a filter that is applied to all of the stations in the subarea.

D. Wavelet-Based Coherence Detection
The wavelet coherence is computed to quantify the presence of coherent structures in multiple TEC time series from different stations observing the same satellite within each subarea.In common practice, a time and scale smoothing (or averaging) process is applied in calculating the wavelet coherence quantities to examine the locally linear coupling of two time series [30], [31].However, in [30], Section 5], it is concluded that the smoothing process would destroy the χ 2 distribution of the wavelet cross spectrum, which is used for significance testing.This makes a smoothed wavelet cross spectrum difficult for significance testing.Instead, we defined a wavelet coherence based on the magnitude-squared coherence (MSC) [32] of multiple TEC time series within a local dense GNSS network.In effect, this replaces time averaging by an averaging over different station pairs, which are assumed to be observing the same disturbance.The distribution of the MSC under white noise can be found [32].
First, we select two subsets of GNSS stations within each subarea.Each subset contains GPS stations distributed within a small area of radius r s = 30 km.We assume that the selected TEC time series within each of the two subsets in each subarea are different realizations of the time series containing common disturbance structures (see Fig. 3).MSC of the N s station-pairs in subsets s1 and s2 within each subarea is computed as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where W h (a, b) is the wavelet coefficient defined in (2), and the subscript (c, si) represents the receiver in the cth station-pair from subset si.MSC measures similarities of wavelet coefficients of the TEC time series, at each localized time and frequency sampling point.In the lower frequency range (i.e., larger scale), the time length of an atom increases.Due to the flexible adjustment of the frequency (or scale) and time length, the wavelet MSC is able to capture TIDs induced by acoustic fields (infrasound signals) and gravity waves simultaneously [23].However, it should be noted that the MSC is not able to measure the similarities of disturbance signals outside of the time length of each atom.A consequence of this is that disturbances that are not within the same subarea and the same time window will be incoherent by definition.

E. Threshold and Significance for Coherence Detections
In wavelet analysis, the cone of influence (COI) sets boundaries in the time-frequency space where coherence spectra are not subject to edge effects produced as a result of finite TEC time series.For the Morlet wavelet in this study, the COI for the singularity at the edge drops by a factor √ 2/f [30].To reduce the likelihood of false detection, a threshold was set based upon a significance level.Following the definition of generalized coherence [32], the generalized coherence detection threshold (λ α ) at the significance level α with using k stations is defined as 1) .A 5% significance level is applied for the detection process, corresponding to a 95% confidence interval in statistics.Furthermore, in [30] and [33], a significance test is used to evaluate the wavelet coherence spectrum against background noise for distinguishing true features from spurious results based on the null hypothesis of the spectrum being red or white noise, which can be simply modeled as a lag-1 autoregressive process [AR (1)] with an AR coefficient of 1 or 0, respectively.In this work, we test against the null hypothesis of red noise.Statistical tests for wavelet coherence using the abovementioned assumptions have been conducted in previous research [33], [34], [35].

F. Filtering Process and Signal Synthesis
After finding coherent structures that are above a desired threshold, λ α , in a subset of GNSS network measurements, a filtering and reconstruction process is applied to reconstruct disturbance signals corresponding to each identified high coherent structure.These structures are interpreted as different types of TEC perturbations in the wavelet domain.According to a study of time-frequency analysis [28], if a signal is essentially concentrated on a limited region in time-frequency space (or phase plane), then only the coefficients of its WT corresponding to the local time-frequency points within or close to the region are needed for an reconstruction of the signal.The waveletbased filtering approach used here is based on this concept.The relevant wavelet coefficients corresponding to time-frequency locations within or close to the high-coherence region on the phase plane are selected and used to reconstruct, via an inverse WT, the filtered disturbances present in each high-coherence region.
We define the selected wavelet coefficients (from the Morlet WT of original time series of receiver h) for the reconstruction of the desired signals as  where R 2 (a, b) is the wavelet MSC from (5) computed on measurements from the subset of N s station-pairs.Filtering [defined by (6)] is applied to all measurements within the subarea.In addition, these wavelet coefficients for an N -point TEC time series and J scales can be efficiently calculated using a fast-Fourier transform with computation complexity O(JN(logN ) 2 ) [36].
After the desired wavelet coefficients are selected, Morlet's wavelet reconstruction, introduced in [30] and [37], is applied.Since we discretized the scale as a j = a 0 2 j/M , the Morlet's Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 6.Solid line represents the averaged waveforms of all TEC disturbance waveforms (blue dot) aligned with the time delay of each TEC time series measurements.The stacked waveforms (the solid line) reconstructed from the wavelet coefficients corresponding to regions 1-5 in Fig. 5 with high coherence quantities above the threshold in phase space at the corresponding times of occurrence.Blue dots represent the quantities of TEC perturbations from each pair of the PRN15 satellite with receivers in subarea 08.reconstruction of the signal Ĩ[k] of [30, Section 3] is given as follows: where R( Wh (a j , kT )) is the real part of the wavelet coefficients at time step k with sampling period T , R(W δ (a j , kT )) are the real parts of wavelet coefficients for a delta function centered at the origin using the same wavelet basis, and (a j ) is discretized scale j = 0, 1, . .., J.
Fig. 4 shows an example of the wavelet coherence spectrum of the GEONET subarea 08.The time-frequency spectrum shows that the energy of the coherent disturbances was concentrated in the two dominant frequency bands spreading locally in time.The time-frequency regions corresponding to these high coherent structures are shown in Fig. 5.Each small blue circle corresponds to a wavelet coefficient.There are five groups classified from the wavelet coefficients.We choose these five groups to demonstrate wavelet filtering.We reconstruct the disturbance signals from TEC time series collected in the subarea 08 based on the five high-frequency band coherent structures, identified as numbers 1-5 in Figs. 4 and 5.The reconstructed waveforms are shown in Fig. 6, in which the reconstructed waveforms are normalized by the maximum magnitude of the TEC perturbation observed in each time series.

G. Horizontal Propagation Speed and Direction Estimation
After detecting coherent ionospheric disturbances through comparing wavelet coherence against a threshold, we apply the method introduced in [15] to estimate the propagation speed and direction of the disturbances.The assumptions implemented in this approach are: planar wave propagation, flat Earth, and a thin-shell ionosphere located at the height of the maximum electron content.The height of a thin-shell ionosphere was set to the height of maximum electron density determined from the IRI model 2 on the event dates.
As shown in Fig. 7, Δ r pp is the relative position vector between r 1 (t 1 ) and r 2 (t 2 ), the ionosphere pierce point (IPP) trajectories of a station pair viewing the same satellite from stations 1 and 2. Station-pairs are chosen for performing wavelet coherence in two clusters s1 and s2 on approximately on opposite corners of the subarea, as shown in Fig. 3.After thresholding (6) and filtering (7), the cross-correlation is done with all the station pairs in a subarea.The time-lag Δt max between the 2 [Online].Available: http://nssdcftp.gsfc.nasa.gov/models/ionospheric/iri/corresponding pairs of filtered TEC time series I 1 (t) and I 2 (t) obtained from cross-correlation is assumed to be equal to the travel time of the propagation wavefront between the two IPP positions r 1 (t 1 ) and r 2 (t 2 ).Each station pair produces one component of the slowness vector.The resulting propagation direction and velocity over the entire subarea can be estimated from the complete set of measurements from all station pairs using a linear least-squares method [15].In principle, this method could fail when the relative pierce point velocity is perpendicular to the propagation direction k of the TID.This method predicts no Doppler effect when the relative pierce point velocity is perpendicular to the propagation direction k, but a 20% mean error can be introduced in the propagation velocity estimation due to Doppler effects [14].When this technique is applied to a large and dense network, such as GEONET in Japan, the full network can be divided into smaller subareas (over which the planar wave propagation assumption holds) with one velocity estimate obtained for each subarea.Variation of propagation direction/speed over a large area is, thus, visualized through comparing results from each subarea, as shown in Figs.8-10.

H. TID Altitude Estimation
This section describes our approach to determine the observation altitude of TIDs in the far field and to estimate their vertical propagation speeds.Observations of TIDs' vertical propagation not only provide an understanding of the evolution of tsunami-induced AGW and TID interactions, but are necessary to understand the inversion mechanism between tsunami wave heights and the amplitude of their corresponding TIDs since the amplitude of a TID depends on altitude [4].Here, we introduce a method to estimate the vertical propagation speeds of TIDs observed in the far-field networks and apply it to the case of the 2011 Tohoku-Oki tsunami and the 2015 Illapel tsunami as viewed by the U.S. Plate-Boundary Observation (PBO) network.In both cases, the tsunami waves propagated over the Pacific ocean and were stopped and reflected by the U.S. west coast.The associated AGWs generated by these tsunamis continued propagating inland from the U.S. west coast.For the far fields in U.S. west coast, we assume that these AGWs propagate as planar waves over each small subarea.Under this assumption, the IPP locations corresponding to the wavefront of TIDs should align as a planar wave within the subarea if TID altitude is determined correctly.Fig. 11 shows, for the 2011 Tohoku-Oki tsunami, the intersection points of IPP trajectories and TID wavefront locations assuming two different TID observation altitudes (300 and 400 km).Significant wavefront distortion is seen when the assumed TID altitude is set to 300 km.We vary the ionospheric height, H i , in subarea i from the bottom to upper of ionospheric altitudes as an input parameter in our velocity estimation post wavelet-based detection, and we minimize a least squares cost function, C(H i ), corresponding to the residuals from the velocity estimation for each value of height.The cost function, derived from [15,Eq. 10], is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where V (H j ) is the estimated velocity for the entire wave in subarea j, and Δ r i (H j ) are the individual estimated IPP displacements corresponding to the lagged observation of the waves with lag Δt i for each of the n IPP pairs within subarea j.We use a grid search method to minimize the cost function over a set of selected heights to estimate the TID altitude.An observed result of this optimization is that the resulting IPPs corresponding to the wavefront become planar, as seen by the improved collinearity in the IPPs (black versus pink) in Fig. 11.By tracking the vertical height of the wave across different subareas as it evolves, the vertical TID velocity of plane waves can be estimated in the far field from differences in estimated heights and observation times of heights corresponding to adjacent subareas.

A. 2011 Japan Tohoku-Oki Earthquake and Tsunami
On 11 March, 2011, an M w = 9.0 earthquake occurred at 05:46:23 UTC near the east coast of Honshu, Japan (38.322 • N, 142.369 • E).The earthquake triggered a tsunami with run-up reaching 10 m in Sendai.The first tsunami wave hit the west coast of the United States about 9 h after the main shock.The earthquake and the resulting tsunami comprised one of most disastrous natural hazard events in recent history.Numerous studies have been conducted to characterize this hazard [38].We applied our wavelet coherence method to detect TEC perturbations in data from large GPS/GNSS networks in Japan, the west coast of the United States, and New Zealand (see Fig. 2).For this event, we considered Japan GEONET as the near-field network, and the other two networks as far field.
1) Near-Field Results: Wavelet coherence results in Fig. 4 show structures in two dominant frequency bands (0.9-2.0 and 2.0-8.3 mHz) after the main shock in all subareas in Japan.
In the near-epicenter fields, the wavelet-based detection method is able to identify and isolate different types of coherent disturbances around the epicenter, as shown in examples in Fig. 4. In each case, the coherence spectrum reveals five common structures numbered as Groups 1-5 in Fig. 4. Short-period TIDs observed in the measurements from PRN15 and PRN18 were found to propagate with speeds in two different ranges; 1.5-2.8km/s (Group 1) and 652-920 m/s (Groups 2 and 3).The resulting velocity vectors for each subarea are shown in Figs.8-10.Group 1 TIDs propagate in the northwest and southwest directions, and their fundamental periods, propagation speeds, and directions are in agreement with the Fig. 11.Alignment of IPPs from TID detection under appropriate choice of TID height (pink) versus unsuitable choice of TID height (black).The corresponding IPP paths are represented with blue and green, respectively.Notice that under the optimal choice of TID height, the IPPs represent a planar wave.motions of the ground surface waves [39].We suggest that the Group 1 TIDs are associated with the seismic Rayleigh waves.
Group 2 TIDs have propagation speeds (≈700 m/s) consistent with the acoustic speeds in the ionosphere at around 350 km.They propagate radially from the main-shock location, as shown in Fig. 9. Group 3 TIDs, as shown in Fig. 10, were observed about 12 min after the aftershock (36.27 • N, 141.14 • E) at 06:15 UTC on 11 March, 2011).Their propagation direction radiates out from the aftershock location at approximately 850 m/s.
The lower frequency band (Groups 4 and 5) consists of disturbances that propagate with speeds ranging between 101 and 265 m/s.These were observed within 30 min of the main shock, which is consistent with the time of the first tsunami arrival (first recording of the tsunami was 25 min after main shock at DART buoy 214183 ).Their propagation speeds and directions are consistent with the recorded tide gauge measurements and tsunami waves reconstructed by MOST4 [40].The propagation velocities and the times of observation for the TIDs resulting from Rayleigh, acoustic, and tsunami waves over Japan are also consistent with the parameters observed in previous TEC analysis, such as [16].The observed TID velocities for the Rayleigh and tsunami waves agree with the corresponding velocity estimates from GPS radio occultation observations on FORMOSAT-3/COSMIC [41].
2) Far-Field Results: Far-field disturbances were observed in data from the PBO network located in the western United States.Fig. 12 shows the TID velocity vectors estimated using PRN21.The estimated averaged speeds of the TIDs are about 200 m/s, consistent with tsunami speeds near the U.S. west coast.Propagation directions are also in agreement with the tsunami waves.The bottom right of Fig. 12 also shows the comparison between the GPS observations of TIDs and the MOST model simulations of the tsunami waves.This comparison indicates that the occurrence times of observed ionospheric disturbances are consistent with the arrival of the tsunami waves.These results are also in agreement with the tsunami propagation speed predicted from bathymetry data5 (≈197 m/s) and far-field arrival times, periods, and velocities published in [20].
For the far-field New Zealand data, we divide the GPS network into the four subareas shown in Fig. 2(d).Disturbances with frequencies in the range 0.6-2.0mHz were detected in data from PRN16 and PRN23.Horizontal velocities of the disturbances were estimated to be in the range 80-328 m/s, having directions consistent with tsunami waves.These results are shown in Fig. 13.
Dependence of the tsunami speed on the square root of the ocean depth may account for differences in propagation direction between the between coastal area and the open oceans.The effect of bathymetry can be seen in the variation propagation directions estimated using the New Zealand Geonet for this event (see Fig. 13).While this network is in the far field where a single planar wave is expected, the ocean bathymetry results in multiple tsunami wavefronts with their corresponding TIDs arriving at different times from different directions.
3) TID Altitude Estimation: Fig. 14 shows the TID tion heights of the 2011 Tohoku-Oki tsunami results estimated by the method introduced in Section II-H from PRN15, PRN21, and PRN22 GPS data using the PBO network in the U.S. west coast.The surface is plotted using MATLAB's surf command with the "FaceColor" setting set to "interp."Ionospheric disturbances were found to propagate more than 2000 km inland after the tsunami hit and stopped at the coast shown in Fig. 14.The first-arrival TIDs were observed between 16:00 and 17:30 UTC after the tsunamis waves propagated to the U.S. west coast.Disturbance height was estimated independently in each subarea in Fig. 14 using the method in Section II-H.

B. 2015 Illapel Earthquake and Tsunami
On 16 September, 2015, an M w = 8.3 earthquake occurred at 22:54:32 UTC 48 km west of Illapel Chile (31.573 • S, 71.674 • W), at a depth of 22.4 km.We applied our wavelet coherence method to detect TEC perturbations in data from GPS/GNSS networks in Japan, the west coast of the United States, Chile and Peru, and New Zealand (see Fig. 2).For this event, we considered the stations in Chile and Peru as the near-field network, and the other two networks as far field.
1) Near-Field Results: We analyzed GPS measurements made at 15-s intervals from the 42 stations in the four (4) subareas in Chile and Peru (see Fig. 15).Significant TEC perturbations were observed in two dominant frequency ranges: 1.0-4.0 and 2.8-8.3 mHz, as shown in Fig. 16.The short-period perturbation is detected about 8-15 min after the main shock with an estimated horizontal propagation speed of 417 m/s (18% error due to the fewer stations).The speed and frequency indicate that the short-period disturbance could be associated with seismicinduced atmospheric acoustic waves.Long-period waves were observed about an hour after the main shock with an averaged horizontal speed ≈118 m/s (19% error).The uncertainty is as calculated using the method in [15].Characteristics of the long-period perturbations are consistent with tsunami-induced gravity wave propagation.
2) Far-Field Results: We present far-field results using the GPS networks in Japan, U.S. west coast, and New Zealand (see Fig. 2).Disturbances within the frequency range of 0.6-2.0mHz were detected.Fig. 17 shows the results for the western United States.Disturbances were detected between 11 and 13 UTC with horizontal velocities from 59 to 319 m/s in directions consistent with arrival Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.    of the tsunami waves.The observed TID heights ranged from 280 to 580 km.Fig. 18 shows the far-field results from the New Zealand GEONET network using data from PRN 24.We detected TIDs in New Zealand at 12.75 UTC on 17 September, 2015.The observed TID propagation speeds of ≈342 m/s and directions toward the northwest are in agreement with the tsunami direction and characteristics.
At the time of tsunami arrival in Japan, the wavelet analysis indicates high-coherent GPS-derived TEC perturbation (PRN26 data) detected with frequencies from 0.6 to 2.1 mHz between 20.5 and 22.5 UTC on 17 September, 2015.The estimated propagation speeds are in the range of 116-462 m/s, as shown in Fig. 19, at ionospheric height from 220 to 300 km.The direction of propagation is to the northwest.
3) TID Altitude Estimate: Results of GPS PRN 9 data analysis indicate that the first-arrival TIDs were detected at 10.5-12 UTC on 17 September, 2015.These TIDs propagated in the direction of tsunami waves through the U.S. west coast regions.Similar to the 2011 Tohoku-Oki analysis, Fig. 20 shows the observed TID altitudes (from Section II-H) vary from 280 to 580 km as TIDs propagate from the south to the north of California.In the case of the New Zealand GEONET data, the TID observation altitudes vary between 260 and 320 km.

IV. CONCLUSION
In this article, we have presented an application of wavelet coherence to the detection of ionospheric disturbances in dualfrequency GPS network data and the estimation of propagation speed and direction estimation.The significance of this advancement is its capability to detect and distinguish multiple different TIDs occurring within a common spatial region and time window.Using this method, we isolated and identified TIDs corresponding to seismic Rayleigh waves, acoustic waves, and tsunami-induced AGWs in both near-field and far-field data for two earthquakes and their resulting tsunamis: the 2011 Japan Tohoku-Oki event and the 2015 Illapel event.In the case of the 2011 Tohoku earthquake and tsunami, the TID detection times, frequencies, and velocities from our wavelet coherence method agree with previous results in both the near-field and far-field cases, as shown in Section III.A summary of the detection and velocity results for both events in near and far-field observations are given in Table I.
We also provided a novel method, using a planar wave model and a least-squares fit, to estimate the TID altitude variation in the far-field network (U.S. PBO network) as the waves propagate inland for both events.Our results show an increase in TID altitude as the wave propagates inland in both cases.The estimated vertical velocities agree with the vertical propagation of tsunami-induced AGWs using the equation for vertical group velocity provided by Occhipinti et al. [42], [43] evaluated at the corresponding TID observation heights.The high-altitude wind, mountain waves, and locally atmospheric temperature variations, which are not studied in this research, could affect AGW propagation.An example of the effects of mountain waves can be found in [44].

Fig. 2 .
Fig. 2. In total, four GNSS networks are used in this work: (from top left, clockwise) U.S. PBO, Japan GEONET, New Zealand Geonet, and Chile network.We define each grid as a subarea shown in the black colored boxes.The blue-colored triangles represent the locations of GPS stations.

Fig. 3 .
Fig. 3. Illustration of GNSS locations, a subarea, and subsets inside the kth subarea.Within subsets 1 and 2, the indices are given as [c, s1] and [c, s2], respectively.c represents the index of c pair.In this research, we choose c = 1, 2, 3.In addition, the center of subsets 1 and 2 are set up at (d 1 /4, d 2 /4) and (3d 1 /4, 3d 2 /4), respectively, where d 1 and d 2 represent the sides of the rectangular subarea.

Fig. 4 .
Fig. 4. Example of the wavelet coherence from GEONET subarea 08 on 11 March, 2011.The joint time-frequency spectrum shows that the energy of the coherent disturbances was concentrated in two dominant frequency bands spreading locally in time.The dashed white curve shows the distributions of COI for the singularity at the edge drops by a factor √ 2/f .Regions above these two curves ensure that the edge effects are insignificant.The black contour lines represent the wavelet coherence above 5% significance level (95% confidence interval).

Fig. 5 .
Fig.5.Plot of the sample point locations (the blue circles) from GEONET subarea 08 (see Fig.4) corresponding to regions with coherence above the threshold in time-frequency space, as shown in Fig.4.In total, five high-coherent regions are numbered for estimating their characteristic independently.

Fig. 7 .
Fig. 7. Geometric relationships among the TID propagation direction k, IPP trajectories, and propagation displacement Δ r pp between two adjacent wavefronts.

Fig. 8 .Fig. 9 .
Fig. 8.Estimated wave propagation directions and speeds (shown as blue vectors), corresponding the first arrived high-frequency coherent disturbance signals (Group 1), as shown in Fig. 4, from GEONET subareas 8, 9, 10, 11, and 12 measurements for the 11 March, 2011, Japan earthquake and tsunami event.The starting point of the disturbance speed vectors are located at the estimated position of the disturbance in the ionosphere at the time of maximum TEC amplitude.The bottomright corner shows the snapshot GPS-derived ground surface motions conducted by JPL ARIA team (http://www.tectonics.caltech.edu/slip_history/2011_taiheiyooki/Displacement/)and NOAA's MOST tsunami waves at 05:48 UTC on 11 March, 2011.

Fig. 10 .
Fig. 10.Estimated wave propagation directions and speeds (shown as blue vectors), corresponding the third arrived high-frequency coherent disturbance signals (Group 3) as shown in Fig. 4, detected in GEONET subareas 8, 9, 10, 11, and 12 measurements for the 11 March, 2011 Japan earthquake and tsunami event.The starting point of the disturbance speed vectors are located at the estimated position of the disturbance in the ionosphere at the time of maximum TEC amplitude.

Fig. 12 .
Fig. 12.Estimated wave propagation directions and speeds (shown as blue vectors), corresponding to the far-field PBO network using PRN 21 for the 2011 Tohoku-Oki tsunami.Inlaid figure: The color-coded background shows tsunami waves modeled by MOST, and the colored dots represent TEC perturbation calculated at each IPP point.The blue dashed lines represent geomagnetic field lines.Comparison indicates that the occurrence times, wavefronts and directions of the observed TIDs are consistent with the arrival of tsunami waves reaching the west coast of the United States.

Fig. 13 .
Fig. 13.Estimated wave propagation directions and speeds (shown as blue vectors), corresponding to PRN 16, 20, and 23 data collected at stations inside four grids of New Zealand GEONET network for the 2011 Tohoku-Oki earthquake.

Fig. 14 .
Fig. 14.Color-coded surface shows the 2-D interpolation, using MATLAB's surf command (with the "FaceColor" setting set to "interp") of TID observing heights for the far-field observations of the 2011 Tohoku event estimated using GPS data in each grid (black boxes) and the least squares method of minimizing cost function of 8.The black stars represent the locations each grid-based TID estimate based on the corresponding IPPs.The selection criteria of each grid are more than GPS 25 stations (blue dots) inside a 2 • × 2 • area.

Fig. 15 .
Fig. 15.Estimated wave propagation directions and speeds (shown as blue vectors), corresponding to PRN 12, 24, 25, and 29 data collected at the 42 stations in Chile and Peru area.We defined the four subarea (black colored grids) for the wavelet detection and velocity analysis.

Fig. 16 .
Fig. 16.Wavelet coherence analysis of the subarea 2 using PRN 25 data in the near field (Chile and Peru network) for the 2015 Illapel earthquake.

Fig. 17 .
Fig. 17.Estimated wave propagation directions and speeds (shown as blue vectors), corresponding to the far-field TID velocities detected using GPS data in USA PBO network for the 2015 Illapel tsunami.The right-hand side plot shows a snapshot of NOAA's tsunami wave propagation over the same region.

Fig. 18 .
Fig.18.Estimated wave propagation directions and speeds (shown as blue vectors), corresponding to the far-field TID velocities detected using GPS data in the New Zealand Geonet network for the 2015 Illapel tsunami.

Fig. 19 .
Fig. 19.Estimated wave propagation directions and speeds (shown as blue vectors), corresponding to the far-field TID velocities (blue colored arrows) detected using PRN 26 data of Japan GEONET for the 2015 Illapel tsunamis.

Fig. 20 .
Fig. 20.TID height estimate determined using the method introduced in Section II-H for the 2015 Illapel tsunami over the U.S. PBO network.The color-coded surface shows the 2-D interpolation, using MATLAB's surf command, of TID observing heights estimated using GPS data in each grid (black boxes) and the least squares method of minimizing the direction error cost function.The black stars represent the locations each grid-based TID estimate based on the corresponding IPPs.The selection criteria of each grid is more than GPS 25 stations (blue dots) inside a 2 • × 2 • area.
Detection of Atmospheric-Ionospheric Disturbances in TEC Time Series From Large GNSS Networks Using Wavelet Coherence Yu-Ming Yang , Member, IEEE, Abi Komanduru , and James Garrison , Fellow, IEEE

TABLE I SUMMARY
OF TID DETECTIONS PRESENTED IN SECTION III