Pre-Seismic Signature Detection using Diurnal GPS-TEC and Kriging Interpolation Maps (ASK-VTEC Technique):
11 May 2011, M9.0 Tohoku Earthquake Case Study

Thammaboribal P.,¹* Tripathi, N. K. ¹ and Lipiloet, S.²

¹Asian Institute of Technology, School of Engineering and Technology, P.O. Box 4, Klong Luang, Pathumthani 12120, Thailand

prapasgnss@gmail.com*

²Department of Civil Engineering, Faculty of Engineering, Rajamangala University of Technology Thanyaburi, Pathumthani 12120, Thailand

*Corresponding Author

Abstract

During earthquake preparation, a seismogenic electric field is generated on the Earth's surface, which then penetrates the atmosphere, reaching several hundred kilometers above the lithosphere and causing perturbations in the ionospheric electron density. These seismo-ionospheric anomalies are typically identified using statistical techniques, such as running averages or inter-quartile range (IQR), with the anomalous zones analyzed from Global Ionospheric Maps (GIMs). This study used ASK-VTEC technique for detecting pre-seismic signatures, based on the daily mean of vertical total electron content (AVTEC) and its standard deviation (SVTEC). The spatial distribution of vertical electron density (VTEC) is mapped using Ordinary Kriging (OrK) interpolation instead of GIMs. To validate and assess the performance of this approach, the M9.0 Tohoku earthquake, which occurred on March 11, 2011, was chosen as a case study. The results show that ionospheric anomalies were observed on May 8, 2011, three days before the earthquake. A significant anomalous zone was detected southwest of the epicenter between 06:00 and 10:00 UTC. These findings are consistent with previous studies, and further, the VTEC spatial distribution maps generated by the OrK interpolation technique outperformed GIMs on a local scale. Additionally, the analytical process is less complex compared to conventional methods. Therefore, the pre-seismic signature detection technique proposed in this study offers a viable alternative for investigating VTEC anomalies prior to earthquake events.

Keywords: ASK-VTEC, Earthquake Precursor, Ionospheric Disturbance, SVTEC, Tohoku Earthquake

1. Introduction

Seismo-ionospheric anomalies prior to earthquake occurrence were first detected in 1964 after the Mw 9.2 Alaska earthquake [1]. The study of ionospheric perturbations before earthquakes has been intensively conducted for many earthquake events. The results show that seismo-ionospheric anomalies typically emerge 0 to 10 days before an earthquake [2][3] and [4]. However, some scientists argue that there is no significant statistical relationship between the lithosphere, atmosphere, and ionosphere (LAI) during the earthquake preparation process [3]. As a result, the topic of ionospheric disturbances before earthquakes remains controversial and is still debated among seismologists. The physical mechanism of LAI coupling is not yet fully understood and continues to be discussed by many researchers. However, there are three hypotheses explaining this mechanism [4] and [5]. The first hypothesis is the chemical channel, based on the emanation of radon gas. The second hypothesis involves acoustic gravity waves (AGWs) generated during the earthquake preparation process, which may propagate upward to the ionosphere and cause irregularities in electron density [6]. The third hypothesis is the electro magnetic channel, which involves positive charge carriers [7]. The electrostatic electric field generated by positively charged ions at the Earth's surface penetrates through the atmosphere to the ionosphere, causing perturbations in electron density [8]. This theory is supported by a laboratory experiment that demonstrated that when rock samples are subjected to high stress, electrical charges are produced, flowing from the stressed part of the rock to the unstressed part. This process ionizes the air at the rock surface, resulting in changes in the conductivity of the air around the rock [9].

On a global scale, when rocks are stressed during the earthquake preparation process, positive charge carriers (holes) are released, traveling to the unstressed part of the rock and ionizing the near-ground atmosphere. This leads to changes in air conductivity due to the production of positive ions, which then penetrate through the atmosphere to the ionosphere, causing variations in electron density [10]. However, some laboratory experiments conducted by other researchers do not support the P-hole theory, as they argue that the Earth's crust is fluid-saturated, making it unlikely for electrical charge buildup to occur during the earthquake preparation process [11].

Seismo-ionospheric coupling is one of the candidates for investigating earthquake precursory signals. Increasingly, studies on this topic have been conducted using both ground-based stations and satellite data to confirm such mechanisms. The most widely used technique for detecting seismo-ionospheric anomalies prior to earthquakes is a statistical-based analytical technique, such as sliding median combined with interquartile range (IQR), which is used to construct upper and lower boundaries for detecting ionospheric anomalies [8][9][10][11][12] and [13]. Other techniques, such as artificial neural networks (ANN) [14][15] and [16], Particle Swarm Optimization (PSO) [17], radon anomalies [18], and very low frequency (VLF) radio waves [19], are also utilized. The most commonly used detection technique is the statistical-based method, where the threshold is set to median ± 2IQR. In cases where the observed VTEC lies very close to the threshold, it becomes questionable whether the value can be considered an anomaly.

The spatial distribution of total electron content (TEC) on a global scale is illustrated by Global Ionospheric Maps (GIMs), which are modeled using mathematical functions, such as spherical harmonic (SH) expansion. The spatial resolution of GIMs is 2.5° × 5° in latitude and longitude, respectively [12], and the temporal resolution is 2 hours. Since 1° at the equator corresponds to approximately 110 km and 95 km in latitude and longitude, respectively [13], the spatial resolution of GIMs at the equator is about 275 km in latitude and 475 km in longitude. For this reason, local ionospheric distribution cannot be accurately reproduced by GIMs [20], especially within the earthquake preparation zone. The distribution of vertical electron content (VTEC) can also be estimated from VTEC values at a set of reference GNSS stations using a "grid-based" technique, such as spatial interpolation. Many studies indicate that the VTEC instantaneous maps generated by the grid-based technique are roughly equivalent to, or even better than, the mathematical function-based technique [21] and [22], and they also show high accuracy when compared to the final product provided by the CODE analysis center [23].

Thus, in this study, the spatial distribution of VTEC maps was created by applying the "ordinary Kriging" interpolation technique to illustrate the variations in AVTEC and SVTEC values over the study area, instead of using GIMs. The performance of the interpolated maps was evaluated by collecting data from the Crustal Movement Observation Network of China (CMONOC) and the International GNSS Service (IGS) from DOY60 to DOY90 of 2015. Ordinary Kriging (OrK), Universal Kriging (UrK), inverse distance weighting (IDW) [24], and planar fit interpolation techniques were used to model regional VTEC maps in China. The results showed that the performance of UrK and OrK were similar, and both interpolation techniques produced better results than IGS GIMs when applied for ionospheric correction in GPS positioning. Moreover, OrK is particularly suitable for creating surface maps on a local scale [18].

To validate the performance of seismo-ionospheric anomaly detection using the ASK-VTEC technique [25] and ordinary Kriging interpolation maps, the M9.0 earthquake, known as the "Tohoku-Oki earthquake," was selected as a case study. This earthquake was chosen because studies on seismo-ionospheric anomaly detection related to this event have been widely reported. The submarine earthquake struck the east coast of Honshu, Japan, on March 11, 2011, at 05:46:23 UT. The epicenter was located at 38.927°N and 142.373°E in latitude and longitude, respectively. The focal depth of the earthquake was 29.0 km, which is considered a shallow earthquake. Previous studies show that anomalies related to this earthquake were observed as early as May 8, 2011, 0 to 4 days prior to the earthquake [26] and [27]. Therefore, the consistency of the anomalous days detected by the alternative techniques proposed in this study can be verified to determine whether the results correspond with those of previous studies.

2. Data and Methodology

The accuracy of GNSS measurements can be improved by eliminating or minimizing sources of error during the measurement process [28]v[29] and [30]. One of the largest sources of error is the electron density in the ionosphere, where the speed of the signals is delayed by particle charges. To minimize this error, it is necessary to measure the electron density in the ionosphere. A dual-frequency GNSS receiver receives L1 and L2 signals transmitted from navigation satellites in space.

The electron density in the ionosphere, known as the "total electron content" (TEC), along the signal path (STEC), is determined by the difference in the travel times of both signals from the satellite to the receiver, which is referred to as the "ionospheric time delay," as shown in Equation 1 [31].

Equation 1

Where f₁ and f₂ are the frequencies of L₁ and L₂ signal which equal to 1575.42 MHz and 1227.60 MHz, respectively, P₁ and P₂ are the pseudo range measurements of L₁ and L₂ signal, respectively. TEC value derived from Equation (1) is the electron density along the ray path between the satellite and receiver so that it is called "slant TEC (STEC)". STEC can be converted to vertical total electron content (VTEC) at the pierce point of the ray path by Equation 2 [32]:

Equation 2

Where z' is zenith angle of satellite which is expressed in Equation 3:

Equation 3

Where R_E is the radius of the earth which approximately equal to 6375 km; α is the GPS satellite elevation angle; and h is height of iono- spheric layer based on thin-shell theory, which is 350 km in this study. The unit of TEC is expressed as TECU (1 TECU = 10¹⁶ electrons/m²). The measured pseudo ranges contain bias values both in satellites and receiver so called "differential code bias (DCB)". The receiver contains larger DCB than that of the satellites, thus, the TEC estimation is affected by receiver DCB than satellite DCB [33]. The DCB values of reference stations can be obtained from ftp://ftp.aiub.unibe.ch/CODE, provided by the university of Bern, thus, equation 2 can be rewritten as in Equation 4 [34].

Equation 4

Where b_r and b^s are DCB values of receiver and satellite, respectively.

The earthquake preparation zone (EPZ) is determined as per Equation 5 [35] :

R = 10^0.43M

Equation 5

Where R is the radius of earthquake preparation area [km] and Mis the moment magnitude of earthquake.

2.1 Data Collection Data Processing and

Methodology

The moment magnitude of the earthquake is M9.0, according to Equation 5 the radius of EPZ is equivalent to 7,413 km. The boundary of EPZ was created by applying buffer operation as illustrates in Figure 1.

Figure 1 Earthquake preparation zone and selected IGS for data collection. Red solid star represents the earthquake epicenter, green solid circle represents IGS reference stations and red circle represents EPZ

Table 1: Selected IGS reference stations for data collection

The data were collected from the International GNSS service (IGS) reference stations that locate within the earthquake preparation zone as shown in Figure 1. In this study, the data from 29 IGS reference stations were collected, the list of the collected data illustrate in Table 1. The Receiver Independent Exchange Observation (RINEX.o) files for 10 days before and 5 days after the event were collected from the stations listed in Table 1. Previous studies have shown that anomalies were observed on March 8, 2011, which is 3 days prior to the earthquake event. Therefore, the performance of the ionospheric anomaly detection method proposed in this study can be validated. The data can be accessed at ftp://cddis.gsfc.nasa.gov/ gnss/data/daily. The daily data were recorded from 00:00:00 to 23:59:30 UTC with a sampling rate of 30 seconds. The data were processed using the open-source software "GPS-TEC," developed by Gopi Seemala, and only the GPS-TEC values were analyzed. The outputs from the software include both Slant Total Electron Content (STEC) and Vertical Total Electron Content (VTEC); however, only the VTEC values were analyzed. The daily mean of VTEC (AVTEC) and the standard deviation of VTEC (SVTEC) at each reference station were determined using the following equations:

The daily mean of VTEC (AVTEC) as well as standard deviation of VTEC (SVTEC) at each reference station were determined from Equations 6 and 7 [25].

Equation 6

Equation 7

Where AVTEC is the average of VTEC on day-to-day basis, SVTEC is the standard deviation of AVTEC, N is the total number of measured VTEC from 00:00:00 to 23:59:30.

The ionospheric anomalies were detected from the variation plots of AVTEC and SVTEC. Both values were calculated on a day-to-day basis to ensure that they were not influenced by the TEC values of previous days. This anomaly detection technique is relatively straightforward and less complicated to implement compared to conventional techniques, such as the sliding median approach.

The ionospheric anomalies were further investigated by analyzing the spatial distribution of AVTEC and SVTEC, which were created using the ordinary Kriging interpolation technique. This technique is advantageous because, unlike deterministic interpolation methods that rely solely on values within the search radius, it also accounts for statistical relationships, such as spatial autocorrelation, in surface generation. In this study, the variation of VTEC over the study area, typically presented by GIMs, was replaced by AVTEC and SVTEC maps. To investigate the day-to-day variations in VTEC prior to the Tohoku earthquake, interpolated maps were created for the 10 days and 5 days before and after the event. Although a running window of 27 days is typically required to account for periodic changes due to the short-period oscillation of the ionosphere, previous studies on this earthquake case showed that anomalies were observed 0-4 days prior to the event. Since the main objective of this study is to propose an alternative VTEC anomaly detection approach based on diurnal VTEC, the 10-day and 5-day periods before and after the event were sufficient to validate whether the results are consistent with previous studies that used conventional detection techniques.

Although the temporal resolution of GIMs is 2 hours [26], the temporal resolution of AVTEC and SVTEC maps was set to 1 day as a rough anomaly scanner. This is because anomalous days can be detected by significantly high variations in AVTEC/SVTEC compared to previous days. With this assumption, the variation of VTEC is not influenced or disturbed by VTEC variations from prior days. The duration of TEC anomalies during the earthquake preparation period is typically 4-6 hours on anomalous days [26] and [36].

Table 2: Sources of solar activity index and geomagnetic storm indices

Based on statistical theory, the distribution characteristics of VTEC values on anomalous days would change due to disturbances in the dataset caused by the seismo-ionospheric mechanism. Therefore, the diurnal VTEC (AVTEC) and its standard deviation (SVTEC) on anomalous days are significantly different from those on other days. According to this statistical assumption, it is possible to investigate VTEC anomalies using AVTEC and SVTEC. Finally, the AVTEC/SVTEC variation maps, with a temporal resolution of 2 hours (similar to GIMs), were analyzed on anomalous days to investigate whether an anomalous zone appears within the EPZ.

2.2 Influence of Geomagnetic Storms and Solar Activity to Ionospheric Perturbation

Ionospheric TEC disturbance is mostly associated with space weathers, so the influence of the space weather on detected anomalous day must be investigated in order to determine whether the detected anomaly is attributed to space weathers or earthquake preparation process. The magnitude of geomagnetic activity is measured by two indicators as disturbance storm time index (Dst) and geomagnetic index (Kp). Solar activity magnitude is presented by solar radio flux at 10.7 cm index (F10.7). These three indices can be achieved from the links given in Table 2. It is difficult to determine whether the anomaly is induced by the influence of an earthquake or space weather [37] and [38]. The three indices listed in Table 2 were used to distinguish whether the anomaly was caused by a seismo-ionospheric event, space activities, or a geomagnetic storm. If the solar-geomagnetic index exceeds the following limits: |Dst| index < 30 nT, Kp index > 3, and F10.7 > 150.

2.3 Ordinary Kriging Interpolation (ORK)

The variation of VTEC was illustrates by AVTEC and SVTEC maps, the maps were created by applying geostatistical interpolation technique so called "Ordinary Kriging (OrK)", this technique is successfully applied to obtain ionospheric information in SBAS, WAAS or GAGAN [39]. The assigned weights in this interpolation technique are not only relied on the distances between unknown points and measured points but spatial autocorrelation is also taken account in the weight determination, so that it gives the best linear unbiased prediction values [40][41] and [42]. The prediction values Z(x₀) at any locations within the study were determined from the known values Z(x_i) , and the weights that derived from variogram as summarized by Equation 8:

Equation 8

Where Z(x₀) is the prediction value at x₀location, Z(x_i) is the known values that obtained at x_ilocation which is AVTEC and SVTEC at x _i location in this study, and is the weight function that depends on the semivariogram, the distance between prediction and known point, and the spatial relationship among the known values around the prediction point. The determination of the weights is based on minimizing the variance of the variable as per Equation 9.

Var(Z(x + h) – z(x)) = 2γ(|h|)

Equation 9

Where var is the variance function, x and x+h are two adjacent locations and γ represents the variogram. The expectation of Z(x+h) and Z(x) is defined in Equation 10.

E(Z(x + h) – z(x)) = 0

Equation 10

Where E is the expectation operator.

In order to satisfy the unbiasness in the prediction, the difference between the predictions and true values should be zero. To ensure that the predictors are unbias for unknown measurement, the sum of the weight equal to 1 is imposed as Equation 11.

Equation 11

The solution to the minimization, constrained by unbiasness is done by applying Lagrange method, and the kriging equation is expressed as a matrix in Equations 12 and 13.

Equation 12

Equation 13

3. Results

The outputs derived from data processing using GPS-TEC software contain slant TEC (STEC) as well as VTEC values. In this study, only VTEC is analyzed and discussed. The AVTEC and SVTEC at each reference station on daily basis were determined.

3.1 Average Diurnal Vertical Total Electron Content (AVTEC)

AVTEC of 10 days before and 5 days after the earthquake events were plotted for investigating the variation of diurnal VTEC. The variation of AVTEC of 29 IGS stations illustrates in Figure 2. The figure shows the variation of AVTEC at each IGS reference station from DOY060 to DOY075. It can be seen that high variations of AVTEC were observed on DOY062, DOY064, and DOY067, which are 8 days, 6 days, and 3 days prior to the earthquake occurrence, respectively. Periodic oscillations were observed at the TCMS and TWTF stations, both located southwest of the epicenter. Although the expected normal level of AVTEC at each station is not known, and it is difficult to define what can be considered an anomaly, Figure 2 displays that anomalous AVTEC values were observed at both stations. After DOY067, it is clearly seen that the AVTEC at both stations decreased, and high variations were no longer observed. The average AVTEC from DOY060 to DOY066 at TCMS and TWTF were 18.98 TECU and 18.63 TECU, respectively, while the maximum AVTEC values observed on DOY067 were 31.36 TECU and 31.38 TECU, respectively. The peak of AVTEC at some stations within the EPZ was also observed on DOY062, DOY064, and DOY067, but the latter day showed a significant high variation of AVTEC. The AVTEC variations in the northern part of the study area showed a flat behavior, and the AVTEC values were relatively low compared to the southern part.

The AVTEC at PIMO also showed a high value compared to the other stations. The AVTEC at PIMO gradually increased from 20.48 TECU on DOY060 and reached its peak of 33.02 TECU on DOY071. The variation of AVTEC at GUUG, GUAM, and CNMR, which are located at the same latitude as PIMO, exhibited the same behavior as at PIMO. However, periodic oscillations were not observed as they were at TCMS and TWTF. To investigate the spatial distribution of AVTEC over the study area, time series AVTEC maps from DOY060 to DOY075 (1–16 May 2011) were created using the ordinary Kriging interpolation technique. The time series map of AVTEC is illustrated in Figure 3. The spatial distribution of AVTEC over the study area, illustrated in Figure 3, clearly shows that the AVTEC at lower latitudes contains higher AVTEC values than those at higher latitudes. This is because the electron density along these latitudes is influenced by the northern crest of the equatorial fountain, which occurs only during the daytime. AVTEC anomalies were significantly observed on DOY064 and DOY067 in the southwest of the epicenter.

Figure 2: TEC anomaly detection from diurnal AVTEC.
Red star represents the day of earthquake event

Figure 3: Time-series map of AVTEC from May 1, 2011 to May 16, 2011 (DOY060 to DOY075).
Red solid star represents the earthquake epicenter

However, it is difficult to determine whether the AVTEC anomalies on these days could be considered pre-seismic signatures of the Tohoku earthquake, as the anomalous zone remained in the southwest of the epicenter even after the earthquake occurred. The interpolated maps reveal that the distribution of AVTEC in the study area was mainly influenced by the AVTEC values at the IGS stations located at lower latitudes.

3.2 Standard Deviation of Diurnal Vertical Total Electron Content (SVTEC)

According to the previous studies, on the anomalous day, the VTEC anomalies occurred between 12:00LT and 18:00LT with the duration of about 2 hours, so that the electron content was likely to be disturbed by the exceeding particle charges generated by earthquake preparation process and caused the variation in data distribution. The assumption of the VTEC anomaly detection derives from the daily standard deviation of VTEC (SVTEC) is that an abnormal peak of SVTEC will be observed on the anomalous day because the standard deviation of VTEC on anomalous day should be higher than the ordinary days. The time series distribution plot of SVTEC 10 days before and 5 days after the earthquake illustrates in Figure 4.

Figure 4 illustrates that the characteristics of SVTEC oscillation are similar to those of AVTEC, with abnormal peaks observed on DOY062, DOY064, and DOY067, which were 8 days, 6 days, and 3 days prior to the earthquake occurrence, respectively. The maximum SVTEC values appeared at TCMS and TWTF on DOY067. Since the distance between these two IGS reference stations is approximately 24.76 km, the difference in SVTEC values between the stations was a maximum of 0.57 TECU. The SVTEC at HKWS and HKSL was also higher compared to the other stations on DOY067, but no significant variation in SVTEC at these stations was observed. The average SVTEC from DOY060 to DOY066 at TCMS and TWTF was 15.17 TECU and 14.98 TECU, respectively, while the maximum SVTEC values observed on DOY067 were 24.31 TECU and 24.15 TECU, respectively.

Figure 4: TEC anomaly detection from diurnal SVTEC.
Red star represents the day of earthquake event

Figure 5: Time-series map of SVTEC from May 1, 2011 to May 16, 2011 (DOY060 to DOY075)
Red solid star represents the earthquake epicenter

Thus, the high variations in SVTEC on DOY067 are likely to be considered as SVTEC anomalies prior to the Tohoku earthquake. To investigate the spatial variations of SVTEC over the study area, a statistical interpolation technique called "Ordinary Kriging" was used to generate SVTEC distribution maps. The maps were generated on a day-to-day basis to identify the anomalous days. The time series of SVTEC spatial distribution from DOY060 to DOY075 (May 1, 2011, to May 16, 2011) is illustrated in Figure 5. The spatial distribution of SVTEC in Figure 5 shows that the characteristics of SVTEC variation were similar to those of AVTEC, with SVTEC concentration inversely proportional to latitude.

Figure 6: Spatial distribution of regional VTEC from 02:00UT to 24:00 UT on March 8, 2008
Red solid star represents the earthquake epicenter

However, the abnormal SVTEC zone was more significant than the AVTEC zone. According to Figure 5, the SVTEC anomalous zone was first observed on DOY064 in the southwest of the epicenter, the same region as the AVTEC anomaly. The anomalous zone then clearly reappeared in the southeast direction of the epicenter on DOY066 and reached its peak value on DOY067. Thus, it is suspected that the SVTEC anomalies on DOY066 and DOY067 were most likely associated with the Tohoku earthquake, as the anomalous zone was located very close to the earthquake epicenter.

To assess the performance of the AVTEC and SVTEC distribution maps modeled by the Ordinary Kriging Interpolation (OrK) technique, the VTEC on DOY067 (the anomalous day) was created with a temporal resolution of two hours, as GIMs. This allows the consistency of the results between the interpolated maps and GIMs to be evaluated. The output maps are illustrated in Figure 6. Basically, the distribution of VTEC is presented by global ionospheric maps (GIMs) generated from IGS-GPS data with spherical harmonic expansion and is utilized for identifying VTEC anomalous zones. However, in this study, GIMs were replaced by OrK-interpolated maps with the same temporal resolution of two hours as the GIMs. Figure 6 shows the variation of VTEC on the anomalous day (May 8, 2011) from 02:00 UT to 24:00 UT. It is clearly seen that VTEC started increasing from 06:00 UT, reached its maximum level from 08:00 UT to 10:00 UT, and gradually decreased afterward. The anomalous zone was located in the southwest of the epicenter, in the same region as the AVTEC and SVTEC anomalous zones on DOY067. In terms of time and space, the results are consistent with previous studies [27] and [43]. Thus, it is reasonable to infer that the VTEC distribution can be represented by OrK-interpolated maps because they provide equivalent or even better output maps compared to GIMs.

3.3 Influence of Solar Activity and Geomagnetic Storm

The electron density in the ionosphere is normally affected by solar activities and geomagnetic storms. Therefore, the levels of both phenomena on the anomalous days were monitored to investigate whether the VTEC anomaly is attributed to such phenomena or to the earthquake preparation process. The magnitude of geomagnetic activity is indicated by two indices: the disturbance storm time index (Dst) and the geomagnetic index (Kp). When the ionosphere is disturbed by a geomagnetic storm, the values of |Dst| < 30 nT and Kp > 3 are observed. The magnitude of solar activity is indicated by the solar radio flux at the 10.7 cm index (F10.7). If the index is greater than 150 SFU, it indicates that the ionosphere is disturbed by solar activities.

Figure 7: Variations of geomagnetic storm indices ( Dst and Kp) and solar activity index (F10.7) from DOY060 to DOY 075 (May 1, 2011 to May 16, 2011)
red sold star represents Tohoku earthquake event day,
and red dashed line represents geomagnetic storm and solar activity thresholds