1. Introduction
The possible relation between thunderstorm activity and seismic activity is an old question. A long time ago some rough correlation between earthquakes and lightning strokes has been shown during summer time in the Kwantô district in Japan [1]. Afterwards only case studies have been mainly done. Unusual variations of the intensity of the sferics received at 6 kHz and 9 kHz at Agartala, India (latitude, 23˚N; longitude, 91.4˚E) have been reported by [2] some days before an earthquake occurring on October 8, 2005 at Muzaffarabad, Pakistan (latitude, 34.53˚N; longitude, 73.58˚E). The number of sferics was highest on the day of the earthquake. A sferic (or atmospheric) is the wave emitted by a lightning stroke which propagates in the Earth-ionosphere waveguide. Prior to an M8.2 earthquake occurring on 24 May 2013 near the Kamchatka peninsula, increase of the sferic amplitudes has been observed [3]. It was noticed in [4] that anomalous electromagnetic signals at frequencies of 3 - 60 kHz prior to earthquakes can be due to local thunderstorms. A change of spectral density in high frequency bands has been found in [5] from several days to several hours before an earthquake. Some powerful thunderstorm activities have been observed prior to major earthquakes [6]. Anomalous thunderstorm activities above epicenters have been observed approximately 5 days before an earthquake event [7]. Concerning the M7.2 Hyogoken-Nanbu earthquake of January 17, 1995 in Japan, strong lightning discharges have been detected in association with intense radio noises one week prior to the main shock [8] [9]. It has been noticed an increase in the amplitude of lightning electromagnetic signals (atmospherics) received at Yakutsk, 12 - 14 days before the M9 Tohoku earthquake in Japan [10]. Recently, abnormal thunderstorm activity has been observed over the epicenter area one month before the M8.3 Illapel earthquake in Chile [11]. But several authors using low resolution data were not able to check the nature of the electromagnetic emissions observed prior to earthquakes and even some authors try to remove waves due to lightning strokes considering that they could not be related to earthquakes [12].
A more global analysis has been done in [13] [14] close to Taiwan. They have performed a statistical correlation between lightning activities and M ≥ 5.0 earthquakes occurring during 1993-2004. Their analyses indicate that the emissions enhance 5 - 7 days before the earthquakes above the epicenters and are proportional to the earthquake magnitude.
The purpose of this paper is to examine if there is likely a connection between the occurrence of earthquakes and the ionospheric signature of thunderstorm activity in the epicenter area using the data of the DEMETER satellite. The DEMETER data and the earthquake selection are presented in Section 2. The method used to process these data is described in Section 3 where results are also shown. They are discussed in Section 4 whereas conclusions are given in Section 5. Note that our study is not related to earthquake lights (see for example [15] [16] [17]) because these lights are faint and diffuse and correspond to a very different phenomenon.
2. The Data
2.1. The Whistlers
Thunderstorm activity is very common in the Earth’s atmosphere and each lightning stroke triggers an electromagnetic wave in a huge frequency range which propagates in the Earth-ionosphere waveguide. Depending on the frequency range and the ionospheric conditions, a part of the wave energy can cross the ionosphere and propagates in the magnetosphere along magnetic field lines with a frequency dispersion. It is called a whistler. Onboard the low altitude satellite DEMETER, RNF is a dedicated experiment based on a neural network to count these whistlers and to classify them by dispersion with a time resolution of 0.1024 s. Details about the RNF experiment can be found in [18] [19], and [20]. Here we will consider the whistlers with low frequency dispersion, i.e., the whistlers which corresponds to waves emitted just below the satellite. After a learning phase, RNF was operational from May 2005 until December 2010 the end of the mission. The lightning activity is naturally function of many parameters as the location and the seasons. But a recent study has shown that a map of whistler activity is very different from a map of lightning activity [20]. The main point is that the whistlers with low dispersion are strongly attenuated close to the magnetic equator.
DEMETER is nearly sun synchronous and measurements are only done at two different local times: 10.30 and 22.30 LT. Only night time whistlers recorded by DEMETER have been considered here because they are more numerous than during day time. This is due to the thunderstorm activity which is lower in the morning, and to the whistler intensity absorption since the ionospheric density is higher.
2.2. The Earthquakes
Earthquakes with magnitude M ≥ 4.8 have been selected in the USGS database (http://www.usgs.gov). It represents a database of 18,752 events with time, location, magnitude M and depth d. Considering that the observation of whistlers on board DEMETER dramatically decreases close to the magnetic equator and that the aim of our study is to investigate a possible relation with the earthquake locations, all earthquakes close to the magnetic equator i.e. with |mag. lat.| ≤ 20˚ have been removed (see Figure 1). It remains N = 6847 events. As we have seen in the previous section that the whistlers propagate along the magnetic field lines we do not consider the position of the epicenters at the Earth’s surface but the magnetically conjugate points of these epicenters at the altitude of the satellite (660 km) in the same hemisphere. It means that, either in the Northern hemisphere or in the Southern hemisphere, the epicenter positions are slightly shifted towards the equator to take into account the whistler propagation (see Figure 2). For example, concerning an earthquake with an epicenter located at (39.22˚N, 41.08˚E) the study of the whistler rate will be done at the location (35.28˚N, 40.73˚E). In the following this location will be called Epicenter Conjugate point (ECP).
Another key problem related to the statistics with earthquakes is the occurrence of aftershocks. When we are looking for possible effects a few hours or
Figure 1. Using the USGS database the map shows the location of the earthquakes with M ≥ 4.8 which are considered in the statistical analysis. Earthquakes with magnetic latitudes of the epicenters less than 20˚ have been removed.
Figure 2. The Earthquake Conjugate Point (ECP) is defined as the intersection between the orbit of the satellite and the magnetic field line (red curve) of which the foot point at the Earth’s surface corresponds to the earthquake epicenter.
days before earthquakes, the existence of aftershocks (which usually enter in the statistical analysis) can alter the results. The mean reason is that, at the time of the shocks, it is known that there are perturbations propagating in the atmosphere and then in the ionosphere. The method we use to remove the effects of the aftershocks is to not handle twice the same data. As the data is sequentially processed it means that the data related to a main shock are always considered in comparison with the data of its aftershocks. It must be noted that this method is not really efficient as thousands of earthquakes with magnitude lower than 4.8 are not considered. One must admit the hypothesis that the magnitude of the earthquakes is proportional to the amplitude of the effects we are looking for. Then, earthquakes with smaller magnitudes must have a small influence and can be removed from the statistics. This will be demonstrated later on.
3. The Data Processing
The whistler activity observed onboard DEMETER depends on the location, the time, the season, and also the solar cycle [20]. Then it is necessary to normalize the whistler occurrence for each earthquake in order to have a homogeneous data set. In a first step, each earthquake has been considered separately and all orbits close to the ECP have been searched during a time interval of 70 days around the earthquake time in order to determine a background value of the whistler occurrence. The distance between the orbits and the ECP has been limited to 1000 km, and 5 concentric rings with a 200 km step around the ECP (0 - 200 km, 200 - 400 km, 400 - 600 km, 600 - 800 km, 800 - 1000 km) have been considered to gradually check the whistler occurrence when the satellite moves away from the ECP (Figure 3). This largest distance of 1000 km is in relation with the estimation of the size of earthquake preparation zones [21].
The data processing used to drive the background whistler rate for a given earthquake is described hereafter. Considering all orbits passing less than 1000 km away from the ECP from 35 days before and up to 35 days after the earthquake occurring time, we can build a 70 × 5 matrix filled by the recorded whistler rates for the 70 days and the 5 rings. Depending on the position and time of the earthquake, all cells are not necessarily filled up but they will be used to define a background, i.e., an average whistler rate which is characteristic of the region where the EPC is, and of the season when the earthquake i occurs. It means that for an earthquake i all values in the 70 × 5 cells are averaged to calculate a background whistler rate
. This background whistler rate
is normalized to take into account the surfaces of the rings around the ECP (Figure 3). It is easy to show that, if S is the surface of the 200 km radius circle around the ECP, the 4 successive rings have a surface of 3 S, 5 S, 7 S, and 9 S. Then
is calculated as
(1)
where vk is the averaged whistler rate measured in the ring k during the 70 days. Among the 70 days we select a time interval around the time of the earthquake i between −15 days and +5 days to count the whistler rate wj,k in each cell (j, k) where j varies from 1 to 20 and k from 1 to 5 (20 days and 5 distances). Then we normalized these wj,k of each cell (j, k) by dividing by
Figure 3. Pattern of rings around the ECP. Their radiuses increase by step of 200 km. The red part of the straight line represents the projection of a satellite orbit where data are recorded. In this example they are no data recorded in the area less than 200 km from the ECP.
(2)
Therefore for each earthquake we have a grid with 100 cells. Finally, all N earthquakes are considered as we will use a superposed epoch method to represent the data. When all earthquakes have been processed we have nj,k normalized values in the cell (j, k) and a total of NV normalized values in the grid.
(3)
In a second step we calculate the average quantity µ of these NV normalized values and the corresponding variance σ in order to obtain standardized values in the entire grid.
(4)
In each cell of the grid of the superposed epoch method we plot the average of these standardized values:
(5)
where njk is always less than N. The results are shown in Figure 4 and Figure 5 top for M ≥ 4.8 and d ≤ 20 km, and in Figure 5 bottom for M ≥ 5.4. They are discussed in the next section.
Figure 4. Representation of the superposed epoch method as function of the days around the earthquake time (all earthquakes occur at the time 0) and as function of the distance from the ECP. It concerns earthquakes with M ≥ 4.8 and d ≤ 20 km. On the right, the color scale gives the average of the standardized values.
4. Discussions
The property of the standardized values
given by Equation (4) is that their mathematical expectation is zero, i.e.
(6)
Another property of the mathematical expectation is that if X and Y are two random variables, then the mathematical expectation of the sum of these two variables is equal to the sum of the mathematical expectation of X and the mathematical expectation of Y, i.e., E [X + Y] = E [X] + E [Y]. Then we can write:
(7)
It means that if we consider that the fluctuations of the standardized values must be identical in all cells of the grid the averaged values in all cells must be zero or very close to zero. This is the case in Figure 4. Looking to the amplitudes of the color scales of Figure 4 and Figure 5 top one can see that the values in the grid are very small and oscillating around zero. Obviously this is not the case in Figure 5 bottom where one averaged value departs from zero and is particularly higher in the cell ([0, - 1 day], [0 - 200 km]). This underlines the unusual behavior of the data in this cell.
Figure 6 shows the result of a further analysis which has been done to check the influence of the magnitude M and the depth d of the earthquakes. It represents the average of the standardized values in the cell ([0, - 1 day], [0 - 200 km]) as function of M and for a depth less than 20 km and less than 50 km. One can see that the effect (high whistler rate) is only significant when the magnitude
Figure 5. In order to underline the variation as function of the magnitude, the top panel is identical to Figure 4 but with the same color scale as the bottom panel which is for earthquakes with M ≥ 5.4 and d ≤ 20 km.
is larger than 5.2. The effect is lower when the depth of the earthquakes increases as it is in the analysis shown in [14].
It is possible to quantify the effect of seismic activity on the whistler occurrence. Figure 7 is very similar to Figure 6 but it represents the percentage of whistler rate increase relatively to the background, i.e. the following quantity
(8)
where
is given by Equation (2) and the indexes j and k corresponds to the cell [0, - 1 day] and [0 - 200 km] respectively.
Figure 6. Variation of the amplitude in the cell [0, - 1 day], [0 - 200 km] as function of the earthquake magnitude for a depth ≤ 20 km and a depth ≤ 50 km.
Figure 7. Percentage of whistler rate increase as function of the earthquake magnitude for a depth ≤ 20 km and a depth ≤ 50 km.
One can see that the whistler rate increases by 11% at less than 200 km from the ECP and one day before earthquakes with M ≥ 5.5 and d ≤ 20 km. It is difficult to check if this percentage again increases for earthquakes with larger magnitudes because their number dramatically decreases and the statistic would not be meaningful anymore. It is also shown in Figure 7 that the percentage decreases for deeper earthquakes as expected.
The same data processing has been done for “random” earthquakes, i.e. the time of each earthquake has been shifted forwards 25 days. The result is shown in Figure 8 for “random” earthquakes with M ≥ 5.4 and d ≤ 20 km and it can be compared with Figure 5 bottom. One can see that there is no link with the time of the “random” earthquakes. The reason is that, with our shift of the time of earthquakes, we have suppressed the time relation between whistler observations onboard DEMETER and large magnitude earthquakes.
Why such small increase of whistler observations onboard DEMETER can occur prior to earthquakes? This is due to important variation of parameters in the atmosphere and the corresponding consequences on the ionospheric absorption of the energy of sferics. Many papers related to atmospheric observations at the time of earthquakes have been published. Change of meteorological parameters and particularly increase of temperature has been noticed by [22] before strong earthquakes in Kamchatka. It has been mentioned that the perturbations of the electric conductivity of the atmosphere caused by the ionized gas release may induce lightning discharges [23]. Possible effects of radon and aerosol injection in the atmosphere prior to earthquakes have been mentioned by [24]. Anomalous cloud formation one week before a M7.8 earthquake in Iran has been reported in [25]. Thermal atmospheric irregularities have been observed in [26] prior to a M7.8 earthquake in Mexico.
Modelling has been also done. A model was discussed by [27] to explain how thunderclouds can be produced in the atmosphere above the sea before earthquakes. It involves atmospheric gradient of temperature and electrical charges induced by gas release on the sea surface. The model described in [26] is related to atmospheric changes due to air ionization by radon, water molecules attachment
Figure 8. Similar to Figure 5 (including the color scale) but when the time of the earthquakes with M ≥ 5.4 and d ≤ 20 km is shifted forwards 25 days. Data processing is identical as before.
to these ions, and at the end, a heat increase. More generally, these phenomena are associated to the so-called Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) described for example in [28] [29] [30] [31] [32], and to the change in the global electric circuit which exists between the Earth’s surface and the bottom of the ionosphere at the time of earthquakes [33] [34].
5. Conclusions
The occurrence of low dispersion whistlers observed by the satellite DEMETER has been correlated with the seismic activity. The number of low dispersion whistlers is linked to the number of lightning strokes occurring below the satellite except in the vicinity of the magnetic equator. Then the study is related to earthquakes located outside this area. Earthquakes with M ≥ 4.8 and various depths have been considered. Our analysis includes a comparison with the background thunderstorm activity and a superposed epoch method to consider all earthquakes. The main result is that an 11% increase of whistlers is observed the day before earthquakes with M ≥ 5.5 and d ≤ 20 km at a distance less than 200 km from the epicenters. No significant relation is observed for earthquakes with M ≤ 5.2.
The major question is: Could this result be useful for earthquake prediction? Unfortunately not because the lightning strokes which produce whistlers are the most common events occurring in the Earth’s atmosphere (about 2000 thunderstorms are active at any time [35]). This study is just a new proof that there is a LAIC mechanism prior to large and shallow earthquakes.
The DEMETER satellite is over an earthquake epicenter only a few minutes per day and only at a given local time. To remove this constraint, it is possible in the future to use the data of the World Wide Lightning Location Network [36] in order to perform a similar analysis.
Acknowledgements
This study is related to data recorded by the RNF experiment of the microsatellite DEMETER which was operated thanks to the French Centre National d’Etudes Spatiales (CNES). The DEMETER data shown in the paper can be downloaded from https://cdpp-archive.cnes.fr/.