Anomalous TEC variations associated with the powerful Tohoku earthquake of 11 March 2011

Abstract. On 11 March 2011 at 14:46:23 LT, the 4th largest earthquake ever recorded with a magnitude of 9.0 occurred near the northeast coast of Honshu in Japan (38.322° N, 142.369° E, Focal depth 29.0 km). In order to acknowledge the capabilities of Total Electron Content (TEC) ionospheric precursor, in this study four methods including mean, median, wavelet transform, and Kalman filter have been applied to detect anomalous TEC variations concerning the Tohoku earthquake. The duration of the TEC time series dataset is 49 days at a time resolution of 2 h. All four methods detected a considerable number of anomalous occurrences during 1 to 10 days prior to the earthquake in a period of high geomagnetic activities. In this study, geomagnetic indices (i.e. Dst, Kp, Ap and F10.7) were used to distinguish pre-earthquake anomalies from the other anomalies related to the geomagnetic and solar activities. A good agreement in results was found between the different applied anomaly detection methods on TEC data.


Introduction
Widespread studies concerning earthquake forecasting during the last decades have resulted in the recognition of many earthquake precursors in the lithosphere, atmosphere and ionosphere.Daily variations of the ionosphere are affected by season, geographic location, thermospheric winds, traveling ionospheric disturbances, acoustic impulses to the atmosphere (such as those disturbances occurring after earthquakes, tsunamis, volcanic explosions and nuclear explosions), tropospheric weather (such as major storms), gravity waves over mountain ranges, and other unknown parameters.In the absence of significant solar and geomagnetic disturbances, ionospheric perturbations can be attributed to earthquakes.These anomalies usually happen in the D-layer, E-layer, and F-layer and may be observed 1 to 10 days prior to the earthquake and continue a few days after it (Hayakawa and Molchanov, 2002;Pulinets and Boyarchuk, 2004;Molchanov and Hayakawa, 2008;Akhoondzadeh, 2011).The effects of the pre-seismic activity on the ionosphere can be investigated using ionospheric electron density variations.
Ionosondes stations are an efficient means for monitoring peak electron density in the F2-layer during seismic activity.But, the spatial and temporal resolutions of ionosondes measurements are rather limited and therefore it is difficult to establish a systematic relationship between ionospheric disturbances and seismic occurrences.Currently, thousands of GPS receivers are used to monitor the Earth's surface deformations.Total Electron Content (TEC) data retrieved from GPS measurements have made a considerable contribution to the understanding of seismo-ionospheric variations (Liu et al., 2004;Devi et al., 2008;Zhao et al., 2008).Liu et al. (2004) statistically described the temporal parameters of the seismo-ionospheric precursors detected during 1-5 days prior to the earthquakes using TEC data for 20 major earthquakes in Taiwan.Hypotheses exist to explain the seismic electromagnetic mechanism based on geophysical and geochemical processes: -Direct electromagnetic wave production in a wide band spectrum by compression of rocks close to earthquake epicenter could be likely related to piezoelectric and tribo-electric effects (Parrot, 1995); -Rising fluids under the ground would lead to the emanation of warm gases (Hayakawa and Molchanov, 2002); -Activation of positive holes that can reach the ground surface (Freund, 2002); -Emission of radioactive gas or noble ions such as radon which changes air conductivity and leads to changes of ionosphere potential (Pulinets et al., 2003); -Penetration of atmospheric gravity waves (AGW), which are induced by the gas-water release from the earthquake preparatory zone into the ionosphere (Molchanov and Hayakawa, 2008).Preseismic electric field and its polarity cause the electrons in the F-layer to penetrate to lower layers and therefore to create anomaly in the ionospheric parameters.The thin layer of particles created before earthquakes due to ions radiation from the earth has a main role in transferring the electric field to above atmosphere and then to the ionosphere.The penetration of this electric field to the ionosphere was first analytically calculated by Park and Dejnakarintra (1973) and then applied to the seismo-electromagnetic process by Kim et al. (1994) and Pulinets et al. (2000).The vertical electric field on the ground surface is transformed into an electric field perpendicular to the geomagnetic field lines.This zonal component leads to plasma density anomalies, which are observed prior to the earthquake.(Namgaladze et al., 2009).In equatorial and low latitudes, TEC measurements indicate that the seismo-ionospheric variations lead to equatorial anomaly intensification.It could be possibly caused by an extra zonal electric field originated from the ground vertical electric (Pulinets, 2009).The actual physical mechanism for generation of an ionospheric precursor to an earthquake is currently unknown, and none of these mechanisms have been demonstrated to be reliable or strong enough to generate a disturbance to ionospheric electron density at the magnitudes observed.In this paper, applied methods determine the timing of ionospheric disturbances relative to the earthquake time, independent of mechanism.In this study, the capabilities of the mean, median, wavelet transform, and Kalman filter methods to detect the ionospheric TEC perturbations before the March 11, 2011 Tohoku earthquake are demonstrated.

TEC data
Recently, the extending network of GPS receivers has generated an increasing amount of data regarding the ionosphere state.TEC is the integrated number of the electrons within the block between the satellite and receiver or between two satellites.The GPS satellites transmit two frequencies of signals (f 1 = 1575.42MHz and f 2 = 1227.60MHz).The received GPS signals in ground stations contain many effects such as ionosphere, troposphere, hardware, and random errors.The ionosphere, unlike the troposphere, is a dispersive medium and its effects can be evaluated with measurement of the modulations on the carrier phases recorded by dualfrequency receivers (Sardon et al., 1994).
In this study, TEC variations have been analyzed using Global Ionospheric Map (GIM) data provided by the NASA Jet Propulsion Laboratory (JPL).The GIM is constructed from a 5 • × 2.5 • (Longitude, Latitude) grid with a time resolution of 2 h.GIM data are generated on a daily basis using data from about 150 GPS sites of the International GNSS Service (IGS) and other institutions.
Studies show that the irregularities in the electron concentration happen when the area on the ground surface occupied by the anomalous field exceeds 200 km in diameter (Pulinets et al., 2003).However, the previous studies show that the location of maximum intensity of the affected area in the ionosphere does not coincide exactly with the vertical projection of the epicenter of the impending earthquake and is shifted towards the equator in high and middle latitudes (Pulinets et al., 2003).Dobrovolsky formulated the radius of the affected area using the formula R = 10 0.43M where R is the radius of the earthquake preparation zone, and M is the earthquake magnitude (Dobrovolsky et al., 1979).
In the case of the Tohoku earthquake, the radius of the anomalous area based on the Dobrovolsky formula is estimated as 7413 km.Using reported geographic latitude and longitude concerning the earthquake's epicenter, TEC values of the nearest grid point in GIM have been analyzed.Since the location of the earthquake epicenter of Tohoku earthquake is 38.322 • N, 142.369 • E, and the nearest grid point to earthquake location in GIM is 37.5 • N, 140 • E, then the distance of the Tohoku earthquake location to the grid point is about 304.47 km, which makes the grid point data valid for the analysis.

Geomagnetic data
The ionospheric parameters are affected by solar geophysical conditions and geomagnetic storms especially in the equatorial and polar regions.Also, auroral activity has an important role in the mid-latitude ionosphere perturbations.The ionosphere current and equatorial storm-time ring current in periods of solar-terrestrial disturbances produce significant geomagnetic field disturbances.The parameters measured in such a way may display variations even in the absence of seismic activity making it difficult to separate preseismic ionospheric phenomena from the ionospheric disturbances due to the solar-terrestrial activities (Ondoh, 2008).Therefore, to discriminate the seismo-ionospheric perturbations from geomagnetic disturbances, the geomagnetic indices D st , K p , A p and F 10.7 accessed through NOAA (http: //spider.ngdc.noaa.gov)have been checked.The K p index monitors the planetary activity on a worldwide scale while the D st index records the equatorial ring current variations (Mayaud, 1980).The F 10.7 index represents a measure of diffuse, nonradiative heating of the coronal plasma trapped by magnetic fields over active regions, and is an excellent indicator of overall solar activity levels.The ionospheric effect of a geomagnetic storm has a global impact which is observable all over the world, while the seismogenic effect is observed only by stations with a distance of less than 2000 km from the potential epicenter (Pulinets et al., 2003).

Methodology
In this section, the mean, median, wavelet, and Kalman filter methods which will be used in detection of preseismic TEC anomalies are explained.

Anomaly detection using mean
Under the assumption of a normal distribution, the mean and the standard deviation of data are utilized to construct their higher and lower bounds in order to separate seismic anomalies from the background of natural variations.In calculation of statistical parameters, the length of the period was selected as about 49 days in order to avoid affects by the seasonal variations.The higher and lower bounds of the mentioned range can be calculated using the following equations: (1) Where x, x high , x low , µ, σ and Dx are the parameter value, higher bound, lower bound, mean value, standard deviation, and differential of x, respectively.According to this, if the absolute value of Dx would be greater than k, (|Dx| > k), the behavior of the relevant parameter (x) is regarded as anomalous.Variations of the ionosphere parameters depend on local time.Therefore, in calculation of TEC the µ and σ values were evaluated over the total time interval of interest for each period of 2 h in local time.

Anomaly detection using median
Daily variations of the ionosphere depend on season, geographic location, thermospheric winds, traveling ionospheric disturbances, and other unknown parameters.The unknown variations preclude the possibility of using methods based on normal distribution of data.As the fluctuation of the ionospheric parameters very often does not follow a Gaussian probability function, some researchers (Liu et al., 2004;Pulinets and Boyarchuk, 2004) use the median value and the interquartile.The median and the interquartile range of data are used to specify higher and lower bounds in order to distinguish seismic anomalies from the background of natural variations.The higher and lower bounds of the mentioned range can be calculated using the following equations: Where x, x high , x low , M, IQR and Dx are the parameter value, higher bound, lower bound, median value, interquartile range, and differential of x, respectively.For a given x, the values of M and IQR have been calculated for the whole period of interest for any interval of 2 h.If an observed TEC falls out of either the associated lower or higher bound, we conclude with a confidence level of about 80-85 % that a lower or higher abnormal signal is detected (Liu et al., 2004).

Anomaly detection using wavelet transformation
In this study, to obtain the time variability of the main periodicities, the wavelet transformation (Eq.7) has been applied on the TEC time series of earthquakes.
Where, s is the scaling factor, b is the location parameter, * is the complex conjugate of continuous wavelet function, and f (x) is the time series under analysis.Due to the variability pattern of data, the Daubechies 1-D wavelet has been applied.In this study, 1-D wavelet transformation has been applied for the whole period of interest for any interval of 2 h in local time (LT) to identify anomalies in the data.The low frequency seasonal components and high frequency noise have been eliminated using the components of the wavelet transform.The high perturbations of TEC are then detectable by wavelet coefficients greater than a pre-defined threshold value.In this study, µ + 2.5 × σ has been selected as an optimum threshold value to detect unusual values of the wavelet coefficients.µ and σ are the mean and the standard deviation parameters respectively.In quiet geomagnetic condition (i.e.K p < 2.5 , D st > −20nt and A p < 20), the wavelet coefficients of TEC values greater than µ + 2.5 × σ are regarded as anomaly values.Continuity in detected anomalies during several hours in each day before earthquake indicates that observed anomalies with longer time periods are unusual and might be related to impending earthquake.

Anomaly detection using Kalman filter
Kalman filter is a recursive solution to optimize the described systems in the state space.This filter is a set of mathematics equations to optimize prediction equations using estimation of state variables and minimization of error covariance.It is suitable for the stationary as well as dynamic and linear processes, and it can be applied for non-linear systems using Taylor expansion equations.The filter has high capabilities in the determination of inner variables and simultaneously solves state and measurement equations in order to reach optimized unobservable states.In other words, this method uses observed variables (y 1 , y 2 , ..., y t ) to estimate state (x i ) with minimum error.Depending on (i), filtering, prediction, and interpolation states are the cases as following definitions (Haykin, 2001): If i = t, i > t or i < t, this method is known as filtering, prediction or interpolation, respectively.Equations ( 8) and ( 9) are state and measurement equations: k t is the Kalman coefficient and must be defined based on the minimum of post-error covariance (Eq.11).
Regarding the mentioned equations, the measurements would be reliable when covariance of measurement error is close to zero.Kalman filter equations are classified into two categories: (1) time update; time retrieval equations update state and covariance matrixes based on the pre measurements (Eqs.12 and 13), (2) measurement update; measurement retrieval equations are for feeding back time update effects in the system and reaching towards an optimum state based on the measurements (Eqs.14, 15 and 16).
Therefore at the beginning, the prediction process is done; then it is corrected based on the observations and again prediction process is repeated.If however, the state and measurement equations are nonlinear (such as time series of earthquake precursors), they could be changed into linear equations using Taylor expansion which is called the extended Kalman filter.This is one of the striking characteristic of the Kalman filter (Haykin, 2001).In this study, the extended Kalman filter has been used to improve its parameters by a training process over the first half of total TEC data.If the difference between the observed TEC value (y t ) and the predicted TEC value is greater than a pre-defined threshold value (i.e.µ±2.5×σ ;µ and σ are the mean and the standard deviation parameters, respectively), the observed TEC value in quiet geomagnetic condition (i.e.K p < 2.5 , D st > −20 nt and A p < 20) is regarded as anomaly.

Results and discussion
This case study is focused on an earthquake which occurred near the northeast coast of Honshu in Japan with a magnitude of M w = 9.0 on 11 March 2011 at 05:46:23 UTC (UTC=LT-9:00).The characteristics of the Tohoku earthquake accompanied by its main foreshocks' information can be found in Table 1.
Figure 1a and b illustrates the variations of K p and A p geomagnetic indices, respectively, during the period of 1 February to 21 March 2011.An asterisk indicates the earthquake time.The X-axis represents the days relative to the earthquake day.The Y-axis represents the local time.The high geomagnetic activities are seen on 34 and 21 days before the earthquake onset.The high K p values between 22:00 and 24:00 LT on 10 days before the earthquake and also between 00:00 and 14:00 LT on 9 days before the event can be interpreted as high geomagnetic activities.The K p value reaches the values of 4, 4.5 and 5 between 02:00 and 06:00 LT, on the earthquake day and increases to a maximum value of 6.0 from 02:00 to 08:00 LT, 1 day after the main shock.The A p value reaches the maximum value of 67, 8 h before the main shock.These unusual variations of K p and A p indices from 1 day before to 1 day after the earthquake can hide pre-and post-seismic ionospheric anomalies.
Figure 1c shows the variations of D st geomagnetic index during the period of 1 February to 21 March 2011.The observed unusual D st values from 9 to 10 days and also 32 to 34 days before the earthquake indicate the high geomagnetic activities.It can be concluded that the detected perturbations on TEC variations during these periods might not be related to a seismic activity.The D st value exceeds the lower boundary value (i.e.−20(nT)) at 16:00 LT, 1 day before the earthquake and then gradually decreases and reaches the minimum value of −76.5 (nT) at earthquake time.The unusual variations of D st values are seen to 1 day after the  earthquake.This study examines the TEC variations during the period of Tohoku earthquake to find the pre-seismic anomalies in low geomagnetic activities (K p < 2.5, D st > −20 nt and A p < 20).
Figure 1d shows the variations of solar radio flux (F 10.7) during the period of 1 February to 21 March 2011.F 10.7 is often expressed in SFU or solar flux units.The F 10.7 value gradually increases from about 14 days before the earthquake and reaches the maximum value of 164.30 (SFU) on 8 March 2011, which is 3 days before the event.High levels of sunspot activity lead to improved signal propagation on higher frequency bands, although they also increase the levels of solar noise and ionospheric disturbances.These effects are caused by impact of the increased level of solar radiation on the ionosphere (Willson and Hudson, 1991).
After four years without any X-flares, the sun has produced two of the powerful blasts in less than one month, one taking place on 19 February 2011, and the other taking place on 9 March 2011.A solar flare is an intense burst of radiation coming from the release of magnetic energy associated with sunspots.March 9th ended with a powerful solar flare.Earth-orbiting satellites detected an X1.5-class explosion from behemoth sunspot 1166 around 23:23 UTC.This continues the recent trend of increasing solar activity associated with our sun's regular 11-yr cycle, and confirms that solar cycle 24 is indeed heating up, as solar experts have expected.Solar activity will continue to increase as the solar cycle progresses toward solar maximum, expected in the 2013 time frame (http://spaceweather.com/).
In addition, on 10 March 2011 around 06:30 UTC, a Coronal Mass Ejection (CME) did strike the magnetic field of the earth.This was a result of an M3 flare that took place on 7 March 2011 which released the fastest CME since September 2005.
Figure 2a shows TEC variations derived from GPS stations close to the epicenter.These TEC data have been downloaded via NASA JPL website (http://iono.jpl.nasa.gov).By visual inspection (without performing analysis), unusual TEC values are clearly seen around the earthquake day, especially between 3 days before to 1 day after the earthquake.
When implementing the mean method, Dx which will be called DTEC here, is calculated using Eq.(3). Figure 2b shows variations of DTEC during the period of 1 February to 21 March 2011 which is around the Tohoku earthquake  date.Figure 2c shows detected TEC anomalies using mean method based on:|DTEC| > 2.5.The noticeable anomalies are clearly seen from 4 days before to 1 day after the earthquake.The TEC anomaly reaches the maximum value (i.e.DTEC=3.0569) at earthquake time.This anomaly coincides with the minimum observed D st value at earthquake time.Then to distinguish pre-earthquake anomalies from the other anomalies related to the geomagnetic activities, the four conditions of |DTEC| > 2.5, K p < 2.5 , D st > −20nt and A p < 20 respectively, are jointly used using AND op-erator to construct the anomaly map (Fig. 2d).The TEC value exceeds the higher bound (µ + 2.5 × σ ), 3 days prior to the earthquake at 12:00 LT with the value of 7.26 % of the higher bound.It had also been reached to its maximum value (i.e.37.81 %) at 14:00 LT on the same date (Fig. 2d).It should be noted that the main foreshock with a magnitude of M w = 7.3 has been happened at ∼12:00 LT, 3 days before the event (Table 1) and it coincides with the highest unusual TEC variations observed prior to earthquake.It is seen that  the other detected anomalies in Fig. 2c have been masked by high geomagnetic activities.
After executing the median method, Fig. 3a shows variations of DTEC calculated using Eq. ( 6), during the period of 1 February to 21 March 2011.Disturbances in DTEC variations from 3 days to 1 day after the earthquake can be clearly seen.Figure 3b illustrates the unusual DTEC values based on: |DT EC| > 2.5.The DTEC value reaches the maximum value (DTEC=2.7547) at earthquake time.
When considering the geomagnetic indices, Fig. 3c shows the TEC anomaly map during quiet geomagnetic conditions.In other words, anomalous TEC values are only depicted at times when |DTEC| > 2.5, K p < 2.5 , D st > −20(nT) and A p < 20. Figure 3c illustrates that an increase (67.29 %) in TEC is clearly observed at 14:00 LT, 3 days before the earthquake.Variations of TEC clearly exceed the higher bound (M + 2.5× IQR) of the order of 32.07 % at 14:00 LT, 2 days before the earthquake.4b clearly shows TEC anomalies detected using wavelet transformation from 8 days before to 7 days after the earthquake.The peak of anomaly reaches the value of 25.80 % above the threshold value at earthquake time.When taking geomagnetic indices to consideration, some detected anomalies are masked by high geomagnetic activities (Fig. 4c).The resulted anomaly map using the wavelet transform represents an anomaly by the order of 19.2 % at 12:00 LT, 3 days before the event and also another anomaly (16 %) at 12:00 LT, 1 day prior to the earthquake.The peak of anomaly reaches to 22 % above the threshold value 3 days prior to the event at 14:00 LT (Fig. 4c).
Figure 5a represents the variations of detail coefficients from 1 February to 21 March 2011.The striking unusual variations are seen on 23, 24, 33, and 34 days before the earthquake that could not be related to the forthcoming earthquake.But the unusual variations during 4 days before to 1 day after the earthquake are considerable.Figure 5b illustrates the disturbances of detail coefficient variations without taking the high geomagnetic activities conditions into account.The anomaly map resulted using detail coefficients of applied wavelet transform in quiet geomagnetic activity represents an anomaly by the order of 12 % at 14:00 LT and also another anomaly (42 %) at 16:00 LT, 3 days prior to earthquake (Fig. 5c).
Figure 6a shows the differences between the predicted TEC values using the Kalman filter method and the observed TEC values during the 14 days before to 10 days after the earthquake.Half of data has been used for training the method to reach the optimum parameters.It is seen that these differences have reached high values during 3 days before to 1 day after the earthquake (Fig. 6a, b).Based on quiet geomagnetic conditions (K p < 2.5 , D st > −20(nT ) and A p < 20 ), an unusual increase of TEC (112.33 % above the threshold value) can be seen at 16:00 LT, 3 days before the earthquake in Fig. 6c.This figure also illustrates an increase of 15.75 % from the normal state 3 days before the earthquake at 12:00 LT.The characteristics of other detected anomalies are presented in Table 2.

Conclusions
In this study, the capabilities of four methods including mean, median, wavelet transform, and Kalman filter have been evaluated for appropriate detection of earthquake anomalies in TEC variations occurring before the Tohoku earthquake.The mean and median methods detect any unusual variations falling outside the predefined bounds.The approximation and detail coefficients of applied wavelet transform significantly show abnormal variations of the studied time series.The capability of the Kalman filter method in anomaly detection depends on the determination of the related parameters such as error covariance matrix in the training process.It can be concluded that the proposed method gets its credibility from the overall capabilities of the four integrated methods.This is done by accepting the major anomalies presented in all methods while neglecting the minor ones presented only by some methods.
The results also indicate that the highest deviations from the normal state that were regarded as anomaly appeared within the time interval 1-3 days before the earthquake.
As mentioned before, positive anomalies were detected in Tohoku earthquake.The sign of an anomaly could be justified by the ionosphere's large day-to-day variability owing to the solar irradiation variability, meteorological influences, and solar wind energy input (Rishbeth and Mendillo, 2001).Pulinets et al. (2003) also indicated that the sign of the anomaly could be related to the sign of the anomalous electric field at the ground surface and to the position of the studied ionospheric regions in relation to the epicenter location (east or west).
Because geomagnetic activity was not quiet during the days around the studied earthquake day, the observed ionospheric disturbances can not be interpreted independent of solar-terrestrial events.Therefore investigation of solarterrestrial coupling mechanism during seismic activity is necessary to understand the complex relationship between the TEC ionospheric anomalies and the occurrence of large earthquakes.However, in this study, only the pre-seismic TEC anomalies in geomagnetic quiet periods have been considered.

0 Fig. 1 .
Fig. 1.Data of geomagnetic indices during the period of 1 February to 21 March 2011 showing variations of (a) K p , (b) A p , (c) D st and (d) F 10.7.A black asterisk indicates the earthquake time.The X-axis represents the days relative to the earthquake day.The Y-axis represents the local time.

Fig. 2 .
Fig. 2. Results of analysis for the Tohoku earthquake showing (a) TEC variations, (b) DTEC variations after implementing the mean method,(c) detected TEC anomalies using mean method without considering the geomagnetic conditions and (c) detected TEC anomalies using mean method during quiet geomagnetic conditions.

Fig. 3 .
Fig. 3. Results of analysis for the Tohoku earthquake using median method, (a) DTEC variations, (b) detected TEC anomalies without considering the geomagnetic conditions and (c) detected TEC anomalies during quiet geomagnetic conditions.

Fig. 4 .
Fig. 4. Results of analysis for the Tohoku earthquake using the approximation coefficients of the wavelet transform method, (a) TEC variations after executing the method, (b) detected TEC anomalies without considering the geomagnetic conditions and (c) detected TEC anomalies during quiet geomagnetic conditions.

Fig. 5 .
Fig. 5. Results of analysis for the Tohoku earthquake using the detail coefficients of the wavelet transform method, (a) TEC variations after executing the method, (b) detected TEC anomalies without considering the geomagnetic conditions and (c) detected TEC anomalies during quiet geomagnetic conditions.

Fig. 6 .
Fig. 6. Results of analysis for the Tohoku earthquake using the Kalman filter method, (a) the differences between observed and predicted TEC values, (b) detected TEC anomalies without considering the geomagnetic conditions, and (c) detected TEC anomalies during quiet geomagnetic conditions.
P and N are probability distribution function and normal distribution function, respectively.Q and R are standard deviation parameters.F is the transition matrix taking the state x t from time t to time t + 1. H is the measurement matrix.If we suppose x t is real state at time t, then we can define pre-estimation error (e − t = x t − ∧x + t (post-estimation of state) using linear integration of ∧ x − t (pre-estimation of state) and measured error (y t − H ∧ x − t ) as Eq.(10).