1. Observations
Kp index is a complex average of the measurements by 13 magnetic observatories around the world. Kp records of almost one hundred years are available. The index goes from 0 to 9 in 1/3 increment reported 8 times a day in 3 hour intervals. The Kp index quantifies the global deviations of the magnetic field from the usual diurnal variations. Here we report on three major seismic events of M ≥ 6 in 2016 and 2017. In all three cases, as Figure 1 attests, the earthquakes struck at times of strong global magnetic disturbances, during periods of elevated or high Kp values. Figure 1(a) shows that the Kumamoto earthquake in Japan on April 14 2016 coincides with the surge of Kp during a SSC (Sudden Storm Commencement), which had started on 2016 04 14 07:35 UT. Figure 1(b) shows the L’Aquila, Italy, earthquake of October 26 2016, coincides with the broad Kp surge with a SSC on 2016 10 25 09:22. Figure 1(c) shows the M8.1 earthquake in Chiapas, Mexico, on September 8 2017, which struck in the subduction zone, at the onset of a magnetic storm caused by the strongest coronal mass ejection
Figure 1. The Kp index variations in the period of the ±28 days around three earthquakes, (a) Kumamoto Japan April 14 2016 M = 6.2, (b) L’Aquila Italy October 26 2016 M = 6.1, (c) Chiapas Mexico September 8 in 2017 M = 8.1, note: top caption means “4.8 hour (04:50) 15.0 latitude, −93.9 longitude, 56.7 km depth”.
(CME) in a decade. In this case two SSCs were recorded at 2017 09 06 23:44 and 09 07 23:00.
Figure 2 shows the plot of the geomagnetic data relative to the Kumamoto earthquake in Figure 1(a), recorded at the geomagnetic observatory in Kanoya, Japan, over four days prior to the event. The storm intensity is shown in detail. It took a sharp downward turn about 2 hrs before the quake. It then turned up rapidly, marked as “First hit” in Figure 2(d), coinciding with the onset of the Kumamoto earthquake on April 14 2016 at 12:26 UT.
Inspired by this coincidence as well as similar cases such as in the 1995 Great Hanshin earthquake of January 17 1995, also known as the Kobe earthquake, and the 2011 Tohoku Earthquake, we analyzed the seismic data of thousands of earthquakes of M ≥ 6 in relation to the variations of the Kp index from 1932 to 2016 (USGS Earthquake Catalogue [1] GFZ Potsdam WDC data, Kp index [2] ).
2. Methods and Results
To visually compare the Kp indices and earthquakes, we stacked intervals of the
Figure 2. (a), (b), (c), (d) for 11, 12, 13, 14 April 2016 The diurnal variations of the geomagnetic fields for the last 4 days before the Kumamoto earthquake on April 14 2016, observed at the Kanoya Observatory, and catalogued at Kakioka Magnetic Observatory data [3] , H, Z, F: horizontal, vertical and total component, respectively, D: magnetic declination (angle of H versus geographic north, east is positive); the magnetic storm starts with a SSC pulse, “First hit” marks the Kumamoto earthquake time.
Kp variations over ±28 days before and after earthquakes. Figure 3(a) (averaged among the earthquakes at each time instance) plot the Kp variations, centered to the time of the earthquake for 4666 events of M ≥ 6, period 1932-2016, in the Pacific Rim region, longitude of 120˚ - 160˚E and 70˚ - 130˚W.
The average Kp in Figure 3(a) has the following features.
・ During the last 12 to 14 days before the earthquakes, “period A” (marked in red), the Kp indices consistently display a broad but large surge of 0.15 (2.15 ~ 2.3), and during the last 2 days before and 2 days after the earthquakes, “period B” (marked in blue), the Kp indices consistently surge by 0.1 units (2.15 ~ 2.25)
・ Higher frequency variations of 0.025 peak to valley
・ The standard deviation of σA = 0.031 over the entire time interval of ±28 days
The plot shows that the double-surge pattern is synchronous to the earthquake onsets. In the following text we call this pattern, “A-B Maxima”. The Kp variations are rather cyclic and consistently present before earthquakes. After the earthquakes, the Kp variations settle back to random characteristics. Notably the
Figure 3. (a) Time Variations of the average of the Kp index for 28 days before and after the earthquakes, 4666 events in the Pacific Rim region. (b) Time Variations of the ap index, (c) Time Variations of the standard deviations of Kp.
Kp variations over the last 12 - 14 days prior to the earthquakes are pronounced and less affected by statistical procedures, as has been verified by changing the number of the earthquakes included in the statistical evaluation.
In addition, since the Kp index is logarithmic, we downloaded the linear ap indices, in nano-Tesla units, from the same observatory [2] and stacked them as shown in Figure 3(b). In both cases the time series are dominated by peaks and valleys that are common to both. In both cases the A-B Maxima patterns stand out above the background. For the statistical analysis, we preferentially use the Kp index because the frequency histograms of the Kp values conform better to a normal distribution, which is a requirement for the subsequent statistical analysis.
The standard deviation of the Kp index is shown in Figure 3(c) (standard deviation among all 4666 earthquakes at each time instance). The mean of the standard deviation is σE = 1.45. Notably this mean is the most critical number which gauges how each Kp index for each earthquake deviates from their average. The figures of Kp plots and the statistical parameters are described in more details in Table 1.
Table 1. Descriptions of the curves in Figure 3(a), Figure 3(c) and the statistical values σA and σE.
Note: Kp values shown in the table are not the actual ones.
3. Discussions
3.1. Sequences of Random Earthquakes―Hypothetical Earthquakes
Figure 4(a) and Figure 4(b) show the Kp variations around the 4,666 earthquakes, based on randomly selected times with random intervals. We thereby create hypothetical earthquakes and look at the Kp signals in the window of ±28 days around the earthquake times. For a visual comparison, the A-B Maxima plot from Figure 3(a) is attached at the bottom of Figure 4. Note that the intervals
Figure 4. (a) and (b) Variations of the average Kp index when the times of earthquake onsets are randomly selected (hypothetical earthquake sequences), comparing with the variations for the actual earthquake sequence, re-plot of Figure 3(a).
were determined so that the total number of hypothetical earthquakes is close to the actual number of earthquakes. Once the earthquake time is randomly determined, the computer brings in the Kp data for that period of ±28 days, whose central time is the earthquake time. Figure 4(a) and Figure 4(b) show clearly a random behavior of Kp, appearing as almost a stochastic stationary time series, while a distinct and large pattern is seen for the real earthquakes. The standard deviations of the time series are 0.0155 and 0.0176 for the sequences of the hypothetical earthquakes, significantly smaller than the standard deviation 0.031 for real earthquakes. The statistical study is described in the subsequent session.
3.2. Statistical Evaluation
First an order estimate is made. The standard deviation of the Kp variation is σE = 1.45 for the Pacific Rim region with n = 4,666 earthquakes. The deviation of the Kp of n earthquakes average is σE/
= 1.45/
= 0.021. The peak-to-valley variation, shown in Figure 3(a) is 0.15. Using the half of the variation 0.075, the cumulative probability at Z factor, Z = 0.075/0.021, is 99.98%, hence this variation is concluded to be outside of the statistical deviation with 99.98% probability.
Next, a conventional statistical study is performed and its result is described in Appendix. The study verifies that the Kp index before the earthquake shows a surge which is statistically significant, implying that the Kp surge is related to and synchronizes with the earthquakes.
For a further analysis, 54 sequences, each consisting of 4,666 hypothetical earthquakes, were generated and a hypothesis test was performed. Hypothesis (null): The Kp variation has nothing to do with the earthquakes occurrence. Hence the sequence of the real earthquakes should show a variation, which is within the variations of the sequences of the hypothetical ones. Figure 5 gives the result. The black solid line is the statistical distribution of σ (standard deviation) of 54 sequences of hypothetical earthquakes (4,666 earthquakes/sequence). The black line is almost same as the red line, which is the normal distribution. The most significant finding is the blue dotted line, σ of the real earthquake sequence, which is clearly separated from the distribution area of the hypothetical earthquake sequences. Thus, the null hypothesis is rejected, and the evidence of Kp surges, relating to the earthquake onsets, is confirmed.
A similar result is shown in Figure 6 for the Mediterranean-Asian region (10˚ - 100˚E, 0˚ - 50˚N). The blue dotted line, showing the σ of the real earthquake sequence, is now closer to the distribution peak of the hypothetical earthquake sequences. Yet, the blue line is still about 2σ away from the peak, indicating 4% ~ 5% probability of no-difference (insignificance). Therefore, we conclude with 4% ~ 5% risk that statistically the real sequence is outside of the hypothetical sequences.
Figure 5. Probability of the standard deviations of the Kp variations of 54 sequences (4,666 hypothetical quakes/sequence). Note: The standard deviation of 0.031 for the real quake sequence in the Pacific Rim region.
Figure 6. Probability of the standard deviations of the Kp variations of 54 sequences (1,200 hypothetical quakes/sequence). Note that the standard deviation is 0.046 for the sequence of the real earthquakes in the Mediterranean-trans-Asiatic zone.
It is to be noted that, if the Kp surge had nothing to do with earthquake onset, the blue dotted line would move to the left and falls on top of the peak area.
3.3. Lack of Precision of Kp Index
Kp index is a historical index of nearly a hundred years of records. Its draw-back is the increments; 0.0, 0.33, 0.67, 1.0, ・・・, 8.0, 8.33, 8.67, 9.0, consisting of only 28 sparse ticks. However, the final Kp variation curve is a stack of the curves of thousands of earthquakes. When we use the effect of
, i.e. central limit theorem, the resolution with the 0.33 step improves to 0.005 for the case of n = 4666, assuring the adequate resolution for detecting the Kp surges. The fine ripples, for example, seen in Figure 3(a), may reflect this sparseness problem.
3.4. Effect of Aftershocks and Mega Earthquakes
Numerous strong aftershocks may falsify the results obtained by the method of stacking the temporal Kp signals. Aftershock activity is solely caused by the main shock and follows this shock by a specific temporal attenuation law. Thus, it cannot be correlated with solar induced Kp values. In order to keep this influence small, we used a seismic data set of larger magnitude earthquakes, M ≥ 6.7, as a threshold which is not often exceeded by aftershocks within a few days after the main shock. We stacked the Kp plots of 1092 earthquakes of M ≥ 6.7 for the entire Earth for the period 1932-2016. The result is shown in Figure 7. A very similar pattern as shown previously in Figure 3(a) appears again.
To prove the statistical significance, the data of thousands of earthquakes are needed. Surprisingly, mega earthquakes of M ≥ 8 show the same A-B Maxima even though only around 50 such earthquakes occurred in the time window considered. These mega earthquakes coincide highly with Kp surges, compared to the M6 level earthquakes.
Figure 7. Kp curve showing the A-B Maxima again, synchronizing with 1092 global earthquakes of magnitude M ≥ 6.7, aftershocks mostly removed.
Figure 8. Kp variations of the single 2011 Tohoku earthquake and Tsunami in Japan.
3.5. The 2011 Tohoku Earthquake and Tsunami
Another significant result was obtained by plotting the Kp sequence for the strongest and most devastating earthquake in the recent decade. Figure 8 shows the ±28 day Kp plot of the Tohoku earthquake, occurred on March 11 2011 at 14:46 JST, 5:46 UTC. The Kp variation is stunningly similar to that previously shown in Figure 3(a), the averaged Kp over many earthquakes in the Pacific Rim region. The implication of the similarity is profound and deserving of the future research that could cast light on the geomagnetic field disturbances being a factor for triggering of mega-earthquakes. The autocorrelation and FFT analysis was also performed. In addition to the dominant 27 ~ 28 day period associated with the Sun rotation, a distinctive 10-day periodicity was found in the Kp variations.
4. Analysis Summary
The distinct patterns found in the Kp fluctuations prior to earthquakes indicate synchronization of the geomagnetic surges and the seismicity. We called this pattern A-B Maxima. Maxima A occurs 8 ~ 10 days before earthquakes and lasts for few days. Maxima B occurs right before or during earthquakes.
An overwhelming statistical significance is obtained from Pacific Rim region data, by applying conventional statistical methods such as t-check, F-check and ANOVA. One factor, requiring a consideration, is that the data points in the time series are usually not independent to each other, which results in less degree of freedom (DOF). Although the reduced DOF is still adequate in proving the statistical significance, a new method is devised, in which the correlation between the Kp surges of hypothetical earthquakes in random time sequences and real earthquake sequences are compared. This comparison confirms that the synchronization is 100% significant in the Pacific Rim region, while other seismic zones provide a weaker, yet significant ~95% correlation. A further research on individual earthquakes may solve a question, whether the A-B Maxima pattern itself effectively triggers earthquakes, or the broad surges trigger earthquakes regardless of the pattern, in such a case the A-B Maxima pattern would just reflect a common pattern of the solar activity.
5. Review of Historical Studies
Aristotle, the Greek philosopher, described that earthquakes occur more frequently during the night than during the day [4] . Starting at the beginning of the 20th century, studies were performed to determine whether the time sequences of earthquakes follow any systematic pattern. After many studies it became evident that the seismicity exhibits distinct diurnal as well as seasonal cycles in many different earthquake zones [5] - [13] .
Such cycles of seismic activity can only be attributed to solar influence. Consequently, studies were done to model such solar-terrestrial effects. The powerful electric current vortices in the ionosphere, generated by solar radiation, and the associated magnetic field variations are assumed to be the essential source for the effect. Due to the penetration of those magnetic field variations into the electrically conductive Earth’s lithosphere, and associated electric “telluric” currents, additional mechanical forces are generated in seismic rupture zones [14] [15] [16] . Another observation that supports the theory of solar influence is the clustering of earthquakes in 11-year cycles in accordance with the solar cycles [17] - [23] .
In this overall context, several recent studies focused on the cycles of the solar polar magnetic field, which manifest themselves in solar flare activity. A high degree of statistical significance was shown for a correlation of those polar field oscillations and the occurrence of strong earthquakes [24] [25] . Similarly, it has been stated that, during the Maunder solar minima, the strongest earthquakes and most violent volcanic eruptions took place during transition phases of the heliospheric magnetic field [26] .
Other studies deal with the question of whether geomagnetic storms play a role in the triggering of earthquakes. These magnetic storms are generated by a sudden surge in the intensities of the solar wind (plasma) streams, which originate on the sun’s surface and travel with high speed toward our planet. Disturbances of the geomagnetic field caused by this solar influence are classified by the so-called 3-hour Kp Index [27] . Several papers confirm this Kp index to be an appropriate indicator for the solar influence, including magnetic storms, on seismicity [28] [29] [30] [31] [32] .
However, another investigation claims to have demonstrated no statistical significance for any solar-terrestrial triggering of earthquakes [33] .
6. Conclusions
The Kp index for the times of earthquakes between 1932-2016 was statistically analyzed. Stacking of thousands of Kp data shows an effect of the geomagnetic field on earthquakes triggering. As described in Analysis Summary, a distinct pattern of the Kp fluctuations prior to earthquakes was found, indicating the synchronization of geomagnetic surges and seismicity. These synchronizations are quite complex, reflecting the regional characteristics and the earthquake magnitude itself. M8 class earthquakes are associated with the Kp surge more than M6 class ones.
The geomagnetic disturbance, typically the magnetic storm, is one of the major factors which synchronize with earthquakes.
This study offers a scientific support to the past numerous researches by the predecessors and the current researchers.
Acknowledgements
We acknowledge a substantial contribution offered by John Scoville during the discussions in the NASA Research Park and express our special thanks to Suzanne Smart, Vancouver Canada, for her devoted editing work. The improvement in wording and phrasing is made by Roger Johnson, an “American” English speaker/researcher, in San Jose California. We thank him for his contribution.
Appendix
Conventional Statistical Assessment
For a further scrutiny, t-test is applied to the stepwise changes in the Kp curves. The step changes were concluded to be significant. Next, One Way ANOVA (Analysis Of Variance) is performed. F number in ANOVA is defined,
,
where SSA is the sum of squared differences of factor A (time factor in our case), i.e. between groups statistically. m is the width of time window. SSE is the sum of squared differences of error (between earthquakes, i.e. inside of each group statistically). σA is sample standard deviation for factor A. In the signal processing world, σA is signal while
is noise. When we use the full span window of ±28 days, m = 28 * 2 * 8 + 1 = 449 (note 8 points/day). For the same Pacific Rim region, σA is 0.031. While,
is 0.021. Thus, F is 2.12. Since F value of 5% significance is F = 1.11 with the large degrees of freedom (DOF), (m = 449, n = 4666), the probability of the Kp variations to be within the statistical variation is none.
Reduced Degree of Freedom
It is known that the use of ANOVA for time series requires an attention to avoid the error because the values at each time are not independent and hence the reduction of the DOF takes places. This reduction is estimated using the relaxation time, derived from the autocorrelation of the time series. We found even with this reduced DOF, the synchronization of the Kp variation is still significant for the Pacific Rim region.