Anomalous Variation in GPS Based TEC Prior to the 5 Earthquakes in 2009 and 2010 ()
1. Introduction
Extensive researches have resulted in the recognition of many earthquake precursors in the field of the lithosphere, atmosphere and ionosphere. The effect of pre-seismic activity on ionosphere can be investigated by using the ionosphere electron density. The variation of electron concentration related to earthquake has been reported by Ondoh [1] and Liperovsky et al. [2] . Large scale ionospheric irregularity in the F-region is observed by many authors [3] - [6] . The electric field triggered by an earthquake preparatory process generates E × B drift. Depletion and enhancement in density profile may be a result of earthquake associated E × B drift [7] . GPS based total electron content (TEC) variation prior to an earthquake was observed by [8] . TEC data can be obtained from dual frequency GPS receiver, which can help to provide a better understanding of the upper atmosphere environment. The TEC is measured in TECU and 1 TECU = 1016 el/m2. Also the atmosphere-ionosphere interaction is studied earlier [9] . Lithosphere can be studied by using ground based sensors such as ultra-low frequency (ULF), very low frequency (VLF), and extreme low frequency (ELF) [10] . Analytical studies show that ionosphere disturbances caused by seismic activity can be observed several days prior to an earthquake. Earlier the relationship between seismic event and unexpected variation in the ionosphere has been reported by [11] and [12] .
The present paper reports the analysis of TEC data obtained from five different IGS Stations. These data are for the period corresponding to five different earthquakes that have occurred in different parts of the world.
The atmospheric refractivity is proportional to the pressure and the temperature gradient of the lower atmosphere. Prior to an earthquake, vertical profile of humidity, pressure and temperature changes [13] . These local parameter variations lead to change the radio refractive index and hence affect the ray path of satellite signals [14] . During the seismic action role of change refractive index variation was studied by [13] [15] . Present paper points out that 2 - 7 days prior to an earthquake, high Vertical TEC (VTEC) values are observed from a large number of satellites indicating an enhancement in the field of view of the receiver. This is only possible when positions of satellites remain below and come within the view of the receiver.
2. Data and Analysis
2.1. Earthquake Data
In this paper, we take 5 earthquakes into consideration: Loyalty Island, Izu Islands, Ryukyu Islands, Haiti and Chile. Table 1 gives the epicentral locations and details of IGS Stations of the selected earthquakes. The earthquake details are collected from http://earthquake.usgs.gov/.
2.2. TEC Data
The RINEX data obtained from GPS receivers, contain the C1 (C/A code pseudo range, in meters, on L1 frequency), P2 (P code pseudo range, in meters, on L2 frequency), L1 (L1 carrier phase, in cycles, on L1 frequency) and L2 (L2 carrier phase, in cycles, on L2 frequency) with a time resolution of the 30 s.
The slant total electron content (STEC) estimated, from an IGS data set in RINEX format, as
(1)
where 
(1227.60 MHz) and 
(1575.42 MHz) are current GPS broadcast frequencies.
STEC converted into VTEC using a suitable mapping function of different ionosphere pierce point (IPP) locations. The mapping function S(E) is defined as,
![]()
Table 1. Details of earthquake and IGS station.
(2)
(3)
where
(4)
= mean Earth radius, 6371 km, 
= elevation angle and 
is zenith angle, and 
= altitude of the IPP = 350 km, is the height of the ionospheric shell above the earth’s surface [16] .
To identify anomaly in ionospheric TEC, we compute in this paper the median 
of the 25 days VTEC and the associated inter-quartile range, (IQR), to construct the upper bound 
+ IQR and lower bound 
− IQR at a certain local time (LT).
The analytical study of 5 earthquakes listed in Table 1 is presented by single figure which have four sub-fi- gures: first figure (a) presents VTEC profiles for a period of 25 days, (b) presents disturbance storm time (Dst) index and (c) Shows solar flux (F-10.7) variation of the analyzed period with a purpose to refer to the geomagnetic and solar condition. The Dst-index data are obtained from the world data center, Kyoto, Japan through the website: http://wdc.kugi.kyoto-u.ac.jp/. The solar F10.7 cm data have been obtained from the website http://spidr.ngdc.noaa.gov/spidr. (d) Represents the plots (No. of satellites recording high VTEC - Elevation plots) plot for the 1 week prior and 1 week after the earthquake. As the magnitudes of the TEC records depend on satellite position when low-elevation satellite passes are considered, the greater the number of satellites seen by the receiver, the greater the scatter in the TEC profiles. Any modification in the number of satellites recording TEC seen by the receiver indicates ray path modification. The plots have been considered, so as to understand the effect of change in ray path caused due to modification in refractivity around the epicenter prior to an earthquake.
3. Results
This paper presents the results of variation in VTEC, prior to 5 earthquakes in 2009 and 2010 as summarized in Table 2.
3.1. Earthquake at Loyalty Island
The epicenter, date and magnitude of the Loyalty Island earthquakes are shown in Table 1. It is observed that prior to the earthquake on 19 January 2009 at Southeast of Loyalty Island there is an increase in VTEC diurnal profile. It is observed from Figure 1(a) that 3, 9, and 12 days prior to the earthquake there is enhancement greater than the upper bound in VTEC (represented by red line). Figure 1(d) represents the plot of total number of satellite’s observations that recorded high TEC > 16 TECU at different elevations: 1 week prior and 1 week after the earthquake. More number of satellites including elevation as low as 10˚ - 30˚ recording high TEC > 16
![]()
Table 2. Represents the summary of the data and results obtained for all the earthquakes included in the paper.
![]()
Figure 1. (a) VTEC profile for the period of 4 - 29 January 2009 of the Noumea France IGS station. VTEC diurnal profile indicating an increase in diurnal VTEC (red line) 3, 9, and 12 days prior to the earthquake. VTEC monthly mean has been shown (Violet line). (b) The figure displays Dst-index and (c) solar F-10.7 cm. The geomagnetic condition is found to quiet with small variation in Dst index. The star symbol represents the earthquake day. (d) Shows a number of satellite passes recorded at different elevation for 7 days prior and 7 days after the earthquake. The appearance of more number of low elevation satellites recording VTEC > 16 TECU 7, 6, 5 and 3 days before the earthquake comparable to those observed after earthquake.
TECU (16 TEC is selected as mean TEC for peak of 25 days) were observed 7, 6, 5 and 2 days prior to earthquake.
3.2. Earthquake at Izu Island
It is seen from Figure 2(a) that there is enhancement greater than the upper bound in TEC 2, 4, 6 and 7 days also there is decrease in VTEC less than the lower bound 9 and 10 days prior to the earthquake (represented by red line). As indicated in Figure 2(d) there are large no of satellite recorded 7, 6, 4 and 2 days prior to the earthquake on 12 August 2009.
3.3. Earthquake at Ryuku Island
As indicated in Figure 3(a) there are no enhancement greater than the upper bound in TEC prior to the earthquake on 17 August 2009 at Southeast of Ryukyu Island Japan. Figure 3(d) shows the large no. of satellite recorded 7 days prior the earthquake.
![]()
Figure 2. (a) VTEC profile for the period of 29 July - 22 August 2009 of the Ksmv, Japan IGS station. VTEC diurnal profile indicating an increase in diurnal VTEC (red line) 2, 4, 6 and 7 days prior to the earthquake. TEC monthly mean has been shown (violet line). (b) The figure displays Dst-index and (c) solar F-10.7 cm. The geomagnetic condition is found to quiet with small variation in Dst index. The star symbol represents the earthquake day. (d) Shows a number of satellite passes recorded at different elevation for 7 days prior and 7 days after the earthquake. The appearance of more number of low elevation satellites recording VTEC > 10 TECU 7, 6, 4 and 2 days before the earthquake comparable to those observed after earthquake.
3.4. Earthquake at Haiti
Figure 4(a) shows the enhancement greater than the upper bound in VTEC 2 and 3 days prior to the Haiti earthquake on 12 January 2010. Figure 4(d) shows the large number of satellite recording high TEC > 14 TECU were observed recorded 2, 3 and 5 days prior to the Haiti earthquake.
3.5. Earthquake at Chile
It is also observed increase in VTEC greater than the upper bound 7 days prior to earthquake on 11 March 2010 at Chile (Figure 5(a)). It is also observed large number of satellite recording high TEC > 26 TECU were recorded 7 day prior to earthquake (Figure 5(d)).
3.6. TEC Variation in the Earthquake Free Period
In order to check the behavior of TEC variation without an earthquake preparation period, the VTEC profiles for
![]()
Figure 3. (a) VTEC profile for the period of 3 - 22 August 2009 of the Tnml, Republic of China IGS station. VTEC diurnal profile indicating no significant enhancement prior to the earthquake. VTEC monthly mean has been shown (Violet line). (b) The figure displays Dst-index and (c) solar F-10.7 cm. The geomagnetic condition is found to quiet with small variation in Dst index. The star symbol represents the earthquake day. (d) Shows a number of satellite passes recorded at different elevation for 7 before the earthquake days prior and 7 days after the earthquake. The appearance of more number of low elevation satellites recording VTEC > 10 TECU 7 days before the earthquake comparable to those observed after earthquake.
the earthquake free period are obtained corresponding to one of the five IGS station SCUB in Cuba. The choice of station is based on the fact that, there was no continuous earthquake free period in other four stations for the year 2009 and 2010. It is observed that the TEC did not show significant variation during this period as shown in Figure 6.
4. Discussion
4.1. TEC Variation
There is no common opinion in scientific society on physical mechanism to explain observed ionospheric variation associated with earthquakes. Variation in ionospheric TEC prior to an earthquake has been reported by many researchers ( [7] [14] [17] - [20] ). Radon exhalation from the ground and consequent formation of large ion clusters as a result of ionization and plasma chemical reactions have been proposed as an agent of seismo-iono- spheric coupling mechanisms. Alpha decay of radon gas released from the ground can also ionize the atmosphere. Electric resistivity of lower atmosphere may change by this ionization of lower atmosphere, which could disturb the global electric circuit and redistribute ionospheric electrons [21] . Parrot [19] suggested that propagation of the direct wave due to compression of rocks close to earthquake epicenter could be more likely related to the piezoelectric and turbo-electric effect. Electric charge at earth's surface and electric current in the
![]()
Figure 4. (a) VTEC profile for the period of 29 December 2009 - 22 January 2010 of the sub, Cuba IGS station. VTEC diurnal profile indicating an increase in diurnal VTEC (red line) 2 and 3 days prior to the earthquake. VTEC monthly mean has been shown (Violet line). (b) The figure displays Dst-index and (c) solar F-10.7 cm. The geomagnetic condition is found to quiet with small variation in Dst index. The star symbol represents the earthquake day. (d) Shows a number of satellite passes recorded at different elevation for 7 days prior and 7 days after the earthquake. The appearance of more number of low elevation satellites recording VTEC > 14 TECU 2, 3 and 5 days before the earthquake January 2010 comparable to those observed after earthquake.
atmosphere-ionosphere system appears due to the stress of the rocks [22] [23] . Then electric field/current in the ionosphere and joule heating could modify or redistribute the electron concentration in height. Rising liquids under the ground would lead to the emanation of warm gases was proposed by [24] . Molchanov and Hayakawa [25] described that primary agent was an upward energy flux of atmospheric gravity waves (AGW) which were induced by gas-water release from earthquake preparatory zone. Penetration of AGW waves into the ionosphere leads to the modification of natural (background) ionospheric turbulence [2] [26] . It is possible that pre-seismic vertical electric field on the ground surface transformed into an electric field perpendicular to geomagnetic field line, which produces a perturbation over the F-region ionosphere. Once the F-region get perturbed within that Zone, it will pre-start to propagate along the conducting magnetic field lines and spread over wider areas discussed by [27] .
4.2. Atmospheric Refractivity
During seismic active periods due to modification in lower atmospheric variables like temperature, pressure and humidity the atmospheric refractive changes, which causes the variation in the number of satellite passes ob-
![]()
Figure 5. (a) VTEC profile for the period of 26 February- 21 March 2010 of VALP Chile IGS station. VTEC diurnal profile indicating an increase in diurnal VTEC (red line) 7 day prior to the earthquake. VTEC monthly mean has been shown (Violet line). (b) The figure displays Dst-index and (c) solar F-10.7 cm. The geomagnetic condition is found to be quiet with small variation in Dst index. The star symbol represents the earthquake day. (d)Shows a number of satellite passes recorded at different elevation for 7 days prior and 7 days after the earthquake. The appearance of more number of low elevation satellites recording VTEC > 26 TECU 7 day before the earthquake comparable to those observed after earthquake.
served thereby changing in the field of view (FOV) of the receiver. To understand the evolution of the FOV of the receiver (or transmitter) of the GPS system, the dynamic situation of the lower atmosphere should be examined. An increase in the propagation horizon of VHF signal has been reported by many researchers [26] [28] [29] . The modification of lower atmospheric parameters can be related to the radio refractive index (RRI) [15] .
As ray curvature C changes with refractive index gradient, its relation to the earth’s curvature does also changes therefore the concept of “effective earth’s radius” ![]()
is introduced and ratio of effective to true earth radius is given by:
![]()
(5)
where a = 6371 km and ![]()
is the rate of change of refractive index with height.
![]()
(6)
The relation of N with temperature, pressure and water vapor is given below
![]()
Figure 6. (a) VTEC profile for the period of 1 - 22 September 2009 of SCUB Cuba IGS station for the earthquake free period. (b) The figure displays Dst-index and (c) solar F-10.7 cm. The geomagnetic condition is found to be quiet with small variation in Dst index.
![]()
(7)
For the normal lapse rate ![]()
and K = 4/3 this factor defines a standard field of view of line of sight propagation links. For any change in ![]()
from the above value results in an increase or decrease of such path link hops [14] . Near to earthquake preparation zone normal lapse rate could reduce ![]()
and so the second term of the equation becomes less negative. This makes ![]()
more negative compared to normal lapse rate. A critical condition occurs when ![]()
making![]()
, i.e. earth appears as a flat surface to the satellite ray path thereby increasing the line of sight of radio horizon. Further, when![]()
, the ray takes guided path. This describes why a large number of satellites (beyond the standard Line of Sight (LOS)) are seen by the receiving antenna for low elevation of zones covering 10˚ - 30˚. With respect to the results shown in Figures 1-5 from the TEC data, it is seen that it may be possible to realize atmospheric changes prior to the Tohoku earthquake.
5. Conclusions
1) It is observed that we find anomalous variation in TEC for all 5 earthquakes included in the paper, which is in agreement with those reported by earlier literature, and hence it can be observed that ionosphere behavior changes, a few days prior to a major seismic shock within or far from the earthquake preparatory zone.
2) It is observed that 2 - 7 days prior an earthquake, high VTEC values are observed from a large number of satellites, indicating an enhancement in the FOV of the receiver. It means that we see satellite passes appear before earthquake which are not seen under normal situation. This is only possible when the ray path of GPS signal changes. At both IGS stations we have observed the change in the number of satellite passes prior to an earthquake. We determine that changes in atmospheric refractivity lead to change in FOV of satellite signal. The results agree with those reported by [14] [15] .
Acknowledgements
The authors would like to thank the data centers, IGS for RINEX data, Kyoto for Dst data and NOAA for solar flux data.