Meteorological Anomalies During Earthquake Preparation: A Case Study for the 1995 Kobe Earthquake (M = 7.3) Based on Statistical and Machine Learning-Based Analyses

Abstract

The purpose of this paper is to discuss the effect of earthquake (EQ) preparation on changes in meteorological parameters. The two physical quantities of temperature (T)/relative humidity (Hum) and atmospheric chemical potential (ACP) have been investigated with the use of the Japanese meteorological “open” data of AMeDAS (Automated Meteorological Data Acquisition System), which is a very dense “ground-based” network of meteorological stations with higher temporal and spatial resolutions than the satellite remote sensing open data. In order to obtain a clearer identification of any seismogenic effect, we have used the AMeDAS station data at local midnight (LT = 01 h) and our initial target EQ was chosen to be the famous 1995 Kobe EQ of 17 January 1995 (M = 7.3). Initially, we performed conventional statistical analysis with confidence bounds and it was found that the Kobe station (very close to the EQ epicenter) exhibited conspicuous anomalies in both physical parameters on 10 January 1995, just one week before the EQ, exceeding m (mean) + 3σ (standard deviation) in T/Hum and well above m + 2σ in ACP within the short-term window of one month before and two weeks after an EQ. When looking at the whole period of over one year including the day of the EQ, in the case of T/Hum only we detected three additional extreme anomalies, except in winter, but with unknown origins. On the other hand, the anomalous peak on 10 January 1995 was the largest for ACP. Further, the spatial distributions of the anomaly intensity of the two quantities have been presented using about 40 stations to provide a further support to the close relationship of this peak with the EQ. The above statistical analysis has been compared with an analysis with recent machine/deep learning methods. We have utilized a combinational use of NARX (Nonlinear Autoregressive model with eXogenous inputs) and Long Short-Term Memory (LSTM) models, which was successful in objectively re-confirming the anomalies in both parameters on the same day prior to the EQ. The combination of these analysis results elucidates that the meteorological anomalies on 10 January 1995 are considered to be a notable precursor to the EQ. Finally, we suggest a joint examination of our two meteorological quantities for their potential use in real short-term EQ prediction, as well as in the future lithosphere–atmosphere–ionosphere coupling (LAIC) studies as the information from the bottom part of LAIC.

1. Introduction

Earthquake (EQ) prediction can be classified into three different categories in terms of the time scale: long-, medium-, and short-term prediction [1]. Among the three, short-term prediction, with a lead time of about one week, is still a very challenging topic in geoscience because there is a lot of demand from society to save human lives and prevent economical losses (e.g., [1]). In order to predict an EQ, it is needless to say that we need to find any kind of EQ precursor, and there has been achieved a huge amount of progress in EQ precursor studies during the last few decades (e.g., [2,3,4]). After extensive studies, we have found that various electromagnetic phenomena do take pace prior to an EQ [5,6]: (i) electromagnetic radio emissions in a wide frequency range from DC (direct current)/ULF (ultra-low frequency), ELF (extremely low frequency), and VLF (very low frequency)/LF (low frequency) to VHF (very high frequency) and (ii) seismogenic perturbations taking place in the atmosphere and ionosphere. The former category of seismogenic effects is obtained with passive observations, while the latter is found from active observations with the use of existing transmitter signals in the VLF/LF and MF (medium frequency)/VHF ranges and also with the data of ionosondes, GPS TEC (total electron content), and satellite observations (e.g., [2,4,7,8]). Based on the extensive investigation of different precursory phenomena, it is becoming our recent consensus that the ionosphere (both upper F region and lower D/E regions) is extremely sensitive to pre-EQ lithospheric activity (e.g., [9,10,11]), which led us to discover the new and attractive concept of lithosphere–atmosphere–ionosphere coupling (LAIC) because of the strong interactions among different layers [12,13,14]. As such, in recent years there have been published a considerable number of papers strongly focused on the elucidation of the mechanism of the LAIC process by making full use of multi-parameter and multi-layer analysis based on satellite- and ground-based measurements [15,16,17,18,19,20,21,22,23,24,25,26,27].

The mechanism of the LAIC process has been investigated extensively during the last 10 years and a few review papers on the LAIC process have been published in recent years (e.g., Conti et al. [28], Picozza et al. [29], and Chen et al. [30]). A few hypotheses of LAIC process have already been proposed in Hayakawa et al. [2,6]; the first is the so-called chemical hypothesis, in which the emanation of radioactive radon, charged aerosols, and/or gases plays the main role, leading to the modification of atmospheric conductivity and the generation of an electric field, thereby driving the variation in ionospheric plasma density [7,31,32]. Additionally, air ionization in this hypothesis leads to the generation of thermal anomalies near the Earth’s surface as the consequence of different physical/chemical processes which seem to be closely related with this paper. The second is the acoustic hypothesis, in which atmospheric oscillations, including atmospheric gravity waves (AGWs) and acoustic waves, are excited by the precursory deformation of ground motion and/or gas emanation or thermal irregularities propagating upwards to the lower and upper ionosphere and leading to perturbations in the ionosphere [33,34,35,36]. The third is the electromagnetic hypothesis, in which electromagnetic waves generated in any frequency range (either in the lithosphere or in the atmosphere) propagate upwards into the ionosphere and magnetosphere, inducing particle precipitation into the upper atmosphere due to wave–particle interactions in the magnetosphere (e.g., [2]). Finally, a fourth electrostatic channel is proposed, based on laboratory experiments, in which positive holes are generated when the ground of interest is stressed by accumulated pressure [37]. These processes have been discussed extensively based on multidisciplinary measurements by different authors (e.g., Ouzounov et al. (Eds), 2018 [3]), but none of the above hypotheses have been evidenced by any definite observational data, necessitating further studies until the process of LAIC is well understood [25].

In this paper we will focus on the simplest meteorological parameters on the Earth’s surface with an emphasis on the use of “open” data. After we recognize that the information on the Earth’s surface and the near-Earth plays a vital role as an essential bottom part of the LAIC, let us review the history of the Earth’s surface remote sensing in the aspect of the use of meteorological parameters in predicting an EQ. Following early works on the application of satellite imaging of the Earth’s surface in the visible and near infra-red (IR) part of spectrum with high resolution, these images were found to be useful to study the lineaments, morphological structures, and tectonic movements in seismically active regions. In comparison to near-IR and visible images, Earth’s surface images as obtained in the thermal IR part of the spectrum are generated due to surface temperature [38,39,40]. Then, paying particular attention to EQ prediction, Tronin et al. (2002) [41] indicated clearly the presence of positive thermal anomalies as based on NOAA/AVHRR satellite thermal images that are associated with the large linear structures and fault systems of the Earth’s crust. Based on long-term (7-year) data, they have indicated temperature anomalies associated with EQs in China and Japan, and the temperature increases range from 3–6 °C about one week before the EQ. This paper provided a potential impact on possible EQ prediction from satellite remote sensing, resulting in the subsequent intensive remote sensing from satellites on various anomalies of mainly meteorological parameters on Earth’s surface, such as air temperature, humidity, OLR (outgoing longwave radiation), and SLHF (surface latent heat flux). As such, many papers have been published with these openly available satellite remote sensing data [42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58].

These meteorological parameters are quite simple to treat, so we want to explore the possibility that such simple and familiar meteorological parameters, such as temperature, humidity, etc., can be effectively utilized for real short-term EQ prediction, preferably using open “ground-based” data with higher temporal and spatial resolutions than the satellite remote sensing data. Those meteorological (thermal) parameters are much easier to treat than the complicated phenomena such as ionospheric perturbations, because they require sophisticated observations and analysis techniques even though they are the main player of LAIC studies [2,4,6,8,15,16,17,18,19,20,21,22,23,24,25,26,27].

Therefore, the present paper aims at proposing the use of simple meteorological parameters and related quantities, as given in Hayakawa et al. (2022) [22] and Schekotov et al. (2023) [59], employing openly available data from a dense network of meteorological observations by Japan Meteorological Agency (JMA). We will check whether such meteorological parameters can be powerful or not for real short-term EQ prediction, first on the basis of the conventional statistical confidence bounds with the use of mean and standard deviations and spatial distributions of the anomaly intensity with more than 40 stations, secondly with the use of machine/deep learning techniques, within the framework of AI (artificial intelligence), and finally by comparing those results. The organization of this paper is as follows. Section 2 deals with the target 1995 Kobe EQ and we discuss again the solar-terrestrial environment around the EQ. In Section 3 we introduce the ground-based meteorological data from AMeDAS (Automated Meteorological Data Acquisition System) of JMA and the two meteorological parameters used in this paper. Section 4 presents the analysis results; initially, those of the conventional statistical analysis, based on the confidence bounds and spatial distributions of both quantities, and secondly, the results of machine/deep learning with a hybrid use of NARX and LSTM, followed by their comparison. Section 5 and Section 6 are the discussion and conclusion, respectively.

2. EQ Studied in This Paper and Solar and Geomagnetic Activity

2.1. Target EQ and AMeDAS Stations

The target EQ of this paper is the famous 1995 Kobe EQ that happened at 5 h 46 m on January (JST) 1995 at the epicenter (geographical coordinates of 34°35.9′ N, 135°02.1′ E) as shown in Figure 1 (red star with notation of EQ) with a M (magnitude) = 7.3 and with a depth of 16 km (e.g., [60]). For this EQ, we discovered the first convincing evidence of ionospheric perturbations four days to one day before the EQ with the subionospheric VLF propagation data from the Omega transmitter at Tsushima to the observatory at Inubo (Hayakawa et al., 1996 [60]); Nagao et al. (2002) [61] also summarized different kinds of electromagnetic EQ precursors for this EQ. Further, we have plotted some AMeDAS meteorological stations close to the EQ epicenter (such Kobe, Himeji, etc.) together with the fault regions possibly related with this EQ.

2.2. Solar-Terrestrial Environment

Figure 2 illustrates the temporal evolutions of geomagnetic activity (Dst and Kp index) and solar radiation flux at the wavelength of 10.7 cm (f10.7) during the whole period of 1 May 1994 through 31 May 1995 (over one year) including the day of the EQ, although we think that, unlike ionospheric parameters, the meteorological parameters in this paper are likely to be much less influenced by solar-terrestrial conditions and perturbations. It is seen from the plot of Dst that we notice some geomagnetic disturbances with Dst below −50 nT only in October and November, 1994, but the geomagnetic activity after the beginning of December 1994 till the day of the EQ was extremely quiet, because |Dst| was less than 20 nT. Further, the solar radiation flux (f10.7), as an indicator of overall solar activity levels, was also found to be less than 100. Taking into account the quiet conditions of the solar-terrestrial conditions just before the EQ, we are ready to investigate short-term EQ precursors within the short-term EQ prediction span of one month before and two weeks after the EQ.

3. Meteorological Parameters from the Japanese AMeDAS Data

In this paper, we want to show the usefulness and importance of the openly available data from ground-based observations of meteorological parameters, even though there have been recently published many papers on the study of these meteorological parameters, such as temperature, humidity, etc., from the “open” data of satellite remote sensing observations [42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58].

3.1. AMeDAS Data Analysis

In Japan, the JMA has established the so-called AMeDAS since 1975 and it has been in operation very regularly. Four meteorological parameters are being measured at each station: (1) precipitation (or rainfall) (in units of 0.5 mm), (2) air temperature (in units of 0.1 °C) and air relative humidity (in %), (3) sunshine duration, and (4) wind direction/speed, and this AMeDAS network includes as many as 1300 stations all over Japan, so those stations are displaced, on average, with an interval of 17 km. These meteorological data are available from the following site of https://www.data.jma.go.jp/risk/obsdl/index.php (accessed on 1 April 2023).

3.2. Meteorological Parameters Used in This Paper

Next, we will study the following two meteorological quantities proposed in Hayakawa et al. (2022) [22] and Schekotov et al. (2023) [59]. The fundamental idea is based on the following supposition. During the EQ preparation phase, radon, Rn may rise to the Earth’s surface more intensely along seismic faults or volcanic fumaroles [62,63,64,65] and the subsequent air ionization leads to a decrease in the air humidity and an increase in its temperature [4]. As such, the first quantity of interest is the simple ratio of temperature (T (in °C) and humidity (relative humidity) (in %)), i.e., T/Hum, because of our expectation of an enhanced seismogenic effect of increasing T and decreasing Hum. The second quantity is ACP (atmospheric chemical potential), which is expressed by the following equation using thermodynamic parameters (Pulinets and Boyarchuk, 2004 [7]):

ACP (in eV) = 5.8 × 10⁻¹⁰ (20 T + 5463)² ln (100/Hum)

(1)

Here, we indicate the reason why we call this the chemical potential: because at the moment of water molecule condensation/evaporation or attachment/detachment to the ions, the latent heat release/absorption is equal to the chemical potential [22]. This ACP value is supposed to show an increase during the process of air ionization by any agents, including galactic rays, etc., and one of them is the increased radon (Rn) release during the pre-EQ seismic activity, as reported by several researchers (e.g., [62,63,64,65,66]).

As mentioned in the Introduction, we have plotted 4 AMeDAS stations near the EQ epicenter in the western part of Japan, Kobe (34°42′ N, 135°13′ E), Osaka (34°41′ N, 135°31′ E), Sumoto (34°19′ N, 134°51′ E), and Himeji (34°54′ N, 134°40′ E), as representative stations close to the EQ epicenter.

4. Analysis Results

4.1. Statistical Analysis Based on the Mean and Standard Deviation

The two meteorological quantities of T/Hum and ACP have been examined for a particular Tokyo EQ with a magnitude of 5.9 (Hayakawa et al., 2022 [22]) and also for Kamchatka EQs (Schekotov et al., 2023 [59]). The first quantity was our proposal, but the second parameter of APC was initially proposed by Pulinets and Boyarchuk (2004) [7] (detailed discussion in [4]) as an integrated parameter of Earth’s surface meteorological changes associated with pre-EQ activity, but there have been very few extensive studies based on ACP.

The basic data of T and Hum are utilized only around midnight, and we choose LT = 01 h, because we have confirmed that daytime is strongly influenced by strong solar radiation and we think that the time period around midnight is considered to be most suitable for our analysis by excluding this solar effect.

First, we will examine the general behavior of these parameters. Figure 3a (upper panel) illustrates the analysis results (or raw data) for T/Hum during our whole period from 1 May 1994 through 31 May 1995, or 6.5 months before and 5.5 months after the EQ (for a total of more than one year). Each vertical black line indicates the current value for each day and the abscissa indicates the date. The day of the EQ (17 January 1995) is indicated by a vertical red line. There are plotted a few colored curves in the figure; the bottom blue curve refers to the mean value (m), green, m + σ (standard deviation), orange, m + 2σ, and red, m + 3σ. Here, the values of m and σ are estimated during 30 days before the current day as in [22]. The same notations are adopted in the upper panel of Figure 4 for the raw data of ACP. The peak on 10 January 1995 is marked with a red circle on its top. Further, we have included the time periods of geomagnetic storms (though not so strong) (light-red boxes) and a typhoon (blue vertical dotted line) for reference.

Because we have a period of more than one year, we can explore the general behavior of the changes of T/Hum (Figure 3a) and ACP (Figure 4a). First, we look at the original data (a) of the two quantities to understand the fundamental characteristics of those quantities. Figure 3a indicates that there exists a strong seasonal dependence, such as enhancements mainly during the summer time of July, August, and September due to the stronger sunshine. Correspondingly, we notice a lot of variability during this period of enhanced values of T/Hum in summer, while, vice versa, less variability during the wintertime. Further, we notice some anomalous and irregular days with enhanced T/Hum, especially in summer (e.g., on 9 September 1994), whose origins are uncertain, such as irregular and nonlinear behavior or a local effect. In sharp contrast to this, Figure 4a suggests that there seems to exist no seasonal variation and this quantity of ACP looks much more stable for the whole-year analysis than the above parameter of T/Hum. Considering the characteristics of these two parameters, it is very desirable for us to use the bottom panel, i.e., detrended values of δ(T/Hum), for the first quantity, by extracting the monthly mean values. However, in the case of ACP, we notice no significant difference in finding statistical anomalies in Figure 4a,b, in which the largest peak is observed on 10 January 1995.

Now we will try to identify any anomalous variations by using the confidence bounds with m and σ. In our previous works it has been natural to take a criterion of m + 2σ when dealing with different physical parameters (e.g., [2,10]), but in this paper we adopt a much stricter criterion of m + 3σ as used conventionally in the field of astrophysics. The use of m and σ is fundamentally based on the assumption that the fluctuation of the relevant parameter follows a Gaussian distribution, so we have checked the distribution of fluctuations of both quantities of T/Hum and ACP for summer (June, July, and August 1994) and winter (December 1994 and January and February 1995) separately because of the presence of a strong seasonal effect only for T/Hum. Figure 5 illustrates the fluctuation distributions of δ(T/Hum) in (a) summer and (b) winter, while the corresponding plots for δ(ACP) are given in Figure 6 for (a) summer and (b) winter. We will explain the simpler case (ACP) of Figure 6a, which illustrates the probability density (ordinate) of fluctuations in blue boxes and the purple thick curve refers to the best-fitted Gaussian distribution. The degree of fitness of the data by Gaussian distribution is estimated by an index of kurtosis (k), and in this case k = 3.25, which means a good Gaussian approximation because a perfect fit can be characterized by k = 3. Figure 6b refers to the same δ(ACP) but for winter and we have again a rather good fit with Gaussian distribution. Therefore, we can conclude that the use of m and σ is very acceptable and reasonable when discussing the behavior of ACP. On the other hand, the situation is considerably different for δ(T/Hum), as seen from Figure 5a,b. Both figures for δ(T/Hum) in summer and winter suggest that the distributions seem to somehow deviate from a Gaussian distribution, as judged from their kurtosis values and some related information on the degree of fitting, and in this case it is preferable to use the median and interquartile. Nevertheless, we have used above and also will continue to still use the m and σ in what follows to keep consistency with the case of ACP, but we have to be cautious with considering the nature of slight deviations from a Gaussian distribution of the fluctuations while discussing T/Hum.

We are interested in short-term EQ prediction, so we pay attention to our time period from about one month before the EQ and two weeks after the EQ. Figure 3b indicates that we can notice a single and a very significant peak exceeding m + 3σ only on 10 January 1995, just one week before the EQ. The value on this anomalous day is indicated by a red circle on the top of the vertical bar. Similarly, we can identify the same anomaly on the same day in Figure 4a,b for ACP, but its value is not exceeding this criterion of m + 3σ, but is well above m + 2σ. As seen from Figure 3a, there are four days exceeding + 3σ; (i) 9 September 1994, (ii) 10 January 1995, (iii) 16 March 1995, and (iv) 24 May 1995, but the origins of those peaks, except for our target of 10 January 1995, are very uncertain, such as temporary climate effects, a local effect, or the effect of non-Gaussian distribution, etc. Nonetheless, the use of this parameter is not so useful as an EQ precursor during summertime, but it will be very powerful during the wintertime (from autumn to spring).

Finally, we can conclude that a significant peak appears in both studied parameters on 10 January 1995, just one week before the EQ, which is marked by the peaks with their tops covered by red circles in Figure 3 and Figure 4. Further, we have checked the meteorological maps around the Kobe area during a few days around 10 January 1995 to try to study whether there existed any meteorological disturbances or not and we have found that there was fair weather without any rainfall not only on 10 January, but also during a few days around the EQ day. Furthermore, the three other anomalous days were also examined regarding δ(T/Hum) and were found to have had fair weather without any significant meteorological disturbances. Finally, it is highly likely that an anomaly on 10 January 1995 is a possible precursor to the EQ. Further, we will apply AI (machine/deep learning) to our data in the following Section, the results of which will be compared with the present statistical results.

4.2. Spatial Distributions of Anomaly Intensities on 10 January 1995

In order to obtain further support to the close relationship of the anomaly on 10 January 1995 to the EQ, we have investigated the spatial distributions of anomaly intensity of both quantities. The anomaly intensity is estimated by the peak value normalized by the standard deviation (σ) at each station and we have used about 40 AMeDAS stations, including the four stations in Figure 1. Figure 7a refers to the spatial distribution of T/Hum, while Figure 7b refers to ACP, in which the studied area is wider than that of Figure 1. The small black dots in the figures represent the locations of AMeDAS stations. The numerical values of 3, 2, etc., are the peak values normalized by the standard deviation (σ), so 3 means 3σ in the figures. It is seen from both figures that the distributions of both quantities are different from each other, and we will describe them one by one. Figure 7a indicates that the area of anomaly intensities exceeding 3σ (3 in the figure) is widely distributed, just encircling the EQ epicenter, and the anomaly intensity decreases at farther distances from the EQ epicenter. For example, the anomaly intensity is very much depleted in the western part of Shikoku Island. As for the case of ACP in Figure 7b, the anomaly intensity exhibits a maximum a little west of the EQ epicenter, but with the same tendency of decreasing in the anomaly intensity at farther distances from the EQ epicenter. As the common feature, the western part of Shikoku island area is weakly perturbed as in both of Figure 7a,b. Probably, this difference might be related with the physical nature of those two quantities, which will be a target of our future work as related to the positions of fault regions.

4.3. Detection of EQ-Related Anomalies Using AI (A Combination of NARX and LSTM Models)

(a)

Combination of NARX and LSTM

We performed basic preprocessing, such as missing value completion, normalization, filtering, creation of delayed data, and data division. In this section, we have conducted a machine/deep learning analysis using the Nonlinear Autoregressive model with eXogenous inputs (NARX) (e.g., [49,67,68,69,70]) and Long Short-Term Memory (LSTM) [71,72,73,74,75,76,77,78,79,80,81] of the same data of the two quantities in the previous subsection. So far, these two models have been widely used, but they have been implemented independently (e.g., [49]), but this paper proposes a combinational use of NARX and LSTM, which is becoming a well-known approach in the AI field, providing a powerful tool for time series forecasting, and we believe that this combination will enable more accurate forecasting by taking advantage of the features of each model. That is, LSTM has an advantage of learning long-term dependencies, which is very effective for capturing the relationship between data that are separated in time. For example, it works effectively when patterns in the distant past have an impact on the future. On the other hand, NARX is a self-regressive model with external inputs and it is good at capturing short-term changes and clear causal relationships by taking advantage of the fact that the current output depends on the past output and external inputs. As a synergy of the LSTM+NARX combination, by combining LSTM and NARX, the strengths of each model can be used in a complementary way. While LSTM learns long-term trends and temporal patterns, NARX can incorporate short-term changes and external factors. NARX explicitly handles past data and external inputs, while LSTM implicitly remembers features. This may improve feature selection and model interpretability and this combination is expected to provide more robust predictions for data with a lot of noise or data that are strongly influenced by external factors. This combination has seen a lot of success in different fields. For example, Massaoudi et al. (2019) [82] have proposed a hybrid model combining a NARX model with an LSTM to be applied to solar power forecasting and they have reported that the proposed model has improved forecasting accuracy compared to other models. Further, Cocianu et al. (2022) [83] have combined a NARX model with an LSTM-based forecasting algorithm and optimized it using an evolutionary algorithm, and they have found that the accuracy of the prediction of financial time series data was improved compared to the standard LSTM model. We can cite some other additional papers on the usefulness of the hybrid use of LSTM and NARX [84,85].

In this paper we obtain the difference between the measured values of T/Hum and ACP and those predicted with the combined model of the NARX/LSTM model, as is illustrated in Figure 8, and perform Bollinger band analysis [86] to evaluate the anomalies related to EQs.

One of the advantages of NARX is that it can directly incorporate external inputs (exogenous variables). As possible inputs, we consider the geomagnetic activity (Dst, Kp) and solar radiation flux (f10.7) and thus we can consider external factors which may affect the forecast target.

(b)

Design of NARX model

We design the structure of NARX and use the past target variables and the lagged exogenous inputs as the features; these are incorporated in the NARX architecture as an input layer of the LSTM model using the NARX architecture. The NARX model is expressed by the following equation.

y(t) = f(y(t − 1), y(t − 2), …, y(t − n), u(t − 1), u(t − 2), …, …, u(t − m)) + ϵ(t)

(2)

where y(t) is the output at a time t to be forecasted (objective variable), y(t − 1), y(t − 2), …, and y(t − n) are the outputs for the past n days (autoregressive part), u(t − 1), u(t − 2), …, and u(t − m) are exogenous variables (input variables) for the past m days, f is a nonlinear function (often using neural networks or other nonlinear models), and ϵ(t) is the error term (noise). The delay in the NARX model is also related to how the model makes use of past information. The number of delays and their settings should be adjusted according to the characteristics of the data and the purpose of the forecast. The selection of appropriate delays is important because the range of information trained by the model and the accuracy of the model change depending on the delay settings. In this analysis, we empirically set the delay = 10 (days).

(c)

Design of LSTM model

The behavior of the LSTM is described by the following Equation (3).

ft = σ_sigmoid(W_f·[h_t−1, x_t] + b_f)

(3)

where ft is the output of the forgetting gate at the current time, W_f is the weight matrix of the forgetting gate, h_t−1 is the hidden state at the previous time, x_t is the current input data, b_f is the bias of the forgetting gate, and σ_sigmoid(.) is the sigmoid function (which limits the gate output between 0 and 1). The forgetting gate determines which information from the previous cell state C_t−1 is retained and which information is forgotten, as shown in Equations (4) and (5).

i_t = σ(W_i · [h_t−1,x_t] + b_i)

(4)

$${\tilde{C}}_{\mathrm{t}}=\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left({\mathrm{W}}_{\mathrm{c}}·[{\mathrm{h}}_{\mathrm{t}-1},{\mathrm{x}}_{\mathrm{t}}+{\mathrm{b}}_{\mathrm{c}}]\right)$$

(5)

where i_t is the output of the input gate at the current time, W_i is the weight matrix of the input gate, b_i is the the bias of the input gate, $\tilde{C}$_t is the new candidate cell state, W_c is the weight matrix of the candidate cell state, and tanh is the hyperbolic tangent function (limits output from −1 to 1). The input gate controls which new information is added to the cell state. It adjusts how much information to add and $\tilde{C}$_t is a new candidate state. The new cell state C_t is updated using the output of the forget gate and the input gate, as in the following Equations.

o_t = σ_sigmoid (W_o · [h_t−1, x_t] + b_o)

(6)

h_t = o_t · tanh (C_t)

(7)

where o_t is the output of the output gate at the current time, h_t is the current hidden state (the output transmitted at the next time), W_o is the weight matrix of the output gate, and b_o is the bias of the output gate. Note that the output gate decides which part of the cell state to use as the output. The final hidden state h_t is computed based on the output gate o_t and the cell state C_t. In the overall LSTM flow, the forgetting gate determines which information is forgotten from the cell state and the input gate determines how much new information is added to the cell state. Next, the cell state is updated and the output gate determines which information is used as the next output. By repeating this sequence of operations and learning long-term dependencies, LSTM becomes a model suitable for forecasting and generating time-series data.

(d)

Model optimization and evaluation

We have attempted to improve the accuracy of the model by optimizing the hyper-parameters of the model (e.g., number of LSTM units and number of delays in NARX) by using a method of Optuna [87]. Optuna is fast compared to other methods such as grid search because it efficiently searches for parameters using Bayesian optimization. The accuracy of the model after training is evaluated in terms of Mean Absolute Error (MAE) and Mean Squared Error (MSE).

(e)

Bollinger Band Analysis

Bollinger band analysis is an analytical technique often used in the financial field [86], but it is also applied in our field to detect anomalous values and fluctuations. By using Bollinger band analysis for geophysical data, it is possible to detect anomalous fluctuations and sudden events. First, some time series data are prepared and Bollinger bands are calculated by moving average. Next, the standard deviation (σ) for the same period is obtained and the upper and lower bands are constructed by adding or subtracting ±2σ or ±3σ from or around the moving average. Finally, the data are plotted together with the data to be analyzed. Applying this method to our data, anomalies can be detected by monitoring periodic or non-periodic data for deviations from the band.

(f)

Analysis Results

To examine short-term variations, we first focused on the period from May 1994 to May 1995 including the data of the Kobe EQ (17 January 1995). Figure 9 illustrates the AI analysis result, in which we plot the deviation of the actual (observed) values from the NARX/LSTM-predicted values. The ordinate refers to the deviation of the observed value from the predicted value and the grey and yellow areas indicate the Bollinger bands of ±2σ and ±3σ. For the sake of comparison, possible time spans of geomagnetic disturbances which might disturb our results are also indicated by light-red boxes. The date of the typhoon that passed near Kobe is also indicated by a blue vertical dotted line. The model was generated using data from the years of 1990 to 1993. As a result, there were four days as mentioned before during the whole period of about one year when the predicted T/Hum exceeded the +3σ Bollinger band, and one of these days (10 January 1995) with a red circle on the top of the value is endorsed as a possible precursor (7 days before the EQ) to the EQ. The abnormality on this day was greater than that on the other three days.

Then, we move on to the corresponding analysis results for ACP in Figure 10, with the same notations as in Figure 7 in the form of the deviation of real values from the AI-predicted values. This figure suggests that an ACP anomaly was confirmed to be observed on the same day of 10 January 1995 during our interest time window of short-term EQ prediction of one month before and two weeks after an EQ. Even though the peak value did not exceed +3σ, we have to emphasize that it was well above +2σ and is approaching +3σ as seen in Figure 10. Another essential conclusion is that this parameter of ACP is an extremely stable indicator of thermal anomalies on the Earth’s surface.

5. Discussion

Two parameters (T/Hum and ACP) have been investigated to understand whether these two meteorological quantities can be used for real short-term EQ prediction. Initially, a combination of two parameters, i.e., the ratio of T over Hum, was suggested in [22,56] in order to have an enhanced response of expected enhanced temperature and decreasing relative humidity [7,22,56]. The anomaly on 10 January 1995 was commonly observed in both quantities. We have to note that there were observed three additional peaks exceeding m + 3σ only in the case of T/Hum, and we have tried to understand their origins, but they are still unknown. However, in the case of ACP, ACP was a stable predictor, such that the peak was the largest on 10 January 1995 during the whole period of over one year. Further studies on the spatial distributions of both quantities have found that the anomaly intensity is enhanced in the vicinity of the EQ epicenter and it decreases at farther distances from the EQ epicenter, which provide further support to the close relationship of the anomaly with the EQ.

We have confirmed by AI-based analyses the presence of the anomalous day of 10 January 1995, just one week before the EQ, within our short-term window of one month before and two weeks after an EQ. The AI-based result has provided an endorsement to the conventional statistical result that this anomaly on 10 January 1995 was a possible precursor to the EQ.

Despite that it was found that this new quantity of T/Hum is very powerful, especially around wintertime (from autumn to spring), this quantity in summertime exhibits a high variability, leading to unknown peaks (or statistical anomalies), probably due to unknown reasons such as the slight deviation of fluctuation distributions from a Gaussian distribution, abnormal natural climate changes, or local effect(s), and this parameter was not so useful in predicting EQs during this period. In reality, the studied Kobe EQ happened in winter, so we were lucky in obtaining a clear anomaly on 10 January 1995, just one week before the EQ. On the other hand, another parameter, the ACP, was found to show a high stability during our whole period of over one year and indicates the presence of a significant precursory peak on the same day of 10 January 1995. Initially, this integrated parameter was proposed in [7], but no thorough investigation of the usefulness of this parameter has been performed so far. However, recently researchers have tried to use it extensively as a predictor of any EQs, but unfortunately only qualitatively [88,89,90]. Principally, the surface latent heat flux (SLHF) created during the abrupt phase transitions of water in the atmosphere can be monitored with the help of the integrated parameter of ACP calculated using the air temperature and humidity. This concept has appeared when considering the process of air ionization by different sources: galactic cosmic rays in the upper atmosphere, thunderstorm discharges in the troposphere, and radon emanating from the Earth’s crust in the near-ground layer of the atmosphere [7]. The present paper might provide a quantitative estimation of the usefulness of ACP for the study of pre-seismic signatures. Our machine/deep learning-based model proved to be very successful in finding the same anomalous day of 10 January 1994 as a possible precursor to the EQ, which provides an endorsement to the result based on the conventional statistical analysis.

Finally, we comment on recent works on different ideas to use meteorological perturbations with major EQs discussed by Freund et al. (2022) [91] and Daneshvar et al. (2023) [92]. Freund et al. (2022) [91] have discussed connections between atmospheric perturbations, e.g., thunderstorm activity and extreme weather conditions, and major EQs along with the LAIC mechanism concerning EQ prediction models, and they have found that a positive correlation exists within 30 days before six major EQs in 2017 in Japan. Further, ref. [92] has provided a systematic assessment of the relationship between climatic variables and major EQs in Iran (2011–2021). These variables include total cloud cover, low cloud cover, total precipitation, SLHF, and total column rainwater, and the combination of a cross-correlation and receiver operating characteristic (ROC) was used to develop the spatial and temporal analytical relationship. They have revealed that an increase in climatic parameters could provide valuable information about impending EQ activity within 8 to 20 days. Considering these recent works together with the present study, we can suggest that different kinds of climatic disturbances can be possible precursors to EQs, despite the mechanism of these anomalies on why and how they are coupled with pre-EQ activity being quite uncertain, but we have to consider them in the context of the LAIC process.

6. Conclusions and Outlook

This paper aimed at suggesting the potential use of simple meteorological parameters based on the Japanese “ground-based” AMeDAS meteorological observations which are available “openly to the public”. The following have emerged from the present study.

(1)

The two parameters of T/Hum and ACP at midnight have been suggested in order to identify any possible precursors to EQs and a case study has been performed for the famous 1995 Kobe EQ (M = 7.3) on 17 January 1995 as an example. The conventional statistical analysis shows that clear precursors are detected on 10 January 1995, just one week before the EQ in both quantities when we pay particular attention to the time window of short-term EQ prediction (one month before and a few weeks after an EQ). However, when we look at a longer period of about one year including the day of the EQ, three additional extreme anomalies appeared, except in winter, only in the case of T/Hum and their origins are still uncertain. However, ACP is a very stable predictor in such a way that the largest anomaly is detected on 10 January during the whole period of over one year. Further studies of the spatial distributions of both quantities (even though they are a little different from each other) have provided a further support to the close relationship of the anomaly with the EQ. Because these quantities may serve as a proxy of pre-EQ radon emanation, the openly available AMeDAS data from JMA with higher temporal and spatial resolutions will be of potential importance as a possible candidate for real short-term EQ prediction in future.

(2)

The above anomaly in both parameters on the same day has been verified with the use of AI machine/deep learning techniques: that is, a hybrid use of NARX and LSTM models, which improves time series prediction accuracy. Further, the AI analysis yielded that the abnormality on 10 January 1995 was greater than that on the other three days. Of course, we recommend the application of these AI techniques to different seismogenic phenomena.

(3)

Further, we suggest a combined use of the above two meteorological quantities (T/Hum and ACP) by taking advantage of each parameter for real short-term EQ prediction.

(4)

The information of meteorological perturbations on the Earth’s surface will, of course, be of potential significance in elaborating future LAIC studies.

(5)

Further statistical research such as more event studies is highly required based on the long-term and multi-station AMeDAS data, leading to the study of the confusion matrix, etc.