Correlation Analysis for Total Electron Content Anomalies on 11th March, 2011
aa r X i v : . [ phy s i c s . g e o - ph ] J un Correlation Analysis for Total ElectronContent Anomalies on 11th March, 2011
Iwata and Umeno
Takuya Iwata, Ken UmenoDepartment of Applied Mathematics and Physics, Graduate school ofInformatics, Kyoto University, Yoshidahon-machi, Sakyo-ku, Kyoto, Japan
Ionosphere is a shell of a large amount of electrons and it is disturbed by vari-ous causes such as volcanic eruptions [
Heki (2006),
Dautermann et al. (2009)],solar flares [
Donnelly (1976)], earthquakes [
Cahyadi and Heki (2015)], and soon. Analyzing Total Electron Content (TEC) data is one of the most pop-ular method to monitor the ionosphere. TEC data contains the number ofelectrons integrated along the line of sight between the GNSS satellites andthe observation stations on the ground. Japan has a dense GNSS observationnetwork (the GNSS earth observation network, GEONET) to collect dailyTEC data in Japan. GEONET is composed of more than 1000 stations andcollects TEC data everyday. The GNSS data from GEONET are availablefreely for everyone.Ionospheric disturbances shortly after large earthquakes (Coseismic iono-spheric disturbances, CID) have been observed and reported for severalcases [
Calais and Minster (1995),
Ducic et al. (2003),
Liu et al. (2010)]. Asfor the 2011 Tohoku-Oki earthquake, various kinds of ionospheric distur-bances following the earthquake have been studied [
Astafyeva et al. (2011),
Tsugawa et al. (2011)].Preseismic ionospheric disturbances have also been studied and some signsbefore large earthquakes have been reported [
Heki (2011),
Heki and Enomoto (2015),
Jin et al. (2015)]. Thus, it is important to clarify whether such preseismic1nomalies really exist or not and establish a valid analysis method to pre-dict large earthquakes if preseismic anomalies exist. Such analysis methodsto detect preseismic ionospheric anomalies with the use of TEC data havenot been established until now. In this paper, we propose the correlationanalysis method to detect TEC anomalies. This method is practical becauseit does not need data after the corresponding earthquakes and can be easilyimplemented as automatic computation of correlations.
TEC changes over time are calculated by analyzing the phase differences be-tween the two carrier waves with the different frequencies (1.2 GHz and 1.5GHz) from GNSS satellites. TEC values obtained by analysing GNSS signalsare largely affected by the elevation angle of the satellites, i.e. the lower ele-vation angle increasing the apparent penetration length of the line of sight,and results in larger slant TEC values. That is, TEC value is susceptible tonot only the condition of the ionosphere but also to the elevation angle of thesatellites. This fact is one of the major obstacles to TEC data analysis. Toovercome this difficulty, vertical TEC (VTEC) is used. We can get VTECdata by computing the vertical components of TEC data, i.e. multiplyingTEC by the cosine of the elevation angle of the satellites and we can getVTEC.Not only TEC changes over time, but also the positions of GNSS satellitesat each time are importnat information. Following the custom for simplifica-tion, we make an assumption that there is a thin layer about 300 kilometersabove us and calculate the intersection of the line of sight with this layer.Such an intersection is called as Ionospheric Pierce Point (IPP), and its pro-jection onto the ground is called as Sub-ionospheric Point (SIP). IPP andSIP give the approximate position where the TEC data are derived.GEONET has a huge amount of TEC data in Japan, but only limitedamount of data are freely available via the Internet, i.e. TEC data for thelast three or four years are freely available. Here, we used TEC data obtainedbefore the 2011 Tohoku-Oki earthquake (Mw 9.0), which is the largest earth-quake in Japan during this period. We used TEC data from 2 months beforethe earthquake to the earthquake day for checking the validity of our pro-posed method to be described in the next section.2
Analysis method
The analysis method we propose here is based on the idea of applications ofthe correlation detection method used in VLBI (Very Long Baseline Interfer-ometry) and spreading spectrum communications technology to the analysisof TEC anomalies and given as follows:STEP 0. Choose the central GNSS station and set up three parameters, i.e. t sample , t test , and M . t sample is the length of data used for regression training and t test is the length of data used for regression test , and the value M de-notes the number of GNSS stations we use.STEP 1. At each station i and at each time epoch t , let SampleData bethe data from t to t + t sample and T estData be the data from t + t sample to t + t sample + t test .STEP 2. Fit a curve to SampleData by the least square method.STEP 3. Calculate a deviation of the
T estData from the model curve, rep-resenting as an ”anomaly”. The anomaly at station i at time t ′ is denotedas x i,t ′ , where −∞ < x i,t ′ < ∞ and we assume < x i,t ′ > = 0.STEP 4. Calculate a summation of correlations between the anomalies atthe central GNSS station and the surrounding stations as follows: C ( T ) = 1 N × M M X i =1 N − X j =0 x i,t + t sample + j ∆ t x ,t + t sample + j ∆ t (1) T = t + t sample + t test Here, N is the number of data in TestData , ∆ t is a sampling interval in TestData , which means ∆ t = t test / ( N − i = 0 means the centralGNSS station, where t sample = 2 . t test = 0 .
25 [hours] and M = 30[stations]. t sample and t test are significant parameters. If t sample is too large, we need a lotof data to calculate C ( T ) and if t sample is too small, we can not obtain reliablemodel curves in STEP 2. If t test is too large, deviations of the T estData fromthe model curves become large and if t test is too small, we can not catch thechanges of TEC data accurately. For these reasons, we set these parametersas written above. In STEP 1, if there are insufficient number of SampleData or T estData , then let C ( T ) be NA (missing value).In STEP 2, we can use arbitrary functions such as polynomial functions,Fourier series, and Gaussian functions for nonlinear regression. In this paper,we use a 7th polynomial function, Fourier series and Gaussian functions fordata extrapolation. In later section, we investigate the effectiveness of using3uch functions for nonlinear regression.The results do not depend on our choice of extrapolating functions, but theysensitively depend on the parameters set up in STEP 0. Generally, fitting a curve to data, i.e. nonlinear regression, needs to considersome points.As a result of regression analysis, we get a function f ( t ) as a model functionfor TEC data. Theoretically, we can choose arbitrary functions for f ( t ).For example, the most simple and often-used functions are polynomial func-tions. A D -th polynomial function is defined by f ( t ) = D X i =0 a i t i . (2)The number of parameters is D +1. Such polynomial functions are simple butnot the best candidate since f ( t ) increase or decrease infinitely as t increaseswhereas the data do not.On the contrary, Fourier series are periodic and bounded functions. A D -thFourier series is defined by f ( t ) = a + D X k =1 ( a k sin πkT t ! + b k cos πkT t !) . (3)The number of parameters is 2 D + 1 and T means the period of data. Inthis case, setting T = 24 [hours] is good since GNSS satellites have a periodof approximately 24 hours.We can also use Gaussian functions for f ( t ). A D -th order Gaussian extrap-olating function is defined by f ( t ) = a + D X i =1 a i exp − ( t − µ i ) σ i ! . (4)The number of parameters is D + 1. Hyper parameters µ i and σ i have to beset in advance. 4 Results
Figure 1 shows the result of correlation analysis on 11th March, 2011 whenthe 2011 Tohoku-Oki earthquake occured. We chose the GPS satellite 26and the 0214 (Kitaibaraki) GNSS station as the central station. The x-axisis time T in the coodinated universal time UTC and the y-axis representsan accumulated correlation, here, briefly wirriten C ( T ) defined in Eq. (1).The black line represents T=05:46 [UT] , i.e., the exact time when the 2011Tohoku-Oki earthquake occured. Here, we chose a 7th polynomial functionas a refference curve.This result indicates that there exists an anomalous trend before the earth-quake.Figure 2 shows the tracks of SIP and the 50 GNSS stations surrounding theKitaibaraki (0214) station used for the correlation analysis to get Fig.1. Figure 3 shows the results of correlation analysis on non-earthquake dayssuch as 2011/02/19 and 2011/03/01. The central GNSS station used for theanalysis is the same 0214 station (Kitaibaraki). In comparison to Fig. 1,correlations C ( T ) are quite small and actually quiet. The maximal value ofthe absolute value of correlation C ( T ) on these days is at most 5, whereas C ( T ) just before the earthquake on March 11, 2011 recorded more than 25,which is five times than the maximal value of absolute value of correlationof these normal days.Figure 4 compares the real TEC data and the results of our correlationanalysis. Correlation analysis detected the anomaly that is difficult to detectby merely looking real TEC data. Figure 5 shows the change of C ( T ) at all stations in Japan on the earthquakeday, 11 th March, 2011. These results suggest that the anomalies can be seennear the epicenter and just on the day of the earthquake occured. In Fig. 5,the anomalies can be seen also in southwest Japan. This anomalous area issmaller than near the epicenter. Although we are not sure of the origin of5his anomalous area, it could be explanation of this phenomenon that thesesouthwest Japan area is near the plate boundary as shown in Fig. 6.Figure 7 shows the change of C ( T ) and their locations at four other sta-tions on the earthquake day. In comparison to Fig. 1, C ( T ) are relativelysmall at the four stations. Note that here C ( T ) at Fukui seems to be par-tially anomalous, which is similar to Fig. 1. We think that such partialanomaly of C ( T ) at Fukui is observed because our correlation analysis candetect ”anomaly” and its SIP of 0580 Fukui station was near the epicenterof the 2011 Tohoku-Oki earthquake.Hence, we used the GPS satellite 26. There were other satellites above Japanat this time, however, adequate enough length of TEC data before the earth-quake is required to observe preseismic ionospheric condition. We think thatthe GPS satellite 26 had such appropriate length of TEC data because thetrack line between the SIP and the station is near the epicenter of the 2011Tohoku-Oki earthquake as shown in Fig.8. Figure 10 shows the comparison of results of correlation analysis with differ-ent extrapolating functions used in STEP 2.We used 7-th polynomial functions, 5-th polynomial functions, 3rd Fourierseries and 7-th Gaussian functions. It is seen that the results largely dependon the fitting functions. However, the anomalous trend before the earthquakecan be seen in every case.Figures 11,12 and 13 show the results of correlation analysis on non-earthquakedays with different fitting functions. The scale of the vertical axis vary withthe fitting functions. However, we can confirm that the correlation values onthe earthquake day are anomalous as compared with non-earthquake dayswhichever the fitting functions we use.
The physical mechanism of the anomaly has been researched so far and sev-eral models which explain the preseismic ionospheric anomalies are intro-duced [
Kuo et al. (2011),
Kuo et al. (2014)]. As for the 2011 Tohoku-Oki case,the preseismic anomaly is simulated by using these models [
Kuo et al. (2015)],where, the electric coupling model is used and the simulation results suggestthat this model can explain the mechanism of the ionospheric anomaliesbefore large earthquakes. To connect such anomaly with the physical mech-anism, we think that more studies to define an ”appropriate anomaly” basedon careful analysis. 6
Conclusion
In this paper, we introduced correlation analysis method into TEC dataanalysis. We detected the TEC anomalies about one hour before the 2011Tohoku-Oki earthquake by applying this analysis method. The TEC anoma-lies just about one hour before the earthquake showed characteristic patternswhich do not depend on our choice of extrapolation functions for nonlinearregression. Although we still do not know the actual physical mechanism re-sponsible for it, the present study based on the correlation analysis could beused for a link between the anomalies observed and the physical mechanism.Further investigation of other large earthquakes such as earthquakes in Chileand understanding of physical mechanism of the anomaly should be pursuedfor practical application of correlation analysis to detect preseismic anoma-lies of potential great earthquakes such as the 2011 Tohoku-Oki earthquakeon 11th March, 2011.
References [ Astafyeva et al. (2011)] Astafyeva, E. L., P. Lognonn´e, and L. M., Rolland(2011), First ionospheric images of the seismic fault slip on the exampleof the Tohoku-Oki earthquake,
Geophys. Res. Lett. , , L22104.[ Cahyadi and Heki (2015)] Cahyadi, M.N. and K. Heki (2015), Coseismicionospheric disturbance of the large strike-slip earthquakes in North Suma-tra in 2012: Mw dependence of the disturbance amplitudes,
Geophys. J.Int. , , 116–129.[ Calais and Minster (1995)] Calais, E., and J. B. Minster (1995), GPS de-tection of ionospheric perturbations following the January 17, 1994,Northridge earthquake,
Geophys. Res. Lett. , (9), 1045–1048.[ Dautermann et al. (2009)] Dautermann, T., E. Calais, and G. S. Mattioli(2009), Global Positioning System detection and energy estimation of theionospheric wave caused by the 13 July 2003 explosion of the Soufrire HillsVolcano, Montserrat,
Geophys. Res. , , B02202, 10.1029/2008JB005722.7 Donnelly (1976)] Donnelly, R. F. (1976), Empirical Models of Solar Flare XRay and EUV Emission for Use in Studying Their E and F Region Effects,
J. Geophys. Res. , (25),4745–4753.[ Ducic et al. (2003)] Ducic, V., J. Artru, and P. Lognonn´e (2003), Ionosphericremote sensing of the Denali Earthquake Rayleigh surface waves,
Geophys.Res. Lett. , (18),1951, 10.1029/2003GL017812.[ Heki (2006)] Heki, K. (2006), Explosion energy of the 2004 eruption of theAsama Volcano, Central Japan, inferred from ionospheric disturbance,
Geophys. Res. Lett. , , L14303, 10.1029/2006GL026249.[ Heki (2011)] Heki, K. (2011), Ionospheric electron enhancement precedingthe 2011 Tohoku-Oki earthquake,
Geophys. Res. Lett. , . L17312.[ Heki and Enomoto (2015)] Heki, K. and Y. Enomoto (2015), Mw depen-dence of preseismic ionospheric electron enhancements,
J. Geophys. Res.Space Phys. , , 7006–7020, 10.1002/2015JA021353.[ Iwata and Umeno (2015)] Iwata, T and K. Umeno (2015), Proposal of corre-lation analysis for detecting ionospheric electron anomalies as the precursorof huge earthquakes,
International Symposium on GNSS 2015 .[ Jin et al. (2015)] Jin, S. G., G. Occhipinti, and R. Jin (2015), GNSS iono-spheric seismology: Recent observation evidences evidences and character-istics,
Earth-Sci. Rev. , , 54–64, 10.1016/j.earscirev.2015.05.003.[ Kuo et al. (2011)] Kuo, C. L., J. D. Huba, G. Joyce, and L. C. Lee(2011), Ionosphere plasma bubbles and density variations induced by pre-earthquake rock currents and associated surface charges,
J. Geophys. Res. , , A10317, 10.1029/2011ja016628.[ Kuo et al. (2014)] Kuo, C. L., L. C. Lee, and J. D. Huba (2014), An im-proved coupling model for the lithosphere-atmosphere-ionosphere system,
J. Geophys. Res. Space Physics , , 3189–3205, 10.1002/2013JA019392.[ Kuo et al. (2015)] Kuo, C. L., L. C. Lee, and K. Heki (2015), Preseis-mic TEC changes for Tohoku-Oki earthquake: Comparisons betweensimulations and observations,
Terr. Atmos. Ocean. Sci. , , 63–72,10.3319/TAO.2014.08.19.06(GRT).[ Liu et al. (2010)] Liu, J. Y., H. F. Tsai, C. H. Lin, M. Kamogawa, Y. I.Chen, C. H. Lin, B. S. Huang, S. B. Yu, and Y. H. Yeh (2010), Coseismicionospheric disturbances triggered by the Chi-Chi earthquake,
J. Geophys.Res. , , A08303, 10.1029/2009JA014943.8 Tsugawa et al. (2011)] Tsugawa, T., A. Saito, Y. Otsuka, M. Nishioka, T.Maruyama, H. Kato, T. Nagatsuma, and K. T. Murata (2011), Ionosphericdisturbances detected by GPS total electoron content observation after the2011 off-the-Pacific coast of Tohoku Earthquake,
Earth Planets Space , ,875–879. 9 ainshock Central Station ID: 0214 time C ( t ) Figure 1: The result of correlation analysis on 11th March, 2011. The verticalaxis shows the correlation C ( T ) and the horizontal one the time t [UTC]. Theblack line indicates the exact time 05:46 [UTC] when the 2011 Tohoku-Okiearthquake occured. 10 IP epicenter Figure 2: The blue line represents the SIP track of the pair of the 0214station (Kitaibaraki) and the GPS satellite 26, and the red points repre-sent the location of surrounding 50 stations around the central station. Thestar represents the epicenter of the 2011 Tohoku-Oki earthquake, and theblack circle on the blue line represents the SIP position at the earthquakeoccurrence time. 11 time C ( t ) time C ( t ) time C ( t ) time C ( t ) Figure 3: The result of correlation analysis on non-earthquake days. Thevertical axis shows the correlation C ( T ) and the horizontal one the time t [UTC]. The day 40, 30, 20, 10 days before the earthquake, respectively. Weused 7th polynomial functions as fitting curves.12 time T E C time C ( t ) −100−5002:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 time T E C time C ( t ) −40−200 02:00 03:00 04:00 05:00 06:00 07:00 08:00 time T E C time C ( t ) Mainshock −40−30−20−100 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 time T E C Mainshock
Central Station ID: 0214 time C ( t ) Figure 4: The comparison between the real TEC (left) and the results ofcorrelation analysis (right) on earthquake day and non-earthquake days.13
30 135 140 145
Longitude La t i t ude −5.000 − −1.250− 2.500− 6.250− 10.000− 13.750− 17.500− 21.250− 25.000
130 135 140 145
Longitude La t i t ude −5.000 − −1.250− 2.500− 6.250− 10.000− 13.750− 17.500− 21.250− 25.000
130 135 140 145
Longitude La t i t ude −5.000 − −1.250− 2.500− 6.250− 10.000− 13.750− 17.500− 21.250− 25.000
130 135 140 145
Longitude La t i t ude −5.000 − −1.250− 2.500− 6.250− 10.000− 13.750− 17.500− 21.250− 25.000 Figure 5: The change of C ( T ) at epochs shortly before the 2011 Tohoku-Oki earthquake, i.e. 1 hour, 25 minutes, 7 minutes, and 4 minutes beforethe actual time 05:46 [UTC] when the main shock of the 2011 Tohoku-Okiearthquate occured on March 11, 2011. The black circle represents the epi-center. 14
30 135 140 145
Longitude La t i t ude (cid:1) (cid:0) (cid:2) (cid:3) (cid:4) (cid:5) (cid:6) (cid:7) (cid:8) (cid:9) Figure 6: The result of correlation analysis on 11th March, 2011. The blackcircle represents the epicenter of the 2011 Tohoku-Oki earthquake. The bluelines present the plate boundaries. 15 ainshock
Central Station ID: 0896 time C ( t )
130 135 140 145
Longitude La t i t ude SIP epicenter
Mainshock
Central Station ID: 0580 time C ( t )
130 135 140 145
Longitude La t i t ude epicenterSIP Mainshock
Central Station ID: 0053 time C ( t )
130 135 140 145
Longitude La t i t ude epicenterSIP Mainshock
Central Station ID: 1028 time C ( t )
130 135 140 145
Longitude La t i t ude epicenterSIP Figure 7: The change of C ( T ) on 11th March, 2011, at Aomori (0896), Fukui(0580), Ishikawa (0053), Okayama (1028) stations, and their positions in themap. The vertical axis shows the correlation C ( T ) and the horizontal one thetime t [UTC]. The red curves represent the tracks of the SIPs from 4:00 [UT]to 6:00 [UT], and the black points represent the SIPs when the earthquakeoccured. The black lines indicate the exact time 05:46 [UTC] when the 2011Tohoku-Oki earthquake occured. 16 IP epicenter Figure 8: The blue line represents the SIP track of the pair of the 0214station (Kitaibaraki) and the GPS satellite 26, and the red points repre-sent the location of surrounding 50 stations around the central station. Thestar represents the epicenter of the 2011 Tohoku-Oki earthquake, and theblack circle on the blue line represents the SIP position at the earthquakeoccurrence time. 17 ainshock
Central Station ID: 0896 time C ( t )
30 135 140 145
Longitude L a t i t u d e SIP epicenter
Mainshock
Central Station ID: 0580 time C ( t )
30 135 140 145
Longitude L a t i t u d e epicenterSIP Mainshock
Central Station ID: 0053 time C ( t )
30 135 140 145
Longitude L a t i t u d e epicenterSIP Mainshock
Central Station ID: 1028 time C ( t )
30 135 140 145
Longitude L a t i t u d e epicenterSIP Figure 9: The change of C ( T ) on 11th March, 2011, at Aomori (0896), Fukui(0580), Ishikawa (0053), Okayama (1028) stations, and their positions in themap. The vertical axis shows the correlation C ( T ) and the horizontal one thetime t [UTC]. The red curves represent the tracks of the SIPs from 4:00 [UT]to 6:00 [UT], and the black points represent the SIPs when the earthquakeoccured. The black lines indicate the exact time 05:46 [UTC] when the 2011Tohoku-Oki earthquake occured. 18a) (b)(c) Mainshock
Central Station ID: 0214 time C ( t ) (d) Mainshock
Central Station ID: 0214
024 04:30 05:00 05:30 06:00 06:30 time C ( t ) Mainshock
Central Station ID: 0214 time C ( t ) Mainshock
Central Station ID: 0214 time C ( t ) Figure 10: The comparison of results of correlation analysis with differentfunctions on the earthquake day. 0214 station and satellite 26 are used. (a)7-th polynomial function (b) 5-th polynomial function (c) 3rd Fourier series(d) 7-th Gaussian function. The black lines indicate the exact time 05:46[UTC] when the 2011 Tohoku-Oki earthquake occured.19
24 07:00 07:30 08:00 08:30 09:00 09:30 time C ( t )
024 06:30 07:00 07:30 08:00 08:30 time C ( t ) time C ( t ) time C ( t ) Figure 11: The result of correlation analysis on non-earthquake days. Thevertical axis shows the correlation C ( T ) and the horizontal one the time t [UTC]. The day 40, 30, 20, 10 days before the earthquake, respectively. Weused 5th polynomial functions as fitting curves.20 .02.55.07.510.0 07:00 07:30 08:00 08:30 09:00 09:30 time C ( t ) time C ( t ) time C ( t ) time C ( t ) Figure 12: The result of correlation analysis on non-earthquake days. Theday 40, 30, 20, 10 days before the earthquake, respectively. We used 3rdFourier series as fitting curves. 21
246 07:00 07:30 08:00 08:30 09:00 09:30 time C ( t ) time C ( t ) time C ( t ) time C ( t ) Figure 13: The result of correlation analysis on non-earthquake days. Theday 40, 30, 20, 10 days before the earthquake, respectively. We used 7thGaussian functions as fitting curves. 22 ainshock
024 4.5 5.0 5.5 6.0 6.5
TIME C O R ainshock
024 4.5 5.0 5.5 6.0 6.5
TIME C O R ainshock
024 5 6 7
TIME C O R ainshock
024 4.5 5.0 5.5 6.0 6.5
TIME C O R
30 135 140 145
Longitude La t i t ude − − − − − − − − − −−
Longitude La t i t ude − − − − − − − − − −−
30 135 140 145
Longitude La t i t ude − − − − − − − − − −−
Longitude La t i t ude − − − − − − − − − −−
30 135 140 145
Longitude La t i t ude − − − − − − − − − −−