Analysis of Meteorological Anomalies in Kamchatka in Connection with the Seismic Process

1. Introduction

In the study of the interaction of processes at the lithosphere-atmosphere boundary, one of the intriguing questions is the study of the relationship between the variability of meteorological conditions in the near-surface conditions of the Earth with the final stages of preparation of strong earthquakes. Consideration of this question is especially relevant in relation to the problem of forecasting earthquakes based on possible meteorological precursors.

A number of publications have demonstrated the effect of increasing temperature and decreasing humidity of the surface air of the atmosphere before strong and moderate earthquakes [1,2,3,4,5,6]. An explanation of this phenomenon is given within the framework of the model of complex relationships in the lithosphere-atmosphere and ionosphere-magnetosphere (LAIMC) system [7]. According to LAIMC, at the final, relatively short-term stage of earthquake preparation, thermal anomalies in the atmosphere at low altitudes may form due to the release of latent heat of evaporation during condensation of water vapor on ions due to the ionization of air molecules by alpha particles emitted by nuclei of radioactive gases, mainly radon, entering the atmosphere from the upper horizons of the earth's crust.

Verification of the LAIMC model requires the presence of indisputable data on the occurrence of anomalous radon flows from the Earth's surface on the eve of strong earthquakes or on establishing a connection between data from long-term meteorological observations in seismically active regions, in particular data on surface air temperature, with earthquakes that have occurred.

Observations of the volumetric activity of radon in near-surface conditions are usually of an experimental nature. In addition to reports of the phenomenon of an increase in the volumetric activity of soil gaseous radon, there are also data on a decrease in its activity before earthquakes [8]. Therefore, the question of anomalous discharges of gaseous radon from the lithosphere into the atmosphere before strong earthquakes remains open to this day.

A promising direction in studying processes at the lithosphere-atmosphere boundary during earthquake preparation is the use of long-term meteorological observations in seismically active regions. Such observations are carried out on the basis of metrologically supported methods for issuing weather forecasts, assessing the state of the natural environment, etc. In the work of the authors [9], using observation data from meteorological stations in the seismically active region of the Kamchatka Peninsula, the relationship between increased and decreased values of air temperature and atmospheric pressure, as well as their contrasting changes, at the final stages of preparation of local strong earthquakes accompanied by noticeable tremors was analyzed. In this work, an empirical method of comparing their average daily values with the behavior of annual average seasonal functions was used to identify anomalous changes in meteorological parameters. In addition, formalized statistical methods were used to identify anomalies in the behavior of three statistics characterizing the variability of time series of air temperature and atmospheric pressure in a sliding time window of 112 days with a step of 1 day. Using empirical and statistical methods for processing time series of meteorological parameters, anomalous effects were identified at time intervals of 7, 30 and 112 days and their seismic forecasting significance in relation to the strong earthquakes that occurred. The detection of the random nature of the manifestation of various types of meteorological anomalies before the earthquakes under consideration, as well as the absence of a systematic manifestation of the effect of an abnormal increase in air temperature before earthquakes, cast doubt on the realism of generating thermal near-surface anomalies of atmospheric air at the short-term stage of preparation of strong seismic events in the Kamchatka region.

It should be noted that in the mentioned work [9], the assignment of time intervals of 7, 30 and 112 days, during which the manifestation of a relationship between surface meteorological anomalies and earthquake preparation processes was assumed, was made on the assumption of a relatively short-term manifestation of thermal effects before strong earthquakes. While the study of various types of earthquake precursors in the Kamchatka Peninsula area shows that the characteristic durations of their manifestation before earthquakes range from days-months to the first years [10,11,12]. This demonstrates the need to develop methods for analyzing the relationship between meteorological parameters and the seismic process, allowing such a relationship to be analyzed in a wide range of periods of its manifestation, including the months and years preceding strong Kamchatka earthquakes.

This paper continues the study of the relationship between anomalies in time series of meteorological parameters and the seismic process in Kamchatka, the first results of which were presented in [9]. In this paper, the Huang method was used to decompose long-term meteorological time series, which allows constructing a sequence of orthogonal oscillations for each series, which are called "internal oscillation modes" in the method. The Hilbert transform of the sequence of these oscillation modes makes it possible to obtain sequences of instantaneous amplitudes and frequencies. The paper studies the relationship between the time points of maxima of instantaneous amplitudes of meteorological parameter series and the time moments of strong earthquakes using a parametric model of interacting point processes. The result of the study is the discovery of the effect of leading the time points of maxima of instantaneous amplitudes of air temperature of the time moments of earthquakes at the low-frequency level of the Huang decomposition.

2. Initial Data

The work used meteorological observation data and data on earthquakes that occurred in the Kamchatka Peninsula region for the period from November 4, 1996 to September 30, 2024 (10193 days). A total of 27 years 11 months (27.9 years).

Meteorological Data

The data of observations of atmospheric pressure, air temperature and precipitation at the Pionerskaya meteorological station of the Kamchatka Territorial Administration for Hydrometeorology and Environmental Monitoring (53.08° N, 158.55° E) from the database of the POLYGON Information System of the Kamchatka Branch of the Geophysical Service of the Russian Academy of Sciences were used [13] (accessed October 1, 2024).

At the Pionerskaya meteorological station, a meteorological mercury thermometer with a division value of 0.2°C is used to measure the air temperature. It is located in a protective wooden booth at a height of 2 m above the earth's surface. Atmospheric pressure was measured with a mercury barometer. Precipitation was measured with a precipitation vessel installed on a horizontal surface at a height of 2 m and protected from gusts of wind. The amount of precipitation is the height of the water layer in mm formed by rain, drizzle, melted snow and hail.

The sampling step of meteorological observations and the corresponding discretization of time series of meteorological parameters is 3 hours (Figure 1).

Earthquake Data

To characterize the seismic process in the Kamchatka Peninsula area during the period from November 4, 1996 to September 30, 2024, a sample of 418 earthquakes with local magnitudes that occurred in the depth range from 0 to 635 km in the responsibility area of the Kamchatka Branch of the Geophysical Survey of the Russian Academy of Sciences (KB GS RAS) was used. The earthquake sampling was carried out from the Kamchatka and Commander Islands Earthquake Catalogue: (http://www.gsras.ru/new/infres/, [ID 85], accessed October 1, 2024) using the Unified Information System of Seismological Data KB GS RAS [14] and the Interactive Map of the Distribution of Earthquake Epicenters from the Kamchatka and Commander Islands Earthquake Catalogue (see https://sdis.emsd.ru/map/).

Figure 2 shows the epicenters of all 418 earthquakes, the boundaries of the KB GS RAS area of responsibility, and the location of the Pionerskaya weather station.

In the catalog of Kamchatka and the Commander Islands compiled by KB GS RAS, the representative energy class of earthquakes (decimal logarithm of seismic energy output) for the entire responsibility zone is $K=11.5$ and corresponds to the local magnitude $M_L=5.0$. To obtain the studied sample of earthquakes, earthquakes with $K=12.5$ for the considered time period were selected from the specified catalog. The magnitudes $M_L$ in the sample of earthquakes were recalculated from the values of the energy classes $K$ using the formula $M_L=0.5K-0.75$ [15]. Thus, setting the minimum magnitude $M_L=5.5$ obviously ensured the representativeness of the resulting sample of 418 earthquakes for the entire territory of the KB GS RAS responsibility zone.

The 418 earthquakes under consideration have magnitudes $M_L=5.5-7.8$. The table presents data on the strongest earthquakes that occurred according to KB GS RAS data and data on their moment magnitudes $M_W$ from the USGS/NEIS catalog.

Table 1. Data on the strongest $M_W\ge 6.5$ earthquakes that occurred in the KB GS RAS area of responsibility according to: https://sdis.emsd.ru/info/earthquakes/catalogue.php , https://earthquake.usgs.gov/earthquakes/search/.

Date yyyymmdd	Time hh:mm:ss	Latitude, φ, N °	Longitude, λ, E °	Depth, km	ML	R, km	MW NEIC	L, km	R/L
20030616	22:08:02	55.30	160.34	190	6.6	331	6.9	50	6.7
20031205	21:26:14	55.78	165.43	29	6.7	533	6.7	41	13.1
20040610	15:19:55	55.68	160.25	208	6.7	370	6.9	50	7.4
20060420	23:24:28	60.98	167.37	1	7.1	1019	7.6	100	10.2
20081124	9:02:52	53.77	154.69	564	6.8	623	7.3	74	8.4
20121116	18:12:39	49.06	155.87	78	6.7	486	6.5	33	14.5
20130228	1:05:48	50.67	157.77	61	6.8	271	6.9	50	5.4
20130301	13:20:49	50.64	157.90	62	6.8	272	6.5	33	8.1
20130524	5:44:47	54.75	153.79	630	7.8	726	8.3	199	3.6
20131001	3:38:19	52.88	153.34	608	6.8	700	6.7	41	17.1
20160130	3:25:08	53.85	159.04	178	7.1	200	7.2	67	3.0
20170329	4:09:22	56.97	163.22	43	6.8	521	6.6	37	14.1
20170602	22:24:47	53.99	170.55	32	6.6	795	6.8	45	17.6
20170717	23:34:08	54.35	168.90	7	7.3	692	7.7	110	6.3
20181013	11:10:20	52.53	153.87	499	7.0	592	6.7	41	14.5
20181220	17:01:54	54.91	164.71	54	7.3	448	7.3	74	6.1
20200325	2:49:20	49.11	158.08	48	7.7	438	7.5	90	4.9
20230403	3:06:56	52.58	158.78	105	6.6	120	6.5	33	3.6
20240817	19:10:25	52.79	160.38	43	7.0	126	7.0	55	2.3

Note: R is the hypocentral distance to the Pionerskaya weather station, $L$ is the size of the earthquake source, calculated using the formula $L={10}^{(0.43M_W-1.27)}$ [17].

Table 1. Data on the strongest $M_W\ge 6.5$ earthquakes that occurred in the KB GS RAS area of responsibility according to: https://sdis.emsd.ru/info/earthquakes/catalogue.php , https://earthquake.usgs.gov/earthquakes/search/.

Date yyyymmdd	Time hh:mm:ss	Latitude, φ, N °	Longitude, λ, E °	Depth, km	ML	R, km	MW NEIC	L, km	R/L
20030616	22:08:02	55.30	160.34	190	6.6	331	6.9	50	6.7
20031205	21:26:14	55.78	165.43	29	6.7	533	6.7	41	13.1
20040610	15:19:55	55.68	160.25	208	6.7	370	6.9	50	7.4
20060420	23:24:28	60.98	167.37	1	7.1	1019	7.6	100	10.2
20081124	9:02:52	53.77	154.69	564	6.8	623	7.3	74	8.4
20121116	18:12:39	49.06	155.87	78	6.7	486	6.5	33	14.5
20130228	1:05:48	50.67	157.77	61	6.8	271	6.9	50	5.4
20130301	13:20:49	50.64	157.90	62	6.8	272	6.5	33	8.1
20130524	5:44:47	54.75	153.79	630	7.8	726	8.3	199	3.6
20131001	3:38:19	52.88	153.34	608	6.8	700	6.7	41	17.1
20160130	3:25:08	53.85	159.04	178	7.1	200	7.2	67	3.0
20170329	4:09:22	56.97	163.22	43	6.8	521	6.6	37	14.1
20170602	22:24:47	53.99	170.55	32	6.6	795	6.8	45	17.6
20170717	23:34:08	54.35	168.90	7	7.3	692	7.7	110	6.3
20181013	11:10:20	52.53	153.87	499	7.0	592	6.7	41	14.5
20181220	17:01:54	54.91	164.71	54	7.3	448	7.3	74	6.1
20200325	2:49:20	49.11	158.08	48	7.7	438	7.5	90	4.9
20230403	3:06:56	52.58	158.78	105	6.6	120	6.5	33	3.6
20240817	19:10:25	52.79	160.38	43	7.0	126	7.0	55	2.3

Note: R is the hypocentral distance to the Pionerskaya weather station, $L$ is the size of the earthquake source, calculated using the formula $L={10}^{(0.43M_W-1.27)}$ [17].

3. Research Method

3.1. Empirical Mode Decomposition

Let $x(t)$ be some analyzed time series with discrete time index.

Empirical mode decomposition (EMD) [18,19] represents the decomposition of the signal into empirical oscillation modes:

$$x(t)={\displaystyle \sum _{j=1}^n{h}_j(t)+r_n(t)}$$

(1)

where $h_j(t)$ is the j-th empirical mode, $r_n(t)$ is the remainder, $n$ is the number of empirical modes.

The algorithm for decomposing into a sequence of empirical modes is iterative for each level $j$. We denote by $m$ the iteration index: $m=0,1,\mathrm{...},M_j$, where $M_j$ is the maximum number of iterations for level $j$. The iterations are described by the formula

$$h_{\hspace{0.33em}\hspace{0.33em}j}^{(m+1)}(t)=h_{\hspace{0.33em}\hspace{0.33em}j}^{(m)}(t)-z_{\hspace{0.33em}j}^{(m)}(t),$$

(2)

Here $z_{\hspace{0.33em}j}^{(m)}(t)=\left(p_{\hspace{0.33em}j}^{(m)}(t)+q_{\hspace{0.33em}j}^{(m)}(t)\right)/2$, where $p_{\hspace{0.33em}j}^{(m)}(t)$ and $q_{\hspace{0.33em}j}^{(m)}(t)$ are the upper and lower envelopes for the signal $h_{\hspace{0.33em}\hspace{0.33em}j}^{(m)}(t)$, which are constructed using spline interpolation (usually a 3rd order spline) over all local maxima and minima of the signal $h_{\hspace{0.33em}\hspace{0.33em}j}^{(m)}(t)$.

Iterations (2) are initialized with a zero step for the first level ($j=1$) of the expansion $h_{\hspace{0.33em}\hspace{0.33em}1}^{(0)}(t)=x(t)$. Then the upper and lower envelopes $p_{\hspace{0.33em}1}^{(0)}(t)$ and $q_{\hspace{0.33em}1}^{(0)}(t)$ are found, the average line $z_{\hspace{0.33em}1}^{(0)}(t)$ is calculated and $h_{\hspace{0.33em}\hspace{0.33em}1}^{(1)}(t)$ is found using formula (2). For $h_{\hspace{0.33em}\hspace{0.33em}1}^{(1)}(t)$ the upper and lower envelopes $p_{\hspace{0.33em}1}^{(1)}(t)$ and $q_{\hspace{0.33em}1}^{(1)}(t)$ are determined and the average line $z_{\hspace{0.33em}1}^{(1)}(t)$, $h_{\hspace{0.33em}\hspace{0.33em}1}^{(2)}(t)$ is found, and so on, up to the last iteration index $M_1$, after which the first empirical mode $h_1(t)$ is considered to be found.

The condition for stopping the iterations is usually chosen in the form of the inequality:

$${\displaystyle {\sum }_k{(h_{\hspace{0.33em}\hspace{0.33em}j}^{(m+1)}(t)-h_{\hspace{0.33em}\hspace{0.33em}j}^{(m)}(t))}^2}/{\displaystyle {\sum }_k{(h_{\hspace{0.33em}\hspace{0.33em}j}^{(m)}(t))}^2}\hspace{0.17em}\le \hspace{0.17em}\hspace{0.33em}\delta$$

(3)

where $\delta$ is some small number, such as 0.01. Once the mode $h_j(t)$ is found, an iterative process is started to determine the empirical mode $h_{j+1}(t)$ of the next level. This process is initialized by the formula for the initial iteration index $m=0$:

$$h_{j+1}^{(0)}(t)=x(t)-h_j(t),$$

(4)

According to formula (4), the high-frequency part is subtracted from the original signal and the resulting, lower-frequency signal is considered as a new signal for subsequent decomposition. The construction of empirical oscillation modes continues until the number of local extrema becomes too small to be used to construct envelopes.

As the empirical mode level number $j$ increases, the signals $h_j(t)$ become increasingly low-frequency and tend to an unchangeable form. The sequence $h_1(t),\hspace{0.33em}{h}_2(t),\mathrm{...},h_n(t)$ is constructed in such a way that its sum gives an approximation to the original signal $x(t)$, which can be represented as (1) [18,19].

Empirical modes of oscillations, known as Intrinsic Mode Functions (IMF), are orthogonal to each other, thus constituting a certain empirical basis for the decomposition of the original signal. In what follows, the decomposition levels will be called IMF levels.

In practical implementation of the method, technical difficulties arise due to edge effects, since the behavior of envelopes beyond the first and last points of local extrema is ambiguous. There are several approaches to overcome this difficulty, in particular, a mirror continuation of the analyzed sample back and forth over a sufficiently long period of time. It was used in this work.

3.2. Ensemble Empirical Mode Decomposition

One of the key drawbacks of the EMD method is the problem of mode mixing, which occurs when one empirical mode includes signals of different scales or when signals of the same scale are distributed among different empirical modes. For example, if the signal exhibits "intermittency", i.e. short-term sections of higher-frequency behavior appear against the background of a smooth signal, then the EMD decomposition results in mixing of behavior modes with different frequencies, since the relatively rare points of local extrema of smooth behavior alternate with significantly more frequent points of local extrema of the high-frequency component.

To overcome this effect, Huang and Wu [19] proposed the ensemble empirical mode decomposition (EEMD) method. It is a regularization of the EMD method in which finite-amplitude white noise is added to the original data. This allows the true empirical modes to be determined as the average over an ensemble of trials, each of which is the sum of the signal and white noise.

The EEMD algorithm includes the following steps:

1. Add a white noise realization to the original data.

2. Decompose the white noise-added data into empirical modes.

3. Repeat steps 1 and 2 a sufficiently large number of times with independent white noise realizations.

4. Obtain an ensemble mean for the corresponding empirical modes.

In this way, numerous “artificial” observations are simulated in the real time series:

$$x^{(i)}(t)=x(t)+{\epsilon }_i(t)$$

(5)

where ${\epsilon }_i(t)$ is the i-th independent realization of white noise.

The true component in the behavior of time series, according to the EEMD definition, for a sequence of all levels of empirical modes is calculated as their average value of the expansions of artificially noisy modes (5).

It is important to note that the use of EEMD largely eliminates the above-mentioned mode mixing problem [Huang and Wu, 2008]. Adding independent white noise to the sample has a regularizing effect, since it simplifies the construction of envelopes (after adding a small white noise, there are immediately many local extrema). The operation of averaging over a sufficiently large number of independent realizations of white noise allows one to get rid of the influence of the noise component and to isolate the true internal modes of oscillations of time series of meteorological parameters.

For each of the time series (Figure 1), EEMD waveforms were calculated by averaging 2000 decompositions of the original signals, to which independent Gaussian white noise with a standard deviation equal to 0.1 of the standard deviation of each decomposed signal was added.

3.3. Hilbert Transform

The Hilbert transform [20] $H_X(t)$ of a signal $X(t)$ is determined by the formula:

$$H_X(t)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{X(u)}{\pi (t-u)}}\hspace{0.17em}du=X(t)\ast \left({\displaystyle \frac{1}{\pi t}}\right)$$

(6)

where $f(t)\ast g(t)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}f(u)\cdot g(t-u)\hspace{0.17em}du}$ – convolution of two functions. If $\tilde{f}(\omega )$ и $\tilde{g}(\omega )$ – Fourier transforms of convoluted functions,$\tilde{f}(\omega )={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}f(t)\cdot e^{-i\omega t}\hspace{0.17em}dt}$, then, as is known, the Fourier transform of a convolution is equal to the product of the Fourier transforms of the functions being convolved. The Fourier transform of $1/(\pi t)$ equals:

$${\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{e^{-i\omega t}}{\pi t}}\hspace{0.17em}dt=-i\cdot sign(\omega )=\left\{\begin{array}{c}-i,\hspace{1em}\omega >0\\ 0,\hspace{1em}\omega =0\\ i,\hspace{1em}\omega <0\end{array}\right.$$

(7)

Thus, if $\tilde{X}(\omega )$ is a Fourier transform of $X(t)$, then

$${\tilde{H}}_X(\omega )=-i\cdot sign(\omega )\cdot \tilde{X}(\omega )=\left\{\begin{array}{c}-i\cdot \tilde{X}(\omega ),\hspace{1em}\omega >0\\ 0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\omega =0\\ i\cdot \tilde{X}(\omega ),\hspace{1em}\omega <0\end{array}\right.$$

(8)

If present $\tilde{X}(\omega )=|\tilde{X}(\omega )|e^{-i\phi (\omega )}$, then

$${\tilde{H}}_X(\omega )=\left\{\begin{array}{c}|\tilde{X}(\omega )|e^{-i(\phi (\omega )+\pi /2)},\hspace{1em}\omega >0\\ 0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\omega =0\\ |\tilde{X}(\omega )|e^{-i(\phi (\omega )-\pi /2)},\hspace{1em}\omega <0\end{array}\right.$$

(9)

Let $x(t)$ be some analyzed time series with discrete time index.

In practice, it is more convenient to calculate the Hilbert transform using the concepts of an analytical signal:

$$Z_X(t)=X(t)+i\cdot H_X(t)=|Z_X(t)|\cdot e^{i\vartheta (t)}\equiv A_X(t)\cdot e^{i\vartheta (t)}$$

(10)

where $A_X(t)=\sqrt{X^2(t)+H_X^2(t)}$ are the amplitudes of the signal $X(t)$ envelope, and $\vartheta (t)$ is the instantaneous phase. The derivative $\nu (t)=d\vartheta (t)/dt$ is called the instantaneous frequency. The Fourier transform of the analytical signal ${\tilde{Z}}_X(\omega )=\tilde{X}(\omega )+i\cdot {\tilde{H}}_X(\omega )=\tilde{X}(\omega )\cdot (1+sign(\omega ))$ can be written as:

$${\tilde{Z}}_X(\omega )=\left\{\begin{array}{c}2\tilde{X}(\omega ),\hspace{1em}\omega >0\\ \tilde{X}(0),\hspace{1em}\hspace{1em}\omega =0\\ 0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{0.33em}\omega <0\end{array}\right.$$

(11)

after which the Hilbert transform is equal to the imaginary part of the result of the inverse Fourier transform of ${\tilde{Z}}_X(\omega )$:

$$H_X(t)=\mathrm{Im}\left({\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\tilde{Z}}_X(\omega )e^{i\omega t}d\omega }\right)$$

(12)

For a discrete-time signal $X(t),\hspace{0.17em}t=0,\mathrm{...},(N-1)$, this transform can be calculated using the discrete Fourier transform:

$$d_X^{(N)}({\omega }_k)={\displaystyle \sum _{t=0}^{N-1}X(t)\cdot \exp (-i{\omega }_{k}t)},\hspace{1em}{\omega }_k={\displaystyle \frac{2\pi }{N}}(k-1),\hspace{1em}k=0,1,\mathrm{...},(N-1)$$

(13)

after which the 2nd part of the Fourier coefficients (corresponding to negative frequencies) should be set to zero $h_{\hspace{0.33em}X}^{(N)}({\omega }_k)=0,\hspace{1em}k=N/2+1,\mathrm{...},(N-1)$: while the 1st part should be doubled $h_{\hspace{0.33em}X}^{(N)}({\omega }_k)=\hspace{0.17em}\hspace{0.17em}2\cdot d_X^{(N)}({\omega }_k),\hspace{1em}k=1,\mathrm{...},N/2$. After this, the Hilbert transform is calculated as the imaginary part of the inverse discrete Fourier transform:

$$H_X(t)=\mathrm{Im}\hspace{0.33em}\left({\displaystyle \sum _{k=0}^{N-1}{h}_{\hspace{0.33em}X}^{(N)}({\omega }_k)\cdot \exp (i{\omega }_{k}t)/N}\right),\hspace{1em}t=0,1,\mathrm{...},(N-1)$$

(14)

Thus, after determining the instantaneous amplitudes and frequencies of the EEMD, the Hilbert-Huang expansion can be represented as follows:

$$x(t)={\displaystyle \sum _{j=1}^n{h}_j(t)+r_n(t)},\hspace{1em}{h}_j(t)=\mathrm{Re}\hspace{0.17em}\left\{A_{\hspace{0.17em}j}(t)\cdot \exp \left(i\cdot {\displaystyle \underset{0}{\overset{t}{\int }}{\nu }_j(s)\hspace{0.17em}ds}\right)\right\}$$

(15)

In the work [21], the Hilbert-Huang expansion was used to search for the predictive properties of the tremor of the earth's surface measured by space geodesy on the Japanese islands.

The further plan of action was to compare the times of the largest local maxima of the correlations between the amplitudes of the envelopes at different levels with the times of sufficiently strong earthquakes. In this case, the number of the largest local maxima of the correlations of the amplitudes of the envelopes is chosen to be equal to the number of seismic events under consideration.

Below, in Figure 3, Figure 4 and Figure 5, the EEMD waveform graphs are presented for levels 1 to 12 of all time series of meteorological parameters, the amplitude graphs of their envelopes and instantaneous frequencies calculated using the Hilbert transform.

The Hilbert-Huang decomposition resembles the wavelet decomposition in its form (functionally): as in the wavelet decomposition, with an increase in the IMF level number, the corresponding component becomes increasingly low-frequency. However, there is a significant difference, which lies in the dependence of both the amplitude and the frequency on time, which is expressed by formula (15). In the orthogonal wavelet decomposition, the frequency band corresponding to the detail level with the number is fixed and corresponds to the frequency interval $[{\Omega }_{\min }^{\hspace{0.17em}(k)},\hspace{0.33em}{\Omega }_{\max }^{\hspace{0.17em}(k)}]$, where ${\Omega }_{\min }^{\hspace{0.17em}(k)}=1/(2^{(k+1)}\cdot \Delta t)$, ${\Omega }_{\max }^{\hspace{0.17em}(k)}=1/(2^k\cdot \Delta t)$, where $\Delta t$ is the discretization step in time [22]. In the Hilbert-Huang decomposition, there is no such one-to-one correspondence. Therefore, to estimate the characteristic frequency of each IMF level, a value equal to the time average of all instantaneous frequency values was used. Figure 6 shows the average values of instantaneous frequencies of meteorological parameters depending on the IMF level number.

3.4. Influence Matrix

The next stage of the analysis is the study of the relationship between two sequences of random events, which in this study represent the sequence of the largest local maxima of the envelopes of the initial series of meteorological parameters determined by the algorithm described above, and the sequence of 418 earthquakes with $M_L\ge 5.5$.

To solve this problem, a parametric model of the intensity of two point processes is used.

Previously, in [23,24], this method was used to test the hypotheses that local extremes of the mean values of certain seismic noise and magnetic field statistics precede the times of strong earthquakes. In [21], this model was used to estimate the seismic forecasting properties of the positions of local maxima of instantaneous amplitudes of the envelopes of ground tremor on the Japanese islands, measured by GPS.

Let's consider the computational algorithm. Let $t_{\hspace{0.17em}j}^{(\alpha )},\hspace{0.33em}j=1,\mathrm{...},N_{\alpha };\hspace{0.33em}\alpha =1,2$ represent the moments of time of 2 streams of events. In our case, these are:

1) a sequence of time moments corresponding to the largest local maxima of the envelope amplitudes at some IMF levels of the EEMD decomposition;

2) the sequence of times of seismic events with a magnitude not less than a given value (in our case, with a magnitude of M_L≥5.5).

Let us represent the intensities of these two streams of events as follows:

$$\hspace{0.33em}{\lambda }^{(\alpha )}(t)=b_{\hspace{0.17em}0}^{(\alpha )}\hspace{0.17em}+\hspace{0.17em}{\displaystyle \sum _{\beta =1}^2{b}_{\beta }^{(\alpha )}\cdot g^{(\beta )}(t)}$$

(16)

where $b_{\hspace{0.17em}0}^{(\alpha )}\ge 0,b_{\beta }^{(\alpha )}\ge 0$ are the parameters, $g^{(\beta )}(t)\hspace{0.17em}$ 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 \hspace{0.17em})}$$

(17)

According to formula (16), the weight of the event with number $j$ becomes non-zero for times $t\hspace{0.17em}>t_j^{(\beta )}$ and decays with characteristic time $\tau$. The parameter $b_{\beta }^{(\alpha )}$ determines the degree of influence of the 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_{\hspace{0.17em}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_{\hspace{0.17em}0}^{(\alpha )},b_{\beta }^{(\alpha )}$.

The log-likelihood function for a non-stationary Poisson process is equal to over the time interval $[0,T]$ [25]:

$$\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}},\hspace{1em}\alpha =1,2$$

(18)

It is necessary to find the maximum of functions (18) with respect to parameters $b_{\hspace{0.17em}0}^{(\alpha )},b_{\beta }^{(\alpha )}$. It is easy to obtain the following expression:

$$b_{\hspace{0.17em}0}^{(\alpha )}{\displaystyle \frac{\partial \ln (L_{\alpha })}{\partial b_{\hspace{0.17em}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}$$

(19)

Since the parameters $b_{\hspace{0.17em}0}^{(\alpha )},b_{\beta }^{(\alpha )}$ must be non-negative, each term on the left side of formula (19) is equal to zero at the point of maximum of function (18) – 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 equal to 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 }}$$

(20)

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

$$b_{\hspace{0.17em}0}^{(\alpha )}+{\displaystyle \sum _{\beta =1}^m{b}_{\beta }^{(\alpha )}\cdot {\overline{g}}^{(\beta )}={\lambda }_{\hspace{0.17em}0}^{(\alpha )}\equiv }{N}_{\alpha }/T$$

(21)

where

$${\overline{g}}^{(\beta )}={\displaystyle {\int }_0^T{g}^{(\beta )}(s)ds}/T$$

(22)

- the average value of the influence function. Substituting $b_{\hspace{0.17em}0}^{(\alpha )}$ from (21) into (18), we obtain the following maximum problem:

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

(23)

where $\Delta g^{(\beta )}(t)=g^{(\beta )}(t)-{\overline{g}}^{(\beta )}$, under restrictions:

$$b_1^{(\alpha )}\ge 0,b_2^{(\alpha )}\ge 0,{\displaystyle \sum _{\beta =1}^2{b}_{\beta }^{(\alpha )}{\overline{g}}^{(\beta )}}\le {\lambda }_{\hspace{0.17em}0}^{(\alpha )}$$

(24)

Function (23) is convex with negative definite Hessian and, therefore, problem (23-24) has a unique solution. Having solved problem (23-24) numerically for a given $\tau$, we can introduce the elements of the influence matrix ${\kappa }_{\beta }^{(\alpha )},\alpha =1,2;\hspace{0.33em}\beta =0,1,2$ according to the formulas:

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

(25)

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

$${\kappa }_0^{(\alpha )}+{\displaystyle \sum _{\beta =1}^2{\kappa }_{\beta }^{(\alpha )}=1},\hspace{1em}\alpha =1,2$$

(26)

As a result, we can determine the influence matrix:

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

(27)

The first column of matrix (27) is composed of Poisson shares of average intensities. The diagonal elements of the right submatrix of size $2\times 2$ consist of self-excited elements of average intensity, while the off-diagonal elements correspond to mutual excitation. The sums of the component rows of the influence matrix (27) are equal to 1.

3.5. Relationship Between Local Maxima of the Instantaneous Amplitudes of Meteorological Parameters and Seismic Events

The further plan is to use the influence matrix apparatus to assess the relationship between the times of maximum amplitudes of the envelopes of meteorological time series and the times of strong earthquakes.

The threshold of local earthquake magnitudes for the Kamchatka Peninsula region in the KB GS RAS control zone was defined as $M_L\ge 5.5$. There were 418 such seismic events during the studied period (average recurrence of 15 events/year or 1.25 events/month).

The working hypothesis was investigated, which consisted in the fact that for certain levels of decomposition, the moments of time of the largest local maxima of the amplitudes of the envelope meteorological parameters "on average" precede the moments of time of earthquakes. In confirming this hypothesis, we assume that the maxima of meteorological parameters precede the occurrence of strong earthquakes, i.e. there is an effect of an advanced connection between meteorological parameters and strong earthquakes.

For a correct comparison of two event flows, their average intensities must be approximately equal. This means that the number of the largest local maxima of the amplitudes of the envelopes of meteorological parameter series must be equal to the number of earthquakes, i.e. 418. It should also be taken into account that with an increase in the number of the decomposition level, both the waveforms of the levels themselves and the amplitudes of their envelopes become increasingly low-frequency. As a result, it is possible to select 418 largest local maxima of the amplitudes of the envelopes only for a certain number of lower decomposition levels.

As a result of sorting through the IMF levels for all three series of the EEMD decomposition (Figure 3, Figure 4 and Figure 5), the 6th IMF level was chosen as the optimal one, for which the “direct” influence of the moments of time of the largest local maxima of the amplitudes of the envelopes on the moments of time of seismic events with a magnitude of at least 5.5 significantly exceeded the “reverse” influence.

This means that in the range of low instantaneous frequencies of the 6th level of decomposition, the amplitudes of the envelopes of meteorological time series precede the moments of time of sufficiently strong earthquakes and one can assume the presence of a prognostic effect.

Another possibility for extracting the amplitudes of the envelopes of narrow-band components of time series is to use wavelet decomposition with subsequent calculation of the amplitudes of the envelopes from each level of detail using the Hilbert transform. However, this method turned out to be less effective for studying the leading influence of meteorological parameters on seismic events than the Hilbert-Huang decomposition, with the exception of the air temperature time series. The behavior of the air temperature time series is quite smooth (Figure 1b), as a result of which the optimal Daubechies wavelet with 10 vanishing moments was determined for this time series from the condition of minimum entropy of the distribution of squares of wavelet coefficients [22]. The prognostic effect of such an analysis of the air temperature time series turned out to be most effective for the 6th detail level using wavelet decompositions. With a sampling time step of 3 hours, this level corresponds to a range of periods of 8-16 days, that is, it is in the low-frequency range of 0.53-0.27 month^-1.

Figure 7(a) shows a sequence of earthquakes with magnitudes $M_L\ge 5.5$. Then, Figure 7(b, c, d) show the graphs of the envelope amplitudes at the 6th IMF level of the EEMD decomposition of the time series of atmospheric pressure, air temperature and precipitation. In addition, Figure 7(e) shows the graph of the envelope amplitude of the wavelet decomposition of the air temperature time series at the 6th level of detail using the Daubechies wavelet with 10 vanishing moments.

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)$. In our calculations, we took $L_{\min }$ equal to 4 years, and $L_{\max }$ - 7 years, $N_L=100$.

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.05 years.

3. For each position of the time window of length $L_k$, the elements of the influence matrix (27) are estimated for a given relaxation time $\tau$ of the model (16-17), corresponding to the mutual influence of the two analyzed processes. We took a value $\tau$ equal to 0.5 years. 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),\hspace{0.33em}{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),\hspace{0.33em}{c}_j(L_k))$, we select elements $(t_j^{*}(L_k),\hspace{0.33em}{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 choose 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 $L_k,\hspace{0.33em}k=1,\mathrm{...},N_L$ of time windows. In our calculations we took the length $\epsilon$ equal to 0.1 years.

Thus, the complete set of free parameters of the influence matrix method is: $\tau$, $L_{\min }$, $L_{\max }$, $N_L$, $\Delta t$, $\epsilon$.

4. Average Values of the Components of the Influence Matrices

Figure 8 (a1, b1, c1) shows the graphs of the average values of the components of the influence matrices of the sequence of the largest local maxima of the envelopes of time series of atmospheric pressure, air temperature and precipitation at the 6^th IMF level of their EEMD decomposition into a sequence of earthquakes with a magnitude of М_L≥5.5 for 100 lengths of time windows from 4 to 7 years.

Figure 8(d1) shows the graph of the average value of the influence matrices of the largest local maxima of the amplitudes of the envelope of the wavelet decomposition of the air temperature time series for the 6^th detail level on earthquakes. Finally, Figure 8 (e1) shows the “secondary” average value of all average components of the influence matrices in Figure 8 (a1, b1, c1, d1).

The graphs in the right column of Figure 8 (a2, b2, c2, d2, e2) show all the average values ​​of the components of the influence matrices, corresponding to the “reverse” advance of the earthquake sequence on the positions of the largest local amplitudes of all time series of meteorological parameters.

The red lines in Figure 8 represent the average values for all time window positions in the range of windows lengths from 4 to 7 years.

The graphs in Figure 8 show a pronounced leading (predictive) effect of the behavior of the amplitude maxima of the envelopes of meteorological time series in relation to seismic events. In this case, the maximum average value of the averaged components of the influence matrices of 0.426 was obtained for the air temperature time series (Figure 8 d1). The inverse influence of seismic events on the behavior of meteorological parameter series is significantly lower, and for the air temperature series it becomes negligibly small (Figure 8 b2, d2).

In all diagrams in the left column of the graphs in Figure 8, systematic bursts of the influence of local maxima of the amplitudes of the envelopes of time series of meteorological parameters on seismic events occur in 2011-2012, which preceded the seismicity activation in the Kamchatka Peninsula area in 2013, which included four strong seismic events, including the Sea of Okhotsk earthquake of 24.05.2013 with $M_W=8.3$ (Table 1), which is the strongest seismic event during the period of detailed seismological observations in Kamchatka since 1962. In individual diagrams in Figure 8 (a1, d1, e1), this maximum can be traced from 2007, i.e., during the period of up to 6 years preceding the sharp increase in the activity of strong earthquakes in Kamchatka.

The next surge in the leading influence of local maxima of meteorological parameters on the seismic process in 2019-2020 (Figure 8) preceded the earthquake of March 25, 2020 with a magnitude of $M_W=7.5$ in the southern part of the considered area and the activation of seismicity immediately near the Pionerskaya meteorological station, including the earthquakes of April 3, 2023, $M_W=6.5$ and August 17, 2024, $M_W=7.0$ (Table 1). While the activation of seismicity in the northern part of the considered area in 2017-2018 was not accompanied by the leading manifestation of local maxima of meteorological parameters. This allows us to assume a certain connection between the manifestations of the discovered seismic prognostic effect of the advanced increase in the amplitudes of the envelope meteorological parameters with the relative proximity of earthquake foci and, possibly, with their seismotectonic localization.

The graph in Figure 8(e1) describes a strong averaging of local maxima of the component of the influence matrices corresponding to the “direct” advance of the set of meteorological parameters of local maxima of the moments of time of earthquakes with magnitudes of at least 5.5. In this case, averaging was performed for all three analyzed meteorological parameters, assuming that their changes occur in an interconnected manner and are controlled mainly by the advection of air masses due to their global movements in the Kamchatka Peninsula area.

At the same time, from the analysis of the left column of the graphs in Figure 8, it is evident that the local maxima of the envelopes of the precipitation series (graphs 8(c1) demonstrate the least prognostic effect, compared to the atmospheric pressure and air temperature series. Therefore, the components of the influence matrices were averaged only for the atmospheric pressure and air temperature series. As a result, a graph was obtained (Figure 9), similar to that shown in Figure 8(e1), but with the exception of the precipitation time series.

5. Discussion

The graph in Figure 9 allows us to compare the behavior of the components of the averaged influence matrices for two time series of meteorological parameters both with the moments of the strongest earthquakes (Figure 9(a))) and taking into account the proximity of the location of the foci of strong earthquakes in relation to the Pionerskaya meteorological station (Figure 9(b)).

To characterize the proximity of individual strong earthquakes in relation to the location of the weather station, we used the $R/L$ ratio, where $R$ is the epicentral distance to the Pionerskaya weather station, km; $L$ is the characteristic size of the earthquake source (in km) depending on its magnitude. The calculation of $L$ was carried out using the formula $L={10}^{(0.43M_W-1.27)}$ [17]. Thus, the $R/L$ value (Table 1) shows the distance of the earthquake from the weather station in the number of earthquake source sizes and can be considered as an indicator of the intensity of the earthquake preparation process in the area of the weather station.

The lower panel of Figure 9 shows only those earthquakes for which the value of $R/L\le 6$, i.e. when these earthquakes occurred, the meteorological stations were located within the middle (intermediate) zone of their sources. There were only 6 such earthquakes. These were the strongest events, which were accompanied in the area of the meteorological station by tremors with an intensity of 4 to 5-6 points on the MSK-64 scale, and for them the ratio $R/L=2.3-5.4$.

Let us consider the manifestations of the maxima of the averaged graph of the components of the influence matrices of the largest local maxima of the amplitudes of the envelopes of the atmospheric pressure and air temperature series with the time distribution of the strongest earthquakes, as well as the closest of the strong earthquakes in relation to the weather station. It is evident that the manifestations of the local maxima of the graph precede the three strongest earthquakes with magnitudes $M_W=8.3$, $M_W=7.2$ and $M_W=7.5$ with a characteristic lead time of such maxima of 1.5 – 2 years. The earthquakes of 2023 ($M_W=6.5$) and 2024 ($M_W=7.0$), which occurred at a relatively short distance from the weather station ($R/L=3.6$ and $R/L=2.3$), were not accompanied by preceding maxima in the changes of the components of the influence matrices. The same applies to the strong seismic events of 2017-2018, which occurred at distances of 450-800 km from the weather station within the junction of the Kuril-Kamchatka and Aleutian seismically active regions.

It should be concluded that the proposed method for identifying meteorological prognostic anomalies is mainly focused on the general seismic activity or, figuratively speaking, on the "seismic temperature" of the region. Therefore, identifying prognostic anomalies before specific strong earthquakes is not yet reliable and requires further research in this direction. As for the assessment of the average measure of the lead time of the largest local maxima of instantaneous amplitudes relative to the time of the flow of earthquakes with a magnitude of at least 5.5, shown in Figure 9, a strong non-stationarity of this measure, the presence of both a strong lead and almost zero, is of particular interest. This non-stationarity can be caused by the impact of atmospheric events in the Pacific Ocean, such as typhoons, and can also be the subject of further research in the problem of interaction of processes in the "atmosphere-lithosphere" system.

6. Conclusions

A method for analyzing the relationship between long-term meteorological observation data and the seismic process is proposed and implemented using the example of the Kamchatka seismically active region using the decomposition of meteorological time series by the Huang method and the construction of a sequence of orthogonal oscillations (intrinsic modes functions) for each series. Using the Hilbert transform, sequences of distribution of instantaneous amplitudes were obtained from the sequence of these oscillation modes, which were compared with the moments of strong earthquakes within the framework of a parametric model of interacting point processes.

The most important result of the study is the detection of the effect of advance of time points of maximum instantaneous amplitudes of air temperature, atmospheric pressure and precipitation of earthquake time moments at the low-frequency (6th) IMF level of Huang decomposition. This indicates the detection of the effect of advance connection between changes in meteorological parameters, primarily air temperature, and the occurrence of strong earthquakes for characteristic periods of variations of 8-16 days.

The manifestation of maxima of average values of the components of the influence matrices of the largest local maxima of the amplitudes of the envelopes of atmospheric pressure and air temperature series was detected with a lead time of 1.5-2 years before three strong earthquakes in the Kuril-Kamchatka seismic active zone with magnitudes of 8.3, 7.2 and 7.5, which occurred at relative distances from the meteorological station $R/L\le 6$. To obtain more substantiated estimates of the seismic forecasting properties of the graphs of average values of the components of the influence matrices of the largest local maxima of the amplitudes of the envelopes of meteorological parameter series in relation to the strongest ($M_W\ge 6$) earthquakes in the region, as well as the closest of such earthquakes in relation to the observation area, it is necessary to conduct additional studies.