The neutral temperature in the ionospheric dynamo region and the ionospheric F region density during Wenchuan and Pingtung Doublet earthquakes

One of the possible candidates which modifies the ionosphere before large earthquake is electric field. We presume that the electric field associated with large earthquakes is generated in the ionosphere dynamo region (100–120 km). This paper tries to identify the evidence of the contribution of the neutral atmosphere in the dynamo region. The relationship between the critical frequency at the F2 peak (foF2) and the height profile of the neutral atmosphere temperature was studied for two large earthquakes: Wenchuan, 2008 and Pingtung Doublet, 2006. It is found that the wave amplitude of the vertical wavelength 20–30 km which is usually superposed on the height profile of the neutral atmosphere temperature enhances when the foF2 increases. The correlation between the wave amplitude and foF2 is found better along a longitudinal direction than along latitude direction.


Introduction
The disturbance of the ionosphere which appears prior to large earthquakes has been repeatedly reported.Among the reports, most of the papers in the early stage discuss on the change in foF2, which is observed by ground-based ionosondes.The number of these papers has rapidly increased in the last 2-3 yr and the quality of the papers seems to be improving; the reports are more convincing.Zhao et al. (2008), Liu et al. (2009), and Pulinets and Ouzounov (2010) report that the enhancement and/or reduction anomalies of the Correspondence to: Y. Y. Sun (yysun@jupiter.ss.ncu.edu.tw)ionospheric electron density, which is measured from the ground-based GPS receiver, ionosonde, and FORMOSAT-3/COSMIC appeared above the epicenter 3-6 days before the 12 May 2008 M = 7.9 Wenchuan earthquake.Le et al. (2011) studied the pre-earthquake ionospheric anomaly statistically by using the total electron content (TEC) data from the global ionosphere map (GIM) for a total of 736 M = 6.0+ earthquakes in the global area during 2002-2010.The results indicate that the anomalous behaviors of TEC within just a few days before the earthquakes are related to the forthcoming earthquakes with high probability.Liu et al. (2011) report for the first time that there is a successive long time (nearly 24 h) enhancement in GPS TEC one day before 12 January 2010 M = 7.0 Haiti earthquake, and the TEC enhancement anomaly appears specifically and persistently in a small region of the northern epicenter area; these two factors suggest that the enhancement anomaly before the earthquake is highly related to the Haiti earthquake.Recently, research by combining satellite data as well as ground-based observations has started (Oyama et al., 2008(Oyama et al., , 2011)).It is becoming clearer that some large earthquakes modify the ionosphere in some way (Depueva and Rotanova, 2001;Singh et al., 2004;Liu et al., 2009).However, we do not know whether all large earthquakes (M > 6.5) disturb the ionosphere, as the morphology of the ionosphere disturbance has not been established yet.
Small disturbance of the ionosphere is even found at the geomagnetic conjugate point of the epicenter for some earthquakes (Zhao et al., 2008;Liu et al., 2009Liu et al., , 2011;;Pulinets and Ouzounov, 2010).A phenomenon which is very similar to the Equatorial Ionization Anomaly (EIA) has been reported to appear over the epicenter recently (Oyama et al., The solid dark lines denote the hourly average.Observations were not available during 10:00 to 14:00 LT.

2011
).These observations imply the role of the electric field (Ruzhin et al., 1998).Provided that the electric field is working, the next question is "Where and how is the electric field generated"?One of the plausible mechanisms is the internal gravity wave.It might be generated on the ground before the earthquake by some mechanisms, propagate upwards and then change the wind pattern around the height of 100-120 km.Once the wind pattern has changed at the dynamo region, the original electric field is modified and the change of the electric field is then propagated to the higher ionosphere along the magnetic field line (England et al., 2006).Rozhnoi et al. (2007) present evidence of atmospheric gravity waves induced by the seismic activity found several days before earthquakes.Their frequency range is 0.28-15 mHz, which corresponds to the period range of 1-60 min.This frequency is still too high to explain the behavior of the foF2 which we mention here.
In order to find the evidence of the suggested mechanism, we tried to study the relationship between the height profile of neutral temperature (Tn), and foF2 for two large earthquakes, the Wenchuan and Pingtung Doublet earthquakes.Wenchuan earthquake (M w = 7.9, Depth = 19 km, Epicenter = 31.021

Data
In this study, the foF2 measurements are provided by the Wuhan (30.5 • N, 114.4 • E) and Pingtung (22.6 • N, 120.6 • E) ionosonde stations for the Wenchuan and Pingtung Doublet earthquakes during 1-14 May 2008 and 19-28 December 2006, respectively.Information on Tn was provided by the SABER (Sounding of the Atmosphere using Broadband Emission Radiometry) instrument on board the TIMED (Thermosphere Ionosphere Mesosphere Energetics and Dynamics) satellite (Zhang et al., 2006).As a first step in our research to find evidence of the coupling between foF2 and the dynamo region, we studied the variation of intensity of Tn in the height range of 15-140 km.Then we examined the irregular structure of height profiles of Tn.Generally, the Tn profile around 100 km shows wavy structures with the vertical wavelength of 20-30 km modified by internal gravity waves.Hereafter, we will describe the data analysis.1a).The amplitude of the diurnal variation is about 40 K. Simultaneously, Tn around Pingtung shows 360 K during nighttime and 390 K during daytime (Fig. 1b).The amplitude diurnal variation is about 30 K, which is about 10 K lower than that for Wenchuan region.By contrast, the local time variation is not clear at the altitude region of 15-20 km (Fig. 1c and d).Even if the local time variation exists, the amplitude is less than 3 K.The average Tn for Wenchuan and Pingtung are about 204 K (Fig. 1c) and 198 K (Fig. 1d), respectively.
Figure 2 depicts the variation of the daily averaged Tn at the two height regions in the years 2002-2009 over the epicenters of the Wenchuan and Pingtung Doublet earthquakes.At the height of 115-120 km, Tn shows two minima in April and October-November above the Wenchuan area (Fig. 2a).
The averaged Tn is about 355 K in both months.Two peaks of the averaged Tn occurred around January and June-July.The averaged Tn is about 400 K.While for Pingtung area (Fig. 2b), two minima of about 360 K were found in April-May and October-November.Two Tn peaks of 400 K were found in January and July-August, respectively.By comparing Tn in the Wenchuan and Pingtung areas, the amplitude of annual variation was larger in Wenchuan than in Pingtung.At the height region of 15-20 km, the Wenchuan area showed a peak of Tn of 210 K in January, and a minimum (Tn = 203 K) around June-July (Fig. 2c), whereas Tn in the Pingtung area showed a peak of Tn (202 K) around August-October, and a minimum in February and March (Tn = 199 K) (Fig. 2d).
Figure 3 illustrates the deviation of the daily averaged Tn by subtracting the reference (the daily average of Tn during 2002-2009).For the Wenchuan earthquake, the two big drops of the Tn deviation appeared in 1-20 April and 1-10 May 2008 at heights of 15-20 km (Fig. 3b).At the height region of 115-120 km, another positive peak existed around 20 April 2008 (Fig. 3a).The minimum at 15-20 km seems to correspond to the peak at the height for 115-120 km.This anti relationship seems to be true.Since another peak was found between 10 October and 1 December 2008 (Fig. 3a), accordingly we cannot conclude that the peak located around 20 April is due to the Wenchuan   earthquake.For the Pingtung Doublet earthquake, at the height of 15-20 km, the largest deviation of Tn was found on 25 December 2006 right before the earthquake (Fig. 3d), and the largest drop existed at the height of 115-120 km around 5 December 2006 (Fig. 3c).However, the deviation of the other years is also large; again we cannot draw a firm conclusion on this special feature of the Tn variation associated with the Pingtung Doublet earthquake with which we can persuade the readers.), which most of the time exists.The irregular structures might be mainly due to the internal gravity wave.To find the various wave componets from the height profile of Tn, the time-frequency analysis method which is named Hilbert-Huang transform (HHT) (Huang et al., 1998) is conducted.
HHT consists of Empirical Mode Decomposition (EMD) and Hilbert spectral analysis.EMD is able to decompose the data into several orthogonal components named Intrinsic Mode Functions (IMFs), and the Hilbert spectral analysis can subtract the instantaneous information of the wavelength and amplitude from these IMFs. Figure 4 illustrates the two sets of IMFs which are decomposed from the individual Tn profile using EMD. Figure 4a displays that the Tn profile (22:38 LT, 8 May 2008, 31.7 • N, 111.3 • E) can be decomposed into 5 components (4 IMFs and 1 DC term).The wavelength of IMF1 and IMF2 is shorter than 20 km, and IMF4 and IMF5 are the long-term curves.However, 6 components (5 IMFs and 1 DC term) are isolated from the Tn profile which measure at 23:08 LT, 9 May 2008, 26.2 • N, 122.3 • E (Fig. 4b).The IMF3 in Fig. 4b indicates the wavy structure with the wavelength of 20-30 km, and the amplitude is greater than 25 K at the altitude of 90-140 km. Figure 5b displays the Hilbert sepctrum which was contructed from the 12 Tn profiles (Fig. 5a) over the epicenter of the Pingtung Doublet earthquake about two hours before the earthquake occurrence.We find that wavy structures with the wavelength longer than 15 km appear at the altitude of 80-140 km (Fig. 5b).To further examine the possible relationship between the earthquake and the Tn wavy structures, we conducted HHT for three areas with a different area over the Wenchuan and Pingtung Doublet earthquakes (Fig. 6).The areas for Wenchuan and Pingtung are marked by red and blue squares, respectively.
Figure 7 depicts the change of the Hilbert spectra for the Pingtung Doublet earthquake at 17:00-20:00 LT.The top three panels (Fig. 7a-c) display the spectra for the smallest area (17-27 • N, 111-131 • E), and unfortunately the number of the Tn profiles is not great enough to show the general picture.The middle three panels (Fig. 7d-f) (12-32 • N, 101-141 • E) show that pronounced wavy structures can be found on 26 December (DOY 360) 2006 (Fig. 7f, same as Fig. 5b), but the 24 and 25 December (DOY 358 and 359) which are indicated in Fig. 7d and e   and the intensity of the wavelength longer than 30 km becomes stronger.
The same data analysis process was conducted for Wenchuan Earthquake on 7-9 May (DOY 128-130) 2008 as shown in Fig. 8.The right panels (Fig. 8c, f, and i) show that the intensity of the wave amplitude, with the wavelength longer than 25 km, at the altitude higher than 60 km, is faded when the area expanding from 26-36 • N, 93-113 • E to 21-41 • N, 83-123 • E, and 16-46 • N, 73-133 • E. However, the continued weakness and strength of the wave amplitude can be found in the left (Fig. 8a, d and g) and middle (Fig. 8b, e  and h) panels, respectively, as the area is expanded.The fading red on the spectra (Figs.7f and i, and 8f and i) represents that the wavy structures with the amplitude of about 20 K around 100 km in altitude concentrated over the epicenters.
Figure 9 shows the daily variation of the foF2 (dark points) and the wave amplitude (dark stars) of the wavy structures during 15 days before and after the occurrence of the Pingtung Doublet and Wenchuan earthquakes.The daily variation in foF2 at Pingtung (Fig. 9a) displays that the increase of the foF2 occurred around noon of 23 December 2006 (−3 day) and the increase gradually reduces toward 25 December (−1 day).In the evening of 25 December, foF2 is as low as 2 MHz.This similar pattern seems to exist during 20-22 December (−6 to −4 day).The foF2 increases during the daytime of 20 December (−6 day) with a small decrease in the afternoon and decreases towards midnight of 22 May (−4 day) and in the early morning of 23 May (−3 day).The daily variation of the foF2 at Wuhan (Fig. 9b) for the Wenchuan earthquake (Zhao et al., 2008) shows that the increase of foF2 can be seen in the afternoon on 9 May (−3 day).The foF2 gradually decreases toward 11 May 2008.The similar feature of foF2 is seen during 5-7 May (−7 to −5 day).The foF2 increases in the morning of 5 May (−7 day), and decreases in the afternoon.This tendency continues for 3 days up to 7 May (−5 day), and finally, it becomes lower than the normal values.Especially during the nighttime, the reduction rate of the foF2 with respect to the average value is faster.This kind of tendency is often found for large earthquakes and therefore seems to be a common feature.The details of these features will be published in the near future.
Note that in Fig. 9, the wave amplitude of the wavy structures is averaged from the amplitude in the Hilbert spectra at the height range of 80-120 km, and the wavelength of 20-30 km (see Fig. 8f).For the Wenchuan earthquake, the wave amplitude and foF2 show a very similar feature during 11 days before to 3 days after the earthquake (Fig. 9b).
For the Pingtung Doublet earthquake, the similarity is not as clear in comparison with the Wenchaun earthquake.However, the wave amplitude becomes larger at about 2 h before the earthquake occurrence (Figs.9a and 7)and the wave of larger amplitude prevails after the earthquake, as one can see on the afternoon of 28 December (2 day).
To further understand the possible relationship between the foF2 and the wave amplitude of wavy structures around 100 km, we compared the daily maximum of the wave amplitude and foF2 in nine different regions over and around the epicenter of the Pingtung Doublet earthquake. Figure 10 shows the relationship between the foF2 and the amplitude of the dominant wave with respect to earthquake days.The correlation between the foF2 and the wave amplitude is better in the panels d, e, and f, while panels a, b, c, g, h, and the correlation is not shown.These panels suggest that the correlation between the foF2 and the wave amplitude is better in longitude, but not in latitude.This seems to be consistent with the fact that the disturbance is more extended in the longitude direction (Oyama et al., 2008(Oyama et al., , 2011)).Figure 11 shows the relationship between the foF2 versus the amplitude of wavy structures.The linear line in each panel is a least square fit to the data.Positive correlation between the foF2 and the wave amplitude is seen in the panels d, f, and i, while in panels a, b, and c, the correlation is rather negative.Figures 12 and 13 provide similar plots for the Wenchuan earthquake. Figure 12 shows the variation of the foF2 and the amplitude of wavy structures from 15 days before and after with respect to the earthquake day.In panels d, e, and f, the correlation between the two values is good, while in other panels the correlation is not identified.This feature again suggests that the disturbance extends wider in the longitude direction.Figure 13 shows the foF2 versus the amplitude of wavy structures and the positive correaltion can be seen in all panels.

Discussion and conclusion
We studied the variation of height profiles of the neutral temperature, Tn, in order to find the relationship between the foF2 and wavy structures of Tn in the ionospheric E-layer dynamo region which is possibly associated with large earthquakes, such as the Wenchuan and Pingtung Doublet earthquakes, assuming that the electric field which modulates the F-layer ionosphere is generated in the ionosphere E-layer dynamo region.The wavy structure of Tn at the altitude of the E-layer, which especially appears in the vertical wavelength of 20-30 km, seems to control the variation of foF2.The variation of the daily maxima foF2 yields a close agreement with the tendency of the wave amplitude over the epicenter from 10 days before the occurrence of the Wenchuan  earthquake (Fig. 12e), although the corresponding correlation coefficient is only 0.43 (Fig. 13e).This moderate correlation means that some other possible mechanisms might also be participating in generating the variation of foF2, or the enhancement of the amplitude of the wavy structures with a vertical wavelength of 20-30 km in Tn height profiles does not always efficiently modify the E-region dynamo.
On the orther hand, the appearance of Es, during the time before the large earthquake, which was reported by Ondoh (2009), might be also explained by the existence of these waves.The vertical wavelength of 20-30 km was observed in the presence of Es (Yoshimura, 2003).However, the mechanism of the wave generation prior to large earthquakes is still in the darkness.Maybe the change of Tn in the lower  atmosphere due to latent heat flux might be one of the possible sources (Pulinets et al., 2010).
The agreement of tendencies of the foF2 and the amplitude of wavy structures (Figs.9-13) is much clearer for the Wenchuan than for the Pingtung Doublet earthquake.One of the possible reasons for this might be the difference in the magnetic latitude of earthquakes.The epicenter of the Pingtung doublet earthgquake (magnetic latitude, 11.90 • N) is closer to the magnetic equator than the Wenchuan earthquake (magnetic latitude, 20.77 • N).Therefore, the E × B drift might be modulated by the internal gravity wave, which is relatively more active for the Wenchuan earthquake than for the Pingtung Doublet earthquake.Moreover, another reason might be related to the energy release of earthquakes.Liu et al. (2006) indicate that the odds of earthquakes with preearthquake ionospheric anomalies increase with the earthquake magnitude but decrease with the distance from the epicenter to the ionosonde station.
In conclusion, we have studied two large earthqaukes, the Wenchuan and Pingtung, by assuming that the ionopshere is

Fig. 1 .
Fig. 1.The diurnal variation of Tn around the epicenters of Wenchuan (left panels) and Pingtung Doublet (right panels) earthquakes during 2002-2009.Gray points are individual Tn measurements at two height regions of 115-120 km (upper panels) and 15-20 km (lower panels).The solid dark lines denote the hourly average.Observations were not available during 10:00 to 14:00 LT.

Fig. 2 .
Fig. 2. The annual variation of Tn around the epicenters of the Wenchuan (left panels) and Pingtung (right panels) earthquakes during 2002-2009.Gray points are the 5 day moving average of the individual Tn measurements at two height regions of 115-120 km (upper panels) and 15-20 km (lower panels).The dark points are the daily average.

3
Figure 1 displays the local time variation of the average Tn at two height regions of 15-20 km and 115-120 km over the epicenters of the Wenchuan and Pingtung Doublet earthquakes.All data were collected and averaged during 2002-2009.Tn measurements were not available between 10:00 and 14:00 LT.At the height range of 115-120 km (Fig.1a and b), the hourly average Tn illustrates a clear local time variation.During nighttime, Tn shows 345 K, and 390 K during daytime over the epicenter of the Wenchuan earthquake (Fig.1a).The amplitude of the diurnal variation is about 40 K. Simultaneously, Tn around Pingtung shows 360 K during nighttime and 390 K during daytime (Fig.1b).The amplitude diurnal variation is about 30 K, which is about 10 K lower than that for Wenchuan region.By contrast, the local time variation is not clear at the altitude region of 15-20 km (Fig.1c and d).Even if the local time variation exists, the amplitude is less than 3 K.The average Tn for Wenchuan and Pingtung are about 204 K (Fig.1c) and 198 K (Fig.1d), respectively.Figure2depicts the variation of the daily averaged Tn at the two height regions in the years 2002-2009 over the epicenters of the Wenchuan and Pingtung Doublet earthquakes.At the height of 115-120 km, Tn shows two minima in April and October-November above the Wenchuan area (Fig.2a).

Fig. 3 .
Fig. 3.The deviation of Tn (the annual vairaion is removed) observed in the years 2006 and 2008 above the epicenters of the Wenchuan (left) and Pingtung Doublet (right) earthquakes at two height regions of 15-20 km (lower panels) and 115-120 km (upper panels).The dashed red squares indicate some significant Tn deviations.

Fig. 4 .
Fig. 4. Examples of the IMFs which decomposed from the individual Tn profiles using EMD.(a) A case without the IMF3 component and (b) a case with the IMF3 component.

4
Small scale variation of individual neutral temperature profileTIMED/SABER provided us the height profiles of Tn within the height of 15 km to 140 km.Generally, at the height of 20 km, Tn is about 200 K, and it gradually increases toward 50 km, where Tn shows about 250 K, and starts reducing at higher altitude.At around 75 km, the height profile of Tn starts to show irregular structures, and the amplitude of irregular structures increases between 80-120 km(Figs.4 and  5a do not show any clear wave activities.As the area is expanded as shown in Fig. 7i (7-37 • N, 91-151 • E), the pronounced wave signals fade out

Fig. 6 .
Fig. 6.The areas where HHT is conducted for the Wenchuan (red squares), and Pingtung Doublet (blue squares) earthquakes.The epicenters are located in the center of the squares and marked by red stars.Hilbert spectra in these three regions for two large earthquakes are displayed in Figs. 10 and 11.For the Pingtung Doublet earthquake, the first area is 17-27 • N in latitude, 111-131 • E in longitude.The second area is 12-32 • N in latitude, 101-141 • E in longitude.The third area is 7-37 • N in latitude, 91-151 • E in longitude.For the Wenchaun earthquake, the first area is limited to 26-36 • N in latitude, 93-113 • E in longitude.The second area is 21-41 • N in latitude, and 83-123 • E in longitude.The third area is 16-46 • N in latitude, and 73-133 • E in longitude.The gray dots denote the location of TIMED/SABER Tn profiles on 12 May 2008.

Fig. 7 .
Fig. 7.The change of the Hilbert spectra in three square areas shown in Fig. 6 for three days (DOY 358, 359, and 360) for the Pingtung Doublet earthquake.The horizontal bar graphs beside the spectra denote the number and local time of the measuremets.Note that the amplitude of wavy structures with wavelength of 20-30 km enhances on DOY 360.

Fig. 8 .
Fig. 8.The change of the Hilbert spectra in three square areas shown in Fig. 6 for three days (DOY 128, 129, and 130) for the Wenchuan earthquake.The horizontal bar graphs beside the spectra denote the number and local time of the measurements.The red square in (f) encloses the region within 20-30 km wavelength and 80-120 km altitude.

Fig. 9 .
Fig. 9.The comparison of the foF2 and the intensity of the internal gravity wave with the wavelength of 20-30 km at 80-120 km in altitude for the Pingtung Doublet and Wenchuan earthquakes during 15 days before and after the earthquakes.The square area studied are 12-32 • N, 101-141 • E, and 21-41 • N, 83-123 • E for the Pingtung Doublet and Wenchuan earthquakes, respectively.The gray lines are the hourly average foF2 of the 14 days observation for the Wuhan ionosonde station and 10 days observation for the Pingtung ionosonde station, respectively.The dark stars indicate the intensity of internal gravity waves within the wavelength range of 20-30 km.The vertical dotted black lines denote the time of occurrences of the two earthquakes.

Fig. 10 .
Fig. 10.The comparison of the foF2 and the intensity of the gravity wave in 9 areas for the Pingtung Doublet earthquake.The star line displays the daily maxima of the amplitude of the gravity wave of 20-30 km wavelength and the point line denotes the daily maxima of the foF2.The vertical dotted dark lines denote the day of the earthquake occurrence.

Fig. 11 .
Fig. 11.The scatter plot of the daily maxima of the foF2 and the amplitude of the gravity wave in 9 areas for the Pingtung Doublet earthquake.

Fig. 12 .
Fig. 12.The same as for Fig. 10 but for the Wenchuan earthquake.

Fig. 13 .
Fig. 13.The same as for Fig. 11 but for the Wenchuan earthquake.