Statistical discrimination of global post-seismic ionosphere effects under geomagnetic quiet and storm conditions

ABSTRACT The retrospective statistical analysis of total electron content (TEC) is carried out using global ionospheric maps (GIM) for 1999–2015. TEC anomalies are analysed for 2670 earthquakes (EQ) from M6.0 to M10.0 classified into 2205 ‘non-storm’ EQs and 465 ‘storm’ EQs during geomagnetic storms. The geomagnetic storms are specified by relevant thresholds of geomagnetic indices AE, aa, ap, ap(τ) and Dst. Using sliding-window statistical analysis, moving daily–hourly TEC median μ for 15 preceding days with estimated variance bounds is obtained for each grid pixel of GIM-TEC maps. The derived ionosphere variability index, Vσ, is expressed in terms of ΔTEC deviation from the median normalized by the standard deviation σ. Vσ index segmentation is introduced specifying TEC anomaly if an instant TEC is outside the bound of μ ± 1σ. Efficiency of EQ impact on the ionosphere (Eσ) is growing with EQ magnitude and depth representing relative density of TEC anomalies within area of 1000 km radius around EQ hypocentre. Positive TEC ‘storm’ anomalies are twice as much as those of non-storm values. This observation supports dominant post-EQ TEC enhancement with Eσ peak decreasing during 12 h for daytime but growing by nighttime during 6 h after EQ followed by gradual recovery afterwards.


Introduction
The effects of earthquakes in the ionosphere are subject of intense studies during recent decades (Davies & Baker 1965;Koshevaya et al. 1997;Liu et al. 2006a, Liu et al. 2006a, 2006bHarrison et al. 2010;Hayakawa & Hobara 2010;Lin 2010;Arikan et al. 2012;Lin 2012;Astafyeva et al. 2013;Komjathy et al. 2013;Pohunkov et al. 2013;Devi et al. 2014;Perevalova et al. 2014). The diversity of pre-earthquake phenomena, such as local magnetic field variations, electromagnetic emissions at the different frequency ranges, excess radon emanation from the ground, changes in water chemistry, water condensation in the atmosphere leading to haze, fog or clouds, and atmospheric gravity waves rising up to the ionosphere, induces changes in the ionospheric total electron content (TEC) and the F2 layer peak electron density (Pulinets et al. 2003;Chen et al. 2004;Pulinets & Boyarchuk 2004;Rishbeth 2006;Liu et al. 2006aLiu et al. , 2006bDepueva et al. 2007;Varotsos et al., 2008Varotsos et al., , 2011Karatay et al. 2010;Le et al. 2011;Namgaladze et al. 2012;Freund 2013;Devi et al. 2014;Akhoondzadeh 2015;Heki & Enomoto 2015). Changes in magnetic field at the time of the earthquakes have been observed and reported in various publications such as Johnston et al. (1981), Yen et al. (2004) and Varotsos et al. (2009). Modification of the electric field and currents due to electric processes in the lithosphere and the lower atmosphere (Varotsos & 1999, to December, 2015 to the availability of GIM-TEC maps and results of analysis are provided in Section 4. The goal of this study is to obtain new evidence on seismic-ionospheric associations which are summarized in the Conclusions in Section 5.

The statistical analysis of TEC data
In this study, statistical data analysis is performed using global ionospheric maps (GIM) of the TEC provided by Jet Propulsion Laboratory (JPL). TEC is defined as the line integral of plasma density in the Earth's atmosphere and it provides an estimate of the total number of free electrons inside a cylinder with 1 m 2 cross-section area in the column from the bottom of the ionosphere (65 km) to the GPS orbit of 20,200 km. The TEC is an important observable in analysis of temporal variability of the ionosphere and the plasmasphere both under quiet and under storm conditions.
The GIM-TEC maps are generated in a continuous operational way by several Data Analysis Centers since 1998, covering the period more than the entire solar cycle (Hernandez-Pajares et al. 2009). The vertical TEC is modelled by JPL in a solar-geomagnetic reference frame using bi-cubic splines on a spherical grid; a Kalman filter is used to solve simultaneously for instrumental biases and vertical TEC on the grid as stochastic parameters (Manucci et al. 1998). GIM-TEC have been initially provided with 2 h time resolution which are linearly interpolated in time to 1 h resolution, and the hourly files are provided by JPL since December 2008. The JPL maps are generated in the denser map grids (¡90:2:90 in latitude, -180:2:180 in longitude), and time specified for 0.5:1.0:23.5 h UT so these maps are preprocessed by linear interpolation into standard IONEX format for 0:1:23 h UT. The IONEX global map consists of 5183 grid values binned in 87.5 S to 87.5 N in step of 2.5 in latitude, 180 W to 180 E in step of 5 in longitude. The similar structure of map grids is applied when the source GIM-TEC maps are converted to geomagnetic coordinates binned in -87.5 N to 87.5 N in steps of 2.5 in geomagnetic latitudes, and 0 E to 360 E in steps of 5 in geomagnetic longitude using the International Geomagnetic Reference Field (IGRF) model.
According to Liu et al. (2006a), the recurrence time of the M 5.0 earthquakes is 14.2 days. Therefore, in order to determine the reference background TEC distribution, we compute the sliding median of every successive 15 days of TEC at each grid point of the map. In the present study, we use TEC sliding median defined by a 15-day moving window, and the median value is assigned to the final day of the window, i.e. to the 15th day of the window. We use such type of 'forward' median approach because it has a potential for development of forecasting model similar to those in Gulyaeva et al. (2013);Muchtarov et al. (2013): m m;d s ¡ d i ðlÞ D medianðx d i ðlÞ:::x d s ðlÞÞ: (1) In the above equation, x denotes the GIM-TEC value at grid point l. d i and d s represent the first day and the final day of the sliding window, respectively. The subscript m indicates the map under investigation. Statistical study of an ionospheric parameter includes determination of median and dispersion, i.e. variability of the parameter around its median. The standard deviation s represents a measure of the dispersion of distribution which can be computed as where N T denotes total number of days in the sliding window which is set to 15 in this case. An interval within one standard deviation around the median accounts for approximately 68% of the dataset, while two and three standard deviations account 95% and 99.7%, respectively.
The measure of TEC variability is further investigated as the TEC deviation (DTEC) from the median m, normalized by the standard deviation s for N t number of days prior to and during day d s : The algorithm is completed by introducing Ds segmentation with thresholds shown in Table 1 to result in the ionosphere variability Vs index with magnitudes from Vsn D -4 (extreme negative TEC anomaly) in step of DVs D 1 to Vsp D C4 (extreme positive TEC anomaly). Vs index represents the integer magnitude of TEC variability regarding quiet reference median in terms of s grades. Here, the ionosphere quiet state is within DTEC < §1s. If the value of instant TEC is outside of pre-defined bounds of m § 1s, the anomaly of TEC is detected (Akhoondzadeh 2015).
The Vs grade segmentation -4:1:4 is similar to the ionospheric weather W-index (Gulyaeva et al. 2013), yet it differs from W-index by the dynamic thresholds expressed through the variable standard deviation, s. An advantage of Vs index is that it is more physically justified showing DTEC in terms of the relevant standard deviation. This scenario can be easily implemented with any physical parameter, such as the ionospheric critical frequency, foF2, the peak electron density, NmF2, and the peak height, hmF2, using the relevant reference value (mean or median) and the standard deviation for the selected parameter.
Efficiency (Es) of the ionosphere response to impact of earthquakes is represented by the relative density of the extreme negative indices Vsn -2 (mVsn) on the specified fragments of a map corresponding to decreased density of electrons DTEC -1s as compared with quiet reference state; or similarly, the extreme positive indices Vsp 2 (mVsp) on the selected fragments corresponding to increased density of electrons DTEC C1s, to the total number of cells in the fragment(s) around the EQ hypocentre(s) on the map or series of EQs on the relevant maps (mtot): In this study, the available GIM-TEC maps from 1999 to 2015 have been processed to produce output global maps of the 15-day sliding median m, standard deviation s, magnitudes Ds and Vs index in IONEX format with spatial resolution of 2.5 £ 5 , in latitude and longitude, respectively. The histogram of annual frequency of occurrence of the specified magnitudes of Ds, in per cent, relative to the total number of about 45 £ 10 6 grid elements per year (5183 grids £ 24 h £ DOY, with days-of-year, DOY, equal to 365 or 366) is plotted in Figure 1 in increments of 0.5s. We note the asymmetry of the TEC enhancement (Ds > 0) and depletion (Ds < 0) occurrence. The sign of Ds depicts DTEC for an instant TEC being either greater than or less than the quiet reference median. An appreciable number of 'quiet' TEC with Ds D 0 is also seen in Figure 1 when TEC is equal to the median value (Equation 3). There are negligible year-to-year (solar cycle) changes in Ds . We note also a certain percentage of Ds occurrences exceeding §1s which are denoted as TEC 'anomalies'.
The TEC data are extracted from GIM-TEC for the regions surrounding earthquake hypocentre at geographic latitude, ue, and geographic longitude, fe, within the radius of 1000 km, determined by u ue § 10 , f fe § 7.5 . The analysis of TEC for a rectangular region defined by (u i , f i ) to (u s , f s ) is provided by Gulyaeva et al. (2013, Appendix A) for the increments in u and in f given as Du and Df, respectively. However, the space around the EQ hypocentre with radius of 1000 km is not a simple rectangular region. It is rather represented by fragments of 24 cells comprised of a square of 4 £ 4 latitude/longitude grids and a rectangle of 8 £ 2 latitude/longitude grids surrounding (ue, fe) as illustrated in Figure 2.  Global instantaneous Vs map in geomagnetic coordinates frame for 1 January 2012, at 06:00 UT is presented in Figure 2 as an example of global variability index distribution. The time for the map, 06:00 UT, is an integer hour just after the Japan's Izu Islands earthquake on 1 January 2012, at 05:27:54 UT (14:40:27 LT) with Mw D 7.0, and at a depth of 348 km (Lin 2012). The hypocentre of the earthquake was at [31.4 N, 138.2 E] in geographic coordinates, and [22.8 N, 208.1 E] in geomagnetic coordinates which is close to the crest of the equatorial ionization anomaly (EIA). The area selected for the analysis is designated by white points on IONEX grids surrounding the earthquake hypocentre (white star) in Figure 2. This earthquake occurred under quiet geomagnetic conditions (see the next Section for the classification criteria) nevertheless there is appreciable negative Vs anomaly southwards of the hypocentre which is detected earlier by Lin (2012) with the nonlinear principal component analysis (NLPCA) while the principal component analysis (PCA) was unable to detect the anomaly. The evaluation of the global distribution of earthquakes of M6.0C under quiet conditions and the geomagnetic storms in the selected fragments of the globe surrounding the EQ hypocentre is provided in the next section.

Spatial distribution of earthquakes under quiet and storm conditions
The aim of the present study is to reveal a novel empirical evidence of the earthquake related TEC anomalies under quiet space weather conditions and geomagnetic storms. We use earthquake data from the global Catalogue of the Advanced National Seismic System (ANSS) provided by the Northern California Earthquake Data Center (NCEDC 2014). The composite Catalogue of earthquakes created by ANSS is a world-wide earthquake catalogue which is generated by merging the master earthquake catalogues from contributing ANSS member institutions and then, removing duplicate events, or non-unique solutions for the same event. We use the monthly and annual data for earthquakes of magnitude M6.0 to M10.0 from the NCEDS Catalogue for a period from January 1999, to December 2015, according to the availability of the hourly GIM-TEC maps during the solar cycles 23 and 24.
Comparison of earthquakes with the equatorial ring current disturbances has shown that the earthquakes occurred during the Disturbance Storm Time (Dst) storms comprise 13% of the total number of more than 79,000 earthquakes M5.0C for 1964-2013 (Gulyaeva 2014). While the severity of a geomagnetic storm is defined by the Dst index which serves as a standard measure of the energy transfer from the solar wind to the ring current within the magnetosphere (Sugiura 1963), there are also other geomagnetic indices specifying impact on the ionosphere under quiet or disturbed conditions (Deminov et al. 2013). In this study, the EQ series and relevant Vs quantities on a map are referred to the 'storm' conditions (Gonzalez et al. 1994) if at least one of the following criteria is satisfied assuming the 'non-storm' conditions otherwise: AE max 500 nT; aa max > 45 nT; ap max > 30 nT; apðtÞ > 18 nT; Dst min ¡ 30 nT: The above conditions should be fulfilled both for the nearest UT hour (or 3 h UT interval) following the earthquake and the nearest pre-earthquake hour (or 3 h UT interval) to capture storm or sub-storm impact at the time of EQ event. Here AE max is the auroral electrojet AE value for two near-EQ hours, aa max is the mid-latitude aa index value for a given and preceding 3 h intervals; ap max is the maximum ap value for a given and preceding 3-h intervals. The ap(t) is the mean weighted value of ap index (Wrenn 1987): with the characteristic time T D 11 h or t D exp(¡3/T) % 0.76; ap 0 , ap ¡1 ,… are ap values at a given time of EQ and preceding 3 h intervals. The Dst min is the minimum disturbance storm time value for 2 h near EQ time. All the above indices are expressed in nanoTeslas (nT). The periods of storms and sub-storms are included by Equation (5) which may occur at all latitudes from the pole (AE index) through the mid-latitudes (aa, ap and ap(t) indices) to equator (Dst index). The global effect is confirmed by the correlation found between the variation in two independent processes occurring at widely separated regions in space, namely, the ring current intensity and the behaviour of ionospheric densities at high latitudes (Yadav & Pallamraju 2015).
From the total number of 2670 earthquakes of M6.0C during 1999-2015, we have found the majority of events happened under quiet geomagnetic conditions (2205 'non-storm' earthquakes) and 465 'storm' earthquakes (17.4% of the total events list). We note that the per cent of M6.0C storm-time earthquakes for a period of observation during 17 recent years exceeds the storm-time percentage (13%) of M5.0C earthquakes for 50 years of observation (Gulyaeva 2014) due to the extended criteria for the 'storm' classification (Equation 5) than the former specification of the geomagnetic state according only to the ring current Dst storm occurrence.
The global spatial distribution of earthquakes is irregular tending to denser earthquake occurrence in the Pacific region (Levin & Sasorova 2012;Gulyaeva 2014). In the present study, we have estimated the spatial percentage distribution of the 'non-storm' earthquakes M6.0C under quiet magnetosphere (Figure 3(a)) and the 'storm' earthquakes ( Figure 3(b)) for 1999-2015. The 'nonstorm' earthquakes distribution (Figure 3(a)) remind that of M5.0C earthquake zones of enhanced seismic activity (Gulyaeva 2014) which are observed along the tectonic plates boundaries at longitudes from 90 to 190 E and magnetic latitudes from 40 S to 40 N, with dominant earthquake occurrence in the sub-equatorial region of the South magnetic hemisphere. The next appreciable zone of enhanced tectonic activity is revealed around the West coast of South America which also corresponds to a tectonic plate boundary. We note that most of the earthquakes are located within the limits of the closed magnetic field lines, which corresponds to L D 4.17 at the magnetic equator for GPS orbit (Lee et al. 2013) so the TEC variability within the low latitude and middle latitude regions represents the area for the co-seismic and post-seismic ionospheric and plasmaspheric effects.

Results
Temporal-latitudinal graphs of TEC (upper panels) and Vs (lower panels) during three days at the meridian of 85 E are intended to illustrate difference between 'storm' type and 'non-storm' states of the ionosphere parameters under consideration ( Figure 5). Figure 5  An erosion and dissipation of TEC EIA is observed during the geomagnetic super-storm on 18 March (Figure 5(a)) which is normally represented by a two-humps-like latitudinal shape with two peaks at the crests of EIA at about §15 in magnetic latitude with a minimum at magnetic equator which is observed in Figure 5 Though the Nepal EQ happened on the quiet day, the peak TEC at the South crest of EIA (the magnetic conjugate region for Nepal EQ hypocentre area) has been diminished, presumably, due to the EQ impact through the ionosphere conjugation. In particular, TEC at the South EIA peak is decreased from 102 TECU on 24 April to 62 TECU on 25 April and further decreased to 47 TECU on 26 April, i.e. day-to-day TEC depletion is observed after the EQ. More drastic differences between the 'storm' and 'non-storm' co-seismic ionosphere are observed with Vs maps in Figure 5(c,d). In particular, most of the Vs values on map are indicators of positive and negative TEC anomalies for We proceed to statistical evaluation of the Vs signatures under the 'storm' and 'non-storm' conditions in the region of interest. Table 2 presents efficiency Es, in per cent (Equation 4) of EQ impact on TEC anomalies at the nearest integer UT hour after the EQ in several ranges of EQ magnitudes from M6.0 to M10.0 in step of DM D 0.5 M units except for the greatest EQ magnitudes M 8.0 to M10.0. Overall Es value for each subset is also provided in the last row of Table 2.
As can be seen in Table 2, the efficiency of EQ impact on TEC anomalies increases as EQ magnitude, M, gets larger for the both negative Vsn occurrence and positive Vsp occurrence around the EQ hypocentres. The total energy emitted by an earthquake (E, in Joules) (Gutenberg & Richter 1956) is in exponential relation with the magnitude (M) represented by the equation: log E D 1.5M C 4.8 which is applied in the present study for calculation of EQ emitted energy for individual EQ events. The mean energy and the standard deviation are provided for each subset in Tables 2 and 4. The increasing efficiency of the EQ impact on TEC anomalies in terms of M (Table 2) is coherent with the amount of energy allocated during an earthquake (Bath 1956;Levin & Sasorova 2012;Swedan 2015) which gets larger with increasing M as presented in Table 2. This result supports numerous studies on seismic-ionospheric associations because it presents straightforward evidence on dependence of co-seismic TEC variability on amount of EQ energy. The EQ allocated energy is the primary reason for Es dependence on M in our results because all EQs of any magnitude (M6.0C) in either subset are analysed with the same algorithm using the derived Vs index in the vicinity of hypocentre under specified level of geomagnetic activity. Also, it follows from Table 2 that the efficiency of positive TEC anomalies is greater than the negative ones which testifies on the dominant EQ-related plasma density enhancements as compared with its depletion. The stormtime efficiency is larger than the non-storm results which bring the evidence that the ionosphericgeomagnetic storms facilitate TEC enhancements or depletion induced by EQs.
We specify Vs results for daytime earthquakes (the solar zenith angle x < 90 ) and nighttime conditions (the solar zenith angle x > 90 ) during 12 h after EQ for the both 'storm' and 'nonstorm' classes. The time variation of efficiency Es (Equation (4)) after EQ is provided in Figure 6 for daytime, nighttime and total diurnal variation. Symbol SC in the plots stands for the positive 'storm' Vsp, QC for quiet 'non-storm' Vsp, S-for the negative 'storm' Vsn, and Q-for the negative quiet Vsn. Points on the 'Total' subplot curves at 0 h are those values that are listed in the last row of Table 2.
In general, all statistical results for the quiet and storm conditions confirm existence of seismic -ionospheric associations since the efficiency of EQ impact on TEC anomalies is not zero in all cases. For some individual EQs, the TEC anomalies in the sense defined in the present study could be missed in the EQ predefined area within 1000 km radius from the hypocentre but we should keep in mind that the most notable ionosphere variability anomalies are specific for the high latitudes while the EQs regions of occurrence belong to the middle and low latitudes. The most important outcome of results in Figure 6 is that efficiency of EQs on positive Table 2. Efficiency of ionosphere response, Es, %, for the TEC enhancement and depletion, Vs, at the different ranges of EQ magnitude (M), at the nearest integer hour (UT) after EQ. The mean energy (J) and standard deviation std for the earthquakes number m are given for each collection during 1999-2015. TEC anomalies under storm condition is twice as large as those under non-storm anomalies. Peak of Es for storm Vsp occurs by daytime (and total diurnal variation) at the nearest (t D 0) hour after EQ. The value decreases after the EQ to the level of other cases as indicated in Figure 6. When compared with the daytime, the results for nighttime storm Vsp anomalies show an enlargement peak by 6 h after the EQ with a value which is twice as large as the other levels and it decreases during the 6 h after the peak. The mean curves of efficiency of EQ impact on TEC anomalies ( Figure 6) are accompanied by Table 3 depicting the ANOVA (Analysis of Variance) statistical results for Vp and Vn occurrence under quiet and disturbed geomagnetic conditions for post-earthquake hours within the 1000 km radius around the hypocentre. Here F implies Fisher's criteria, and p is the probability of the result assuming the null hypothesis. Analysis of variance (ANOVA) is a collection of statistical models used to analyse the differences among group means and their associated procedures (such as 'variation' among and between groups). In the ANOVA setting, the observed variance in a particular variable is partitioned into components attributable to different sources of variation. In its simplest form, ANOVA provides a statistical test of whether or not the means of several groups are equal, and therefore generalizes the t-test to more than two groups. ANOVA is applied here for comparing Figure 6. Efficiency of the seismic impact on the ionosphere for 12 h after earthquakes with Vs index anomalies for nighttime, daytime and the total data-set under quiet conditions and during the geomagnetic storms.  Table 4 show that the selected algorithm of Vp and Vn estimates is meaningful according to the variables.
To determine the dependence of ionosphere variability on the depth of the EQ hypocentre, the relations of the different magnitudes of EQs with their depth are evaluated. The EQs occurrence for the different ranges of the hypocentre depth in the Pacific region is provided in detail by Levin and Sasorova (2012). The results of evaluation of the earthquake energy and standard deviation for three categories of depths for daytime and nighttime under geomagnetic quiet and storm conditions for 1999-2015 are given in Table 4. Hypocentre depth, D, is grouped into three classes: the shallow depth, D1 70 km; the descent depth, 70 < D2 300 km; the deep depth, 300 < D3 800 km. The occurrence of EQs decreases with increasing depth both for geomagnetic quiet conditions and storms. While the magnitude M is introduced by Gutenberg and Richter (1956) as a measure of energy emitted by EQ, the specification of energy distribution in terms of the depth categories shows the dependence of EQs energy on depth so that the energy of EQs gets larger as the depth increases.

Conclusion
In this study, the structural changes of ionosphere are investigated with respect to disturbances in the ionization levels and geomagnetic field due to storms and earthquakes using a novel Vs index, which is derived using the variability of GIM-TEC. The seismic-ionospheric associations are analysed during 12 h after each of 2670 earthquakes of Richter magnitude from M6.0 to M10.0 separated to 'storm' class of 465 EQs and 'quiet' or 'non-storm' class of 2205 EQs worldwide from January 1999 to December 2015. The median, m, of 15 days prior to the current day at each cell of GIM-TEC map in 2.5 £ 5 of latitude / longitude grids is computed for each hour UT (0, 1, …, 23 h) as a reference value. The standard deviation s from the median represents a measure of the dispersion of distribution. The deviation of instant TEC from the median normalized by the standard deviation, Ds, is converted into an index, Vs, varying from ¡4 to C4, that corresponds to extreme negative or positive deviations, respectively.
Efficiency (Es) of the ionosphere response to impact of earthquakes is estimated as a relative density of the negative indices Vsn ¡2 on the specified fragments of a map (DTEC -1s), or the positive indices Vsp 2 (DTEC C1s), regarding the total number of cells in the fragment(s) of 1000 km radius around the EQ hypocentre(s) on the map or series of EQs on the relevant maps.
It is found that the efficiency of EQ impact on the ionosphere is growing with EQ magnitude M at the nearest integer hour UT after EQ both for the storm and non-storm classes. The positive TEC anomalies are more effective than the negative ones for both storm and non-storm subsets which indicate on the EQ post-effects producing rather increased plasma variability in the ionosphere than its decreasing process.
The Vs values grouped with respect to storm-time earthquakes and quiet-time earthquakes for nighttime (solar zenith angle x > 90 ) and daytime (x < 90 ) occurrences during 12 h after EQ show that post-seismic TEC positive anomalies occur almost twice as much as compared to the negative anomalies under storm conditions. Twice as many positive TEC anomalies during geomagnetic storm in the near-hypocentre region are observed at the first integer hour in UT after EQ with a subsequent decrease during 12 h afterwards for daytime. The increase of TEC positive anomalies by nighttime is observed during 6 h after EQ followed by a gradual recovery after the peak.
Analysis of the EQs energy for three classes of the depth (D 70 km, 70:300 km, 300:800 km) brought an evidence of its dependence on the depth of the tectonic events. While the magnitude M is introduced by Gutenberg and Richter (1956) as a measure of energy E emitted by EQ, M»M(E), the specification of energy distribution in terms of the depth categories shows the energy of EQs growing with the greater depth D, in other words, the EQ magnitude should be represented in a function of two variables: M»M(E,D).
The present results suggest that there is a challenge for more sophisticated techniques to be developed in order to distinguish the earthquake effects on the ionosphere happened on the background of geomagnetic activity. The results of this study will be used as a basis for observing and grouping the disturbances in the ionosphere and geomagnetic field and Vs index can be developed further as a storm and/or earthquake precursor.