Correlations between Earthquake Properties and Characteristics of Possible ULF Geomagnetic Precursor over Multiple Earthquakes

In this study, we improved and adapted existing signal processing methods on vast geomagnetic field data to investigate the correlations between various earthquake properties and characteristics of possible geomagnetic precursors. The data from 10 magnetometer stations were utilized to detect precursory ultra-low frequency emission and estimate the source direction for 34 earthquakes occurring between the year 2007–2016 in Southeast Asia, East Asia, and South America regions. As a result, possible precursors of 20 earthquakes were identified (58.82% detection rate). Weak correlations were obtained when all precursors were considered. However, statistically significant and strong linear correlations (r ≥ 0.60, p < 0.05) were found when the precursors from two closely located stations in Japan (Onagawa (ONW) and Tohno (TNO)) were exclusively investigated. For these stations, it was found that the lead time of the precursor is strongly (or very strongly) correlated with the earthquake magnitude, the local seismicity index, and the hypocentral depth. In addition, the error percentage of the estimated direction showed a strong correlation with the hypocentral depth. It is concluded that, when the study area is restricted to a specific location, the earthquake properties are more likely to have correlations with several characteristics of the possible precursors.


Introduction
A useful earthquake prediction must comprise three elements: (i) precursor detection, (ii) estimation of epicenter location, and (iii) determination of earthquake properties (e.g., hypocentral depth, epicentral distance, and magnitude). Researchers have demonstrated various non-seismic approaches to detect any kind of possible precursory proxy appearing before earthquakes. One of the most popular approaches is by studying anomalous lowfrequency seismo-electromagnetic (seismo-EM) phenomena, i.e., direct lithospheric [1,2], seismo-ionospheric and seismo-atmospheric emissions [3]. In these phenomena, preearthquake seismo-EM emission is assumed to either directly or indirectly originate from rection estimation. The PRA method's effectiveness was improved by introducing a new normalization process, which was then validated through a statistical analysis. A total of 34 earthquakes in Southeast Asia, East Asia and South America occurring between the years 2007-2016 were included in this study.

Instrumentation and Data Acquisition
The geomagnetic field data were acquired from the Magnetic Data Acquisition System (MAGDAS) magnetometer network through the International Center for Space Weather Science and Education (ICSWSE), Kyushu University website (http://magdas2.serc.kyushuu.ac.jp/datausage/index.html). The network comprises a wide array of ring-core fluxgate type magnetometers located around the world [23]. The magnetometers measure parallel to three magnetic components, namely magnetic northward (H), magnetic eastward (D) and downward vertical (Z) [24]. Data from 10 MAGDAS stations (each given a code name) located in Southeast Asia, East Asia and South America were utilized in this study and their details are listed in Table 1. Refer to geographic latitude and longitude respectively throughout this paper. 2 Refers to the year(s) for which data are available.
A total of 34 M ≥ 5.0 earthquakes that occurred between the years 2007-2016 where the epicenters appeared within the range of 180 km from the nearest MAGDAS station were included in this study (data source: www.emsc-csem.org). The minimum magnitude, M and maximum epicentral distance, d (i.e., the distance between the epicenter and the observatory station) were selected based on the empirical threshold of 0.025d ≤ M − 4.5 [25]. The threshold estimates that an earthquake of a magnitude as low as 4.5 potentially generates detectable ULF emission. Hence, M ≥ 5.0 was a reasonable choice since lower magnitude earthquakes occur much more frequently (Gutenberg-Richter law) and typically in swarms that complicate spatial and temporal separation of possible impacts by individual earthquakes. The empirical threshold also predicts that ULF emission generated by an M9.0 (the greatest magnitude in our study) earthquake located 180 km away is detectable. It was decided to use such a distance for all studied earthquakes to maximize the study scope.
In cases where the earthquakes occurred in a sequence, only the mainshock was included. Fundamental earthquake properties such as magnitude, M, hypocentral depth, h, epicentral distance, d, and azimuthal angle of the earthquake epicenter, ϑ EQ (i.e., the angle measured clockwise from the true North at the corresponding station) are listed in Table 2. Additionally, the local seismicity index, K LS [26] was calculated for each earthquake to quantitatively indicate seismicity condition at the station using the following equation: The location of the MAGDAS stations and earthquakes utilized in this study are illustrated in maps in Figure 1. In the maps, the earthquakes are represented as red circles and labeled with numbers, corresponding to the ID numbers as assigned in Table 2. The unbalanced distribution of studied earthquakes across the three regions was due to different frequencies of earthquake occurrence, as well as different magnetometer station densities and data availabilities located near earthquakes.
In order to distinguish possible earthquake precursors from other external sources of variation in the geomagnetic field, ap and Dst geomagnetic indices (unit: nT; sampling period: 1 h) were observed (data source: www.omniweb.gsfc.nasa.gov). Since ap and Dst are recorded by several geomagnetic observatories in mid-latitude/sub-auroral and dip-equatorial respectively [27], they were relevant for our regions of interest. Geomagnetic indices generally represent the global geomagnetic activity level which is affected by solarterrestrial interaction. Values of ap > 27 nT (equals Kp > 4) indicate planetary geomagnetic disturbances [28]; meanwhile, Dst < −30 nT indicates the occurrence of at least weak geomagnetic storms [29]. Therefore, anomalies appearing during these hours must be discarded as they are likely to be due to the global geomagnetic condition [30,31]. That being said, it should be stressed that the geomagnetic indices are moderately correlated only with the G component of the raw geomagnetic field [19]. In contrast, the correlation of the indices with the precursory parameter in this study (i.e., P Z/G ) is greatly diminished due to the processing sequence in the PRA method (refer Section 3.1). Therefore, discarding the anomalies during disturbed periods can be considered as an additional precautionary measure to avoid false precursors.

Precusor Detection
The polarization ratio analysis (PRA) method was first introduced by Hayakawa et al. [32]. The analysis was performed by calculating the ratio of power spectral density (PSD) of vertical (Z) to the PSD of total horizontal (G) field components by using the following formula: where P Z/G is the daily polarization ratio which is taken as the precursory parameter. In the equation, S Z and S G are the PSDs of vertical and total horizontal (G = H 2 + D 2 ) components respectively, averaged over a ULF frequency range, ∆f and over a local nighttime period, ∆T. This step was to minimize human-made noise during the day. Frequency ranges, ∆f in prior studies were arbitrarily chosen [33], hence nine ranges were tested to determine the most effective one for each earthquake: ∆f = 0.01-0.02, 0.02-0.03, . . . , 0.09-0.10 Hz. Similarly, nighttime period, ∆T, also differed across studies; thus, five different periods (∆T = 22-02, 23-03, 00-04, 01-05 and 02-06 LT) were examined for each earthquake to identify the least geomagnetically disturbed period at which precursors appeared. An anomaly is identified when the P Z/G value of a particular day exceeds the mean plus two times the standard deviation (µ + 2σ) of the entire period of observation [15]. If the anomaly appeared during an undisturbed period, then the anomaly was assumed to be a possible precursor. Some earthquakes were preceded by several anomalies in multiple ∆f or/and ∆T; for these cases, the anomaly having the prominently highest amplitude was regarded as the precursor to the corresponding earthquake. Additionally, it is imperative to check that any increment of P Z/G is due to the significant increase of S Z and not the decrease of S G ; the observations of individual S Z and S G are useful for this purpose.
The PRA method has been improved over the years by employing various signal processing techniques such as standardization [34], and fractal dimension [35]. In this study, the method was further improved by introducing a new range-normalization process since means and standard deviations of both components are not directly comparable. After evaluating several range combinations heuristically, it was decided to normalize S G into the range [1,2] (between 1 and 2), and S Z into [1,3] (between 1 and 3) to give a greater weightage to the Z component, since it is the primary indicator of pre-earthquake anomalies. In addition to the µ + 2σ threshold value, the value of P Z/G must also exceed 1.5 to be considered as an anomaly. This is to ensure the reliability of the P Z/G parameter which requires the normalized value of S Z to be the maximum (S Z,max = 3) or one of the maximal values so that, when it is divided by any value of S G , the resultant P Z/G value is maintained at 1.5 or higher. It was found that the combination of the normalization ranges and the newly proposed threshold value was an effective modification to the original PRA method in detecting possible precursors of multiple earthquakes in this study.

Direction Estimation
The polarization ellipse method used in this study was mainly based on the description by Schekotov et al. [3] and Ohta et al. [20]. The major angle, θ, is the angle the major axis of polarization ellipse (red line marked by a in Figure 2) makes with the eastward axis (E) with θ < 0 revolving counterclockwise; its tangent is given by: The H (northward) and D (eastward) geomagnetic components were first bandpass filtered in a frequency band ∆f (determined previously, specific to each earthquake) to obtain quasi-monochromatic signals in the ULF range, U h and U d , respectively. The signals then underwent Hilbert transform to acquire complex signals. From the real and imaginary parts of the complex signals, the instantaneous amplitudes (A h , A d ) and phase angles (ϕ h , ϕ d ) were calculated to be used in Equation (3) to obtain θ.
The seismo-atmospheric emission is assumed to propagate from two possible directions towards the center of the ellipse as shown by the blue arrows in Figure 2. The azimuthal direction of the emission source (wiggly blue lines in Figure 2) is perpendicular to the major axis, hence, the angle between the direction and the northward axis (N) in the clockwise direction can be obtained by Equation (4): where α 1 is the first possible azimuthal direction that exists in the interval of [π/2, 3π/2]. The second possible direction, α 2 , in the interval of [3π/2, π/2] was found by calculating the opposite direction of every α 1 using Equation (5). Both α 1 and α 2 were combined to obtain the complete α in the whole interval of [0, 2π].
In order to ensure that the obtained azimuthal angles are reliable in estimating the source direction, only α(i) which satisfy the following condition were included in the final plotting of azimuth distribution [20]: In Equation (6), the left-hand side is the instantaneous total intensity of horizontal magnetic component while the right-hand side (excluding K coefficient) is its average intensity in a given day. K acts as the minimum requirement coefficient of signal-to-noise ratio, where K = 5 was adopted by prior studies. On the other hand, we found in this study through a heuristic technique that K = 3.5 provides reasonable direction estimations for multiple earthquakes in different regions.

Precursor Detection and Direction Estimation
In this section, the investigation of earthquake ID no. 1 is elaborated and shown in Figures 3 and 4 to exemplify the results of precursor detection and direction estimation. For precursor detection, Figure 3b illustrates a typical observation of P Z/G temporal evolution (black solid line). Values of P Z/G from 45 days before until 15 days after the earthquake were observed to identify statistically significant anomalies. It is to be noted that some data gaps existed during the period of observation for several earthquakes; these gaps are shown by disconnections of line, for example on 17th January 2012. As mentioned previously, the P Z/G parameter value must simultaneously exceed (i) µ + 2σ (black dashed line), and (ii) 1.5, to be considered as an anomaly. In this example, P Z/G (∆f = 0.01-0.02 Hz, ∆T = 01-05 LT) at CEB station satisfied both conditions on 20th January 2012 in the absence of geomagnetic storm and disturbance, where ap and Dst indices in Figure 3a are less than 27 nT and more than −30 nT respectively. During the date, S Z (blue line in Figure 3b) reached the maximum value of 3; meanwhile, S G (red line) remained at a moderate value of 1.718, thus verifying that the anomaly was due to a severe disturbance in the vertical component. This observation aligns with the premise of the PRA method which assumes that seismo-EM emission is more dominant in the vertical component. Therefore, the anomalous deviation can be assumed as the precursor for the 6th February 2012 M6.7 earthquake; the earthquake is represented by the filled magenta triangle symbol in Figure 3a.   Figure 3). The circle sizes indicate the magnitudes of the earthquake. In (c), the angles between the true North (N) and the earthquake of interest, the first ϑ, and the second ϑ are written in magenta, blue and red, respectively.
The direction of the ULF emission source on the date of the precursor appearance was estimated for each earthquake. The 'azimuthal beams' at Cebu (CEB) station on 20 January 2012 is shown in Figure 4c as an example. The so-called azimuthal beam was obtained by classifying the azimuthal angle, α, into 36 sectors around the station; each sector is 10 • wide. In Figure 4c, the darkness of the azimuthal beam represents the number of α falling into that beam; hence, the maximal beams (the darkest) are the most probable directions of the emission source. For each beam, there is an opposite beam, therefore, the maximal beam nearer to the earthquake of interest (outlined magenta circle in Figure 4c) is considered the final estimated direction. Also, the angle which the estimated direction makes with the true North in the clockwise direction is called the estimated azimuthal angle, ϑ. As illustrated in Figure 4c, there were two possible ϑ's: 75 • (blue) and 255 • (red); since ϑ = 255 • is nearer to the earthquake of interest, it is taken as the only ϑ.
The azimuthal beams on the day before (19 January 2012) and the day after (21 January 2012) the date of precursor appearance are also shown in Figures 4a and 4b, respectively, to illustrate the background polarization angle, which is influenced by the local electromagnetic condition and underground electrical composition. The maximal beams on 19 January 2012 pointed to 105 • /285 • ; meanwhile the beams pointed to 95 • /275 • on 21 January 2012. This observation suggests that 95 • -105 • /275 • -285 • were the predominant azimuths of the normal background disturbance, thus the departure from these ranges on 20 January 2012 (Figure 4c) was possibly contributed by external influences such as local lithospheric process related to earthquake.
In order to evaluate the performance of the direction estimation, the error in degree and percentage (∆ϑ and ∆ϑ % , respectively) were calculated using the following formulae: where ϑ EQ is the actual azimuthal angle of the earthquake epicenter as listed in Table 2 (ϑ EQ = 244.59 in this example). By using Equation (7), we obtained that ∆ϑ and ∆ϑ % of earthquake no. 1 equals −10.41 • and 2.892%, respectively. The same calculations were computed for other earthquakes. The histogram of estimated direction error, ∆ϑ, is illustrated in Figure 5 (the values are listed in Table 3), which shows that the majority of ∆ϑs fall within ±30 • indicating that the estimations are sufficiently good with minimal error. The results might also support the association of the ULF emission with the corresponding earthquake. A total of 20 out of 34 earthquakes (58.82%) were preceded by possible precursors. Four precursor characteristics were recorded in Table 3, which are the amplitude of precursor (A), frequency range (∆f ), lead time (τ, i.e., the number of days before the earthquake when the precursor appeared) and percentage error of estimated direction (∆ϑ % ). Additionally, the values of ap and Dst indices during the appearances of the precursor did not exceed the thresholds mentioned in Section 2, indicating that the global geomagnetic activities were low (see Table 3). Earthquake no. 18 which is the 2011 M9.0 Tohoku earthquake was studied previously by Takla et al. [36] in a similar manner in terms of the magnetometer station and method used. The results obtained in this study are in accord with those results regarding the lead time and frequency range. Meanwhile, the results for earthquake no. 10 are rather interesting since τ = 0 days was observed, indicating that the precursor appeared on the same day as the earthquake occurrence. Upon closer inspection, it was revealed that the precursor was caused by a disturbance appearing in filtered (∆ f = 0.03-0.04 Hz) Z and G components simultaneously 36 s after the earthquake, suggesting that the precursor is a result of co-seismic effect.  3 frequency range, 4 lead time, 5 nighttime period, 6 estimated azimuthal angle, 7 estimated direction error, and 8 percentage error of the estimated direction. 9 Maximum ap and minimum Dst during the respective nighttime periods. 10 The precursor appeared on the same day of the earthquake event.
The modified PRA method was validated by performing the superposed epoch analysis (SEA) on the geomagnetic field data during 'no-earthquake' periods to form the control group. The SEA is a statistical tool that effectively superposes multiple binary (true/false) time series observations to reveal periodicities or a temporal correlation in the data [37,38]. In order to perform the analysis, we identified all 61-day periods of data measured by the 10 magnetometer stations (Table 1) throughout the years when no M ≥ 5.0 earthquakes occurred within 180 km from the stations. It was determined that 93 periods were available for the analysis, for which the data were subjected to the modified PRA method as described in Section 3.1. Since the nine frequency ranges, ∆ f and five local nighttime periods, ∆T are characteristically unique, a total of 4185 (=93 × 9 × 5) P Z/G temporal evolutions were obtained. The number of anomalies during geomagnetically undisturbed days were counted in each evolution and then added together. Figure 6 shows the results as a percentage, divided by the total number of observations (n = 4185) so that it can be compared with the percentage of precursor detection during the presence of earthquakes (i.e., the test group). The anomalies which appeared in the control group were less than 0.90% per day throughout the 61-day period with no apparent trend, which suggests that the anomalies were probably random. In total, 5.87% of the observation periods in the control group contained at least one anomaly, which is significantly lower than 58.82% that was observed in the test group. It is concluded that the SEA results provide statistical validity of the modified PRA method in detecting possible earthquake precursors.

Statistics of Precursor Presences
A statistical analysis was performed to observe any dependence of precursor presences on different earthquake properties. Grouping earthquakes based on equal-width class intervals (e.g., M5.0-5.9, M6.0-6.9, etc.) will result in groups of unequal size and may cause biases when calculating percentage of precursor presences (PP % ). Therefore, grouping based on percentile ranks, where earthquakes were divided into three equally sized groups based on their value ranks in terms of all four properties, yields a more accurate statistical representation. Figure 7 shows the distribution of earthquakes that were and were not preceded by precursors (red and white circles respectively) by plotting their epicentral distance, d, against magnitude, M in Figure 7a, and their hypocentral depth, h, against local seismicity index, K LS in Figure 7b. The blue vertical and magenta horizontal lines are the boundaries between the groups and are labelled with corresponding percentile values. The three groups are (i) the 'bottom' group-below the 33rd percentile (P <33 ), (ii) the 'middle' groupbetween the 33rd and 66th percentiles (P 33-66 ), and (iii) the 'top' group-above the 66th percentile (P >66 ).
From Figure 7a, it is apparent that PP % is the highest for both M and d in their respective top groups, where PP % = 83.33% and 90.91%, respectively. The observation regarding M aligns with the assumption that greater-magnitude earthquakes (which also have greater seismicity energy, indicated by higher K LS values) emit higher intensity of electromagnetic emission [3]. This high-intensity emission has a higher chance of being detected by magnetometer stations in the form of geomagnetic variation, thus increasing PP % . On the other hand, the high PP % for more distant earthquakes (indicated by higher d) seems to not agree with previous studies that have been summarized by Hayakawa [25]. However, this observation can be explained by the fact that most earthquakes having higher d in this study coincidentally have higher M, as illustrated in Figure 7a. The same explanation is applicable for the K LS against h plot in Figure 7b, where the majority of deeper earthquakes (higher h) have higher K LS (K LS grows exponentially with M, see Equation (1)). From this statistical analysis, it is deduced that M and K LS are pivotal properties in determining the presence of precursors. For an earthquake having M > 6.7 or/and K LS > ∼ 600, there was an approximately 80% chance of it being preceded by a precursor. Meanwhile, in order to infer the influences of h and d on PP % , a dataset containing earthquakes having similar M and K LS will be required. Figure 7. Distributions of precursor presences in terms of (a) epicentral distance, d against magnitude, M, and (b) hypocentral depth, h against K LS index. Earthquakes were divided into 3 groups based on their percentile ranks, i.e., P <33 , P 33-66 and P >66 . The percentage of precursor presences, PP % is shown for each group.

Correlation Analysis
Since precursors show varying characteristics across different earthquakes, they might possess correlations with the earthquake properties. Therefore, scatter plots of four precursor characteristics (A, ∆f, τ and ∆ϑ % ) against four earthquake properties (M, h, d and K LS ) are shown in Figures 8 and 9. The Pearson correlation coefficient, r, was calculated and written in black in the respective subfigure to identify which characteristic-property pairs have strong correlations, together with corresponding best-fit lines. The values of r were interpreted by following the characterizations for correlation strength suggested by Evans [39]: 0.00-0.19 (very weak), 0.20-0.39 (weak), 0.40-0.59 (moderate), 0.60-0.79 (strong) and 0.80-1.00 (very strong). In addition, the corresponding p-value less than 0.05 (p < 0.05) is considered statistically significant for the null hypothesis (that there is no correlation) to be rejected.  The r values (written in black) were obtained when data points from all stations were included, ranging from |r|= 0.0033 to 0.4480, which show very weak to moderate correlations. These generally weak correlations suggest that there could be a local dependence that affects the characteristics of found precursors. This dependence might be caused by the influences of the observatory station locations on the geomagnetic field components. For instance, a station near the magnetic equator would record a weaker intensity vertical component and stronger horizontal component while a station near the magnetic poles would record contrasting observations. Since the ONW station has the most data points, it was used to see if stronger correlations are exhibited when using data points from this station only; these data points are shown by the blue upright triangles in Figures 8 and 9. Similarly, r values are included in each subfigure, written in blue. We paid attention to the |r| ≥ 0.60 to observe strong or very strong correlations.
The value of r = +0.9562 (very strong) was observed for the τ-M pair (Figure 8c). The positive correlation value suggests that, for greater magnitude earthquakes, the lead time is longer (precursors appear earlier), which corresponds to the study by Ahadi et al. [22]. Similarly, K LS also showed a positive correlation with τ and had r = +0.9709 (very strong), which is the highest, where τ was shown to be longer for higher K LS (Figure 9g). These observations might be due to the greater strain energy accumulated at the focal region, where even if a fraction of this energy is converted into electromagnetic energy, it is sufficiently high to be detected by a near station at any point in time [21]. It is to be noted that the linear correlation between τ-K LS was observed when K LS was scaled logarithmically (due to the positively skewed, heavy-tailed data distribution), which means that the correlation actually corresponds to τ-log 10 K LS . Hypocentral depth, h exhibits two correlations, namely with τ and ∆ϑ % , with r = − 0.7949 (strong) and +0.8378 (very strong) respectively. The τ-h pair (Figure 8g) exhibits a negative correlation which agrees with Ahadi et al. [22], i.e., for earthquakes with hypocenter located deeper underground, the precursor appears later than shallower ones. Additionally, it was shown from ∆ϑ % -h correlation (Figure 8h) that the estimated azimuthal direction tends to have a larger error when the hypocentral depth is greater. This finding also means that shallower earthquakes produce higher reliability of direction estimation, thus offering advantages for the potential victims since shallower earthquakes are typically more destructive. There are other pairs which showed |r| ≥ 0.60; however, their respective p-values did not satisfy the p < 0.05 condition. Figure 10 shows a heatmap visualization of how r (the first number in each box) and p-value (the second number) differ across the properties and characteristics. In addition to all stations and ONW station data points, we also analyzed ONW+TNO stations combination, since the stations are just 100 km apart and located on the same tectonic plate, i.e., Okhotsk plate, hence the earthquake precursors detected in these stations might have similar characteristics. It was found that the correlations between τ-M and τ-log 10 K LS continued to be strong, with r = +0.7177 and +0.7417 respectively (the last row in Figure 10). Besides, lead time, τ is more frequently correlated with earthquake properties than other precursor characteristics, which suggests that it is a dependent characteristic and rendering it useful in predicting the properties of upcoming earthquakes. It is worth mentioning that a finding by Chauhan et al. [21] which mentions a high correlation between A-d, and another finding regarding A-K LS by Schekotov et al. [3], were not observed in this study. It is to note that Schekotov et al. [3] analyzed a different parameter (i.e., ULF depression), therefore their results might be incomparable with those obtained in this study.

Conclusions
In this paper, the correlations between four characteristics of possible geomagnetic precursor (amplitude, A, frequency, ∆f, lead time, τ and percentage error of estimated direction, ∆ϑ % ) and four earthquake properties (magnitude, M, hypocentral depth, h, epicentral distance, d, and local seismicity index, K LS ) have been examined by utilizing vast MAGDAS geomagnetic field data on 34 earthquakes. In achieving the main objective, this study has improved the polarization ratio analysis (PRA) method (for precursor detection) by introducing a new normalization process with the following range combination: S Z = [1,3] and S G = [1,2], while the polarization ellipse method (for direction estimation) has been successfully applied on a ULF range (0.01-0.10 Hz). As a result, 20 possible precursors were found (58.82% detection rate) and characterized. The detection rate is significantly higher than the percentage of anomalies observed during no-earthquake periods (5.87%), suggesting that the detected precursors were not by chance. The percentage of precursor presences was around 80% for earthquakes having M > 6.7 and/or K LS > ∼ 600. The correlation strengths for all stations' data points were found to be very weak to moderate. In contrast, the analysis of the Onagawa, Japan (ONW) station data points showed four pairs with strong to very strong linear correlations: (i) τ-log 10 K LS (r = +0.9709), (ii) τ-M (r = +0.9562), (iii) ∆ϑ % -h (r = +0.8378), and (iv) τ-h (r = −0.7949). Additionally, strong correlations were maintained for τ-log 10 K LS (r = +0.7417) and τ-M (r = +0.7177) when the data points from ONW were combined with those from a near station, i.e., Tohno, Japan (TNO). In conclusion, some earthquake properties are strongly correlated with several precursor characteristics. However, the characteristics are locally dependent and exhibit heterogeneity, which is evidenced by the weak correlations for all stations data, but strong correlations for a single station data. It is recommended that more earthquakes having uniform distribution of properties be studied, in order to possibly discover more appreciable correlations. Data Availability Statement: MAGDAS geomagnetic field data are available on request from Akimasa Yoshikawa (yoshikawa.akimasa.254@m.kyushu-u.ac.jp). Global geomagnetic indices and earthquake data are publicly available at www.omniweb.gsfc.nasa.gov and www.emsc-csem.org, respectively.