Natural Hazards and Earth System Sciences Earthquake forecasting : a possible solution considering the GPS ionospheric delay

The recent earthquakes in L'Aquila (Italy) and in Japan have dramatically emphasized the problem of natural disasters and their correct forecasting. One of the aims of the research community is to find a possible and reliable forecasting method, considering all the available technologies and tools. Starting from the recently developed research concerning this topic and considering that the number of GPS reference stations around the world is continuously increasing, this study is an attempt to investigate whether it is possible to use GPS data in order to enhance earthquake forecasting. In some cases, ionospheric activity level increases just before to an earthquake event and shows a different behaviour 5-10 days before the event, when the seismic event has a magnitude greater than 4-4.5 degrees. Considering the GPS data from the reference stations located around the L'Aquila area (Italy), an analysis of the daily variations of the ionospheric signal delay has been carried out in order to evaluate a possible correlation between seismic events and unexpected variations of ionospheric activities. Many different scenarios have been tested, in particular considering the elevation angles, the visibility lengths and the time of day (morning, afternoon or night) of the satellites. In this paper, the contribution of the ionospheric impact has been shown: a realistic correlation between ionospheric delay and earthquake can be seen about one week before the seismic event. in the friction, generation of electric charges and chemical reactions. The Earth's crust is the rigid external shell of our planet, and it consists of the continental and the oceanic crusts. The slow movement of these plates over the asthenosphere layer and the ocean floor extension is called plate tectonics. The movements of these plates, such as a collision of these plates one under another (convergent or subduction boundaries) or a plate separation (divergent boundaries), lead to a strain accumulation within the Earth's crust, to mechanical deformations and crust rupture when the


Introduction
Earthquake physics is a very complex and broad topic.It involves many scales of the Earth's crustal structure, ranging from tectonic plates to the microscopic processes involved Correspondence to: M. De Agostino (mattia.deagostino@polito.it) in the friction, generation of electric charges and chemical reactions.
The Earth's crust is the rigid external shell of our planet, and it consists of the continental and the oceanic crusts.The slow movement of these plates over the asthenosphere layer and the ocean floor extension is called plate tectonics.The movements of these plates, such as a collision of these plates one under another (convergent or subduction boundaries) or a plate separation (divergent boundaries), lead to a strain accumulation within the Earth's crust, to mechanical deformations and crust rupture when the deformation exceeds the mechanical strength limit.
The intensity of an earthquake is estimated from the oscillations that are created by different kinds of seismic waves (usually, surface and body waves).
The mechanism which controls the generation of an earthquake is still not well known.Many theories have been suggested, starting with the early "Elastic Rebound Theory" (Reid, 1911), SOC theories (Main 1995), up to the latest one which states that an earthquake is "a frictional phenomenon rather than a fracture one" (Scholz, 1998).
The term "earthquake forecasting" refers to the knowledge of the earthquake prognostic parameters -the location, the time of occurrence and its magnitude -for some time before it takes place.The prognostic time window can be distinguished as: long-term, referring to a time window of some decades of years; medium-term, referring to a time window of a few years (2-3); and short-term, referring to a time window on the order of up to a couple of months; while the term "immediate" is sometimes used when the time window is on the order of a few days.
The scientific literature which concerns earthquake forecasting is composed of several numbers of papers.Since the mid 1970s, in fact, seismologists have been confident that earthquake prediction could be achieved within a short period of time (Cicerone et al., 2009).This assumption is the result of the first successful prediction of a large earthquake in Haicheng -China (1975), with a the magnitude equal to 7.4.Over the last half century, some ionospheric variability, observed for hours to days before the earthquake, has been suggested as an ionospheric precursor (Pulinets et al., 2004;Zakharenkova et al., 2006;Xu et al., 2011).However, even today the possibility to identify ionospheric precursors is controversial.One of the key points is whether it is possible to distinguish the solar, geomagnetic (Duma and Ruzhin, 2003) and even tropospheric contributions from the observed ionospheric variability in order to identify likely precursor ionospheric signatures.
The purpose of this paper is to evaluate the efficacy of the ionospheric Total Electron Content (TEC) as an additional earthquakes precursor, considering a dramatic Italian case study: the L'Aquila earthquake (6 April 2009, magnitude 6.3).In this paper, Global Positioning System (GPS) data from some reference stations close to the earthquake epicentre are used in an "one-way" positioning algorithm.This procedure allows to separate the GPS signal ionospheric propagation delay (directly correlated with the TEC) from the geometric errors (i.e.multipaths, random errors) and the propagation delay due to gases in the troposphere.
In the following, after a short review about some possible earthquake precursors, the "one way" positioning algorithm used to compute the GPS ionospheric delay is described in detail.The case study (L'Aquila earthquake) and the computation method are presented in the second part.A comparison between the obtained results and the effective earthquake map concludes the paper.

Earthquake forecasting: some possible precursors
In order to predict an earthquake and in particular to identify an effective earthquake precursor, many different physical quantities and mechanisms have been studied (Eftaxias et al., 2011).Some of the observations realized before a strong earthquake are: the seismic gap, that is the absence of normal seismicity in a seismogenic area for a long period (e.g.Papadopoulos et al., 2010), seismic quiescence, that is the drop in seismicity below its normal level, the change in Earth resistivity (e.g.Dologlou, 2008), the emission of electro-magnetic waves (e.g.Gladychev et al., 2001), the increase in Radon emission from the ground (e.g.Richon et al., 2003), geodetic variations (abnormal ground elevations, e.g.Sobolev, 2011), changes in the chemical composition of underground water (e.g.Ryabinin et al., 2011), change in temperature of the aquifers, changes in the Earth's magnetic field, changes in the Earth's electric field (e.g.Konstantaras et al., 2008), observed changes in the electron density of the ionosphere, and strange animal behaviour, to name some of the most important and well-known.
First, the precursors should be ranked in two categories: the time of appearance before the earthquake and their confidential merit.All available seismic information should be used, starting from the seismic regioning, calculation of the seismic risk, and finishing with the most recent techniques based on the self-organized criticality.The cumulative principle should be used, adding each new indication that appears to the expert alarm system.The short-term precursors, including all types of physical, geochemical, electromagnetic and biological monitoring, should be finally processed.

Radon emissions
Radon emission monitoring is one of the most widely used geological study techniques.As it is known, Radon is the product of radium decay.Radon is an ideal indicator in geological research because it is generated continuously in any geological structure.Its concentration loss due to decay and due to migration into the atmosphere is always compensated by a new production.The loading-unloading process during earthquake preparation is the reason for Radon concentration variations before an earthquake.Many authors, such as Teng (1980), Hauksson (1981), Talwani et al. (1981) and Richon et al. (2003), studied Radon variations before earthquakes and used them as a short-term earthquake precursor.
Most of the Radon anomalies began within 30 days of the earthquake, even though there does not appear to be any diagnostic behaviour of either the beginning or the end of a gas anomaly that gives a consistent clue about when an earthquake is going to happen.

Groundwater level change observations
Changes in groundwater levels have been observed before certain earthquakes and are believed to be in response to a volumetric strain in the Earth's crust (Cicerone et al., 2009).Many of the characteristics of the groundwater change precursors documented in this study, such as the time of the initialization of the anomalies, the dependence of the amplitude of the anomaly and epicentral distance, seem to parallel the same characteristics as in the Radon gas anomalies.

Surface deformations
There has been longstanding interest in looking for surface deformations (uplifts, downdrops, tilts, strains, strain rate changes, etc.) before earthquakes (Rikitake, 1976).Many crustal earthquakes, with magnitudes equal to or greater than 6, have been associated with Earth surface deformations.
Surface levelling and laser-ranging geodetic measurements are the most accurate way to document ground deformations over regions that are tens of kilometres in dimension.However, such measurements are time consuming and expensive, and the feasible time between individual measurements is from months to years.However, modern GPS and satellite-based SAR interferometry measurements are now available to produce geodetic position changes with individual measurements separated by minutes to days.The reported deformations took place months to days before the earthquakes, and the largest amplitude strains and tilts seem to be associated with the largest earthquakes.

Electromagnetic phenomena
Earthquake preparation is usually accompanied by electromagnetic phenomena in different frequency bands (Uyeda et al., 2009).Summarizing, we can state that before strong earthquakes in fair weather conditions, anomalies of the atmospheric electric field have been observed in the form of electric field increases.According to a previous study (Hao et al., 2000), the time of the anomalous electric field appearance could be more than one month before the earthquake.An objection often put forward against seismic electric precursory signals is that these signals may be generated by ionospheric and magnetic anomalies of the Earth's normal field.Consequently, the question that arises (Koulouras et al., 2009) is: is it possible to discriminate electrical signals induced by magnetic and ionospheric anomalies from the electrical signals that are generated at the focal area?

Ionospheric electron content
The ionosphere is a part of the upper atmosphere where there are enough electrons and ions to effectively interact with electromagnetic fields.The electron and neutral densities of the upper ionosphere are shown in the Fig. 1 (Pulinets and Boyarchuk, 2004).
The most complete and developed model has been created for the observed variations of the electron density in the ionosphere associated with an earthquake preparation process.A complete version of its description can be found in Pulinets and Boyarchuk (2004).Ionospheric precursors are part of the more general physical process of earthquake preparation and it is probably the youngest precursor method.The set of plasma-chemical reactions that start as a result of ionization by Radon involves changes in the water vapour and of the electron content in the near ground layer of the atmosphere.
However, statistical ionospheric precursors lead us to say that (Pulinets, 2006): 1. Ionopheric variations reveal the local time dependence (in the form of the sign of the critical frequency deviation from the undisturbed level), as well as "day-before-the-shock" dependence.
2. There exists an earthquake magnitude threshold (equal to 5 degrees) from which the ionosphere starts to "feel" the earthquake preparation process.
3. The anomalous variations in the ionosphere appear, on average, within a time interval of 5 to 1 days before the seismic shock.
4. The size of the modified area in the ionosphere increases with the magnitude of the earthquake.
The ionospheric effect can be estimated considering different data sources, such as the Very Low Frequency (VLF) emissions (Parrot, 1994), ground based measurements of the ionospheric precursors of earthquakes (Chen et al., 1999;Pulinets et al., 2002), and local plasma parameter measurements (Afonin et al., 1999).Different satellite techniques (topside sounding, ionospheric tomography, radio occultation) have been realized taking into account the measurements of the TEC distribution.
In recent years, the use of GPS data derived from reference stations is becoming a useful way of computing a very detailed TEC map to process the signal ionospheric propagation delay from satellite to receiver.Many important and dramatic earthquakes have been analyzed considering the TEC variation (Plotkin, 2003;Gahalaut et al., 2006;Huang et al., 2009;Borghi et al., 2009;Singh et al., 2009;Liu et al., 2010;Hasbi et al., 2011).This has led to interesting results and to confirm that ionospheric activity as an earthquake precursor can be considered.
the propagation velocity and the refractive index for GPS signals are frequency-dependent.The ionospheric effect has the same absolute value for pseudorange and carrier phase measurements, but with opposite signs.
Irregularities in the ionosphere produce short-term signal variations.These scintillation effects may cause a large number of cycle slips because the receiver cannot follow the short-term signal variations and fading periods.Scintillation effects mainly occur in a belt along the Earth's geomagnetic equator and in the polar auroral zone.
The primary purpose of the second frequency in the GPS satellite constellation is to eliminate the effect of the ionosphere on signal propagation.The errors associated with ionospheric refraction, in fact, are dependent on the frequency signal, and thus are different for the L1 and L2 frequencies.They can be eliminated through linear combinations of them.
The ionosphere may be characterized as that part of the upper atmosphere where a sufficient number of electrons and ions are present to affect the propagation of radio waves.The spatial distribution of electrons and ions is mainly determined by photo-chemical and transportation processes.The state of the ionosphere may be described by the electron density in units of electrons per cubic meter.The TEC is an important descriptive quantity for the Earth's ionosphere.TEC is the total number of electrons present along a path between two points, with units of electrons per square meter, where 10 16 electrons m −2 = 1 TEC unit (TECU).Therefore, TEC is directly correlated to the signal ionospheric propagation delay, which can be computed, for example, using a differential network or a precise point positioning.
GPS-derived ionosphere models that describe the deterministic component of the ionosphere are usually based on the so-called Single-Layer Model (SLM), as outlined in Fig. 2.This model assumes that all free electrons are concentrated in a shell of infinitesimal thickness.
The SLM mapping function F I may be written using the equation: where: z,z are the zenith distances at the height of the station and the single layer, respectively; R is the mean radius of the Earth; H is the height of the single layer above the Earth's surface.
In order to realize a TEC map, the so-called geometry-free linear combination, which principally contains ionospheric information, is considered.The equations for un-differenced pseudorange and carrier-phase observables can be written: where: P 4 , 4 are the geometry-free pseudorange and carrier phase observables (in meters); a is a constant equal to 4.03 × 10 17 m s −2 TECU −1 ; f 1 ,f 2 are the frequencies associated with the two carrier phases; F I (z) is the mapping function evaluated at the zenith distance z (see Eq. 1); E (β,s) is the vertical TEC (in TECU) as a function of geographic or geomagnetic latitude β and sun-fixed longitude s; B 4 ,b 4 are constant biases (in meters) due to the initial phase ambiguities and to the pseudorange differences, respectively.
"One-way" positioning, or Precise Point Positioning (PPP), is an absolute positioning technique, which is made without differential techniques (zero-difference positioning), using data collected by only one receiver, accurate orbital and satellite clock data and error models of atmospheric delays.The previous equations, referring to a single receiver and a single satellite, from metric point of view are equal to: where: ρ is the range between the receiver and the satellite; c is the speed of light (in meters per second); dt is the satellite clock error; dT is the receiver clock error; Tr is the tropospheric propagation delay; α is the quadratic ratio between the two GPS carrier-phase frequencies f 1 ,f 2 ; I is the ionospheric propagation delay; λ 1 ,λ 2 are the wavelenghts associated with the two carrier phases; N is the integer carrier phase ambiguity; m P is the pseudorange multipath delay; m is the carrier phase multipath delay; e P are the random errors on pseudorange observables; e are the random errors on carrier phase observables.In matrix notation, the previous equations become (Leick, 2004): It should be noted that this expression is independent of clocks and geometry (receiver and satellite coordinates and tropospheric delay), and thus gives origin to the name geometry-free.The objective of this method is to use this observable in order to estimate the ionospheric signal propagation delay over the carrier phase and pseudorange measurements.In this equation, the effects due to the tropospheric propagation delay and the signal multipath are separated from the ionospheric contribution.Therefore, the state vector devoted to solve could be defined considering the system at each epoch ("step by step"), or using a Kalman filter procedure.In the first case, a raw value of ambiguity is evaluated, whereas in the second case, it is possible to obtain a filtered value of the phase ambiguity that converges, more or less quickly, to an integer value (the so-called fixed ambiguity).The simplest solution in the Kalman Filter is to use a transition matrix F equal to a four-by-four identity matrix, utilizing the following measurement and system noise covariance matrices (meters): (5) 50 2 0 0 0 0 10 2 0 0 0 0 0.1 2 0 0 0 0 0.1 2 After the Kalman filtering, an observation smoothing could be performed (called Kalman smoothing).In this way, at each epoch, all measurements are used, and, ultimately, this solution looks like the solution obtained with the Least Square method.
The ionospheric propagation delay, in units of distance, computed using the "one-way" positioning for satellite PRN29, using the data from the reference station of L'Aquila (AQUI) is shown in Fig. 3.
The filtering and the smoothing procedures markedly reduce the RMS values of the ionospheric propagation delay.In addition, the observations can be weighted with the satellite elevation angle in order to reduce the influence of noise at low elevation angles; this calculation is very important for a single receiver positioning.The elevation-dependent weight proposed by Huber ( 2003) is used: where: z is the satellite zenith distance; a is a coefficient 0a 1.In this case a = 0.3.

Case study: the L'Aquila earthquake
On 6 April 2009, the middle of Italy was subjected to a catastrophic event: a violent earthquake seriously damaged the Regione Abruzzo (Fig. 4).The main event occurred at 03:32 local time (01:32 UTC), and was rated 6.3 on the moment magnitude scale; its epicentre was near L'Aquila, Abruzzo capital, which, together with the surrounding villages, suffered the most damage.There had been several thousand foreshocks and aftershocks since December 2008, more than thirty of which had a Richter magnitude above 3.5.After the main shock, 256 other tremors have been registered.
Following the L'Aquila earthquake, many researchers and research groups have analyzed the available data (GPS measurements, InSAR data) in order to describe the displacements and deformations that occurred during the earthquake as clearly as possible.
In addition, some researchers, such as Fidani (2010), Eftaxias (2010), Contoyiannis et al. (2010), Plastino et  2010), have attempted to analyze the phenomenon with the purpose of discriminating possible earthquake precursors (earthquake lights, electromagnetic anomalies, uranium groundwater anomalies, ionospheric precursors) that could have foreseen the seismic event.Papadopoulos et al. (2010), in particular, used the available earthquake catalogue extending from 1 January 2006 to 30 June 2009 to detect significant changes before and after the mainshock in the seismicity rate.According to their analysis, in the last 10 days before the mainshock, strong foreshock signals became evident in space (dense epicentre concentration in the hanging-wall of the Paganica fault), in time and in size.
Tsolis and Xenos (2010) use the Cross Correlation analysis method (Pulinets et al., 2004) with their Empirical Mode Decomposition method to analyze foF2 signals collected from three ionospheric stations, in order to verify the existence of seismo-ionospheric precursors before the earthquake.Their work shows that precursors may appear as early as 22 days prior to the event, and 2 days before the mainshock.
The authors have attempted to contribute to this topic by making an analysis of the ionospheric propagation delay, where several GPS reference stations located within 100 kms of the earthquake epicentre have been considered.The main purpose of this analysis was to analyze whether the ionospheric delay computed with a "one way" geometry-free approach, using a Kalman filter procedure to model the daily variability, which may be used as an additional earthquake precursor.

Tests
Data from GPS Continuous Operating Reference Stations (CORSs) placed near the earthquake epicentre have been used to study the trend of the ionospheric propagation delay near the L'Aquila seismic event.The CORSs that were considered are reported in Table 1 and in Fig. 5.These stations were chosen according to the earthquake preparation zone radius with respect to the distance from the epicentre, as described in Table 2 (Dobrovolsky et al., 1979).
GPS data from 15 March (Julian Day 74/2009) to 17 April (JD 107/2009) (15 days before and 15 days after the earthquake) were considered for each station.Concerning L'Aquila reference station, close to the mainshock epicentre, a larger time window, from between 9 February (JD 40/2009) and 17 April (JD 107/2009), was considered."One-way" positioning for each station and every day were performed, considering three different satellites (PRNs 8, 11 and 29) of the GPS constellation.The choice of these three satellites was necessary to monitor the behaviour of the ionospheric propagation delay during three different parts of the day (morning, afternoon and night), as shown in Fig. 6.For this reason, the highest satellite for each part of the day was considered.The figure indicates the position of the satellites on the day of the greatest earthquake event (JD 96/2009), but it can be considered roughly equivalent for the other days of the analysis, as the temporal displacement of the satellites is equal to 4 min per day.The three satellites used in this analysis are highlighted in Fig. 6.
The "one-way" positioning in Eq. ( 4) has been implemented in a Kalman filter, using a dedicated software developed in FORTRAN language.
The results were analyzed by means of some toolboxes that have been developed in MATLAB ® language (Fig. 7).The aim of these software is to provide an early-warning tool prior to the earthquake phenomena.Broadcasted ephemeris and observation data from RINEX files with a sampling interval equal to 30 s were used.The results of the performed analysis are shown in the next section.

Results
The ionospheric propagation delay for each reference station using the "one-way" positioning described in Sect. 3 considering the three GPS satellites PRNs 8, 11 and 29 was computed.The AQUI CORS (L'Aquila) was analyzed for a longer time period than the one mentioned in the previous section, in order to have a longer and more significant time series.
A time window from between 9 February (JD 40/2009) and 17 April (JD 107/2009) was considered.The trends of the daily average value of the ionospheric propagation delay for the L'Aquila GPS station are shown in the following figures for each chosen satellite.The week before and the week after the earthquake of 6 April (JD 96/2009) are drawn in red.
Analyzing the figures above, it is possible to identify an anomalous variation in the ionospheric propagation delay patterns about a week before the seismic event (JD 88/2009), especially for satellite PRN 29 (see Fig. 8).The ionospheric propagation delays computed for satellites PRNs 11 and 8  suffer from noise due to solar activity during the morning and the afternoon, and therefore do not allow the possible disturbance caused by the earthquake to be detected.The median and RMS values computed for the AQUI station, considering the period before and after the JD 88/2009, are described in Table 3.
A comparison with the average values of the ionospheric propagation delay computed for the other GPS CORSs provides additional information.Figure 9 shows the time series computed for the Castel Del Monte GPS station (CDRA) in the JDs 74 to 107 range.The median and RMS values computed for each station, considering satellite PRN 29, are summarized in Table 4.
Analyzing the results reported in Table 4, it can be seen that the RMS of the stations can be used even more than the median value as a parameter of interest for the detection of an earthquake precursor.In fact, the RMSs obtained for the L'Aquila station, located a short distance from the epicentre of the earthquake, are not comparable with those obtained for the other ones.
In general, the noise increase is not only seen for the AQUI station, but also for all the stations in the proximity of the  Lophaven et al., 2002).
The spatial variations of the ionospheric propagation delay, with respect to the median value of the period before the JD 88/2009, are shown in Fig. 10.
The behaviour of the ionospheric propagation delay variation on the JDs 88/2009 and 94/2009 are shown in the following figures.

Conclusions and future tests
The developed tests and the relative analysis demonstrate that the ionospheric delay could be an additional precursor of earthquakes.
As demonstrate above, before L'Aquila's earthquake the study of the ionospheric anomalies could be significant, especially around 7 days before the event.Ionospheric variations and RMS values have been effective precursors of the dramatic event of 6 April.A few days before the event, strange variations were observed around the epicentre.
The results are consistent with the previous work of Tsolis and Xenos (2010), showing an anomalous behaviour of the ionosphere one week before the mainshock.The possibility to use GPS and GLONASS CORSs is really interesting, since doing so allows densifying the number of observations that can be analysed in order to obtain a complete tomography of the ionosphere.
The earthquake preparation zone has been correctly detected as an area with a radius equal to 100 kms with respect to the epicentre.The precursor's characteristics are different in different seismically active regions.These regional peculiarities should therefore be studied.
The GPS ionospheric propagation delay analysis can be considered an appreciable precursor but its limits can be observed when the value of magnitude is not so high (<5.5).The tremors that occurred before and after the main tremors are not detected using this approach.
The last conclusion concerns the precursor variation with the season and solar cyclephase.It is well known that there are strong variations in the ionosphere parameters according to the season and solar cycle.It is quite possible that ionosphere sensitivity to the seismic events and the overall precursor characteristics may also change according to the season and the cycle phase.The seasonal dependencies could not only be due to ionospheric variations, but also to changes in the weather conditions.Winds, rain, snow, and fog may all contribute to the electric field generation mechanism.In order to answer this question, it is necessary to conduct further investigations, to make other tests and consider other factors.

Fig. 4 .
Fig. 4. Map of Abruzzo (Italy) and surrounding regions showing the locations of the epicentres of the earthquakes (M ≥ 4.0) from 6 April to 13 April 2009.The black star indicates the main event (M = 6.3) of 6 April 2009.

Fig. 5 .
Fig. 5. Map of the GPS CORSs in the middle of Italy.

Table 1 .
Coordinates of the GPS CORSs.

Table 2 .
Earthquake preparation zone radius for different magnitudes.

Table 3 .
Values of ionospheric delay in AQUI, expressed in meters.

Table 4 .
The ionospheric propagation delay for satellite PRN29 expressed in meters.
earthquake epicentre.A representation of this phenomenon can be obtained interpolating the values computed daily for each station (e.g. through a simple Kriging, using a MATLAB® toolbox developed by