Ionospheric Responses to the Earthquake in the Gulf of Alaska and the Kusatsu-Shiranesan Volcanic Eruption on 23 January 2018

Abstract

Numerous research revealed a strong association between the ionospheric perturbations and various natural hazards. The ionospheric measurements from Global Navigation Satellite System (GNSS) observations provide the state of electron contents in the ionosphere that contributes to investigate the source events. In this study, two geophysical events occurred on 23 January 2018, the 7.9 Mw earthquake in Alaska and Kusatsu-Shiranesan volcanic eruption in Japan, are examined to characterize the fingerprint of each event in the ionosphere. Firstly, we extracted the Total Electron Content (TEC) from GNSS measurements, then isolated disturbed wave signatures from the TEC measurements that is referred to as a traveling ionospheric disturbance (TID). As TIDs are short-term ionospheric variations, the major trend of GNSS TEC measurements should be properly removed. We applied a natural neighbor interpolation method together with a leave-one-out cross validation technique for detrending. After detrending the TEC, the remaining signals are further enhanced by applying a band-pass filter and TIDs are detected from them. Finally, the detected TIDs are verified as the response of the ionosphere to Kusatsu-Shiranesan volcanic eruption and Gulf of Alaska earthquake which propagated through the ionosphere with an average velocity of 530 m/s and 724 m/s, respectively. In addition, a coherence analysis is conducted to discriminate between the signatures from a volcanic explosion and an earthquake. The analysis reveals the TID waveforms from each single event are highly correlated, while a low correlation is found between the TIDs from the earthquake and explosion. This study supports the claim that different geophysical events induce the distinctive characteristics of TIDs that are detectable by the ionospheric measurements of GNSS.

1. INTRODUCTION

The theoretical and experimental studies have shown that the natural and man-made hazards can trigger abnormal disturbances in the electron concentration of the ionosphere that are known as the traveling ionospheric disturbance (TID). The mechanism of hazard-induced anomalies has been described by the wave structures, which propagate through the ionosphere, couple with its ionized or neutral contents, and result in generation of TIDs. Depending on the characteristics of the induced perturbations in the ionosphere, the TIDs can be distinguished as large-scale, medium-scale, and short-scale (Jacobson et al. 1995, Rieger & Leitinger 2002, Vlasov et al. 2011). The large-scale TIDs (LSTID) propagate with high velocity about 1 – 2 km/s and duration of more than 1 hour. The velocity of medium-scale TIDs (MSTID) varies between 200 m/s – 1 km/s with a period of about 10 – 60 minutes. The propagation velocity and duration of short-scale TIDs (SSTID) is usually lower than 200 m/s which lasts for a few minutes.

To characterize TIDs excited by natural or man-made hazards, a variety of ionospheric monitoring techniques can be applied including high-frequency Doppler sounding, ionosonde, Faraday rotation measurements extracted from polarized electromagnetic signals, and more recently the Global Navigation Satellite System (GNSS). Among all the above-mentioned techniques, GNSS has become highly affordable (e.g., Mannucci et al. 1998, Occhipinti et al. 2013, Park et al. 2014, Reddy & Seemala 2015, Savastano et al. 2017) due to its worldwide availability, high sampling rate of measurements, and improved capability in modernized GNSS offered by multiple constellations.

The ionospheric perturbations caused by strong earthquakes has been investigated for several decades (Calais & Minster 1995, 1998, Afraimovich et al. 2001, Artru et al. 2004, Astafyeva & Afraimovich 2006, Shinagawa et al. 2007, Choosakul et al. 2009, Jin et al. 2010, Heki 2011, Astafyeva et al. 2013, Occhipinti et al. 2013, Park et al. 2014, Catherine et al. 2017). Calais & Minster (1995) applied the Global Positioning System (GPS) derived TECs to investigate the ionospheric impacts due to 1994 Northbridge earthquake. They found the generated shock-acoustic waves propagated upward from the focal area in a narrow cone of zenith angles. Shinagawa et al. (2007) explained the earthquake-induced ionospheric variations by resonance of acoustic waves in the vicinity of the epicenter. Choosakul et al. (2009) reported a periodic TEC fluctuation after the 2004 Sumatra earthquake. Astafyeva et al. (2013) related the intensity of TEC perturbations to the magnitude of earthquakes where the greater the magnitude, the more intense TEC disturbance expected. Perevalova et al. (2014) stated that the earthquakes with magnitude of greater than 6.5 Mw can induce significant disturbances in the ionosphere based on their findings from the earthquakes between the years 1965-2013. Catherine et al. (2017) found the ionospheric disturbances affected by the propagation direction of 7.8 Mw 2015 Gorkha earthquake and demonstrated a consistency between the TEC waveforms and the focal mechanism of the earthquake.

Besides the earthquakes, volcanic eruptions can generate noticeable variations of the electron density in the ionosphere. The ionospheric response to volcanic eruptions was first revealed by Cheng & Huang (1992) for 1991 Mount Pinatubo explosion. Cheng & Huang (1992) used ionospheric sounding data to detect volcano-induced perturbations which were characterized by a propagation velocity of 131 – 259 m/s and period of 16 – 32 minutes. Since then, several studies also detected the ionospheric fluctuations caused by volcanic eruptions (de Ragone et al. 2004, Heki 2006, Dautermann et al. 2009, Komjathy et al. 2012, Lin 2013). de Ragone et al. (2004) studied 13 different volcanic eruptions of level 2 – 5 occurred within the years of 1960 – 1980. They identified a decrement in hourly F2 critical frequency (foF2) and an increment in F layer virtual height (h’F) data for 60% of the cases and concluded that the detection of eruption-induced wave structures is highly dependent on the volcanic intensity. By utilizing GPS-derived TECs, Heki (2006) detected ionospheric disturbances about 12 minutes after the 2004 Asama Volcano eruption where the TIDs propagated with a period of 1.25 minutes and velocity of 1.1 km/s. Dautermann et al. (2009) observed the atmospheric perturbations induced by 2003 Soufrière Hills volcano about 18 minutes after the primary eruption with spectral peaks at 1 and 4 mHz and horizontal velocity of 624 m/s. Lin (2013) applied two-dimensional principal component analysis to detect the TEC anomalies induced by 2013 Mexico's Colima volcanic eruption.

In addition to natural hazards, man-made events also lead the ionospheric perturbations that are shown in (Park et al. 2011, Yang et al. 2012, Zhang & Tang 2015, Savastano et al. 2019). Park et al. (2011) investigated the response of ionosphere to the 2009 North Korean underground nuclear explosion (UNE) and suggested the GPS-derived TEC as a reliable tool for determining its epicenter and detecting the ionospheric disturbances. Yang et al. (2012) applied wavelet analysis to localize the TIDs triggered by the 2006 and 2009 North Korean UNE. Their experiments revealed that only the 2009 explosion induced high-frequency perturbations with characteristics of 2 – 5 minutes in duration and 297 – 1322 m/s propagation velocity, while both explosions induced low-frequency disturbances with period of 3 – 12 minutes and velocity of 75 – 453 m/s. Zhang & Tang (2015) also verified two different types of TIDs induced by the 2009 North Korean UNE, namely high-frequency TIDs with a period of 1 – 2 minutes and velocity of 95 – 586 m/s and low-frequency TIDs with a period of 2 – 5 minutes and the velocity of 92 – 483 m/s.

TIDs can be generated by various events besides the abovementioned ones that add challenges to distinguish the TIDs associated with a particular event. Indeed, the distinctive mechanism of geophysical hazards releases different energy and wave structures and the ionospheric response to different events may vary. Therefore, it is critical to discriminate TIDs induced from different types of events. In this study, we selected two events occurred on the same day. The Gulf of Alaska earthquake and Kusatsu-Shirane volcanic eruption, both occurred on 23 January 2018. Because they occurred on the same day, the background ionospheric condition is nearly same in a global scale, and this allows a fair comparison between TIDs occurred by two events. The following sections describe the signal processing methods of TID detection and characterization. For an effective TID detection, we propose a new approach of detrending TEC based on natural neighborhood interpolation and leave-one-out technique. The properties of detected TIDs induced by each event are analyzed and compared their similarity and dissimilarity between the TIDs from two events.

2. TID DETECTION

GNSS-derived TEC is one of most commonly used observables for monitoring the ionospheric responses to natural hazards. The TEC measurements involve the wave signatures as the consequence of many geophysical activities which are more pertinent to small scale TEC variations on top of the background ionospheric condition. The small scale TEC variation related to a particular geophysical event can be considered as the TIDs generated by acoustic, acoustic-gravity, or gravity waves depending on a source event. By analyzing the small-scale TEC variations, the TIDs can be detected through a proper signal processing. One of most critical steps is the removal of dominant trend of TEC time series. A trend of TEC measurements on a continuous data span are determined in different ways, which is challenging due to the dynamic nature of the ionosphere. In this study, we applied an optimal detrending method to each line-of-sight TEC measurement for detecting TID waveforms from the remaining TEC, which is referred to as residual TEC. The TID waveforms are then determined by examining the fluctuations of residual TEC.

2.1 TEC Extraction from GNSS Observations

TEC is precisely measured from the dual-frequency carrier phase GPS GNSS observations by forming a geometry-free (GF) (or ionospheric) linear combination of GNSS dual-frequency carrier phase observables (Klobuchar 1985, Hoffman-Wellenhof et al. 2008). The carrier phase GF combination is given by

\(\Phi_{i,GF}^k=\Phi_{i,1}^k-\Phi_{i,2}^k=\left({1-f_1^2 \over f_2^2 }\right) I_{i,1}^k+\beta+\varepsilon_{i,GF}^k\)                                                       (1)

where \(\Phi_{i,GF}^k\) is a carrier phase GF observation transmitted from a satellite \(k\) to a receiver \(i, f_1\) and \(f_2\) are the frequency values of dual frequency GNSS signal, \(I_{i,1}^k\) is the ionospheric delay of the first frequency, e.g., GPS L1, \(\beta\) is a constant bias throughout a single continuous data span consists of the ambiguities of dual frequency signals, and inter-frequency bias between dual frequency signals. The noise term, \(\varepsilon_{i,GF}^k\), is computed by the law of error propagation from the measurement noises of two carrier phase observations (Park et al. 2014). \(I_{i,1}^k\) can be obtained from Eq. (1) and the TEC is computed from \(I_{i,1}^k\) where \(I_{i,1}^k=-40.3\cdot TEC f_1^{-2}\) with the TEC in TEC units, which is \(10^{16}\) electrons per m². In this case, the TEC is the total number of electron content along the slant line of sight that is referred to as slant TEC (sTEC) that is shown in Eq. (2).

\(sTEC_i^k={1 \over 40.3}  {f_1^2 f_2^2 \over (f_1^2-f_2^2 ) } \Phi_{i,GF}^k-\beta' - \varepsilon '\)                                                                     (2)

\(\beta'\) and \(\varepsilon '\) are the bias and error in (1) scaled by \({1 \over 40.3}  {f_1^2 f_2^2 \over (f_1^2-f_2^2 ) }\). Many studies (e.g., Heki 2006, 2011, Komjathy et al. 2012) derive vertical TEC from the sTEC measurement after properly determining the bias term. This process may increase uncertainty due to the error from the bias estimation or the single layer assumption of the vertical projection. Because the event derived TIDs are associated with the TEC residual rather than the absolute TEC itself, and the relative variation within a certain time is more critical than absolute magnitude of TEC, studies by Galvan et al. (2011), Park et al. (2011, 2014), Yang et al. (2012) analyzed the sTEC for detecting small scale fluctuation of TEC, which is a TID candidate.

In this study, we focus on the sTEC without correcting the bias term which stays as constant over the time period of observations that is referred to as the relative slant TEC (rsTEC). Eq. (3) is the rsTEC at a particular epoch, t. Note that the bias term is not time dependent.

\(rsTEC(t)= sTEC_i^k (t)+\beta'={1 \over 40.3}  {f_1^2 f_2^2 \over (f_1^2-f_2^2 ) } \Phi_{i,GF}^k (t)+\varepsilon '' (t)\)                                         (3)

2.2 Detrending TEC

The event-driven TID waveforms can be isolated by extracting the high order variation from continuous TEC observations recorded at a GNSS receiver. The high order variation on rsTEC can be separated by subtracting the primary trend of TEC from the TEC measurement that is the residual rsTEC. Because isolating notable peaks of residual rsTECs, which are the TID candidates, is complicated due to the dynamic background ionosphere and numerous unknown causes of TEC disturbances, various approaches have been introduced.

TID candidates can be isolated from the rsTEC by computing numerical horizontal derivatives of the signal for removing trends (Park et al. 2011). The first order numerical derivative can eliminate the dominant trend of TECs, as well as the biases existing in the signal. In this approach, the numerical derivatives are obtained by using three sequential data points as follows:

\(\delta s_i= s_i-{(s_{i-1}+s_{i+1} ) \over 2},~~~i=2,…,N\)                                                                     (4)

where \(\delta s_i\) is the first order numerical derivative, \(s_i\) is the i-th data point, and \(N\) is the total number of observations. The higher order of derivatives can be obtained by taking numerical derivative from the previous order. Park et al. (2011) empirically identified that the third order numerical derivative is sufficient to amplify the potential TIDs and no significant changes appear on the TEC derivatives after that by increasing the order. Fig. 1 illustrates an example of rsTEC measurements (top) and the third order TEC derivatives (bottom) showing potential TID peaks.

Fig. 1. The measured rsSTECs from dual-frequency GNSS carrier phase observables (a), and its third order numerical derivative (b). The rsTECs range between [0.36,7.13] TECU and their numerical derivatives vary between [-0.09,0.08] TECU. Supposing the time of event as the source (solid green line), the traveling time is considered as the time before and after the event specifying by negative and positive spots, respectively. According to the numerical derivative signatures, a TID was revealed about 7.5 minutes after the event.

Although the numerical derivative method was confirmed as an effective approach to detect the MSTIDs by several case studies (Park et al. 2011, 2014), this approach can be affected by a data sampling rate. The high sampling rate data may consider the high frequency waves, such as TIDs, as a regular wave and eliminate it as if it is a trend. The low sampling rate data may miss the high frequency TID signatures. To avoid this weakness, this paper investigated another approach to eliminate the trend of rsTEC by adopting an interpolation approach.

Kotulak et al. (2017) compared three types of interpolation methods for computing regional ionospheric maps based on scattered vertical TEC (vTEC) measurements. The authors utilized these five interpolation methods – inverse distance weighting (IDW), polynomial interpolation, and three Voronoi diagram methods that includes natural neighbor interpolation, Non-Sibsonian interpolation, and Quasi-natural interpolation. For a cross-validation of the interpolated value, the leave-one-out technique (Cawley 2006) was applied. In this study, we adopted the natural neighborhood interpolation as a detrending strategy for computing the residual rsTEC.

The natural neighbor algorithm was first introduced by Sibson (1981) to interpolate the value of a query point as the weighted sum of its natural neighbors formed based on the Voronoi diagrams. In three-dimensional space, R³, the natural neighbor interpolation can be equated as follows (Parsania & Virparia 2016):

\(f(q)=\sum\limits_{i=1}^N w_i (q) f(p_i)\)                                                                                 (5)

where \(f(q)\) is the interpolated value at the query point \(q∈R^3\), \(N\) is the number of neighboring points, \(f(p_i)\) is the value of i-th natural neighbor \(p_i∈R^3\), and \(w_i (q)∈[0,1]\) are the weight coefficients. The weight coefficients are dependent on the volume of the created Voronoi cells based on the considered set of data points. Inserting a new point into an existing Voronoi diagram, a new cell will be created which includes the overlapping segments of the initial cells. The proportion of the overlapping volume between the new and the initial cells to the volume of the new cell defines the weight coefficient as (Parsania & Virparia 2016):

\(w_i (q)={Vol(V_(p_i )∩V_q^+ ) \over Vol(V_q^+ )}\)                                                                                 (6)

where \(V_{p_i}\) denotes the Voronoi cell corresponding to i-th natural neighbor point before inserting the query point into the diagram, and \(V_q^+\) is the Voronoi cell related to the query point q. \(V_{p_i }∩V_q^+\) indicates the overlapping segment between \(V_{p_i}\) and \(V_q^+\).

In this study, we adopted this method to detect TIDs that may be a small amplitude fluctuation on the STEC measurement of the continuous data span between one satellite and one receiver on the ground. The interpolation is utilized to predict the TEC value at a certain time by using the TEC values at neighboring data points. The number of neighboring data points differs based on the sampling rate of data and the background ionospheric condition. For the experiments, we experimentally determined the data points. By excluding one point at a time and interpolating the expected TEC value for the time, the trend can be determined. The residual TEC is derived by subtracting the trend from the TEC measurements. In a quiet ionospheric condition, less fluctuations on TECs are observed. If i-th observation is taken out from the vector of observations, the remained data is still capable of estimating a reliable value for the missing observation via the natural neighbor interpolation. Conversely, the estimated value deviates from its corresponding observation if an aberrant change appears in one or some of the natural neighbors. This deviation deteriorates when unaffected natural neighbors are applied to create the interpolant for estimation of a disturbed observation. Although the deviation of estimated TEC from its true value is not a desirable condition in detrending procedure, it can provide the apt measure to isolate the potential TIDs by defining the residual TEC as the difference between the true TEC and its estimated trend. It should be noted that two different methods in Figs. 1 and 2 detected same peaks considerably TIDs that shows the proposed method successfully retrieved the TID detectable from a confirmed method. In addition, the residual TEC values in Figs. 1 and 2 involve many small fluctuations coming from measurement noise and other effects that may increase the uncertainty to confirm as TIDs. Therefore, we applied the Butterworth band pass filter to mitigate these effects.

Fig. 2. The estimated dominant trend of rsTEC signal used in Fig. 1 by applying cross-validated natural neighbor interpolation (a), and its detrended signature (b). The trend varies between [0.36,7.12] TECU and its residuals range between [-0.03,0.03] TECU.

2.3 TID Detection and Analysis

As Fig. 2 presents, the detrending process properly eliminates the trend and the majority of remaining TEC residuals is random noise of the GNSS observations. We consider the peaks greater than 3σ of the TEC residual data span as the potential TID wave packet. The properties of these potential TIDs are examined in terms of the time of arrival and the 3D coordinates of the waves on the ionospheric thin shell of 350 km altitude above the ground. When the properties of peaks are the medium scale TIDs for the event, we identify them as the TID for the source event. In case of Fig. 2, a TID appeared about 7.5 minutes after the event that can be identified based on the known event time and the time of the peak appeared.

Detected TIDs are further analyzed by observing the morphology of wave packet, their duration and the propagation velocity. The morphology is determined based on the coherent analysis for a pair of TID waves. The correlation coefficient between two TID waves represents the level of similarity between the TID waves that can be utilized for discriminating the TID waves from different source events. The other properties, the duration and the propagation velocity of TIDs also indicate the type of the trigger events.

3. GEOPHYSICAL EVENTS ON 23 JANUARY 2018

On 23 January 2018, i.e., the day of year (DOY) 23 in 2018, a volcanic eruption and several earthquakes with magnitudes of 5.3-7.9 Mw occurred on the ring of fire. To ensure the ionospheric disturbances coming from a certain external source rather than a global scale event, we explored the background ionospheric condition. On 23 January, the \(k_P\) index ranged between 0 and 2 which inferred the probable ionospheric disturbances on DOY 23 predominately originated from any sources other than the geomagnetic activities. In addition to the level of geomagnetic activities, the level of solar activities was examined. The NOAA weekly report of SWPC PRF 2213 specified the level of solar activity during 22-28 January in very low degrees, where the solar radiation storms in NOAA Space Weather Scales was S1 which is “minor”, demonstrating that the ionospheric variations on DOY 23 could not be excited due to any solar activities.

Verifying the quiet ionospheric conditions on 23 January 2018, the 7.9 Mw Gulf of Alaska earthquake and Kusatsu-Shirane volcanic eruption were chosen among the geophysical hazards on this day not only to monitor their impacts on the ionosphere, but also to distinguish between the induced signatures with respect to the type of hazard. Fig. 3 shows the epicenter of the Gulf of Alaska earthquake and Kusatsu-Shirane volcanic eruption, as well as the other geophysical events on DOY 23 associating with the ring of fire.

3.1 Kusatsu-Shirane Volcanic Eruption

The Mount Kusatsu-Shirane is an active stratovolcano in Japan and consists of a series of pyroclastic cones and three crater lakes. The summit of this mount is located on 36.618°N in latitude, 138.528°E in longitude, and 2165 m in elevation. On 23 January 2018, the Japan Meteorological Agency (JMA) reported an unexpected eruptive activity for Mount Kusatsu-Shirane, beginning at 01:59 UTC, and issued an alert level of 3 on a scale of 1-5 to prohibit the residents from approaching the volcano. The eruption was coincident with the onset of volcanic tremors, an avalanche, ejection of tephra, and numerous volcanic earthquakes.

The Mount Kusatsu-Shirane is located nearby to four permanent sites from International GNSS Service (IGS) network within the distance of 150 – 900 km from the volcano summit, which are shown by yellow diamonds in Fig. 3b. As described in Section 2, the dual-frequency carrier phases were converted to relative STEC in Eq (3) for each IGS site. The dominant trend of the rsTEC was estimated based on the combination of the natural neighbor interpolation and the leave-one-out cross validation technique.

Fig. 3. (a) Geophysical hazards on DOY023 in 2018 including earthquakes (circles) and volcanic eruptions (triangles) on the ring of fire (shaded red boundary); (b) Kusatsu-Shirane volcanic eruption; (c) 7.9 Mw Gulf of Alaska earthquake.

After the dominant trend of rsTEC estimated, the residual rsTEC were derived by subtracting the trend from the rsSTEC using a proposed method in the previous section. It should be noted that we applied one-hour of window for determining the interpolating point. Fig. 4 shows three examples of rsTEC (top), the residual rsTEC (middle), and the filtered residual rsTEC (bottom). Since the short-term ionospheric signatures of the volcanic eruption were of the major interest, the signals are presented within a 70 minutes time window, from 10 minutes before to 60 minutes after the eruption. For the residual rsTEC signals, the natural neighbor interpolation successfully estimated the dominant trend of TEC with standard deviation of 0.075 and an average of -0.001 TECU. The investigation in the residual rsTEC signals revealed the excitation of TIDs by the eruptive activity of Mount Kusatsu-Shirane during a period of 7.5 – 32.5 minutes after the eruption which were considered as the immediate response of the ionosphere to the explosion.

Fig. 4. The observed rsTEC (top), residual rsTEC derived by the detrending method of cross-validated natural neighbor interpolation (middle), and smoothed residual rsTEC filtered by Butterworth bandpass filter; the black dashed line in all plots indicates the time of the Kusatsu-Shirane eruption.

Fig. 5. Locations of the volcanic eruption (red star), GNSS ground stations (blue triangles), and TIDs at a 350 km altitude thin shell (pink circles).

Fig. 5 depicts the locations of the volcanic eruption, the TID locations based on the ionospheric pierce points (IPP) and the ground stations.

To estimate the propagation velocity of the detected TIDs, their 3D location should be determined based on the IPP on a single layer model (SLM) (Hofmann-Wellenhof et al. 2008). The detailed mathematical description for TID locations can be found at Park et al. (2011). The 3D coordinates of IPP are employed in computation of the slant distance to the source of disturbances, and eventually, calculation of the propagation velocity. Based on the multiple IPP locations and the time of arrival of TIDs, the experiments revealed that Kusatsu-Shirane eruption excited TIDs with an average velocity of 530 m s-1. TID-looking waveforms may or may not be induced by the source event as they could come from other sources. To confirm the similarity among the detected waveforms, we tested the coherency by computing the correlation coefficient between two waveforms. The correlation coefficient of 0.8 between the detected TIDs supported the correspondence between the ionospheric anomalies and the volcanic eruption.

3.2 Gulf of Alaska Earthquake

On 23 January 2018 at 09:31 UTC, an immense earthquake with magnitude of 7.9 Mw rocked Alaskan with an epicenter about 280 km away from the southeast of Kodiak Island and 560 km away from the southwest of Anchorage. The hypocenter of the earthquake was located in Gulf of Alaska on 56.046°N in latitude, 149.073°W in longitude, and 14 km in depth. The earthquake was due to a strike-slip faulting in which the quake is parallel to the compressional force and its energy is released by the rock displacement in a horizontal direction. Since the strike-slip faulting does not generate any notable uplift or subsidence, the Gulf of Alaska earthquake did not significantly change the water-level in Pacific Ocean, and consequently, only a small tsunami was observed several hours after the earthquake with less than 30 cm depths.

To investigate the ionospheric disturbance triggered by the Gulf of Alaska earthquake, the rsTEC were derived from the GNSS measurements at four Continuously Operating Reference Stations (CORS) within the radius of 1000 km from the epicenter. The distribution of the CORS relative to the epicenter of Gulf of Alaska earthquake is shown by yellow diamonds in Fig. 3c. The rsTEC signals were detrended by the natural neighbor interpolation and the leave-one-out. As stated in the previous section for the volcanic event case study, we applied one-hour of window for determining the interpolating point. The signal trends were estimated with a standard deviation of 0.087 and an average of -0.003 TECU. Fig. 6 presents the detrended rsTECs for different pairs of satellite-CORS. The average amplitude of the perturbations on the residual rsTEC was about 0.3 TECU which is comparable with the reported amplitude by Astafyeva et al. (2013) for shallow earthquakes with magnitude higher than 7. In Fig. 6, a few additional peaks beside the most dominant peak appeared with a smaller intensity. These disturbances can be contributed by aftershocks of this earthquake. According to USGS Earthquake Hazards Program, 95 aftershocks with magnitude greater than 4.0 occurred in the first 48 hours after the primary shock of Alaskan earthquake. The first aftershock occurred at 09:47 UTC, about 16 minutes after the primary shock, with magnitude of 5.0 Mw which could originate the smaller TEC perturbations.

Fig. 6. The observed rsTEC (top), residual rsTEC derived by the detrending method of cross-validated natural neighbor interpolation (middle), and smoothed residual STEC filtered by Butterworth bandpassed filter; the black dashed line in all plots indicates the time of the Gulf of Alaska earthquake.

Fig. 7. Locations of the earthquake (red star), GNSS ground stations (blue triangles), and TIDs at a 350 km altitude thin shell (pink circles).

Fig. 7 depicts the locations of the volcanic eruption, the TID locations based on the ionospheric pierce points (IPP) and the ground stations.

The propagation velocity of TID also was computed using the IPP on the SLM and the average velocity was about 724 m/s. To compare the property of the detected TIDs with other case studies, we examined the duration of the TID waves. We defined the wave packet when the rsTEC residual peaks are greater than 3σ for the 1 hour window of continuous data span. This test was applied before applying the band pass filter which may blend the TID wave and random noise. Based on this approach, the computed duration of the wave packet was about 3 minutes that is slightly shorter than the duration of 4-8 minutes for shallow earthquakes with magnitude between 7.2 –7.8 reported by Astafyeva et al. (2013). However, the mechanism of horizontal rock displacement in strike-slip faulting can explain the lower duration of triggered TIDs by Gulf of Alaska earthquake.

Krabbenhoeft et al. (2018) reported that the Gulf of Alaska earthquake initially propagated from the epicenter to the north, and then its direction changed to the eastward. Spotting the detected TIDs on the ionospheric thin shell at altitude of 400 km, a meaningful relation between the rupture propagation and TID propagation can be identified. The distribution of the detected TIDs indicates two different propagation directions, directed in north and northeast of the epicenter which is consistent with the rupture propagation. This consistency not only validates the correctness of detected TIDs, but also confirms the impact of seismic waves propagation on the location of TIDs.

4. DISCRIMINATION OF TID INDUCED BY TWO EVENTS

TID waveforms can be analyzed to discriminate the response of ionosphere between the volcanic eruption and earthquake events. We applied the 3-minute of time window to define the TID signatures as the duration of earthquake-triggered TID waves was revealed as 3 minutes in the experiment. The correlation coefficients were calculated for sets of residual rsTEC waveforms before and after bandpass filtering.

For the case of the Kusatsu-Shirane volcanic eruption, we computed the correlation coefficient (CC) between the TID signatures and the average of the CC was 0.797. To exclude the noise and other effects, we applied a bandpass filter and the CC was even increased to be 0.824. In the case of the Gulf of Alaska earthquake, the detected TID signatures had a correlation of 0.705 and 0.635 for the waveforms before and after the bandpass filter.

To obtain a better insight on the differences between the volcanic eruption and earthquake TIDs, the same coherence analysis across two events was performed. The earthquake TIDs were correlated to the eruption TIDs with a coefficient of 0.378 and 0.302 for the noisy and denoised residual rsTEC waveforms, respectively. The result indicates that the waveforms of TIDs from a single natural hazard, such as earthquake or volcanic explosion, involve similar characteristics among them showing high correlation up to 0.8. In contrast, characteristics across different types of events are distinctive where the highest correlation coefficient was only 0.4.

While both approaches result the stronger correlation between the TIDs from one event and the lower correlation between the TIDs generated across two events. It is seen that the natural neighbor method showed the higher correlation between TIDs from one event. Although the denoising process using a band-pass filter is useful for isolating the TID waveform, the coherent analysis result shows it possibly loses the characteristic as the denoised signals in both methods lead increased correlation between two events.

In addition, we also investigated the frequency of the representative TID waves induced by the earthquake (top) and the volcanic eruption (bottom) as shown in Fig. 8. The left panel presents the TID waves after filtering. The right panel shows the spectral of representative TIDs. The frequency of earthquake-driven TID is ranging from 2 mHz to 8 mHz where the frequency of the volcanic eruption-driven TID is ranging from 4 mHz to 9 mHz.

Fig. 8. Frequency response of residual STECs after band-pass filtering by Butherworth method following the Gulf of Alaska earthquake (top) and Kusatsu-Shirane eruption (bottom).

5. SUMMARY AND DISCUSSION

This study explored the ionospheric impacts from the Gulf of Alaska earthquake and Kusatsu-Shirane volcanic eruption by analyzing the triggered TIDs from both events. As the events occurred on the same day, the background ionospheric behavior in global scale was comparable for both events which makes the comparison of different types of events valid. The GNSS-derived TEC was processed to detect the TID signatures in the ionosphere for two events. The TIDs were isolated by cross-validating the natural neighbor interpolation and the leave-one-out technique for estimation of STEC trends. The detrended STECs, i.e. the residual rsTEC, were post-processed by applying a bandpass filter for mitigating observational noise. Investigating the denoised residual rsTEC, the short-term ionospheric perturbations after the hazards were isolated. The verification of the proposed detrending method in isolation of TIDs was made based on the numerical derivatives of rsSTEC. The experiments revealed that the TIDs propagated through the ionosphere with an average velocity of 529.894 m/s after the Kusatsu-Shirane volcanic eruption, while the averaged velocity of TIDs was estimated 723.655 m/s after Gulf of Alaska earthquake. A coherence analysis showed that the TID waveforms of the volcanic eruption were highly correlated up to 82.4%. The correlation of TID waveforms of the earthquake was also high with a coefficient of 0.705. However, there was not an extreme correlation between the TID waveforms of earthquake and explosion which was about 37.8% of similarity. This low correlation between the TIDs of different sources may have caused by disparate mechanism of TID generation which affects the spectral characteristics of the waveforms. This discrimination can be further improved by investigation of the mechanism wave propagation through the ionosphere after each type of geophysical event.

AUTHOR CONTRIBUTIONS

Ms. Shahbazi contributed to developing methodology, collecting and processing GNSS data, conducting experiments, and writing. Dr. Park contributed to bringing this research topic, suggesting analysis/validation methods, and writing.

CONFLICTS OF INTEREST

The authors declare no conflict of interest.