Capability of TEC correlation Analysis and Deceleration at Propagation Velocities of Medium-Scale Traveling Ionospheric Disturbances: Preseismic Anomalies before the Large Earthquakes

Ken Umeno^*, Ryo Nakabayashi, Takuya Iwata, Minghui Kao

Department of Applied Mathematics and Physics, Graduate School of Informatics, Kyoto University, Kyoto, Japan.

DOI: 10.4236/ojer.2021.104008                

Abstract

Data analysis method (CRA, hereafter) to correlate multiple TEC anomaly signals has detected pre-seismic anomalies before the 2011 Tohoku-Oki earthquake (Iwata & Umeno 2016), the 2016 Kumamoto earthquake (Iwata & Umeno 2017) and the 2016 Tainan earthquake (Goto et al. 2019). However, a critical argument said that those anomalies detected by CRA would not be pre-seismic anomalies published by Journal of Geophysical Research-Space Physics (126), 2021 (JGR-SP (126), hereafter). In this paper, we would point out its incorrect use of statistical anomalies in evaluating CRA as the following points: CRA is shown to increase the signal-to-noise ratio (SNR) to amplify pre-seismic TEC’s small anomaly signals with synchronizing and correlating multiple GNSS receivers’ data. We proved again that pre-seismic anomalies certainly exist before the 2011 Tohoku-Oki earthquake and the 2016 Kumamoto earthquake with additional data analysis. In particular, as a temporal anomaly, deceleration at propagation velocities of medium-scale traveling ionospheric disturbances (MSTID, hereafter) before the 2016 Kumamoto earthquake captured by CRA (Iwata & Umeno 2017) is elucidated as pre-seismic anomalies. Furthermore, we proposed a physical model to predict that 35 m/s change at MSTID propagation velocities estimated by TEC’s CRA requires 0.58 × 10^-3 V/m electric field in the F Layer ionosphere. Contrary to the claim with the incorrect use of statistical anomalies in JGR-SP (126), TEC’s correlation anomalies detected by CRA (Iwata & Umeno 2016 and Iwata & Umeno 2017) clearly provided supporting evidence that physical pre-seismic anomalies really exist.

Keywords

CRA, TEC, MSTID, Preseismic Ionospheric Anomaly, The 2016 Kumamoto Earthquake, The 2011 Tohoku-Oki Earthquake

1. Introduction

CorRelation Analysis (CRA, hereafter) is a general method to extract signal from complicated noise in diverse kinds of signal processing. It can be distant to merge radio signals of Quasars to lock and unlock digital communication as an encryption tool, or is near to extract Wi-Fi signal from noise of home appliances around people’s daily living. CRA to detect total electron content (TEC) anomalies before large earthquakes is based on the very long baseline interferometry’s concept and spreading spectrum communications technology. It has been implemented to report in the references: 2016, Iwata & Umeno (hereafter I & U16) [1], 2017, Iwata & Umeno (hereafter I & U17) [2], and 2019, Goto, et al., (hereafter Goto et al. 2019) [3].

The existence of pre-seismic TEC anomalies before large earthquakes has been debated until today (Heki 2011 [4], Kamogawa & Kakinami 2013 [5], Heki & Enomoto 2013 [6], Masci et al. 2015 [7], Kelley et al. 2017 [8], Muafiry & Heki 2020 [9], Eisenbeis & Occipinti 2021 [10] ). Such debate for a decade is caused by lacking of conclusive physical models to explain pre-seismic TEC anomalies. Under such circumstances, CRA is sequentially targeted to extract pre-seismic anomalies from the 2011 Tohoku Oki earthquake (Mw 9.0, depth 24 km), the 2016 Kumamoto earthquake (Mw 7.3, depth 12 km) and the 2016 Tainan earthquake (Mw 6.4, depth 14.6 km), respectively. Due to the energy preservation law, every earthquake must have a preparatory physical process. Assume one factor of its mechanism makes electromagnetic process be affected by Earth’s magnetic field penetrating to the ionosphere, the signature can be detected with increasing the signal-to-noise ratio in principle.

Recently, Ikuta, et al. [11] examined the results of I & U16 [1] and I & U17 [2] by the statistical analysis and posed a question on CRA capability at detecting a pre-seismic anomaly. One of the purposes of the present paper is to add conclusive evidence to support that TEC correlation anomalies detected in I & U16 and I & U17 are really pre-seismic anomalies with proposing a new physical model of such pre-seismic TEC anomalies. Furthermore, we show that their analysis method of CRA by [11] has fatal errors for judgement on CRA capability. Another purpose of the present paper is to identify a new physical phenomenon “deceleration at the propagation velocities of MSTID” as peculiar characteristics of pre-seismic TEC anomaly phenomena, which is the first time observed by CRA (I & U17).

The general characteristics of CRA are introduced and why [11] has fatal errors are explained in Section 2. Deceleration at the propagation velocities of MSTID as a new phenomenon of pre-seismic ionospheric anomaly is discussed in Section 3. Three physical models showing this phenomenon are also presented. Model 1 predicts that 35 m/s change in deceleration at the propagation velocities of MSTID detected requires 0.58 mV/m electric field in the F-Layer of ionosphere. In Section 4, supportive data analysis of CRA to identify the deceleration at the propagation velocities of MSTID before the 2016 Kumamoto earthquake as a pre-seismic phenomenon is given. Discussion about the physical model validity for deceleration at the propagation velocities of MSTID is then presented in connection with CRA data analysis in Section 5. Finally, conclusion is given in Section 6.

2 General Characteristics of CoRrelation Analysis (CRA)

To detect anomalies from GNSS (Global Navigation Satellite Systems) stations before large earthquakes, CRA computes a correlation among abnormalities in TEC (Total Electron Contents) observed at GNSS stations. The first step of CRA is to choose a GNSS station to correlate [1]. Once we choose a central station, $M\left(\ge 1\right)$ surrounding stations, which are the nearest to the central station, can be selected. One can number the central station and each surrounding station from 0 to M, where the number 0 means the central station and the numbers 1 to M are allocated to the surrounding stations. Let $X_{i\mathrm{,}t}$ be abnormalities of the station i at time t such as prediction errors computed from sample data at the station i. Let $t_s$ be the time length of sample data for learning to predict which were set to 2.0 hours in the CRA in [1] [2] [3].

The crux of the CoRrelation Analysis (CRA) [1] is to compute a correlation given by

$C\left(T\right)=\frac{1}{NM}{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\displaystyle \underset{j=0}{\overset{N-1}{\sum }}}{X}_{0,t+t_S+j\Delta t}\cdot X_{i,t+t_S+j\Delta t}$ (1)

$T=t+t_S+t_{test},$

where $N\left(>1\right)$ is the number of data in a Test Data during the time $t+t_S$ to $t+t_S+t_{test}$ , $\Delta t$ is a sampling interval in the Test Data (usually 30 seconds for TEC data), $t_S$ is the time length of the Sample Data (Learning period) and $t_{test}$ is the time length of the Test Data (Prediction Period). I & U16 and I & U17 set up that $t_S=2.0\left[\text{hours}\right]$ and $t_{test}=0.25\left[\text{hours}\right]$ . The correlation value C(T) can be rewritten as:

$C\left(T\right)=\frac{1}{N}{\displaystyle \underset{j=0}{\overset{N-1}{\sum }}}{X}_{0,t+t_S+j\Delta t}\cdot \left(\frac{1}{M}{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{X}_{i,t+t_S+j\Delta t}\right)=\frac{1}{N}{\displaystyle \underset{j=0}{\overset{N-1}{\sum }}}{X}_{0,t+t_S+j\Delta t}\cdot {\tilde{X}}_{0,t+t_S+j\Delta t},$ (2)

where

${\tilde{X}}_{0,t+t_S+j\Delta t}=\frac{1}{M}{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{X}_{i,t+t_S+j\Delta t}.$

Note that if $M=1$ , C(T) becomes just a normal correlation between $X_0$ and $X_1$ :

$C\left(T\right)=\frac{1}{N}{\displaystyle \underset{j=0}{\overset{N-1}{\sum }}}{X}_{0,t+t_S+j\Delta t}\cdot X_{1,t+t_S+j\Delta t}.$

Thus, from Equation (2), one can see that C(T) can capture a synchronized temporal anomaly patterns correlated between $X_0$ (a value at the central station) and ${\tilde{X}}_0$ (a mean value of the values $X_i$ observed at the surrounding stations). If anomaly patterns of observational points are coherently periodic such as medium-scale traveling disturbances (MSTIDs), C(T) also shows periodic patterns with the same period. On the contrary, if anomaly patterns are coherently non-periodic irregular patterns, C(T) also shows a certain irregular pattern. Thus, not only its value C(T), but also a temporal characteristics of C(T) are vitally important to elucidate anomaly alert.

If N is large, the following relation

${\displaystyle \underset{j=0}{\overset{N-1}{\sum }}}{X}_{0,t+t_S+j\Delta t}\cdot {\tilde{X}}_{0,t+t_S+j\Delta t}=O(\; N\; )$

holds for non-correlated noisy signals $X_0$ and $X_i$ from the central limit theorem (CLT). Thus,

$C\left(T\right)=\frac{1}{N}{\displaystyle \underset{j=0}{\overset{N-1}{\sum }}}{X}_{0,t+t_S+j\Delta t}\cdot {\tilde{X}}_{0,t+t_S+j\Delta t}=O\left(\frac{1}{\sqrt{N}}\right)\to 0\text{for}N\to \infty .$ (3)

On the contrary, for some coherent synchronized signals $X_0$ and $X_i$ due to some anomaly phenomena, it is evident that

${\displaystyle \underset{j=0}{\overset{N-1}{\sum }}}{X}_{0,t+t_S+j\Delta t}\cdot {\tilde{X}}_{0,t+t_S+j\Delta t}=O\left(N\right).$

Thus we can expect a higher C(T) such that

$\left|C\left(T\right)\right|=\left|\frac{1}{N}{\displaystyle \underset{j=0}{\overset{N-1}{\sum }}}{X}_{0,t+t_S+j\Delta t}\cdot {\tilde{X}}_{0,t+t_S+j\Delta t}\right|=O\left(1\right)>0\text{for}N\to \infty ,$ (4)

which clearly distinguishes a signal from noisy signals when N is sufficiently large. An SNR or signal-to-noise ratio at this abnormality detector C(T) can be measured by the ratio between the variances of signal and noise; thus the following general relation holds:

$\text{SNR}=\frac{{\left(O\left(1\right)\right)}^2}{{\left(O\left(\frac{1}{\sqrt{N}}\right)\right)}^2}=O\left(N\right).$

Thus, N is a key parameter of CRA to measure temporal correlations with each temporal abnormalities, where N is regarded as the spreading factor in spread spectrum technology [12] [13].

A pre-seismic ionospheric anomaly, if it exists, should be distinguished from other space weather phenomena such as MSTID and high geomagnetic activity. For the issue on distinction between ionospheric anomaly and MSTID by CRA, the reference [11] argues that the 65 - 168 m/s MSTID propagation velocity range in I & U17 as TEC pre-seismic anomaly is not abnormally low, as compared to the statistics on the propagation velocities reported in the past. They also concluded that TEC anomaly detected for the 2016 Kumamoto earthquake day reported in I & U17 was not a pre-seismic one. Referring [1] [2] [3], this kind of anomaly on I & U17 is Not a statistical anomaly.

We would define it single event anomaly (focus on one time, one location, one satellite and one sub-ionospheric track) of C(T) in CRA as follows. Suppose anomalies occurred around specific area as earthquake’s characteristic. A single event anomaly means an anomaly which is dependent on a specific set of event characterized by time, location, satellite and sub-ionospheric point tracks of the satellite. In Fig. 2 of Goto et al. 19 [3], the TEC correlation anomalies observed by different satellites (GPS17 and GPS28) before the 2016 Tainan earthquake, are sharply different. C(T) of GPS 17 shows a strong anomaly (C(T) higher than 100) while C(T) of GPS28 shows no anomaly (C(T) less than 10) even at the same GNSS stations (gais-wanc of Taiwan) at the same time period (17:00-20:00 UTC) on the earthquake day (Feb. 5., 2016). This means whether TEC pre-seismic correlation anomalies can be detected or not are dependent on the SIP (sub-ionospheric) track of the satellite. In this respect, the reference [11] is incorrect because they calculated based on the statistics of C(T) with all kinds of GNSS satellites used in their analysis. Such kind of statistical analysis confused their judgement and rough use of CRA in [11] to argue the existence of pre-seismic anomaly is inconclusive. TEC correlation anomalies detection depends on a choice of satellite [3].

The reference [11] also performed CRA analysis towards the 2011 Tohoku-Oki earthquake on March 11, 2011 and the foreshock on March 9, 2011. They also re-examined I & U16 based on the statistical anomaly distribution of C(T). They reproduced CRA’s high correlation value on March 11 of I & U16 and further argued that the correlation values C(T) were not so abnormally high compared to the statistic of high C(T) values such that $C\left(T\right)\ge 25$ . (Fig. 2 of [11] ). Again, their logic of the argument in [11] is based on the simple criteria of statistical anomaly values of all over Japan. Earth’s geomagnetic field strength on Tohoku area (higher latitude) is higher than Kumamoto area, ionospheric anomalies computed by C(T) of Kumamoto (lower latitude) tend to be higher than Tohoku area (higher latitude). It is natural to provide different C(T) thresholds at different latitude so the difference of the abnormality criteria of C(T) naturally exist between the case of 2011 Tohoku-Oki earthquake and the case of the 2016 Kumamoto earthquake.

Considering such inconsistency, a threshold of C(T) can be computed by the mean and the variance of its preceding non-earthquake days such as 12 days [14]. The reference [11] also argued that the high values of C(T) of the 2011 Tohoku-Oki earthquake may be attributed to the large Kp index and thus the anomaly detected (high C(T)) before the Tohoku-Oki earthquake) by CRA in I & U16 may be due to high-geomagnetic activity (Kp = 5). We object the argument of [11] by giving a counter example on the non-earthquake days with low C(T) value and large Kp index. Such days with low C(T) and large Kp (Kp = 5) can be illustrated as March 1, 2011 and April 8, 2016, both of which are the non-earthquake days (See Table 1). The days with (Kp = 5) have no abnormality in C(T) as compared to the earthquake days (March 11, 2011 and April 15, 2016). In the data analysis for computing a mean value and the standard deviation (sd) of C(T), the 12 consecutive days before the target date were used for each day. In Table 1, data with low elevation angle (one hour from the beginning and one hour to the end of TEC data observed) were discarded for CRA to avoid high C(T) values due to the artificial anomaly by the low elevation angle. Contrary to the claim of [11], a signature of large Kp index has no apparent relation with high C(T) of CRA which can detect synchronous anomaly with multiple GNSS stations while the high C(T) on 2011/03/11 and 2016/04/15 seems to be related to the two large earthquakes (the 2011 Tohoku-Oki earthquake and the 2016 Kumamoto earthquake), respectively. The high C(T) could be considered as ionospheric pre-seismic anomalies as presented in I & U16 and [14]. Through the above descriptions, we argue that the reference [11] has several fatal errors with incorrect use of statistical anomalies for evaluating CRA.

3. Deceleration of Propagation Velocities of MSTID and the Physical Mechanism

In this section, a general relation between the deceleration at propagation velocities in MSTID found in I & U17 and a predicted change of electric field strength in the ionosphere is derived to provide a physical basis to cause the anomaly patterns detected by CRA. Physical behavior of MSTID can be understood in terms of plasma physics (physics for ionized gases) [15]. Its physical origin of MSTID is known to be related to the Perkins instability of the ionosphere [16].

The equations of motion for electrons of mass $m_e$ and ions of mass $m_i$ in the ionosphere are given by

$\begin{array}{l}{n}_e{m}_e\left(\frac{\partial v_e}{\partial t}+\left(v_e\cdot \nabla \right)v_e\right)\\ =n_e{m}_{e}g-e{n}_e\left(E+v_e\times B\right)-\nabla p_e-n_e{m}_e{\nu }_{en}\left(v_e-v_n\right)+{\displaystyle \underset{i}{\sum }}{R}_{ie}\end{array}$ (5)

Table 1. Maximum values of C(T) of some days with Large Kp index in 2011 and 2016.

$\begin{array}{l}{n}_i{m}_i\left(\frac{\partial v_i}{\partial t}+\left(v_i\cdot \nabla \right)v_i\right)\\ =n_i{m}_{i}g+e{Z}_i{n}_i\left(E+v_i\times B\right)-\nabla p_i-n_i{m}_i{\nu }_{in}\left(v_i-v_n\right)-R_{ie}-{\displaystyle \underset{j\ne i}{\sum }}{R}_{ij}\mathrm{,}\end{array}$ (6)

where $v_e\left(v_i\right)$ is the velocity of an electron (an ion i), $Z_i$ is the ion charge number (multiples of e) of ion i, $n_e\left(n_i\right)$ is the number density of electrons (ions i), ${\nu }_{en}\left({\nu }_{in}\right)$ is the frequency of collisions between an electron (an ion i) and neutral particles, $\nabla p_e\left(\nabla p_i\right)$ is the gradient of pressure acting on electrons (ions i), $R_{ie}$ is the force per unit volume affected by collisions between electrons and ions i, $R_{ij}$ is the force per unit volume affected by collisions between ions i and another kind of ions j and $g$ the gravity force affected by the earth is a vector per unit mass per unit volume. After summing Equation (5) and ${\displaystyle {\sum }_i}$ Equation (6), with ${\displaystyle {\sum }_i}{R}_{ie}$ (in Equation (5)) + ${\displaystyle {\sum }_i}-R_{ie}$ (in Equation (6)) = 0 and ${\displaystyle {\sum }_i}{\displaystyle {\sum }_{j\ne i}}{R}_{ij}=0$ , one can derive the plasma equation:

$\begin{array}{l}\rho \frac{\partial v}{\partial t}+n_e{m}_e\left(v_e\cdot \nabla \right)v_e+{\displaystyle \underset{i}{\sum }}{n}_i{m}_i\left(v_i\cdot \nabla \right)v_i\\ =\rho g+j\times B-\nabla p-{\displaystyle \underset{i}{\sum }}{n}_i{m}_i{\nu }_{in}\left(v_i-v_n\right)\mathrm{,}\end{array}$

where $\rho \left(\equiv n_e{m}_e+{\displaystyle {\sum }_i}{n}_i{m}_i\right)$ is the mass density, $v\left(\equiv \left(n_e{m}_e{v}_e+{\displaystyle {\sum }_i}{n}_i{m}_i{v}_i\right)/\rho \right)$ is the center of mass velocity, $j\left(\equiv -e\left(n_e{v}_e-{\displaystyle {\sum }_i}{Z}_i{n}_i{v}_i\right)\right)$ is the current density, and $p\left(\equiv p_e+p_i\right)$ is the pressure (plasma pressure). Here, by electrical neutrality of plasma, ${\displaystyle {\sum }_i}{n}_i{Z}_i=n_e$ and we neglect the term $-n_e{m}_e{\nu }_{en}\left(v_e-v_n\right)$ because the

electron cyclotron frequency ${\Omega }_e=\frac{eB}{m_e}$ is much greater than the collision frequency ${\nu }_{en}$ and $m_e\ll m_i$ . The layer of the ionosphere for considering MSTID

is the F-Layer with the 300 km height above the ground. One can assume that ions in that layer are of one type, O⁺ for simplicity. Thus, $\rho \simeq n_i{m}_i$ and $v_i\simeq v$ hold because $m_e\ll m_i$ . Furthermore,

$n_e{m}_e\left(v_e\cdot \nabla \right)v_e+{\displaystyle \underset{i}{\sum }}{n}_i{m}_i\left(v_i\cdot \nabla \right)v_i\simeq \rho \frac{Dv}{Dt}\mathrm{,}$

where $\frac{Dv}{Dt}\equiv \frac{\partial v}{\partial t}+\left(v\cdot \nabla \right)v$ . Accordingly, the final form of the equation motion for ionized gas (ionosphere) above the 300 km is:

$\rho \frac{Dv}{Dt}=\rho g+j\times B-\nabla p-n_i{m}_i{\nu }_{in}\left(v_i-v_n\right)\mathrm{.}$ (7)

An ionospheric current $j_{\perp }$ perpendicular to Earth’s magnetic field line $B$ penetrating the ionosphere is given by

$j_{\perp }={\sigma }_P\left(E_{\perp }+v_n\times B\right)+{\sigma }_H\frac{B}{B}\times \left(E_{\perp }+v_n\times B\right)\mathrm{,}$

where $E_{\perp }$ is an electric field vector perpendicular to $B$ , $v_n$ is the mean velocity vector of a gas of neutral particles, ${\sigma }_P$ is the Pedersen conductivity

computed by ${\sigma }_P=\frac{n_{i}e}{B}\left(\frac{{\nu }_{in}{\Omega }_i}{{\nu }_{in}^2+{\Omega }_i^2}+\frac{{\nu }_{en}{\Omega }_e}{{\nu }_{en}^2+{\Omega }_e^2}\right)$ and ${\sigma }_H$ is the Hall current conductivity computed by ${\sigma }_H=\frac{n_{i}e}{B}\left(\frac{{\Omega }_i^2}{{\nu }_{in}^2+{\Omega }_i^2}-\frac{{\Omega }_e^2}{{\nu }_{en}^2+{\Omega }_e^2}\right)$ [17]. In the F-Layer

ionosphere 300 km over the earth, ${\sigma }_P\gg {\sigma }_H$ . Thus one can assume that $j_{\perp }={\sigma }_P\left(E_{\perp }+v_n\times B\right)$ . The obtained equation of motion for a velocity $v_{\perp }$ perpendicular to the geomagnetic field $B$ is:

$\begin{array}{l}\frac{D{v}_{\perp }}{Dt}=g_{\perp }+\frac{e}{m_{i}B}\left(\frac{{\nu }_{in}{\Omega }_i}{{\nu }_{in}^2+{\Omega }_i^2}+\frac{{\nu }_{en}{\Omega }_e}{{\nu }_{en}^2+{\Omega }_e^2}\right)\left(E_{\perp }+v_n\times B\right)\times B\\ -\frac{{\left(\nabla p\right)}_{\perp }}{n_i{m}_i}-{\nu }_{in}\left(v_{i\perp }-v_{n\perp }\right)\mathrm{.}\end{array}$

Propagation $v_{\perp }$ of MSTID is essentially a macroscopically stationary drift motion of an ionized gas (not electrons). Thus, the propagation velocity of MSTID satisfies the continuity equation for an incompressible fluid:

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho v\right)=0\text{with}\nabla \cdot v=0\Rightarrow \frac{Dv}{Dt}=\frac{D{v}_{\perp }}{Dt}=0.$

Therefore, we get the propagation velocity of MSTID which is the perpendicular to $B$ by the following formula:

$v_{\perp }=v_{n\perp }+\frac{g_{\perp }}{{\nu }_{in}}+\frac{e}{m_{i}B}\left(\frac{{\Omega }_i}{{\nu }_{in}^2+{\Omega }_i^2}+\frac{{\nu }_{en}{\Omega }_e}{{\nu }_{in}\left({\nu }_{en}^2+{\Omega }_e^2\right)}\right)\left(E_{\perp }+v_n\times B\right)\times B-\frac{{\left(\nabla p\right)}_{\perp }}{{\nu }_{in}{n}_i{m}_i}\mathrm{.}$ (8)

Suppose an electric field $E_{\perp }$ is changed as $E_{\perp }\to E_{\perp }-\Delta E_{\perp }$ . Such a change in $E_{\perp }$ also changes the propagation velocity of MSTID as $v_{\perp }\to v_{\perp }-\Delta v_{\perp }$ where

$\Delta v_{\perp }=\frac{e}{m_i}\left(\frac{{\Omega }_i}{{\nu }_{in}^2+{\Omega }_i^2}+\frac{{\nu }_{en}{\Omega }_e}{{\nu }_{in}\left({\nu }_{en}^2+{\Omega }_e^2\right)}\right)\Delta E_{\perp }\times \frac{B}{B}\mathrm{.}$ (9)

Finally, one can obtain:

$\Delta v_{\perp }=\frac{{\sigma }_{P}B}{n_i{m}_i{\nu }_{in}}\Delta E_{\perp }=\frac{e}{m_i}\left(\frac{{\Omega }_i}{{\nu }_{in}^2+{\Omega }_i^2}+\frac{{\nu }_{en}{\Omega }_e}{{\nu }_{in}\left({\nu }_{en}^2+{\Omega }_e^2\right)}\right)\Delta E_{\perp }\simeq \frac{e}{m_i{\Omega }_i}\Delta E_{\perp }$ (10)

which is the causal relation between the deceleration at propagation velocities of MSTID ( $\Delta v_{\perp }$ ) and $\Delta E_{\perp }$ , a sudden change of electric field in the ionosphere.

Namely, a sudden change in the opposite direction (in the eastward direction at the midnight, see Figure 1.) causes deceleration at MSTID’s propagation velocities (Model 1). Here assume that ${\Omega }_i\simeq 100\text{rad}\cdot {\text{s}}^{-1}$ for quantitative validation of the model. Note that $e=1.602\times {10}^{-19}\text{C}$ , $m_p=1.673\times {10}^{-27}\text{kg}$ , $m_i=16m_p$ . In this case, $\Delta v_{\perp }=35\text{m}\cdot {\text{s}}^{-1}$ change with the deceleration at propagation velocities in MSTID requires $\Delta E_{\perp }=0.58\times {10}^{-3}\text{N}\cdot {\text{C}}^{-1}=0.58\text{mV}/\text{m}$ change in the F Layer of the ionosphere. Thus, even a small change in electric field in the F layer can be measured by a macroscopic data estimation of deceleration at propagation velocities of MSTID. In other words, even a fairly small change in the electric field strength can be measured by amplification effect with the propagation velocity of MSTID by the following formula:

$a\equiv \frac{\Delta v_{\perp }}{\Delta E_{\perp }}=\text{Const}\mathrm{.}=\frac{eZ}{m_i{\Omega }_i}=5.9848\times {10}^4\text{m}\cdot \text{C}\cdot {\text{s}}^{-1}\cdot {\text{N}}^{-1}\simeq 6\times {10}^4{\text{T}}^{-1}\mathrm{,}$

where a is an amplification factor between $\Delta v_{\perp }$ and $\Delta E_{\perp }$ and can be regarded as a constant parameter. Interestingly, 35 m∙s⁻¹ change in the propagation velocities of MSTID requires 0.58 × 10⁻³ N∙C⁻¹ electric field lines at 300km height is almost consistent with Kelley et al. 2017’s estimation [8] such that an $E\times B/B^2$ drift of 12 m∙s⁻¹ for the dislocation of electrons observed with TEC and its 3D-tomography analysis by Heki et al. [4] [9] [18] requires 0.5 × 10⁻³ N∙C⁻¹ electric field lines at the base of ionosphere, although the above two estimation methods are totally different. Responsible components of plasma in MSTID propagation are ions as $\rho \simeq m_i{n}_i$ in the F region while [4] and [8] consider a model of electron dislocations due to an $E\times B/B^2$ drift in the E region.

Other physical models responsible for MSTID’s deceleration at propagation velocities can be attributed to a reduction of Pedersen conductivity ${\sigma }_P$ such as ${\sigma }_P\to {\sigma }_P-\Delta {\sigma }_P$ by

$\Delta v_{\perp }=\frac{\Delta {\sigma }_{P}B}{n_i{m}_i{\nu }_{in}}{E}_{\perp }$ (Model 2) (11)

or an increase in ion density as $n_i\to n_i+\Delta n_i$ with $\Delta n_i>0$ by

$\Delta v_{\perp }=-\frac{{\sigma }_P\Delta n_{i}B}{n_i^2{m}_i{\nu }_{in}}{E}_{\perp }+\frac{{\left(\nabla p\right)}_{\perp }\Delta n_i}{n_i^2{m}_i{\nu }_{in}}\simeq -\frac{{\sigma }_P\Delta n_{i}B}{n_i^2{m}_i{\nu }_{in}}{E}_{\perp }$ (Model 3) (12)

where we discard the term of the gradient of pressure during the time scale of preservation of the MSTID periodic stripe structure. Oyama et al. 2019 [19] observed the reduction of the Pedersen conductivity prior to the 2011 Tohoku-oki earthquake, which is consistent with Model 2 (the former theory on a reduction of Perdersen conductivity of the F region). They observed the enhancement of O⁺ by DMSP satellites prior to the 2011 Tohoku-Oki earthquake, which is also consistent with Model 3 (the latter theory on an increase of ion density) [19]. Figure 1 summarizes the three physical models presented here where MSTID in the nighttime hour of the mid-latitude northern hemisphere is assumed.

It should be pointed out that among the three models presented here, Model 1 is most suitable for explaining the deceleration of propagation velocities in MSTID because the time scale of a sudden change in electric field in the ionosphere in Model 1 is comparable to the time scale (1 - 2 hours) for the observed deceleration of propagation velocities in MSTID by CRA while other models (Model 2 and Model 3) must have longer time scales to cause change in the reduction of the Pedersen conductivity in the ionosphere (Model 2) or the increase of ion density in the ionosphere (Model 3) to compare with Model 1 (the electric field change model).

4. Evidence for Deceleration of Propagation Velocities at MSTID before the 2016 Kumamoto Earthquake

We analyzed GNSS data obtained by GEONET and then converted them to get

Figure 1. Physical models for deceleration at propagation velocities of MSTID. Three physical models of explain deceleration at propagation velocities with MSTID are depicted. What are depicted here is the case of the mid-latitude northern hemisphere in the nighttime. The southward direction of MSTID propagation is crude approximation for simplicity of presentation. In reality, the nighttime propagation of MSTID on the mid-latitude northern hemisphere is southwestward direction while the daytime propagation of MSTID on the mid-latitude northern hemisphere is southward. The three physical models have their different time-scale in changing corresponding physical variables.

Slant TEC data to perform CRA as I & U16 and I & U17. We selected the 15 GNSS stations located in Kyushu Island in Japan as the central stations (See FigureS1.) and set the same parameter as $M=30$ as I & U16 and I & U17. Figure2 and Figure3 show that MSTID deceleration at propagation velocities is clearly seen. On the earthquake day, the half periods $\Delta T_1$ and $\Delta T_2$ of the MSTID one cyclic period became widen as $\Delta T_1<\Delta T_2$ while the MSTID maintains the spatial periodic stripe structure with the wavelength $\Lambda$ . See FigureS2 and FigureS3 for the MSTID spatial structures on the corresponding days.

Thus, the averaged values of Table 2 over 15 stations depicted in FigureS1 are obtained as:

$\overline{\Delta T_1}=0.203\text{hour},\overline{\Delta T_2}=0.302\text{hour},\bar{\gamma }\equiv \overline{\left(\frac{\Delta T_1}{\Delta T_2}\right)}=\mathrm{0.617.}$

The wavelength $\Lambda$ of MSTID around 15:50 (UTC) on April 15, 2016 is estimated as $\Lambda =157400\text{m}$ by CRA with all the GNSS stations in Japan (See FigureS2). Thus, we obtain the propagation velocities of MSTID:

$v\left(\text{before}\right)=\frac{\Lambda }{2\overline{\Delta T_1}}=107.69\text{m}\cdot {\text{s}}^{-1},v\left(\text{after}\right)=\frac{\Lambda }{2\overline{\Delta T_2}}=72.39\text{m}\cdot {\text{s}}^{-1}.$

A deceleration $\Delta v$ at MSTID propagation velocities is finally obtained by

$\Delta v=v\left(\text{before}\right)-v\left(\text{after}\right)=\frac{\Lambda }{2}\left(\frac{1}{\overline{\Delta T_1}}-\frac{1}{\overline{\Delta T_2}}\right)=35.3\text{m}\cdot {\text{s}}^{-1}.$

On the contrary, the data on April 13, 2016, where the usual MSTID was identified by I & U17 shows the opposite sign: no deceleration at propagation velocities in MSTID is observed. As can be seen in TableS1, Figure 4 & Figure 5, the half

periods $\Delta T_1$ and $\Delta T_2$ on April 13, 2016 are almost same as $\gamma =\frac{\Delta T_1}{\Delta T_2}\simeq 1$ .

Figure 2. Correlation values (0087) before the 2016 Kumamoto earthquake The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 0087 (Koga, Fukuoka Prefecture) is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure 3. Correlation values (0089) before the 2016 Kumamoto earthquake on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 0089 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Table 2. Half periods of MSTID on April 15, 2016 estimated by CRA.

Figure 4. Correlation values (0087) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 0087 as a central station and GPS satellite PRN17.

Figure 5. Correlation values (0089) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 0089 as a central station and GPS satellite PRN17.

Through the remarkable difference between the deceleration of MSTID on the earthquake day (April 15, 2016) and non-deceleration of MSTID on April 13, 2016 as is also seen in FigureS2 and FigureS3, one can consider a deceleration at propagation velocities at MSTID is the characteristics of preseismic phenomena because it is extremely difficult to find such a deceleration at MSTID propagation velocity on the usual MSTIDs [20]. The other twenty six figures, Figures S4-S29 also support that a deceleration of propagation velocities of MSTID occurred on the earthquake day while such a deceleration was not observed on the non-earthquake day. We conclude here that the TEC’s correlation analysis presented here shows the deceleration at propagation velocities of MSTID and its physical existence on the deceleration of MSTID propagation velocities before the 2016 Kumamoto earthquake is conclusive.

5. Discussion

We have shown that MSTID, a typical space weather phenomenon is deeply related to pre-seismic ionospheric anomaly identified as a deceleration in propagation velocity of MSTID before large earthquakes. Similar interrelation between space weather phenomena such as geomagnetic storms and preseismic atmospheric and ionospheric anomalies is also reported [21].

For the issue on distinction between ionospheric anomaly and MSTID, the reference [11] argues that the 65 - 168 m/s MSTID propagation velocity range firstly reported in I & U17 as TEC preseismic anomaly is not abnormally low as compared to the statistics on the propagation velocities reported in the past [20] and Ikuta et al. concluded that TEC anomaly detected for the 2016 earthquake day reported in I & U17 is not a preseismic one. We argue that this kind of anomaly reported on I & U17 is Not a statistical anomaly but a single event anomaly (focus on both space and time duration) whose definition is given in Section 2. Here, with additional data analysis with the half periods of MSTID obtained by CRA in the preceding section, we have shown that a deceleration in MSTID propagation velocities before the 2016 Kumamoto earthquake on April 15, 2016 has certainly occurred as pre-seismic anomaly behavior as reported by I & U17 (See Figure 2 of I & U17) and that the reduction of propagation velocities of MSTID as originally reported by I & U17 has been further clarified in comparison with the normal propagation velocity case of MSTID on April 13, 2016 (See TableS1 and Figures S17-S29). To explain this phenomenon, we have provided three physical models (Models 1-3) of abnormal deceleration in MSTID propagation velocities before large earthquakes. Interestingly, our estimation of 0.58 mV/m electric field requirement in the F-Layer ionosphere for 35 m/s deceleration of MSTID propagation velocities is almost consistent with Kelley’s estimation of 0.5 mV/m electric field requirement at the base of ionosphere for dislocations of electrons firstly claimed by Heki [4] [8] [9] [18]. The 10⁴ times amplification effect with a measurement of MSITD propagation velocities elucidated in Section 3 is comparable with the amplification effect of CRA in increasing signal-to-noise-ratio introduced in Section 2. An electric field of 0.58 mV/m of the F-Layer ionosphere is not detectable in practice, which means a high capability potential of ionospheric anomaly detection with TEC’s CRA. There are also another two physical models (Models 2-3) explaining deceleration of MSTID propagation velocities. Models 2-3 (decrease in the Pedersen conductivity and increase in ion densities) are also consistent with DMSP satellite data of direct observations on O⁺ prior to the 2011 Tohoku-Oki earthquake by [19]. With these physical models, we clarify that abnormality as deceleration at MSTID propagation velocities detected on the 2016 Kumamoto earthquake day is a conclusive physical process caused by a sudden change of some physical parameters before the earthquake. Note here that there is still missing link such as LAI (Lithosphere-Atmosphere-Ionosphere) coupling models [22] [23] to fill a gap of the whole preparatory process of earthquake. Concerning the 2011 Tohoku-Oki earthquake, various ionospheric anomaly phenomena have been reported so far [1] [4] [18] [19] [24] [25]. Mizuno & Takashima 2013 [24], and Igarashi et al., 2020 [25] observed some physical anomalies before the earthquake by direct measurement of physical parameters such as current in air and oblique ionograms between Wakkanai and Kokubunji in Japan, respectively. They indicated supportive and significant physical evidence on the existence of certain abnormal pre-seismic phenomena before the 2011 Tohoku-Oki earthquake. How to prove the pre-seismic anomalies related to high C(T)? To provide a physical model causing ionospheric anomaly such as the deceleration in MSTID propagation velocity before the large earthquakes may be one of the answers.

6. Conclusions

We have shown that contrary to the claim by [11] which is shown to have multiple fatal errors in the judgement of capability of detecting TEC pre-seismic anomalies with the incorrect use of statistical anomaly, TEC’s correlation anomalies detected by I & U16 and I & U17 clearly provided supporting evidence that physical pre-seismic anomalies really exist before the two large earthquakes (the 2011 Tohoku-Oki earthquake and the 2016 Kumamoto earthquake), respectively.

Here we declare that the abnormality criteria should be taken by an AND operation of various abnormality sensing detectors such as the low propagation velocity of MSTID and the low anomalous area rates as discussed in I & U17 while [11] considered these abnormality conditions separately. One important finding of the present study is that space weather phenomena such as MSTID are related to pre-seismic anomalies as confirmed by the CRA results. We also deduced a significant amplification effect in the method of CRA for detecting TEC pre-seismic anomalies and give a plausible physical model of observed deceleration in MSTID propagation velocity detected before the 2016 Kumamoto earthquake by CRA. The model validity of proposed physical models (Models 1-3) and the universality of the new phenomena of deceleration in MSTID propagation velocity towards other large earthquakes elucidated as an ionospheric pre-seismic anomaly with the physical models are left as a future study to be investigated.

Acknowledgements

We thank Prof. Akira Mizuno, Mr. Hiroki Tanaka and Mr. Yuta Tsusaka for useful discussions.

Supporting Information

Supporting Information Section are composed of twenty-nine Figures and one Tablethat help readers understand the manuscript better.

FigureS1 shows the locations of the fifteen GNSS stations used for TEC CoRelation Analysis (CRA) on April 13, 2021, and April 15, 2021 (the day of main shock of the 2016 Kumamoto earthquake).

FigureS2 and FigureS3 show the correlation values at all the GNSS stations in Japan before the 2016 Kumamoto earthquake on 2016/04/15 (UTC) and on 2016/0413 (UTC), respectively.

Figures S4-S16 show the correlation values for each different central station before the 2016 Kumamoto earthquake on April 15, 2016 (UTC).

Figures S17-S29 show the correlation values for each different central station on April 13, 2016 (UTC).

TableS1 shows the list of the half periods of MSTID on April 13, 2016, estimated by CRA.

Figure S1. Location of the 15 selected GNSS stations for CRA.

Figure S2. Correlation values at all the GNSS stations in Japan before the 2016 Kumamoto earthquake on 2016/04/15. We used every GNSS station as a central station and mapped the results into the Japan map. The GPS satellite PRN 17 is used here. The black x marks represents the epicenter. The earthquake occurrence time is 16:25 UTC on April 15, 2016 and the time 15:50 in the figure corresponds to 35 minutes before the main shock.

Figure S3. Correlation values at all the GNSS stations in Japan at 16:00 (UTC) on 2016/04/13. No earthquakes occurred on the day while MSTID was observed. We used every GNSS station as a central station and mapped the results onto the Japan map. The GPS satellite PRN17 is used here.

Figure S4. Correlation values before the 2016 Kumamoto earthquake (0451) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 0451 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S5. Correlation values before the 2016 Kumamoto earthquake (0452) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 0452 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S6. Correlation values before the 2016 Kumamoto earthquake (0453) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 0453 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S7. Correlation values before the 2016 Kumamoto earthquake (0685) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the timet (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 0685 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S8. Correlation values before the 2016 Kumamoto earthquake (0687) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 0687 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S9. Correlation values before the 2016 Kumamoto earthquake (0688) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 0688 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S10. Correlation values before the 2016 Kumamoto earthquake (0710) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 0710 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S11. Correlation values before the 2016 Kumamoto earthquake (0711) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 0711 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S12. Correlation values before the 2016 Kumamoto earthquake (1060) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 1060 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S13. Correlation values before the 2016 Kumamoto earthquake (1062) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 1062 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S14. Correlation values before the 2016 Kumamoto earthquake (1063) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 1063 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S15. Correlation values before the 2016 Kumamoto earthquake (1064) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 1064 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S16. Correlation values before the 2016 Kumamoto earthquake (1069) on April 15, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The black line indicates the exact time 16:25 (UTC) when the 2016 Kumamoto earthquake occurred. The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $0<\Delta T_1\equiv t_2-t_1<\Delta T_2\equiv t_3-t_2$ , a deceleration at propagation velocity of MSTID is clarified. The GNSS station 1069 is used as the central station and the GPS satellite PRN17 is selected for the analysis.

Figure S17. Correlation values (0451) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 0451 as a central station and GPS satellite PRN17.

Figure S18. Correlation values (0452) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 0452 as a central station and GPS satellite PRN17.

Figure S19. Correlation values (0453) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the timet (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 0453 as a central station and GPS satellite PRN17.

Figure S20. Correlation values (0685) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 0685 as a central station and GPS satellite PRN17.

Figure S21. Correlation values (0687) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 0687 as a central station and GPS satellite PRN17.

Figure S22. Correlation values (0688) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 0688 as a central station and GPS satellite PRN17.

Figure S23. Correlation values (0710) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 0710 as a central station and GPS satellite PRN17.

Figure S24. Correlation values (0771) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 0771 as a central station and GPS satellite PRN17.

Figure S25. Correlation values (1060) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the timet (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 1060 as a central station and GPS satellite PRN17.

Figure S26. Correlation values (1062) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 1062 as a central station and GPS satellite PRN17.

Figure S27. Correlation values (1063) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 1063 as a central station and GPS satellite PRN17.

Figure S28. Correlation values (1064) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 1064 as a central station and GPS satellite PRN17.

Figure S29. Correlation values (1069) on April 13, 2016. The vertical axis shows the correlation C(T) and the horizontal one the time t (UTC). The blue lines indicate the times $t_1,t_2,t_3,\left(t_1<t_2<t_3\right)$ when C(T) has extremal values. Because $\Delta T_1\equiv t_2-t_1\simeq \Delta T_2\equiv t_3-t_2$ , a deceleration of propagation velocity of MSTID is not detected. We used the pair of the GNSS station 1069 as a central station and GPS satellite PRN17.

Table S1. Half periods of MSTIDs on April 13, 2016 estimated by CRA.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Iwata, T. and Umeno, K. (2016) Correlation Analysis for Preseismic Total Electron Content Anomalies around the 2011 Tohoku-Oki Earthquake. Journal of Geophysical Research: Space Physics, 121, 8969-8984. https://doi.org/10.1002/2016JA023036
[2]	Iwata, T. and Umeno, K. (2017) Preseismic Ionospheric Anomalies Detected before the 2016 Kumamoto Earthquake. Journal of Geophysical Research: Space Physics, 122, 3602-3616. https://doi.org/10.1002/2017JA023921
[3]	Goto, S., Uchida, R., Igarashi, K., Chen, C.-H., Kao, M.H. and Umeno, K. (2019) Preseismic Ionospheric Anomalies Detected before the 2016 Taiwan Earthquake. Journal of Geophysical Research: Space Physics, 124, 9239-9252. https://doi.org/10.1029/2019JA026640
[4]	Heki, K. (2011) Ionospheric Electron Enhancement Preceding the 2011 Tohoku-Oki Earthquake. Geophysical Research Letters, 38, L17312. https://doi.org/10.1029/2011GL047908
[5]	Kamogawa, M. and Kakinami, Y. (2013) Is an Ionospheric Electron Enhancement Preceding the 2011 Tohoku-Oki Earthquake a Precursor? Journal of Geophysical Research: Space Physics, 118, 1751-1754. https://doi.org/10.1002/jgra.50118
[6]	Heki, K. and Enomoto, Y. (2013) Preseismic Ionospheric Electron Enhancements Revisited. Journal of Geophysical Research: Space Physics, 120, 7006-7020. https://doi.org/10.1002/2015JA021353
[7]	Masci, F., Thomas, J.N., Villani, F., Secan, J.A. and Rivera, N. (2015) On the Onset of Ionospheric Precursors 40 min before Strong Earthquakes. Journal of Geophysical Research: Space Physics, 120, 1383-1393. https://doi.org/10.1002/2014JA020822
[8]	Kelley, M.-C., Swartz, W.-E. and Heki, K. (2017) Apparent Ionospheric Total Electron Content Variations Prior to Major Earthquakes Due to Electric Fields Created by Tectonic Stresses. Journal of Geophysical Research: Space Physics, 122, 6689-6695. https://doi.org/10.1002/2016JA023601
[9]	Muafiry, I.-N. and Heki, K. (2020) 3-D Tomography of the Ionospheric Anomalies Immediately before and after the 2011 Tohoku-Oki (Mw 9.0) Earthquake. Journal of Geophysical Research: Space Physics, 125, e2020JA027993. https://doi.org/10.1029/2020JA027993
[10]	Eisenbeis, J. and Occhipinti, G. (2021) The TEC Enhancement before Seismic Events Is an Artifact. Journal of Geophysical Research: Space Physics, 126, e2020JA028733. https://doi.org/10.1029/2020JA028733
[11]	Ikuta, R., Oba, R., Kiguchi, D. and Hisada, T. (2021) Reanalysis of the Ionospheric Total Electron Content Anomalies around the 2011 Tohoku-Oki and 2016 Kumamoto Earthquakes: Lack of a Clear Precursor of Large Earthquakes. Journal of Geophysical Research: Space Physics. https://doi.org/10.1002/essoar.10506732.1
[12]	Umeno, K. (2000) Chaotic Monte Carlo Computation: A Dynamical Effect of Random Number Generations. Japanese Journal of Applied Physics, 39, 1442-1456. https://doi.org/10.1143/JJAP.39.1442
[13]	Umeno, K. (2001) SNR Analysis for Orthogonal Chaotic Spreading Sequences. Nonlinear Analysis, 47, 5753-5763. https://doi.org/10.1016/S0362-546X(01)00714-3
[14]	Umeno, K. and Iwata, T. (2021) Abnormality Detection Apparatus, Communication Apparatus, Abnormality Detection Method, and Recording Medium. United States Patent, No. 11,016,206 (Issued on May 25, 2021).
[15]	Spitzer, L. (1962) Physics of Fully Ionized Gases. John Wiley and Sons, Hoboken.
[16]	Perkins, F. (1973) Spread F and Ionospheric Currents. Journal of Geophysical Research, 78, 218-226. https://doi.org/10.1029/JA078i001p00218
[17]	Maeda, K.-I. (1977) Conductivity and Drifts in the Ionosphere. Journal of Atmospheric and Terrestrial Physics, 39, 1041-1053. https://doi.org/10.1016/0021-9169(77)90013-7
[18]	Heki, K. (2021) Ionospheric Disturbances Related to Earthquakes. In: Advances in Ionospheric Research: Current Understanding and Challenges, Wiley, Hoboken, 511-526. https://doi.org/10.1002/9781119815617.ch21
[19]	Oyama, K.I., Chen, C.H., Bankov, L., Minakshi, D., Ryu, K., Liu, J.H. and Liu, H. (2019) Precursor Effect of March 11, 2011 off the Coast of Tohoku Earthquake on High and Low Latitude Ionospheres and Its Possible Disturbing Mechanism. Advances in Space Research, 63, 2623-2637. https://doi.org/10.1016/j.asr.2018.12.042
[20]	Otsuka, Y., Kotake, N., Shiokawa, K., Ogawa, T., Tsugawa, T. and Saito, A. (2011) Statistical Study of Medium-Scale Traveling Ionospheric Disturbances Observed with a GPS Receiver Network in Japan. In: Abdu, M.A. and Pancheva, D., Eds., Aeronomy of the Earth’s Atmosphere and Ionosphere, IAGA Special Sopron Book Series 2, Springer, Berlin, 291-299. https://doi.org/10.1007/978-94-007-0326-1_21
[21]	Biswas, S., Kundum, S., Ghosh, S., Chowdhury, S., Yang, S.-S., Hayakawa, M., Chakraborty, S., Chakrabarti, S.K. and Sasmal, S. (2020) Contaminated Effect of Geomagnetic Storms on Pre-Seismic Atmospheric and Ionospheric Anomalies during Imphal Earthquake. Open Journal of Earthquake Research, 9, 383-402. https://doi.org/10.4236/ojer.2020.95022
[22]	Pulinets, S. and Ouozounov, D. (2011) Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) Model—A Unified Concept for Earthquake Precursors Validation. Journal of Asian Earth Sciences, 41, 371-382. https://doi.org/10.1016/j.jseaes.2010.03.005
[23]	Kuo, C.L., Lee, L.C. and Huba, J.D. (2014) An Improved Coupling Model for the Lithosphere-Atmosphere-Ionosphere System. Journal of Geophysical Research: Space Physics, 119, 3189-3205. https://doi.org/10.1002/2013JA019392
[24]	Mizuno, A. and Takashima, K. (2013) Continuous Measurement of Current in Air and Possible Relation with Intense Earthquake. Journal of Electrostatics, 71, 529-532. https://doi.org/10.1016/j.elstat.2012.11.015
[25]	Igarashi, K., Tsuchiya, T. and Umeno, K. (2020) Characteristics of Anomalous Radio Propagation before and after the 2011 Tohoku-Oki Earthquake as Seen by Oblique Ionograms. Open Journal of Earthquake Research, 9, 100-112. https://doi.org/10.4236/ojer.2020.92007