Ultra low frequency (ULF) electromagnetic anomalies associated with large earthquakes in Java Island, Indonesia by using wavelet transform and detrended fluctuation analysis

Indonesia is one of the most seismically active regions in the world and mitigation of seismic hazard is important. It is reported that Ultra low frequency (ULF) geomagnetic anomalies are one of the most convincing phenomena preceding large earthquakes (EQs). In this paper we have analyzed geomagnetic data at Pelabuhan Ratu (PLR) (7.01 S, 106.56 E), Sukabumi, West Java, Indonesia, from 1 September 2008 to 31 October 2010. There are twelve moderate–large EQs ( M ≥ 5) within 160 km from the station during the analyzed period. The largest one is the M = 7.5 EQ (depth= 57 km, epicentral distance = 135 km, 2 September 2009) based on EQ catalog of Indonesian Meteorological, Climatological and Geophysical Agency (BMKG). To investigate the ULF geomagnetic anomalous variations preceding all the EQs, spectral density ratio at the frequency range of 0.01± 0.003 Hz based on wavelet transform (WT) and detrended fluctuation analysis (DFA) have been carried out. The spectral density ratio results show the enhancements a few weeks before the largest EQ. The enhancement persists about one week and reaches a maximum on 16 August 2009. At the same time, the result of the DFA presents the decrease of α value. For other EQs, there are no clear increases of the spectral density ratio with simultaneous decrease of α value. When these phenomena occur, the value of Dst index shows that there are no peculiar global geomagnetic activities at the low latitude region. The above results are suggestive of the relation between the detected anomalies and the largest EQ.


Introduction
Recently, electromagnetic phenomena have been considered as a promising candidate for the short-term prediction of large earthquakes (EQs).There are three types of measurement for electromagnetic phenomena.They are (1) passive ground-based observations for lithospheric emissions in a wide frequency range (from DC to microwave frequency); (2) ground-based observations with the use of transmitter signals as active monitoring of seismo-atmospheric and seismoionospheric perturbations; and (3) satellite observations of plasma perturbations, thermal anomalies and radio emissions associated with EQs in the upper atmosphere (e.g., Molchanov and Hayakawa, 2008;Hayakawa, 2009Hayakawa, , 2012)).
Ultra low frequency (ULF) is the frequency range of electromagnetic waves less than 300 Hz.However, in this study, an alternative definition of ULF is usually given as less than 1 Hz.The usage of the ULF range gives an advantage of detecting precursory EQ signatures due to deep skin depth.
Published by Copernicus Publications on behalf of the European Geosciences Union.

F. Febriani et al.: Ultra low frequency (ULF) electromagnetic anomalies
Scientists all over the world have reported that EQ-related ULF electromagnetic phenomena provide a promising tool to understand EQ preparation processing.Some ULF monitoring systems all over the world have registered anomalous ULF electromagnetic signals preceding moderate-large EQs (e.g., Kopytenko et al., 1990;Fraser-Smith et al., 1990;Hayakawa et al., 1996Hayakawa et al., , 2007;;Nagao et al., 2000;Uyeda et al., 2002;Hattori, 2004;Hattori et al., 2004aHattori et al., , b, 2013;;Schekotov et al., 2006;Hirano and Hattori, 2011;Chen et al., 2010;Han et al., 2011;Wen et al., 2012).Hattori (2004) found that the appearance of ULF electromagnetic anomalies prior to an EQ is related to its magnitude and the epicentral distance.He showed that the detectable distance, R, should be 0.025R < M − 4.5, where M is the magnitude of an EQ.According to this result, the detectable distance of ULF anomaly preceding the EQ with magnitude around 7.0 is about 100 km.Though, an exception was observed by Ohta et al. (2007) for the 2004 Sumatra EQ.
The tectonic settings of Indonesia are very complex due to a meeting point of several active tectonic plates.They are the Eurasian Plate, the Australian Plate, the Indian Plate, the Sunda Plate, the Caroline Plate, the Philippine Sea Plate, and the Pacific Plate as shown in Fig. 1 (McCaffrey, 2009).
The existence of the Eurasian lithospheric plate forms the southeastern extremity of Indonesia.The subduction zone around the Eurasian plate is called the Sunda trench.The Indian and Australian Plates move toward the northeast at a rate of about 50-60 mm year −1 relative to the Eurasian Plate and this results in oblique convergence at the Sunda trench.The oblique motion is partitioned into the thrust fault, which occurs on the plate interface and involves slip directed perpendicular to the trench, and the strike-slip faulting which occurs several hundred kilometers to the east of the trench and involves slip directed parallel to the trench.Such a condition makes Indonesia become one of the most seismically and volcanically active regions in the world.The Java Island has about 128.5 million people (as of 2005) and it is the most densely populated area in Indonesia.Therefore, the Java Island is one of the higher risk areas in Indonesia.
One of the greatest EQs in the Java Island is the 2009 Tasikmalaya EQ, which occurred on 2 September 2009.According to the EQ catalog issued by Indonesian Meteorological, Climatological and Geophysical Agency (BMKG), the EQ magnitude was 7.5.This EQ has caused hundreds of casualties, destruction and damage to infrastructure and thousands of buildings.
Mitigating EQ disasters and establishing the short-term EQ prediction in the Java Island are of importance.In this paper, we would like to evaluate whether there is an ULF electromagnetic anomaly preceding large EQs in Indonesia or not, especially in the Java Island.The concept of EQ prediction research in Indonesia is illustrated in Fig. 2.

Methodology
In this paper, the wavelet transform (WT) and detrended fluctuation analysis (DFA) have been applied to the observed data.In this section, we describe the methodology in brief.

Spectral density analysis based on wavelet analysis (WT)
WT becomes a common tool for analyzing a localized variation of power within a time series.By decomposing the time series into time-frequency space, one is able to determine both the dominant modes of variability and how those modes vary with time (Torrence and Compo, 1998).
The WT is useful in analyzing the time series data which might contain nonstationary power at many different frequencies (Daubechies, 1990).In WT definition, it is firstly assumed that there is a time series, x n , which has equal time spacing, δ t , and n = 0 . . .N − 1, where N is the length of the time series.It is also assumed that the time series has a wavelet function, ψ o (η), that depends on a nondimensional "time" parameter, η.To be "admissible" as a wavelet function, this function must have zero mean and be localized in both time and frequency space (Farge, 1992).
The continuous WT of a time series data x is defined as the convolution of x n with a scaled and translated version of ψ o (η) as shown in the following equation.
where * indicates the complex conjugate, s and n are the wavelet scale and the localized time index, n, respectively.Because the wavelet function, ψ o (η), is in general complex, the WT W n (s) is also complex.The transform can be divided into the real part, {W n (s)}, and imaginary part, {W n (s)}, or amplitude, |W n (s)|, and phase, θ; θ = tan −1 {W n (s)} / {W n (s)} .Finally, the definition of the wavelet power spectrum can be described as|W n (s)| 2 .
Since the WT is a bandpass filter with a known response function (the wavelet function), it is possible to reconstruct the original time series using either deconvolution or the inverse filter.This is straightforward for the orthogonal WT, but for the continuous WT it is complicated by the redundancy in time and scale.However, this redundancy also makes it possible to reconstruct the time series using a completely different wavelet function.The easiest way is by using a delta (δ) function (Farge, 1992).
In this case, the reconstructed time series is just the sum of the real part of the WT over all scales as described in the following form: W n s j s 1/2 j .
(2) The factor ψ o (0) removes energy scaling, while the s 1/2 j converts the WT to an energy density.The factor C δ is a constant for each wavelet function and comes from the reconstruction of a δ function from its WT using the function ψ o (η).
We adopt the Morlet wavelet in this study because of its similarity to the Fourier transform (Morlet et al., 1982).The Morlet wavelet consists of a plane wave modulated by a Gaussian, as defined in the equation below, where η is a nondimensional time parameter, ψ 0 is the nondimensional frequency and here taken to be 6 to satisfy the admissibility condition (Farge, 1992).
Because of the finite-length time series, errors will occur at the beginning and end of WT, which will define the cone of influence (COI).By definition, COI is the region of the wavelet spectrum in which edge-effects become important and is defined as the e-folding time for the autocorrelation of wavelet power at each scale.This e-folding time is chosen so that the wavelet power for a discontinuity at the edge dropped by a factor e −1/2 and ensures that the edge effects are negligible beyond this point.For the Morlet wavelet, this e-folding time is defined by e = √ 2s (s is scale).The similarity between the Morlet wavelet and the Fourier transform can be explained by the relation between the wavelet scale, s, and Fourier frequency, f , given in the following formula: For the Morlet wavelet with ψ o = 6, it will give a value of t = 1.03 s, where t is the Fourier period (t = 1/f ).It indicates that the wavelet scale of the Morlet wavelet is almost equal to the period of the FFT.
An analysis method of spectral density ratio is known to be useful in discriminating EQ-related ULF emissions from other noises.This spectral density ratio analyzed was proposed by Hayakawa et al. (1996), who analyzed the 1993 Guam EQ with the use of the ratio (S Z /S H ), for detecting EQ-related ULF emission by three components fluxgate magnetometer.Here, S Z and S H indicate the spectral intensities of vertical and horizontal magnetic field components.

Detrended fluctuation analysis (DFA)
In previous study, some researcher have applied fractal analysis in the ULF geomagnetic studies (e.g., Ida et al., 2006, Ida andHayakawa, 2006;Hayakawa et al., 2008;Telesca and Hattori, 2007;Telesca et al., 2008).In this study, we use DFA.DFA is a method for determining the scaling behavior of data in the presence of possible trends without knowing their origin and shapes.DFA was proposed by Peng et al. (1994) to avoid spurious detection of scaling and correlations that could be artifact of trends and nonstationary.Such trends have to be well distinguished from intrinsic fluctuation of the system in order to find the correct scaling behavior of the fluctuation.In the recent ULF geomagnetic studies, DFA method has been already used by some scientists to know the characteristic of DFA in relation to the ULF anomalous changes preceding EQs (Telesca and Hattori, 2007;Telesca et al., 2008).
This methodology operates on the series of the interevent times x(i), where i = 1, 2, . . ., N and N is the length of the series.We indicate x ave as the mean interevent time.The series is first integrated As the next step, the integrated time series is divided into boxes of equal length, n.In each box a least-squares polynomial y n (k) of degree p is fit to the data, representing the trend order p in that box.Next we detrended the integrated time series y(k) in each box.The root mean square fluctuation of this integrated and detrended time series is calculated by Repeating this calculation over all box sizes, we obtain a relationship between F (n), which represents the average fluctuation as a function of box size, and the box size, n.If F (n)behaves as a power-law function of n, the data present scaling: Under these condition the fluctuation can be described by the scaling exponent α, that represents the slope of the line fitting log F (n) to log n.The estimate of the α value, performed by a least square method, furnishes information about the type of correlations in the interevent series.Furthermore, it quantifies the growth of the root mean square fluctuations F (n).For an uncorrelated seismic sequence, α = 0.5.If there are only short-range correlations, the initial slope may be different from 0.5 but will approach 0.5 for large window sizes.α > 0.5 indicates the presence of persistent long-range correlations, meaning that a large (compared to the average) interevent interval is more likely to be followed by a large one and vice versa.α < 0.5 indicates the anti persistent long-range correlations, meaning that a large (compared to the average) interevent interval is more likely to be followed by a small one and vice versa.α = 1 indicates flicker noise dynamics.Finally, α = 1.5 characterizes processes with Brownian-like dynamics.

Data analysis and EQs
We have installed a three components fluxgate magnetometer at Pelabuhan Ratu (PLR) (7.01 • S, 106.56 • E), Sukabumi, West Java, Indonesia.The place we installed the sensor is the geophysical, climatological and meteorological station which belongs to BMKG.It is about 150 km and 75 km from Jakarta and Bandung, the capital of West Java Province, respectively.
The system has been in operation since September 2007, as the result of mutual international collaboration among Chiba University (Japan), BMKG, and Indonesian Institute of Science (LIPI).The system records three components of geomagnetic field data (the north-south component (X), the east-west component (Y ) and the vertical component of geomagnetic (Z)) and two horizontal components of electric field.The sampling rate is 10 Hz.
We analyzed geomagnetic data, from 1 September 2008 to 31 October 2010.During the analyzed period, twelve EQs (M ≥ 5) took place within 160 km from the PLR station.Table 1 lists EQs based on the EQ data catalogs issued by BMKG. Figure 3 shows a map of the PLR station and EQs in Table 1.
The largest EQ (as shown by the red star in Fig. 3) occurred on 2 September 2009, in Tasikmalaya, West Java, at 07:55 UT or 14:55 LT.The epicenter is located at 8.07 • S, 107.28 • E; the magnitude is 7.5 and the depth is 57 km.Based on the US Geological Survey (USGS) report, the mechanism is a reverse fault.The epicentral distance is about 195 km from Jakarta, and about 100 km from Bandung (USGS, 2009).The EQ killed at least 70 people and gave severe damages to the western Java Island.
Regarding the threshold chart proposed by Hattori (2004), the largest EQ (M = 7.5) on 2 September 2009 is located in the vicinity of the empirical threshold because the largest EQ was situated on about 135 km from the PLR ULF station.After the largest EQ, the seismic activity is like a swarm.We can not discriminate the effect of the mainshock and aftershocks clearly.Therefore, we hypothesized that this largest EQ was most likely the one to give ULF anomalous changes.Then, we have analyzed the three components of geomagnetic data using the spectral density ratio analysis based on WT and DFA.
Theoretically, due to the global ionospheric or magnetospheric activity, the spectral density ratio between vertical and horizontal components of geomagnetic field is very small because upper atmospheric waves are assumed plane waves with vertical penetration into the ground.Then, this ratio is assumed to be close to zero.However, if in a region where there are induction currents or underground activities such as local seismic activities, these activities can induce the horizontal current.According to the Biot-Savart law, this current can generate a certain value of vertical geomagnetic field.If we apply the spectral density ratio analysis for the detection of ULF geomagnetic precursors, the value of spectral density ratio could increase.However, it is essentially difficult to analyze the magnetic data by using a conventional FFT approach under the existence of the intensive transient noises due to man-made or global solar effects.Recently, simultaneous uses of WT and reference data have been effective in dealing with these difficulties (Harada et al., 2003(Harada et al., , 2004)).Instead of FFT, we use WT in this paper.
In this analysis, we have used data only around midnight (16:00-21:00 UT) because at this time there is less contamination of global solar signals and artificial noises.Actually, the data during 16:30-20:30 UT are available for the investigation because of the edge effect of WT.First of all, we have re-sampled the data down to 1 Hz.To ensure the contamination of natural signal as well as artificial noise, we have checked our data day by day.We found that the geomagnetic data are very much disturbed even in nighttime.Compared to the background data, there are many artificial noises with transient and high intensity.Therefore, if there are intense transient signals defined by mean ± 5 σ , where σ is the standard deviation, we remove these signals as noises.
After removing the intense transient changes, we have applied the spectral density analysis based on WT.In the next step, the new daily average and σ values were calculated.The corresponding mean + 2 σ have been chosen as the threshold of an anomaly, because the data are stable.
However, the data with possibly reduced noises still have transient signals so that the result based on previous procedure is not sufficient.Then, to obtain more convincing result, we have applied analysis of the minimum energy of spectral density ratio analysis based on WT. Figure 4 presents how to pick up the minimum energy of spectral density ratio.We divided four hours of data each day into 30-minute segments and compute the energy for each segment.Since Z components are much affected by artificial noises, we have picked up the minimum energy of Z.Then, we have analyzed the spectral density ratio (S z /S x , S z /S y , and S z /S g , where g means the horizontal component of geomagnetic data) by using the same segment of the minimum energy data.For DFA, data analysis was based on the procedure explained in Sect.2.2.Then, we have chosen mean -2 σ as the threshold, where mean and σ are the daily average and its standard deviation of α values.Finally, the results on spectral density ratio and DFA have been compared to the geomagnetic activity of Dst index, which reflects the global magnetic storm at the low-latitude area.Figure 6 is an example of the WT spectrograms of X, Y and Z components on 16 September 2008.Due to the edge effect, the first and the last 30 min of data are not used in the next analysis.We just analyze the data in the rectangle area in Fig. 6.
Figure 7 illustrates the variation of the daily spectral intensity of S x , S y , and S z based on WT at the frequency of 0.01 ± 0.003 Hz from 1 September 2008 to 31 October 2010.It was found that in this frequency band, the daily averaged value of spectral intensity do not show any clear changes associated with EQs.The bottom panel in Fig. 7 is Dst index, which reflects the global magnetic storm level at the lowlatitude area.
Figure 8 shows results of the spectral density ratio analysis at the frequency range of 0.01 ± 0.003 Hz.The top and second panels indicate the results with and without the minimum energy application, respectively.We clearly found broad anomalies (grey-colored regions in Fig. 8) of S z /S y in time, which appear before and after the 2009 largest EQ.Furthermore, it is found that the value of S z /S y (for example see mean + 2 σ in Fig. 8) decreases due to the reduction of artificial noises.The bottom panel indicates Dst value and there are no significant disturbances during the period of enhancement of S z /S y .
In the next step, DFA has been performed to investigate the simultaneous changes with the results of the spectral density ratio analysis.We focus on 30 days before and after the main shock (M = 7.5 EQ), as shown in Fig. 9.The top, second, and third panels in Fig. 9 illustrate the variations of α value, S z /S y , and the minimum energy application of S z /S y , respectively.In each panel, the horizontal line indicate a statistical threshold of the anomaly.Whereas, the bottom panel shows the variation of Dst.DFA result shows that only the anomaly appeared about a few weeks before the largest EQ with the decrease of α value.At the same time, the result of the daily averaged or that of the minimum energy of spectral density ratio show the enhancements, as presented by the second and the third panels of Fig. 9, respectively.This enhancement appeared about a few weeks before the M = 7.5 EQ and reached a maximum on 16 August 2009.These enhancements persist about a few days.When these phenomena occur, the value of Dst index shows no peculiar global geomagnetic activity.For other EQs, there are no clear enhancements with simultaneous decrease of α value for the analyzed period (for ex., EQ 6 and EQ 7 in Fig. 9).The above results could be suggestive of the relation between the detected anomalies and the largest EQ (M = 7.5, depth = 57 km).To compare with the previous studies, the magnitude-distance relationship of the Tasikmalaya EQ is plotted in Fig. 10.  1.The bottom panel is Dst index.

Discussion
In the previous ULF geomagnetic studies, the analysis of spectral density ratio is popular among scientists (Akinaga et al., 2001;Hattori et al., 2002;Hobara et al., 2002;Hirano and Hattori, 2011).Hattori et al. (2002) has performed the spectral density analysis based on FFT to analyze ULF magnetic anomaly for the 1997 Kagoshima EQ, who found the spectral density ratio increased about one month before the first EQ on 26 March 1997.Hirano and Hattori (2011) have also performed the spectral density analysis to investigate ULF magnetic anomaly associated with the 2008 Iwate-Miyagi Nairiku EQ.Instead of FFT, they performed WT with the Morlet wavelet as a mother wavelet.Based on their results, the spectral density ratio increased some weeks preceding the EQ.In this study, our result shows that the spectral density ratio enhancement took place some months before and after the largest EQ (Fig. 8).In order to reduce the global and nonseismic effects, recent analyses have been used with adequate reference data (Hattori et al., 2002;Hirano and Hattori, 2011).In this study, we do not have any remote stations.Therefore, we adopt DFA together with the spectral density ratio analysis for obtaining convincing results.Telesca et al. (2008) have performed the DFA to investigate ULF geomagnetic signals for the 2000 Izu EQ swarm.Based on their results, the relation of log F (n) to log n has some fluctuations for all components (N-S, E-W, and vertical components).Their results suggest a significant nonuniform scaling behavior in ULF geomagnetic data in relation to the occurrence of intense seismic cluster of the swarm.In this study the relation of log F (n) to log n of the vertical component is uniform, which means that α value is constant for each day.Figure 9 shows that the maximum decrease of α value simultaneously occurred with the enhancement of spectral density ratio, either from the result of the daily average or that with the use of the minimum energy.This feature has appeared only in relation to the largest EQ for the analyzed period.That is, simultaneous significant changes in spectral density ratio S z /S y and α value are of importance.The simultaneous use of spectral density ratio analysis and DFA is found to be effective and useful for a single station analysis to detect EQ-related ULF geomagnetic anomalous changes.
The physical mechanism is not well understood at present.However, three possible generation mechanisms have been proposed for ULF electromagnetic anomalies, which are discussed in Molchanov and Hayakawa (2008).The possible mechanisms to generate such ULF variations are electrokinetic effect (Mizutani et al., 1976;Fittermen, 1978Fittermen, , 1979Fittermen, , 1981;;Ishido and Mizutani, 1981;Dobrovolsky et al., 1981;Fenoglio, 1991), microfracturing (Molchanov andHayakawa, 1995, 1998), and conductivity change (Scholz et al., 1973;Draganov et al., 1991;Surkov and Pulipenko, 1999).Our results in Fig. 9 shows that the anomalies associated with the largest EQ did not last for a long time.Changes for a few days are not acceptable by the mechanism of conductivity change, so that either electrokinetic effect or microfracturing is the preferable mechanisms in this case.1.The blue horizontal line in the first panel gives mean −2 σ for a threshold of anomaly in α value, whereas in the middle two panels indicates mean + 2 σ for a threshold of anomaly in S z /S y .The grey-colored areas are corresponding to the anomalies in each panel.
However, this research is the first EQ-related ULF variation research which took place near at the PLR observatory, West Java, Indonesia.The further investigation is definitely needed to know the exact mechanism of ULF generation for the largest EQ in Indonesia.Analyzing much more case studies of the Indonesian EQs by using longer data and some reference stations as well as performing other methods such as direction finding must be definitely done in the future to provide more convincing results.

Conclusions
The purpose of this paper is to evaluate whether there are ULF electromagnetic anomalous variations associated with large EQs occurred in Indonesia, especially in Java Island region or not.The analysis of spectral density ratio based on WT and DFA have been carried out for a single station in this research to detect anomalous variations.We have analyzed twelve moderate-large EQs from 1 September 2008 to 31 October 2010.The largest EQ (M = 7.5, depth = 57 km) during the analyzed period occurred on 2 September 2009 in Tasikmalaya, Indonesia, with aftershocks.Based on the empirical relationship between magnitude and detectable distance, the largest EQ could be preceded by ULF geomagnetic anomalies.By using the spectral density ratio analysis and DFA, the results show the simultaneous anomalous changes a few weeks before the largest EQ.The simultaneous use of spectral density ratio analysis and DFA is found to be effective and useful for a single station analysis to detect ULF geomagnetic anomalies associated with EQs.

Fig. 2 .
Fig. 2. The concept of EQ prediction research in Indonesia.

Fig. 3 .
Fig. 3. Map of the ULF station (PLR) and EQ epicenters occurred during the analyzed period.The red star indicates the main shock (2 September 2009) and orange stars present other EQs in Table1.

Fig. 4 .
Fig. 4.An example how to perform the minimum energy application for the spectral density analysis based on WT.The first, second, and third panels indicate energy of the X, Y , and Z components on 18 September 2009, respectively.The black square corresponds to the period of the minimum energy of Z component.The spectral density ratio with the minimum application is computed from the energy values of X, Y, and Z components in the square.

Fig. 5 .
Fig. 5.An example of the daily geomagnetic variations.The first, second, and third panels indicate X, Y, and Z components on 16 September 2008, respectively.

Fig. 6 .
Fig. 6.An example of the wavelet spectrogram.The first, second, and third panels indicate X, Y, and Z components from 16:00 to 21:00 UT on 16 September 2009, respectively.

Fig. 7 .
Fig. 7. Temporal variations of the spectral intensity based on WT from 1 September 2008 to 31 October 2010.The first, second, and third panels illustrate the spectral intensity of X, Y , and Z components at the frequency range of 0.01 ± 0.003 Hz, respectively.A red dot presents the daily averaged value.A black triangle indicates the occurrence time of the largest EQ (M = 7.5) and other EQs listed in Table 1.The bottom panel is Dst index.

Fig. 8 .
Fig. 8.The temporal variations of the spectral density ratio from 1 September 2008 to 31 October 2010.The first and second panels indicate the daily averaged spectral density ratio S z /S y without and with the use of the minimum energy for the frequency range of 0.01 ± 0.003 Hz.The bottom panel shows the Dst index.The black triangle indicates the occurrence time of EQs listed in Table1.The vertical broken line shows the occurrence time of the largest EQ (M = 7.5).The horizontal blue line in the upper two panels indicates mean ± 2 σ for a threshold of anomaly in S z /S y .

Fig. 9 .
Fig. 9.The temporal variations of the spectral density ratio for ±30 days around the largest EQ.The zero day is the day of the largest EQ.The first panel shows the DFA results (α) for the vertical component.The second and the third panels indicate the spectral density ratio S z /S y without and with the use of the minimum energy at the frequency range of 0.01 ± 0.003 Hz, respectively.The bottom panel indicates the Dst index.The black triangle indicates the occurrence time of EQs listed in Table1.The blue horizontal line in the first panel gives mean −2 σ for a threshold of anomaly in α value, whereas in the middle two panels indicates mean + 2 σ for a threshold of anomaly in S z /S y .The grey-colored areas are corresponding to the anomalies in each panel.

Fig. 10 .
Fig. 10.The empirical relationship between magnitude and epicentral distance (after Hirano and Hattori, 2011).The Tasikmalaya EQ on 2 September 2009 is added with a red open circle.

Table 1 .
List of EQs analyzed in this study.