Relationship of Low-Frequency Variations of Meteorological Parameters Measured at a Network of Stations in Japan to Seismic Regime

1. Introduction

Changes in meteorological parameters associated with the preparation of earthquakes or occurring after seismic events have long been the subject of intensive research. These changes involve not only the atmosphere but also the electromagnetic fields in the ionosphere. In recent decades, these studies have led to the emergence of "new non-seismic precursors" that differ from traditional precursors associated with changes in the seismic regime, such as variations in the slope of the recurrence graph and oriented towards extracting information almost exclusively from seismic catalogues [1]. These precursors are considered to be the result of the interaction of many processes included in the global electric circuit and affecting air ionization and changes in the thermodynamic parameters of the atmosphere [2,3]. The general concept of the lithosphere-atmosphere-ionosphere interaction (LAIC) of the formation of short-term precursors as a result of the interaction of processes in the Earth's atmosphere, ionosphere and magnetosphere in the final stages of earthquake preparation is considered in [4,5]. From the point of view of the formation of precursors in the atmosphere, the key factor in these works is proposed to be considered as the ionization of air and the formation of water condensation nuclei by ions as a result of the action of radioactive radon on atmospheric gases. Radon emanation from the earth's interior occurs through the boundaries of the earth's crust blocks, including through active underwater tectonic faults. The application of this concept to the analysis of earthquake precursors is considered in articles [6,7]. In works [8,9,10], methods for analyzing radon emanation data (time series and spatial distribution) before strong earthquakes and for assessing potentially dangerous seismic areas are proposed based on the use of fractal statistics and power distribution laws.

Changes in the atmosphere are closely related to processes in the ionosphere. The works [11,12,13] study the electron density of the ionosphere, the formation of the distribution of electric potentials in the ionosphere and the features of electromagnetic ULF/ELF radiation before strong earthquakes, including the Tohoku mega-earthquake of March 11, 2011 in Japan. In the work [14] multiparameter observations before two strong earthquakes in Japan were analyzed from the point of view of the occurrence of low-frequency changes in the magnetic field, total electron content, as well as synchronous anomalous changes in meteorological parameters (temperature and humidity), which could be interpreted as precursors of seismic events. Similar studies are presented in the articles [15,16] before a series of earthquakes in New Zealand and Greece. In the work [17] using spectral and wavelet analysis the occurrence of an atmospheric gravity wave before an earthquake (M = 6.7) in India was studied. The occurrence of short-term precursors in China in the form of an anomalous vertical electric field of the atmosphere, due to the ionization of decay products of radioactive radon at different altitudes, was considered in the article [18]. In the work [19] it is shown that in the preparation zone of a large earthquake in Mexico (M=7.4), variations in temperature and relative humidity of the atmospheric air reach peak values within 10 days before the upcoming earthquake due to later disturbances of the total electron content in the ionosphere.

The occurrence of electrical impulses in the atmosphere in the vicinity of seismically active faults with subsequent air ionization and the formation of aerosols and ascending air currents was studied in [20]. Variations in the total electron content in the ionosphere before the earthquake in Haiti (M=7) were analyzed in [21] using parametric time series models and neural networks. The article [22] shows that some climatological anomalies may appear in the epicentral region weeks before large earthquakes. In particular, for several earthquakes in Central Italy, anomalies in temperature, water vapor and ozone content occurred within 2 months before the seismic events. The article [23] presents the results of detecting anomalies in air temperature, relative humidity and pressure that occurred one day before the seismic event in China. The article [24] shows that disturbances in atmospheric circulation in the form of interruption of climatic cyclone trajectories are a trigger for the occurrence of earthquake swarms in Iran. Similar studies on atmospheric processes such as abnormally low pressure, showers, storms and thunderstorms as triggers that provoke earthquakes in Japan are presented in [25].

The influence of the groundwater level and observed precipitation in the area of swarming earthquakes in southeastern Germany on the migration of hypocenters using the pore pressure diffusion model is presented in [26]. The problem of determining the most informative atmospheric and geochemical precursors in Japan for earthquakes with a magnitude of at least 6 is considered in [27]. The possibility of constructing a short-term forecast based on a model of interaction between the lithosphere, atmosphere and ionosphere was investigated in [28] using the example of a retrospective analysis of atmospheric parameters before a strong earthquake (M=7.8) in Mexico. The selection of short-term precursors based on the anomalous behavior of temperature, relative humidity and atmospheric chemical potential was studied in [29,30], including using the natural time apparatus to select critical behavior [29]. Similar studies are presented in [31] using the example of analyzing the processes of preparation for a strong seismic event in Altai. The use of satellite-based surface latent heat flux observations to extract the precursor behavior of large earthquakes (M=7.8, 7.3, 7.0) in Nepal and Japan using the natural time approach is presented in [32]. The analysis of surface latent heat flux together with atmospheric chemical potential for retrospective prediction of the Crete, Greece earthquake is presented in [33]. Multiparameter precursors of different physical nature describing the ionospheric and atmospheric states associated with the seismic event in Japan (M=7.1) are analyzed in [34] using several types of neural networks.

This work continues the study carried out in work [35] on the analysis of the relationships between anomalies of meteorological parameters (pressure, temperature and precipitation), recorded over a long period of time, from November 4, 1993 to September 30, 2024, with a time step of 3 hours, at the Pionerskaya meteorological station in Kamchatka with the seismic regime on the peninsula. In this study, the analysis was performed not for one station, but for a network of meteorological stations on the Japanese islands for a recording time of 52.5 years, from the beginning of 1973 to mid-2025.

2. Initial Data

The time series of temperature, humidity, barometric pressure and wind speed are freely available from https://www.ncei.noaa.gov/data/global-hourly/archive/csv/ . For meteorological stations belonging to Japan, the number of stations with available temperature data is 386, for humidity data the number of stations is 379, for barometric pressure 251 and for wind speed 387. The earliest start of the time series is dated back to 1946, but there are significant gaps in the data before 1973. Since 1973, the number of stations with available data increases sharply. Therefore, time series starting from 1973 with the longest gap duration not exceeding 0.5 year were further selected. The gaps were filled based on the behavior of the data to the left and right of the gap time interval in time segments whose length is equal to the gap length. The time series were then reduced to a uniform time step of 3 hours by averaging the data within successive 3-hour time intervals. As a result, time series of humidity for 111 stations, temperature for 112 stations, wind speed for 112 stations, and atmospheric pressure for 53 stations were obtained. The length of the obtained synchronous time series is 153452 readings with a step of 3 hours from the beginning of 1973 to July 7, 2025 (52.5 years).

For each meteorological parameter, the time series of the weighted average was calculated using the principal component method [36], in which the squares of the components of the eigenvector corresponding to the maximum eigenvalue of the correlation matrix of the multidimensional time series were taken as weights. By construction, the sum of the weights is equal to one. The station weight reflects the contribution of measurements made at the station to the overall average. In Figure 1, the positions of the meteorological stations for which the time series for humidity, temperature, wind speed and atmospheric pressure were used are shown in colored dots. The color of the dots reflects the rank of the station weight, which contributes to the calculation of the weighted average. A total of 6 ranks were used, for the calculation of which the range of change in weights for each parameter was divided into 6 intervals containing the same number of elements. On the maps shown in Figure 1, the remote station MINAMI TORISHIMA (Lat = 24.289697, Lon = 153.979119) stands out; its weight for all parameters belongs to the lowest rank.

Figure 2 shows the graphs of the average weighted time series of humidity, atmospheric pressure, temperature and wind speed. Once again, we note that these graphs represent one of the variants of the first principal component of the original multivariate time series.

3. Local Extremum Points of Envelopes After Wavelet Filtering

Each of the time series presented in Figure 2 was then subjected to decomposition into orthogonal components using a smooth 20th order symlet with 10 vanishing moments [37]. The choice of the optimal orthogonal wavelet was made based on the principle of minimum entropy of the distribution of the squared wavelet coefficients of all orthogonal Daubechies wavelets with a number of vanishing moments from 1 to 10.

The minimum period corresponding to the level of detail of the wavelet decomposition with number $k$ is $\Delta t\cdot 2^k$, where $\Delta t$ is the time step of the time series. The maximum period of the level with number $k$ is $\Delta t\cdot 2^{k+1}$. Since the time step is 3 hours, the range of periods for the 7th level varies from 384 to 768 hours, i.e. from 16 to 32 days.

For wavelet decompositions at level of detail number 7, the graphs of which are presented in Figure 3, their envelopes were calculated using the Hilbert transform [38].

To clarify these issues, we will apply the influence matrix method, which was previously used in [35,39,40,41] to analyze the relationships between the time moments of seismic events and the time points of local extremes of various statistics of meteorological time series, magnetic field fluctuations, earth surface tremors, and the influence of proton flux on the seismic process.

Figure 4. The graphs are presented: (a) - moments of time of 213 earthquakes with a magnitude of at least 6.5; (b)-(e) - amplitudes of the envelopes of the wavelet components of the 7th level of detail of the weighted averages of humidity, atmospheric pressure, temperature and wind speed. The red dots on the graphs (b)-(e) mark the 213 largest local maxima of the envelopes of the 7th level of detail; the blue dots on the graphs (b)-(e) mark the 213 smallest local minima of the envelopes of the 7th level of detail.

Figure 4. The graphs are presented: (a) - moments of time of 213 earthquakes with a magnitude of at least 6.5; (b)-(e) - amplitudes of the envelopes of the wavelet components of the 7th level of detail of the weighted averages of humidity, atmospheric pressure, temperature and wind speed. The red dots on the graphs (b)-(e) mark the 213 largest local maxima of the envelopes of the 7th level of detail; the blue dots on the graphs (b)-(e) mark the 213 smallest local minima of the envelopes of the 7th level of detail.

Figure 4(a) shows the time sequence of seismic events in the rectangular vicinity of the Japanese islands shown in Figure 1, with magnitudes of at least 6.5 for the time interval from the beginning of 1973 to July 7, 2025. Graphs 4(a, b, c, d) show the envelopes of the 7th level of the wavelet decomposition of the weighted average time series of humidity, pressure, temperature and wind speed and indicate the positions of the 213 most expressive local extrema of the envelopes: red dots - local maxima, blue dots - local minima.

4. Measures of Advance of Local Extrema of Envelope Time Moments with Respect to Times of Earthquakes

In what follows, we will be interested in the question of whether there is an advance of the times of the largest local maxima of the envelopes or the times of the smallest local minima of the envelopes of the 7th level of wavelet decompositions in Figure 4 with respect to time moments of earthquakes. In this case, we will separately consider the pairs of event sequences in Figure 4(a) and the largest local maxima of the envelopes of the 7th level of wavelet decompositions in Figure 4(b, c, d, e), red dots, and the smallest local minima of the envelopes in Figure 4(b, c, d, e), blue dots. In total, we will consider 4 pairs of time sequences, in which one of the sequences will always be the sequence of earthquakes with a minimum magnitude of 6.5, shown in Figure 4(a). Clarification of these issues requires also estimating the “reverse” advance and calculating their difference. If the average value of this difference is positive, then there is a trigger effect of the proton flux on seismicity. In addition, the value of the average difference of the lead measures will give a measure of the trigger effect.

To clarify these issues, we will apply the influence matrix method, which was previously used in [35,39,40,41] to analyze the relationships between the time moments of seismic events and the time points of local extremes of various statistics of meteorological time series, magnetic field fluctuations, earth surface tremors, and the influence of proton flux on the seismic process.

Let represent $t_j^{(\alpha )},j=1,\mathrm{...},N_{\alpha };\alpha =1,2$ the moments of time of 2 event streams. In our case, these are:

1) a sequence of moments of time corresponding to the points of the most "expressive" local extrema of the envelopes of wavelet decompositions at the 7th level of detail;

2) a sequence of moments of time of earthquakes with a magnitude of at least 6.5.

Let us represent their intensities as follows:

$${\lambda }^{(\alpha )}(t)=b_0^{(\alpha )}+b_1^{(\alpha )}{g}^{(1)}(t)+b_2^{(\alpha )}{g}^{(2)}(t)$$

(1)

where $b_0^{(\alpha )}\ge 0,b_{\beta }^{(\alpha )}\ge 0$ are the parameters, $g^{(\beta )}(t)$ is the influence function of events $t_j^{(\beta )}$ of the flow with number $\beta$:

$$g^{(\beta )}(t)={\displaystyle {\sum }_{t_j^{(\beta )}<t}\exp (-(t-t_j^{(\beta )})/\tau )}$$

(2)

According to formula (2), the weight of an event with number $j$ becomes non-zero for times $t>t_j^{(\beta )}$ and decays with characteristic time $\tau$. The parameter $b_{\beta }^{(\alpha )}$ determines the degree of influence of the event flow $\beta$ on the flow $\alpha$. The parameter $b_{\alpha }^{(\alpha )}$ determines the degree of influence of the flow $\alpha$ on itself (self-excitation), and the parameter $b_0^{(\alpha )}$ reflects a purely random (Poisson) component of intensity. Let us fix the parameter $\tau$ and consider the problem of determining the parameters $b_0^{(\alpha )},b_{\beta }^{(\alpha )}.$The log-likelihood function for a non-stationary Poisson process is equal to over the time interval [42]:

$$\ln (L_{\alpha })={\displaystyle \sum _{j=1}^{N_{\alpha }}\ln ({\lambda }^{(\alpha )}(t_j^{(\alpha )}))-{\displaystyle \underset{0}{\overset{T}{\int }}{\lambda }^{(\alpha )}(s)ds}},\alpha =1,2$$

(3)

It is necessary to find the maximum of functions (3) with respect to parameters $b_0^{(\alpha )},b_{\beta }^{(\alpha )}$. It can be shown that at the point of maximum of function (3) the condition is satisfied (for the derivation details see [35]):

$$b_0^{(\alpha )}{\displaystyle \frac{\partial \ln (L_{\alpha })}{\partial b_0^{(\alpha )}}}+{\displaystyle \sum _{\beta =1}^2{b}_{\beta }^{(\alpha )}\frac{\partial \ln (L_{\alpha })}{\partial b_{\beta }^{(\alpha )}}}=N_{\alpha }-{\displaystyle \underset{0}{\overset{T}{\int }}{\lambda }^{(\alpha )}(s)ds}$$

(4)

Since the parameters $b_0^{(\alpha )},b_{\beta }^{(\alpha )}$ must be non-negative, each term on the left-most part of this formula is zero at the point of maximum of function (3) – either due to the necessary conditions of the extremum (if the parameters are positive), or, if the maximum is reached at the boundary, then the parameters themselves are zero. Consequently, at the point of maximum of the likelihood function, the equality is satisfied:

$${\displaystyle \underset{0}{\overset{T}{\int }}{\lambda }^{(\alpha )}(s)ds=N_{\alpha }}$$

(5)

Let's substitute the expression $g^{(\beta )}(t)$ from (2) into (5) and divide by $T$. Then we get another form of formula (5):

$$b_0^{(\alpha )}+{\displaystyle \sum _{\beta =1}^m{b}_{\beta }^{(\alpha )}\cdot {\overline{g}}^{(\beta )}={\lambda }_0^{(\alpha )}},{\overline{g}}^{(\beta )}={\displaystyle {\int }_0^T{g}^{(\beta )}(s)ds}/T$$

(6)

where ${\lambda }_0^{(\alpha )}=N_{\alpha }/T$. Substituting $b_0^{(\alpha )}$ from (6) into (3), we obtain the following maximum problem:

$${\Phi }^{(\alpha )}(b_1^{(\alpha )},b_2^{(\alpha )})={\displaystyle \sum _{j=1}^{N_{\alpha }}\ln ({\lambda }_0^{(\alpha )}+b_1^{(\alpha )}\Delta g^{(1)}(t_j^{(\alpha )})+b_2^{(\alpha )}\Delta g^{(2)}(t_j^{(\alpha )}))\to \underset{b_1^{(\alpha )},\hspace{0.17em}\hspace{0.17em}{b}_2^{(\alpha )}}{\max }}$$

(7)

where $\Delta g^{(\beta )}(t)=g^{(\beta )}(t)-{\overline{g}}^{(\beta )}$, subject to the following restrictions:

$$b_1^{(\alpha )}\ge 0,b_2^{(\alpha )}\ge 0,b_1^{(\alpha )}{\overline{g}}^{(1)}+b_2^{(\alpha )}{\overline{g}}^{(2)}\le {\lambda }_{\hspace{0.17em}0}^{(\alpha )}$$

(8)

Having solved the problem (7-8) numerically for a given $\mathit{\tau }$, we can enter the elements of the influence matrix ${\mathit{\kappa }}_{\beta }^{(\alpha )},\alpha =1,2;\beta =0,1,2$ according to the formulas:

$${\mathit{\kappa }}_0^{(\alpha )}={\displaystyle \frac{b_{\hspace{0.17em}0}^{(\alpha )}}{{\lambda }_{\hspace{0.17em}0}^{(\alpha )}}}\ge 0,{\mathit{\kappa }}_{\beta }^{(\alpha )}={\displaystyle \frac{b_{\beta }^{(\alpha )}\cdot {\overline{g}}^{(\beta )}}{{\lambda }_{\hspace{0.17em}0}^{(\alpha )}}}\ge 0$$

(9)

The quantity ${\mathit{\kappa }}_0^{(\alpha )}$ is a part of the average intensity ${\lambda }_{\hspace{0.17em}0}^{(\alpha )}$ of the process with number $\alpha$, which is purely stochastic, the part ${\mathit{\kappa }}_{\alpha }^{(\alpha )}$ is caused by the influence of self-excitation $\alpha \to \alpha$ and ${\mathit{\kappa }}_{\beta }^{(\alpha )},\beta \ne \alpha$ is determined by the external influence $\beta \to \alpha$. The normalization condition follows from formula (9):

$${\mathit{\kappa }}_0^{(\alpha )}+{\mathit{\kappa }}_1^{(\alpha )}+{\mathit{\kappa }}_2^{(\alpha )}=1,\alpha =1,2$$

(10)

As a result, we can determine the influence matrix:

$$\left(\begin{array}{l}{\mathit{\kappa }}_{\hspace{0.17em}0}^{(1)}{\mathit{\kappa }}_{\hspace{0.17em}\hspace{0.17em}1}^{(1)}{\mathit{\kappa }}_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}\\ {\mathit{\kappa }}_{\hspace{0.17em}\hspace{0.17em}0}^{(2)}{\mathit{\kappa }}_1^{(2)}{\mathit{\kappa }}_2^{(2)}\end{array}\right)$$

(11)

The first column of matrix (11) is composed of Poisson fractions of mean intensities. The diagonal elements of the right submatrix of size 2×2 consist of self-excited elements of mean intensity, while the off-diagonal elements correspond to mutual excitation. The sums of the component rows of the influence matrix (11) are equal to 1. The influence matrices are estimated in a certain sliding time window of length $L$ with offset $\Delta L$ and with a given value of the attenuation parameter $\mathit{\tau }$.

When analyzing variations of the components of influence matrices in sliding time windows corresponding to the mutual influence of the analyzed time sequences, the main attention is paid to their local maxima with their subsequent averaging. Below is described a method for processing these local maxima and averaging them for a set of time window lengths changing within specified limits.

1. The minimum $L_{\min }$ and maximum $L_{\max }$ lengths of time windows are selected and $N_L$ - the number of lengths of time windows in this interval. Thus, the lengths of time windows took the values $L_k=L_{\min }+(k-1)\Delta L$, $k=1,\mathrm{...},N_L$, $\Delta L=(L_{\max }-L_{\min })/(N_L-1)$, $N_L=50$.

2. Each time window of length $L_k$ shifted from left to right along the time axis with some offset $\Delta t$. Let us denote by $t_j(L_k)$, $j=1,\mathrm{...},M(L_k)$ the sequence of time moments of the positions of the right windows of length $L_k$. The number $M(L_k)$ of time windows of length $L_k$ is determined by their time offset $\Delta t$. We used a time window offset $\Delta t$ of 0.2 years.

3. For each position of the time window of length $L_k$, the elements of the influence matrix (11) are estimated for a given relaxation time $\mathit{\tau }$ of the model (1), (2), corresponding to the mutual influence of the two analyzed processes. For definiteness, we will consider some one influence, for example, of the first process on the second. As a result of such estimates, we obtain their values in the form $(t_j(L_k),c_j(L_k))$, where $c_j(L_k)$ is the corresponding element of the influence matrix for the position with the number $j$ of the time window of length $L_k$.

4. In the sequence $(t_j(L_k),c_j(L_k))$, we select elements $(t_j^{*}(L_k),c_j^{*}(L_k))$ corresponding to local maxima of values $c_j(L_k)$, that is, from the condition $c_{j-1}(L_k)<c_j^{*}(L_k)<c_{j+1}(L_k)$.

5. We select a "small" time interval of length $\epsilon$ and for a sequence of such time fragments we calculate the average value of local maxima $c_j^{*}(L_k)$ for which their time marks $t_j^{*}(L_k)$ belong to these fragments. Averaging is performed over all lengths of time windows $L_k,k=1,\mathrm{...},N_L$. In our calculations we took the length $\epsilon$ equal to 0.1 years.

5. Optimal Selection of Parameters and Two Advance Mechanisms

The complete set of free parameters of the influence matrix method is: $\mathit{\tau }$, $L_{\min }$, $L_{\max }$, $N_L$, $\Delta t$, $\epsilon$. The most critical for the calculation results is the choice of parameters $\mathit{\tau }$, $L_{\min }$, $L_{\max }$. For each set of values $(\mathit{\tau },\hspace{0.17em}{L}_{\min },\hspace{0.17em}{L}_{\max })$, it is possible to calculate the average value of the lead measure $G(T_{extr}\to EQ)$ between the moments of time $T_{extr}$ of the most expressive local extrema of the envelopes relative to the moments of time of earthquakes and the average value $G(EQ\to T_{extr})$ of the inverse lead measure. We are interested in the difference of the average measures: $\Delta G(\mathit{\tau },\hspace{0.17em}{L}_{\min },\hspace{0.17em}{L}_{\max })=G(T_{extr}\to EQ)-G(EQ\to T_{extr})$. The choice of parameters $(\mathit{\tau },\hspace{0.17em}{L}_{\min },\hspace{0.17em}{L}_{\max })$ was carried out from the condition:

$$\Delta G(\mathit{\tau },\hspace{0.17em}{L}_{\min },\hspace{0.17em}{L}_{\max })\to \max ,{\mathit{\tau }}_{\min }\le \mathit{\tau }\le {\mathit{\tau }}_{\max },L_{\min }^{(\min )}\le L_{\min }\le L_{\min }^{(\max )},L_{\max }^{(\min )}\le L_{\max }\le L_{\max }^{(\max )}$$

(12)

The a priori limits of parameter variation in problem (12) were as follows (in years): ${\mathit{\tau }}_{\min }=0.05$, ${\mathit{\tau }}_{\max }=10$, $L_{\min }^{(\min )}=1$, $L_{\min }^{(\max )}=3$, $L_{\max }^{(\min )}=4$, $L_{\max }^{(\max )}=10$`. Problem (12) was solved numerically, using the random enumeration method, allowing parallel implementation of a large number of independent search flows (in total 3000 random trials for each variant) with the final selection of the best result.

As a result of calculating all the differences in the average lead measures $\Delta G$ for the optimal parameters found from the solution of problem (12), for the variants of sets of pairs of time moment sequences, those sequences of extremum points were selected for which the difference is maximum. The results of such a choice are presented in Table 1, which presents the optimal values of the parameters $(\mathit{\tau },\hspace{0.17em}{L}_{\min },\hspace{0.17em}{L}_{\max })$ in years for all variants of the considered pairs of data and the corresponding maximum values of the difference $\Delta G$.

The results of solving the problem (12) can be divided into 2 groups. The first group is a set of parameters of the model (1), (2), for which the value $\Delta G$ is maximum for each pair of local extremum points corresponding to each of the 4 meteorological quantities. This group corresponds to the points of local minimums of humidity, atmospheric pressure and temperature and the points of local maximums of wind speed. In Table 1, this group of parameters is highlighted in blue. The second group of parameters is made up of the remaining rows of Table 1, which are marked in black. The second group corresponds to the points of local maximums of humidity, pressure, temperature and the points of local minimums of wind speed.

The main difference between the first group and the second is the difference in the optimal values of the relaxation parameter $\mathit{\tau }$ found from the solution of problem (12) for humidity, pressure and temperature. In the first group, these values are 0.98, 0.98 and 0.78, while in the second group they are an order of magnitude greater: 9.74, 9.30 and 9.78.

Such a strong difference suggests that there are two physically different mechanisms of advance. We hypothesize that the first advance with small relaxation times (blue lines in Table 1) is related to the trigger mechanism of cyclone influence on the seismic process, since cyclones are characterized by low pressure, a drop in temperature, and increased wind. As for the second mechanism with large relaxation times (black lines in Table 1), we assume that it is related to the formation of atmospheric earthquake precursors [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34], which are characterized by an increase in humidity and temperature as a result of the effect of radioactive radon on atmospheric gases and air ionization.

In the following, we will agree to call these two mechanisms the trigger and precursor advance. As can be seen from Table 1, the trigger atmospheric effect exceeds the precursor one, which is one of the reasons for the difficulties in detecting atmospheric precursors in practice.

The following Figure 5 shows the graphs of the change in the lead measures for all variants from Table 1. The left column of the graphs (blue) corresponds to the set of parameters from Table 1 that relate to the trigger mechanism of the lead, and the right column of the graphs (black) relates to the mechanism of the occurrence of atmospheric precursors.

6. Average Lead Measures

To highlight the general properties of the change in the lead measures separately for the trigger and precursor mechanisms, we average the curves in the left and right columns of Figure 5 in successive time fragments 0.1 years long and compare them with the graph of the logarithm of the released seismic energy calculated in a sliding time window of 1 year length also with a shift of 0.1 years for all earthquakes. The results of such calculations are presented in Figure 6.

It is of interest to estimate the extent to which the curves of the change in the average lead measures in Figure 6(b, c) themselves additionally lead the seismic energy emissions shown in Figure 6(a). To do this, we calculate the correlation coefficients between the values $LgE(t)$ of the logarithm of the discharged seismic energy (Figure 6(a)) and the average lead measures $C(t)$ (Figure 6(b, c)), where $t$ is the discrete time index. We calculate the correlation coefficients between $C(t)$ and $LgE(t+k)$, where the time shift $k$ changes within certain limits $-k_{\max }\le k\le k_{\max }$, that is, the correlation function. Due to the fact that the correlated values in Figure 6 differ greatly from the Gaussian ones, we will use the Spearman rank correlation coefficient [43]. Figure 6(a) shows the graphs of the Spearman correlation functions $<C(t)\cdot LgE(t+k)>$, the angle brackets mean the calculation of the correlation coefficient. The limits of the change in the time shift $k$ were taken from -3 to +3 years. From the graphs in Figure 7(a) it can be seen that the correlation functions are strongly asymmetric and shifted towards negative $k$ values, which means that the variations in the average measures of advancement of changes in the release of seismic energy are ahead.

In addition, Figure 7(b) presents the graphs of the estimates of the power spectra of the curves of the change in the average lead measures and indicates the periods corresponding to the spectral peaks. The spectral peak with a period of 12.8 years, which is close to the period of solar activity, is noteworthy. This spectral peak is most strongly expressed for the average lead measure in Figure 6(c), corresponding to the second lead mechanism with large values of the relaxation parameter $\mathit{\tau }$, presumably associated with the formation of atmospheric precursors.

7. Discussion

This section is not mandatory but may be added if there are patents resulting from the work reported in this manuscript.

A method is proposed for analyzing a large number of long-term meteorological time series measured at a network of stations to compare their anomalies with the seismic regime. Data processing is based on the sequential application of the principal component method, wavelet decomposition, calculation of envelopes using the Hilbert transform, and application of a parametric model of coupled point processes to the points of local extrema of the envelope amplitudes and the sequence of seismic events.

The method is illustrated by the example of the analysis of data from a network of meteorological stations in Japan for the time period from early 1973 to mid-2025. As a result of the analysis, two mechanisms were discovered for the advance of local extremum points of the amplitudes of the envelopes of the wavelet decomposition of time series of humidity, pressure, temperature and wind speed at the low-frequency 7th level of detail (periods from 16 to 32 days) relative to the time moments of earthquakes with a magnitude of at least 6.5. The mechanisms differ in the values of the optimal relaxation time of the model of interacting point processes by an order of magnitude.

The first mechanism, which was classified as triggering, is characterized by the advance of local minimum humidity, pressure and temperature points and local maximum wind speed points relative to the time of seismic events. Its cause is the triggering effect of cyclones on the occurrence of earthquakes.

The second mechanism, which was classified as a precursor, is characterized by the advance of local maximums of humidity, pressure and temperature and local minimums of wind speed relative to the time of seismic events. Its cause is the occurrence of atmospheric precursors as a result of the impact of radioactive radon, released from faults in the earth's crust before earthquakes, on atmospheric gases and the subsequent ionization of air.

The averaged lead curves for both mechanisms also lead the seismic energy release calculated for all earthquakes, as follows from the asymmetry of the Spearman correlation function.

For the averaged curve of the precursor-type lead measure, a strong periodicity with a period close to the period of solar activity was discovered.