Fractal analysis of subionospheric LF propagation data and consideration of the lithosphere-atmosphere-ionosphere coupling

Fractal analysis has been applied to the local nighttime data of subionospheric LF propagation, and the fractal dimension is estimated every day in the two distinct frequency ranges (AW: acoustic wave and AGW: atmospheric gravity wave). The data during several years are analyzed for the propagation paths from the Japanese transmitter of JJY to Moshiri (Hokkaido) and to Kochi. As the result of analysis, we come to the conclusion that when we pay attention to the period just around the earthquake, we sometimes detect some significant increases in the fractal dimension either in AW or AGW range. This indicates that the self – organization effect prior to an earthquake in the lithosphere, might be seen even in the lower ionosphere, probably in terms of atmospheric oscillation effect.


Introduction
During the last decade there have been accumulated a lot of evidence on the presence of precursory electromagnetic signature of earthquakes (EQs) (e.g., Hayakawa and Molchanov, 2002;Molchanov and Hayakawa, 2008;Hayakawa, 2009).The fundamental idea on how to extract a precursory EQ effect from any time series data, is based on the simple statistical approach based on how the observed behavior in the time series record is abnormal as compared to the regular variation.For example, we estimate the average value and its standard deviation, or the median and the quantile, which are conventionally used to find an anomaly (e.g., Hayakawa et al., 1996a;Liu et al., 2006).However, these kinds of statistical analyses are not sufficient enough to convince the readers that the observed anomaly is really seismogenic.So that, we have to adopt any other methods Correspondence to: M. Hayakawa (hayakawa@whistler.ee.uec.ac.jp) based on the consideration of any physical process closely related with the nonlinear process taking place in the lithosphere prior to an EQ.
We can show a few examples in this direction especially for seismogenic ULF emissions.In 1996 a new proposal on the use of polarization method (based on the ratio of the vertical magnetic field component to the horizontal ones) has been presented by Hayakawa et al. (1996b) in order to distinguish the seismogenic ULF emission from the conventional magnetospheric geomagnetic variation (+geomagnetic pulsations), which is one kind of physical analyses.Then, the first attempt of fractal analysis on the ULF emissions has been made in Hayakawa et al. (1999), who found that the fractal dimension of ULF emissions was closely related to the nonlinear effect of the lithosphere, or so-called self-organized criticality (SOC) of the lithosphere.Extensive works on mono-and multi-fractal analysis have been performed later not only for ULF emissions (Smirnova et al., 2001(Smirnova et al., , 2004;;Gotoh et al., 2003Gotoh et al., , 2004;;Ida et al., 2005;Ida and Hayakawa, 2006;Hayakawa and Ida, 2008), but also for DC seismic signals (Varotsos et al., 2002) and VHF radio emissions (Eftaxias et al., 2004;Yonaiguchi et al., 2007).
The same situation holds for the seismogenic ionospheric perturbations.Based on the conventional statistical approach and by using the average value and standard deviation, Rozhnoi et al. ( 2004), Maekawa et al. (2006) and Kasahara et al. (2008Kasahara et al. ( , 2010) ) have found the significant correlation between VLF/LF propagation anomalies and EQs (with magnitude greater than 6.0 and with shallow depth).As is mentioned in the previous paragraph, the fundamental agent of possible ionospheric seismogenic perturbation is manifested itself in the lithosphere.
Recently, Hayakawa et al. (2009) have analyzed the effect of Earth's tides in the seismogenic precursory phenomena, in order to study the condition of the lithosphere (that is, whether the EQ focal region is under the critical or supercritical condition).The effect of Earth's tides is found to be present very clearly in the lithosphere as in the case of ULF emissions, while it is also seen in the subionospheric VLF/LF data, but not so clearly as in the case of ULF emissions.
So, the purpose of this paper is to examine whether such a kind of SOC phenomenon is manifested itself or not in the subionospheric VLF/LF propagation data.We use higher temporal resolution data to study the fractal characteristics of VLF/LF data.

Subionospheric LF data and EQs freated
As is already known in the previous papers (Rozhnoi et al., 2004(Rozhnoi et al., , 2007;;Muto et al., 2008Muto et al., , 2009;;Kasahara et al., 2010), the most probable candidate of the lithosphere-ionosphere coupling is the atmospheric oscillation channel (Hayakawa, 2004(Hayakawa, , 2009;;Molchanov and Hayakawa, 2008), in which the atmospheric gravity waves (AGWs) are excited by the variation in the geochemical quantities (such as pressure, temperature, etc.) due to the lithospheric activity near the EQ hypocenter or epicenter.The sampling of subionospheric LF propagaton data of 120 s(=2 m) was used in the previous our conventional analysis (e.g., Maekawa et al., 2006;Kasahara et al., 2008;Hayakawa et al., 2010a, b).However we want to analyze the wider frequency range, so that we try to use the higher temporal resolution data.Sometimes we obtain the data with sampling of 20 s.
Although we have abundant data during over ten years, such higher temporal resolution data with sampling of 20 s are available only for a limited number of propagation paths.
Figure 1 illustrates the relative locations of JJY-MSR (Moshiri, Hokkaido) and JJY-KCH (Kochi).The location of our Japanese LF transmitter with a call sign of JJY is indicated with a blue diamond in Fig. 1 in Fukushima prefecture.And, the VLF/LF receiving stations are indicated by red stars (MSR and KCH).The details of VLF/LF observations are given in Hayakawa (2004Hayakawa ( , 2007Hayakawa ( , 2009)).
The EQs treated for the above propagation paths are also plotted in Fig. 1, with magnitude greater than 5.5.The center of a circle is the location of an EQ, and EQ depth is indicated in color.The EQs within or close to the wave sensitive area defined by 5th Fresnel zone are expected to show their effect on the VLF/LF data at the receiver.

Fractal analysis of LF propagation data
Only the amplitude data (no phase data) are used, and Fig. 2a shows an example of raw LF data for the JJY signal observed at KCH on a particular day (14 June 2007).The nighttime data are only used for fractal analysis; that is, UT=12∼17 h  Moshiri) as red stars.The wave sensitive area for each path of JJY-MSR or JJY-KCH is indicated by an elliptic area.EQs (with magnitude greater than 5.5) treated in this paper are plotted.The center of each EQ is its epicenter, and its size is corresponding to its magnitude, EQ depth is indicated in color.
(LT=21∼26 h), because nighttime data are considered to be apparently and strongly influenced by seismogenic effects.
Fractal analysis is applied to the above nighttime LF data, in which we adopt the Higuchi method seeming to be the best method among several methods (Smirnova et al., 2001;Gotoh et al., 2004;Ida and Hayakawa, 2006).The details of this method is given in Ida et al. (2006) or Hayakawa and Ida (2008).
Figure 2b illustrates the plot of log L(k) versus log k for a particular night in Fig. 2a in which 900 data points are used.L(k) is the length of the curve (the amplitude data) divided by k segments (see the details in Gotoh et al., 2004;Ida et al., 2006;Hayakawa et al., 2008).As you see from Fig. 2b, we cannot fit the whole range by a single straight line, but we find that there is a significant difference in the range below and above k=10 (k=10 means 400 s( 7 m)).This behavior is true for all the data.So that, we divide the whole k range below k=10 and above k=10.Accidentally, this k=10 (∼7 m) is found to be the period of boundary between the frequency range of AW(acoustic wave, less than ∼10 m) and AGW (10 m ∼ a few hours).Further requirement to estimate the fractal dimension, is discussed.That is, we can approximate L(k) as follows: L(k)∝k −D , in which D is fractal dimension.In the following we can estimate this D-value, but we adopt the D-value only when the correlation between L(k) and k −D is larger than 0.99 (which is the condition for  fractality).Also, when we have no data, no plot is given to that day and this is given by a thin vertical line in the figures.Further, when it is definite that some artificial noise is present in the data, we disregard the plot, and these plots are indicated by a succession of thin vertical lines just like a hatching in the following figures.

Analysis results of fractal dimensions
Figure 3 illustrates the one-year result of fractal dimension (D) during one year of 2004 for the path from JJY to MSR.(a) refers the frequency range of AW, while (b) refers to that of AGW.Blue lines are the daily plots and the red curve is the running average of D over ±5 days.There are six EQs within the wave sensitive area, with magnitude greater than 5.5.The depth is indicated in color.
Three EQs happened in May, June, and August, but we notice no significant changes in the fractal dimension in Fig. 3a, but there seem to exist some noticeable changes in the D-value for last 3 EQs in November and December, with all nearly the same magnitude of 5.7 to 5.8 as the previous 3 EQs.Blue lines of D are found to exhibit a significant increase in D before an EQ in November, and another peak in D might be associated with the 2nd or 3rd EQ.The increase in D means an increase in the non-stationary noise.
A very different behavior from Fig. 3a is seen in Fig. 3b for AGW range.We can say that at least just around the EQ we notice any peaks in D-value in this AGW frequency range.Corresponding to two EQs in May and June, there exist two significant peaks in D in the end of May and in the middle May.These might be precursors to these EQs.Significant two peaks are observed in the middle of August, which might be highly likely to be related with an EQ in August.Finally, the two maxima in D are seen to appear in late November and in early December, which seem to be closely related to the later 3 EQs.As seen from Fig. 3b, there are present other peaks in early part of the year, but it seems likely that some changes in D take place in possible association with the EQs (mainly as a precursor).First we look at the result for 2004 in Fig. 4. In October 2004 there was a famous Niigata-Chuetsu EQ with magnitude of 6.7 (Hayakawa et al., 2006 have already summarized different electromagnetic phenomena associated with this EQ and it is found that the ionosphere is perturbed before the EQ).It is definite that the fractal dimension (D) for AW range in Fig. 4a shows an increase before this Chuetsu EQ.The behavior of D after Chuetsu EQ cannot be discussed because of inaccurate data.During other time period there is a significant change in D. On the other hand, we discuss the D for AGW in Fig. 4b.When we look at the period before the Chuetsu EQ, we find these peaks in the end of September and in early October.However, we are not sure that these are related with the EQ.
We move on to Fig. 5 for the year of 2005.There are two EQs for this year; In the middle of February and in the middle of October (19 October).At least before the 1st EQ we can find two peaks in D-change.On the other hand, it seems clear that the fractal dimension (D) in the AW range shows a gradual increase before the 19 October EQ in the off-shore of Ibaraki prefecture.On the other hand, Fig. 5b for the AGW range is seen to exhibit no clear effect during the whole one year.The results for the next year of 2006 are not presented here because about one half of the year is characterized by data gaps.
The last year of 2007 is discussed in Fig. 6 for which year there were no EQs and no seismic activity.Of course, there are seen a few peaks in Fig. 6a for the AW range, but it seems that there is no significant change in D. Similarly, Fig. 6b for the AGW range is also found to show no conspicuous changes in D.

Discussion and conclusions
Fractal analysis has been performed, for the first time, to the subionospheric LF amplitude data (nighttime data).Several years of data are used, and the fractal analysis has indicated that the fractal dimension should be investigated separately for the AW and AGW range.As the result of analysis, we come to the conclusion.
1.When the curve fitting is performed (logL(k) vs. logk), we have found that the whole range in k cannot be approximated by one line, but that there is accidentally a boundary with k=10(∼7 min) in the curve fitting.This boundary is physically just as the boundary between AW and AGW.
2. Though there exist several peaks in the fractal dimension in the AW and AGW ranges, we often find some significant changes in the fractal dimension either in the AW or AGW range when we look at the behaviors just around the EQ.But, these peaks tend to appear prior to an EQ.
We here have to point out that the results obtained in this paper are not very distinct regarding the EQ effect in the dynamics of fractal dimension of subionospheric VLF/LF propagation time series.Nevertheless the abnormal behavior is sometimes clearly observed around the dates of EQs (for example, Niigata Chuetsu EQ).By using these results, we try to answer our initial question whether any SOC effect is taking place even in the lower ionosphere as detected by subionosphereic LF propagation.Hayakawa et al. (2009) have investigated the effect of Earth's tides in seismo-electromagnetic effects and have found that this tidal effect is clearly seen in the lithospheric phenomenon such as ULF emissions, but is also (but not so definitely) in the subionospheric VLF/LF data.One of the two plausible candidate of lithosphereionosphere coupling mechanism is due to atmospheric oscillation (AW, AGW) channel (e.g., Molchanov and Hayakawa, 2008).Any fluctuation in some surface parameters (such as pressure, temperature, etc.) induced due to the pre-seismic lithospheric activity under the SOC situation near the EQ hypocenter or near the EQ epicenter, excites the corresponding atmospheric oscillations (AW, AGW), which propagate upward to couple to plasma parameters in the lower ionosphere as detected by subionospheric LF data.The selforganization in the lithosphere might be seen sometimes in the LF data.The results in the present paper might suggest that either AW or AGW is influenced in the form of change in fractal dimension, which lends support to the idea that the ionosphere is, to some extent, influenced by the SOC effect of the lithosphere.
The similar kind of ideas was applied to the in-situ data of electron density and electric field observed by different satellites traversing in the upper ionosphere (Cosmos-900 and Aureol-3) (Molchanov et al., 2004;Hobara et al., 2005).That is, the fluctuations in the plasma density and electric field in the ionosphere have been investigated in terms of the fractal exponent.Then they have come to the conclusion that the connection of the fractal index with seismicity is not so definite, but the weak connection can be explained in teams of the change in AGW forcing during the EQ preparation process.
As the general conclusion, the connection of our fractal analysis for the subionospheric VLF/LF data with seismicity is found to be not very pronounced, but is suggesting some weak connection.However, this does not deny the necessity of performing more physical signal analyses, and the time is now ripe for further application of physical signal analyses to any seismogenic phenomena in order to convince others on the presence of seismogenic effect and to obtain further understanding of the mechanism.

Fig. 1 .
Fig. 1.Relative location of the LF transmitter, JJY (indicated by a blue diamond) and two receiving stations (KCH: Kochi and MSR:Moshiri) as red stars.The wave sensitive area for each path of JJY-MSR or JJY-KCH is indicated by an elliptic area.EQs (with magnitude greater than 5.5) treated in this paper are plotted.The center of each EQ is its epicenter, and its size is corresponding to its magnitude, EQ depth is indicated in color. Fig2(a)

Fig. 2 .
Fig. 2. (a) A typical example of amplitude data (JJY-KCH path) observed on a particular day (14 June 2007).The local nighttime is indicated as well.(b) A plot on the relationship between logL(k) and logk.The whole area in k can be reasonably divided into two parts with a boundary of k=10 (period ∼ 7 min), and each region can be expressed by an appropriate regression line (red and blue area).

Figure
Figure 4a and b are the fractal results for the path from JJY to KCH for the year of 2004.(a) refers to the AW range, and (b) the AGW range.The corresponding results for the same path but different years are plotted in Figs. 5 and 6 for the years 2005 and 2007.First we look at the result for 2004 in Fig.4.In October 2004 there was a famous Niigata-Chuetsu EQ with magnitude of 6.7(Hayakawa et al., 2006 have already summarized different electromagnetic phenomena associated with this EQ and it is found that the ionosphere is perturbed before the EQ).It is definite that the fractal dimension (D) for AW range in Fig.4ashows an increase before this Chuetsu EQ.The behavior of D after Chuetsu EQ cannot be discussed because of inaccurate data.During other time period there is a significant change in D. On the other hand, we discuss the D for AGW in Fig.4b.When we look at the period before the Chuetsu EQ, we find these peaks in the end of September and in early October.However, we are not sure that these are related with the EQ.We move on to Fig.5for the year of 2005.There are two EQs for this year; In the middle of February and in the middle of October (19 October).At least before the 1st EQ we can find two peaks in D-change.On the other hand, it seems clear that the fractal dimension (D) in the AW range shows a gradual increase before the 19 October EQ in the off-shore of Ibaraki prefecture.On the other hand, Fig.5bfor the AGW range is seen to exhibit no clear effect during the whole one year.