Critical analysis of the electrostatic turbulence enhancements observed by DEMETER over the Sichuan region during the earthquake preparation

Abstract. In this paper, we report initial results from a detailed statistical study of the plasma waves observed by the DEMETER satellite over the Sichuan region during a period of 20 days encompassing the large earthquake of magnitude M =7.9 that occurred on 12 May 2008. The main objective of this paper is to present a statistical method to process and analyze plasma wave data and assist in detecting possible earthquake precursors among larger irregular disturbances arising from the natural variability of the ionized environment of the Earth. This method, presently used for dayside observations, involves two stages. First, VLF wave spectra are processed to recognize the various types of plasma waves usually observed at mid and low latitudes and derive a reduced number of parameters that fully characterize these emissions and may be conveniently used for a detailed statistical study. In a second stage, we perform a statistical analysis of the results by taking into account two "reference zones" displaced respectively 30° eastward and westward from the "epicentre zone". Plasma and wave disturbances possibly induced by earthquakes in preparation are likely to maximize close to the "epicentre zone", while natural disturbances associated, in particular, with the varying magnetic activity are rather uniform over a wider longitude sector, thus enabling the use of observations over the reference zones as a base line. The initial results of this study show a deviation of the power spectrum of electrostatic turbulence in the epicentre zone about 6 days prior to the earthquake but no significant anomalous variations can be observed on other characteristics of plasma waves. From the analysis of the data over the two reference~zones and using recently produced sector magnetic activity indices, we conclude that the enhancement of electrostatic turbulence is associated with magnetospheric processes rather than with pre-seismic activity.


Introduction
An earthquake of magnitude M = 7.9 occurred at 14:28:01 CST (China Standard Time) with the epicentre located at 103 • longitude and 31 • latitude and at a depth of 19 km close to the town of Chengdu in the Sichuan region of China. This event triggered a number of studies to investigate the possible existence of seismic precursors. The search for pre-seismic electromagnetic phenomena is not a new field since a first scientific paper may be traced back to the end of 19th century (Milne, 1890).
Over the last 30 years, numerous studies have been devoted to the search for ionospheric disturbances possibly triggered by pre-seismic processes in the lithosphere. Preseismic atmospheric-ionospheric anomalies as well as telluric current and ultra-low frequency electromagnetic waves have been reported (see Pulinets and Boyarchuk, 2004;Kamogawa, 2006;. These pre-seismic ionospheric and sub-ionospheric anomalies were reported initially by Antselevich (1971) and Gohkberg et al. (1989) using transmission observation of VLF electromagnetic waves propagating trough the earth-ionosphere waveguide.  showed statistically significant VLF/LF sub-ionospheric propagation anomalies before Published by Copernicus Publications on behalf of the European Geosciences Union. M ≥ 6.0 shallow earthquakes in their 7-year observation. From the VHF transmission observation, the atmospheric anomalies before earthquakes were statistically verified by Fujiwara et al. (2004). Liu et al. (2000) investigated the relationship between large earthquakes and ionospheric anomalies in the region of Taiwan using ionospheric critical plasma frequency, foF2, obtained from an ionosonde measurement. It was observed that foF2, corresponding to the electron density in the F2 layer, significantly decreased in the afternoon during a few days before M ≥ 6 earthquakes, such as during 3 ∼ 4 days before M = 7.7 Taiwan Chi-Chi earthquake on 21 September 1999. Similar electron density depression above Taiwan Island was observed by the GPS total electron contents (TEC) (Liu et al., 2001). From such observations, Liu et al. (2006) constructed quantitative criteria for ionospheric anomalies and established the statistical correlation between Taiwan M ≥ 5 earthquakes and ionospheric anomalies. The statistical results indicated that anomalies clearly appeared within 5 days before earthquakes. Recent statistical achievements seem to be highly encouraging for promoting further studies on the pre-seismic lithosphere-atmosphereionosphere (LAI) coupling (see Kamogawa, 2006). The occurrence rate of M7-class earthquakes is about 20 times per year all over the world and this renders advantageous satellite observations because of the capability of covering all seismic regions for the statistical investigation of pre-seismic LAI coupling. Many satellites have been launched in search of pre-seismic ionospheric anomalies (see Pulinets and Boyarchuk, 2004).
DEMETER is the first micro-satellite dedicated to searching and identifying disturbances of plasma and waves in the ionosphere that may be triggered by seismic activity . Launched on 29 June 2004, on a 98 • inclination, nearly helio-synchronous orbit at an altitude of 715 km, DEMETER has provided a nearly continuous survey of the ionospheric plasma and plasma waves over two local time sectors 09:00-10:30 LT and 21:00-22:30 LT since September 2004 over a period of more than 6 years until the end of the mission in December 2010. In our study we have used frequency power spectra of the ELF/VLF (∼20 Hz-20 kHz) electric field component perpendicular to orbit plane computed on-board from the measurements of the ICE instrument .
In the first section we briefly describe the automatic method used to process the VLF spectra developed to enable an accurate characterization of plasma waves and their statistical study. Then we present the data related to the period encompassing the Sichuan EQ and in the final section we describe the improved statistical method set-up to analyze the observations using results from 6 years of nearly continuous survey by DEMETER over a wide zone encompassing the Sichuan region. The electrostatic turbulence displays an anomalous behaviour over the epicentre region 6 days before the EQ but other plasma waves do not display any significant deviation from their average variation. In the last section of our conclusion, the possible coupling between seismic activity and ionospheric turbulence is discussed and we briefly indicate our ongoing studies in this area.

VLF plasma wave automatic characterisation
Several recent studies based on DEMETER data have focused on the analysis of plasma waves, e.g. Nemec et al. (2008), Molchanov et al. (2006) and Boudjada et al. (2008). Due to the variable conditions of occurrence of earthquakes and to the small amplitude of seismoelectromagnetic effects, it was not possible up to now to recognize specific wave disturbances associated with individual events.
Results from Nemec et al. (2008) appear the most convincing at the moment and were based on a thorough statistical analysis of ∼3000 EQ by sorting the intensity of natural waves in 10 frequency bands. However, a significant improvement in such statistical studies and in their interpretation would certainly be achieved if natural emissions (electrostatic turbulence, electromagnetic hiss, etc) could first be recognized and the statistical study performed on a reduced set of representative parameters that fully characterize each of them, such as the frequency band, the frequency of maximum power, the total power, etc.
Indeed, sorting wave data by frequency only and making statistics irrespective of the position of the EQ result in the mixing of several types of natural emissions that may have totally different origin and propagation conditions. As an example, in the frequency interval of a few hundred Hz to a few kHz, two main emissions are observed at mid and low latitudes (i) ELF plasmaspheric hiss, originating from waveparticle interactions near the equatorial plasmapause at 4-6 Earth radii altitudes, and (ii) whistlers, emitted by lightning currents in tropospheric thunderstorms close to the satellite ground track or in the conjugate area. It cannot be expected that a physical process associated with seismic activity may have a similar effect on hiss and whistlers, owing to their different sources and propagation paths. In addition, hiss emissions in the ordinary mode are controlled by the local H + gyro-frequency f H+ with their intensity peaking close to f H+ and decreasing very sharply until the cross-over frequency f co about 60 to 100 Hz below. Below f co and down to ∼150 Hz and depending on the local ionospheric plasma conditions, hiss can propagate in the extraordinary mode, most often with much less intensity than in the ordinary mode. At DEMETER altitudes, f H+ varies from ∼400 Hz at equator to ∼600 Hz at 60 • latitude. Consequently, waves observed close to an epicentre at 500 Hz may be either above or below f co and propagate in the ordinary or extraordinary modes, depending on the EQ latitude and longitude. Since statistical studies have to use a superposed epoch method on many earthquakes (EQ) spread over ±40 • in latitude, one should sort wave emissions by their origin (whistlers, hiss) and propagation mode before performing a statistical study.
In our study we have used the on-board calculated (20 Hz-20 kHz) VLF power spectra averaged over ∼2 s since they are available both in Burst and Survey modes. This analysis has first been performed for dayside data since natural emissions such as ordinary and extraordinary hiss, whistlers etc. are more clearly distinguishable. The same work is being performed for night data. In the first stage we have developed a method to automatically identify plasma wave emissions and characterize them by a reduced set of parameters such as the frequency band, the frequency of maximum power, the total power, etc. The method has been described in detail by Onishi and Berthelier (2010) and it will be briefly summarized here. Basically, the shape of each VLF power spectrum is described by a number of characteristic points representing local maxima, minima and discontinuities in the shape of the spectrum detected through local extrema of the second derivative of the power vs. frequency curve. A further processing of the characteristic points enables us to recognize and sort out the various natural emissions as well as the signals from the ground based VLF transmitters. In addition, the power spectra can be approximated with very good accuracy by a simple linear interpolation between successive characteristic points. The total information on wave emissions is thus reduced to a small set of parameters that are used, in the second stage, to perform the statistical study. An example is shown in Fig. 1.
The following information can be simply extracted from the characteristic points. In order to define a frequency range of each type of emission, characteristic points of a positive 2nd derivative (red points in Fig. 1) are used with additional conditions specific to the emission.
(i) Ordinary hiss (O-hiss) defined by -frequency and intensity of the peak hiss emission close to f H+ at point labelled a -cross-over frequency f co at point 3 and power index (slope in frequency) above f co and below the peak emission in the frequency range indicated by the slope a-b -frequency range above the peak emission and average power spectral index (slope in frequency). When calculating the upper frequency limit of ordinary hiss the threshold on power was set at 0.3 µV 2 m −2 Hz, the background level of ICE.
(ii) Extraordinary hiss (Ex-hiss) observed between points 2 and 3 -frequency range, peak frequency, total intensity (iii) Electrostatic (ES) turbulence observed in the frequency range labelled c below the extraordinary hiss and down to the lowest frequency of the VLF spectrum at 19.7 Hz and characterized by a negative power index (negative slope in frequency).
(iv) Whistler-type emissions are detected independently below and above f H+ the occurrence of sporadic time variations of the power intensity.
(v) VLF transmitter signals at 19.8 kHz from NWC with peak frequency, frequency band and total power, enabling us in particular to study the signal's broadening. Fig. 2 is an overview of the results for an entire orbit. In the top panel of the figure, characteristic points obtained from the automatic method are shown. Black and red dots correspond to negative and positive 2nd derivatives. In the bottom panel, peak power of the ordinary hiss indicated by red dots is determined from the characteristic points by selecting a local maximum nearest to the proton gyrofrequency f H+ . The cross-over frequency is indicated by pink dots and the proton gyro-frequency by a white line. In the same vein, the peak power of extraordinary hiss and electrostatic turbulence are recognized and indicated by blue and black dots, respectively. Whistler-type emissions are shown in white and treated separately from other emissions, even if some characteristic points relative to hiss emissions are detected inside the white parts of spectrogram.

Method of statistical analysis
One of the main difficulties in deciphering spacecraft observations of plasma waves to detect possible pre-seismic effects stems from the presence of natural emissions with large and very highly variable intensity. Except for whistlers generated by thunderstorms, natural plasma waves are mainly controlled by ionospheric and magnetospheric processes associated with auroral activity. The mid to low latitude effects of auroral activity are fairly well represented by magnetic activity indices such as Am (Menvielle and Berthelier, 1991), which may be used to evaluate, for each type of natural plasma waves and for a given period and magnetic activity, the statistical average intensity and the variability usually measured by the standard deviation. In principle, with the help of a large enough database and a thorough determination of the statistical parameters, one should be able to assess the statistical significance of a deviation from the mean value that may be observed at time of seismic activity. However the variety and complexity of ionospheric, magnetospheric and auroral phenomena and, more importantly, their prolonged effects over hours and even days not taken into account when using "instantaneous" indices of magnetic activity, induce a noticeable uncertainty in such statistical results. To attempt to overcome this difficulty, or at least to lessen its effect, we have developed an improved method of analysis by performing the same statistical study over two zones, West and East of the EQ zone, close enough to it but still far enough not to

Overview
Figure 4 displays an overview of the main parameters that have characterized the various types of plasma waves in the ELF frequency range for 6 years of DEMETER operation, from 2005 to 2010, between 3 and 22 May in each year and for the western and eastern reference zones and the EQ zone. The Sichuan earthquake occurred on 12 May 2008 and is indicated by a vertical red line in each panel. To take into account the dependence of magnetospheric phenomena on the Earth's magnetic field, geomagnetic latitude is used instead of geographic latitude and, for each orbit, average values of the parameters are calculated over the geomagnetic latitudes (11 • ∼ 31 • ) that encompass the 21 • geomagnetic latitude of the epicentre. Average and standard deviation, σ , are calculated in each zone for the whole set of data (i.e. for all years and all magnetic activities) and in each diagram the deviation from the mean value of the indicated parameter is plotted after normalization by σ . Data of the western, Sichuan and eastern zones are plotted in blue, red and green, respectively. When the deviation from the mean is greater than 4σ , a square symbol is placed above the data in the corresponding colour.
Several cases with deviations larger than 4σ are observed. In 2005 this was the case of the peak frequency in (Fig.4 1a) and, more distinctly, of the upper frequency limit of the ordinary hiss in (Fig.4 1e). These deviations are associated with large surges in magnetic activity associated with the magnetic storm in the second week of May 2005 and an extended period of sporadic substorm activity observed during the third week. They fit well to the known variations of plasmaspheric hiss at time of large auroral activity. Similarly the surge of magnetic activity corresponding to the beginning of the new solar cycle in 2010 is marked by the increase of the upper frequency limit of the ordinary hiss and of the level of electrostatic turbulence as seen both on the slope index and on the power at 20 Hz.
Early in May 2008, a large increase in the total power of whistlers is observed. A closer examination of this event as a function of time and position shows an intensified power level at the high latitude border of the Sichuan and East reference zones. This is better shown on the global map of the average total power of whistlers shown in Fig. 5. The power intensity above the NWC VLF transmitter on the northwest coast of Australia (slightly off to the north of the location of the transmitter due to the inclination of the Earth's magnetic field) is very large and concentrated. It is known that this region is subjected to heavy thunderstorm and lightning activity. In the northern conjugate region an increased intensity is also observed but distributed over a larger area. Since the total whistler power is calculated at frequencies below 10 kHz, these observations cannot stem from a parasitic influence of the NWC emission at 19.8 kHz. Following Parrot et al. (2009), this observation can be explained by an enhanced transmission of whistler waves through the ionosphere above the NWC transmitter through filamentary structures generated by the interaction of high power VLF waves with the ionospheric plasma.
Finally, on 3 and 6 May 2008, large deviations of two parameters related to electrostatic turbulence were observed, a positive deviation > 4σ for the power index, thus a steepening of the frequency spectrum and a positive deviation > 3σ for the power at 20 Hz, the two effects resulting in fact from an increase of the electrostatic turbulence intensity at 20 Hz. One of the physical mechanisms could be suggested to explain the occurrence of seismic precursors in the ionosphere as the development of a low intensity background of gravity or acoustic waves generated during the EQ preparation either through lithospheric processes in the fault or through surface effects. When propagating in the upper atmosphere, these waves interact with the ionospheric plasma and produce field-aligned irregularities in the topside ionosphere, as was observed in the case of MSTID's of meteorological origin by Onishi et al. (2009). Such irregularities may be driven unstable by the convective motion of the plasma, leading to the development of electrostatic turbulence (ET) which may therefore be considered as a possible precursor signal in the ionosphere. We have paid a great attention to the ET disturbances observed on 6 May.

Electrostatic Turbulence (ET) variations
Since magnetic activity was slightly increased on 3 and 6 May compared to the neighbouring days, we investigated the statistical distribution of the values of the power index as a function of magnetic activity independently for the three zones. The results are plotted in Fig. 6.
In each zone, data from all 6 years are binned and counted in cells of 0.5σ of the variation and 1nT of the αλ longitude sector index and the histograms are normalized for each value of sector index to give a distribution function in percentage at a given magnetic activity. Events on 3 and 6 May are indicated. Because of the reduced set of days with high enough magnetic activity, the average variations of the ET power index for magnetic activity above 20 nT are characterized with a lower accuracy. Nevertheless, the power index appears to increase when magnetic activity is ≥20 nT in the Sichuan and eastern zones and ≥40 nT in the western zone. The linear fits also illustrate the increasing trend of the parameter.
The large deviations observed in the Sichuan and western zones on 6 May occur at local magnetic activity indices of 15 nT in the western zone, 15 nT and 30 nT in the Sichuan zone, and 30 nT in the eastern zone. The two activity levels for the Sichuan zone arise from the different position of the satellite orbit with respect to the boundary between the αλ index sectors N2 and N3 in Fig. 3. The main information is that, on 6 May, the magnetic activity increased regularly from the western reference zone to the eastern reference Nat. Hazards Earth Syst. Sci., 11, 561-570, 2011 www.nat-hazards-earth-syst-sci.net/11/561/2011/  zone. Similarly, on 3 May the large deviation corresponds to a magnetic activity of 10 nT in the western zone, 10 to 15 nT in the Sichuan zone and 15 nT in the eastern zone, and again magnetic activity was regularly increasing from west to east but at a slower pace than on 6 May. To summarize, the event on 6 May is the largest one in the Sichuan zone with the deviation of ET > 4σ from its average level in the Sichuan and western zone, but insignificant (< 1σ ) in the eastern zone. On 3 May, the deviation is > 4σ in the western zone, ∼ 2.5σ in the Sichuan zone and, again, insignificant on the eastern zone. One must be aware that the small number of observations at high magnetic activity (≥ 30 nT) induces a significant uncertainty in evaluating the amplitude of the deviation in the eastern zone where magnetic activity is always the highest (and ≥ 30 nT on 6 May). We focus in the following on the observations on 6 May. As a first conclusion, the deviation of ET is as large in the western zone as in the Sichuan zone, contrary to what would be expected for a disturbance induced by seismic activity in the EQ region.
The study was then extended to the southern latitudes of the three zones as indicated in Fig. 3. Plotted in Fig. 7 is a synthetic colour coded representation of the whole set of ET power index measurements from 2005 to 2010 in the Sichuan zone and the two reference zones in the West and the East. The event of 6 May is clearly visible slightly southward of the epicentre and one can also observe, along the same orbit, a similar variation of the ET parameters at southern latitudes in the conjugate region of the ET enhancement in Sichuan zone. The intensity of the ET disturbance in the South is actually larger than in the Northern Hemisphere over the Sichuan region. The western reference zone displays a similar situation with a large deviation of the ET power index from its average value in the southern conjugate region, more intense and spread over a wider range of latitudes. The existence of conjugate ET disturbances can be explained by the strong electrical coupling between the northern and southern ionospheres. Nevertheless, the existence of disturbances with a maximum intensity in the Southern Hemisphere cannot be expected if the source is in the opposite hemisphere, associated with the Sichuan EQ.
The presence of ET disturbances of similar intensity in the western reference zone compared to the Sichuan zone and the observation of higher intensity ET disturbances in the southern conjugate area of the western and Sichuan regions suggest that the ET event are due to ionospheric and magnetospheric processes associated with a surge of magnetic activity. This is supported by data shown in Fig. 8, where the variations of the αλ longitude sector indices from 4 to 14 May have been plotted. The surge in magnetic activity on 5 and 6 May is clearly visible in both hemispheres. More importantly, the Southern Hemisphere diagram definitely shows the occurrence of a very large spike of magnetic activity in the 00:00-03:00 UT time interval on 6 May, coincident with the DEMETER orbits over the Sichuan and West zones at 01:00 UT and 03:00 UT respectively. This is a rather unusual dissymmetrical enhancement of magnetic activity that naturally explains the higher intensity and wider domain of occurrence of the ET disturbance in the Southern Hemisphere.

Conclusions
We have performed a detailed statistical study to search for possible anomalous variations in the VLF natural emissions observed onboard DEMETER over the Sichuan region during a 20 day period encompassing the Great Sichuan earthquake on 12 May 2008. In the first stage, the raw spectral The data are sorted along three zones, the Sichuan zone (central column) and the two reference zones west (left) and east (right) to the Sichuan zone. In each column, geomagnetic latitudes are on the x-axis, positive in the Northern Hemisphere, negative in the Southern Hemisphere. For each year the time in days varies along the y-axis between 3 and 22 May. The time and geomagnetic latitude of the Sichuan EQ are indicated by horizontal and vertical red lines. The EQ geomagnetic latitude is indicated by a dotted red line in the reference zones. All earthquakes with M ≥ 5 that occurred in these zones and periods are marked with circles, whose diameter is proportional to the magnitude. data were processed to detect and characterize the different types of plasma waves usually observed at mid and low latitudes. Electrostatic turbulence, through its power spectral index and its intensity at 20 Hz, has shown an anomalous variation in time and space 6 days prior to the earthquake. In the second stage an improved statistical method has been developed to search for a likely origin of these disturbances. Using nearly simultaneous observations over two reference zones west and east of the Sichuan zone and a finer representation of magnetic activity by means of longitude sector indices, we opt for the conclusion that all results are consistent with the observed enhancement of electrostatic turbulence being generated by magnetospheric processes that are at the origin of the surge of magnetic activity. These findings were obtained thanks to the new and improved statistical method of analysis that rely on reference zones and longitude sector magnetic indices and also to the very large database acquired by DEME-TER over 6 years of operation. They point out the significant difficulty and possible errors that can be made when assessing with less precaution the origin of disturbances observed during the preparatory phase of large earthquakes. Fig. 8. αλ indices of magnetic activity in N2 (blue triangle) and N3 (green square) in the Northern Hemisphere (top) and in S1 (blue triangle), S2 (green square) and S3 (purple inverse triangle) in the Southern hemisphere (bottom) in the period from 4 to 14 May. The ET disturbance recorded over the Sichuan region and its conjugate area occurs on 6 May. Range of indices are indicated in grey by maximum and minimum indices. Global am-index is plotted in red circles. Vertical dotted lines correspond to 00:00:00 UT of each day.
We are currently extending our analysis to the ionospheric plasma measurements provided by DEMETER during the same Sichuan period and to other large earthquakes observed during the DEMETER observation phase. The results will be reported in forthcoming papers.