*Article* **Variability in the Statistical Properties of Continuous Seismic Records on a Network of Stations and Strong Earthquakes: A Case Study from the Kamchatka Peninsula, 2011–2021**

**Galina Kopylova 1,\*, Victoriya Kasimova 1, Alexey Lyubushin <sup>2</sup> and Svetlana Boldina <sup>1</sup>**


**Abstract:** A study of spatiotemporal variability and synchronization effects in continuous seismic records (seismic noise) on a network of 21 broadband seismic stations on the Kamchatka Peninsula was carried out in connection with the occurrence of strong earthquakes, *M* = 7.2–8.3. Data of 1-min registrations of the vertical movements velocity Earth's surface were used for constructing time series of daily values of the generalized Hurst exponent *α\**, singularity spectrum support width Δ*α*, wavelet-based spectral exponent *β*, and minimum normalized entropy of squared orthogonal wavelet coefficients *En* for all stations during the observation period 2011–2021. Averaged maps and time-frequency diagrams of the spectral measure of four noise parameters' coherent behavior were constructed using data from the entire network of stations and by groups of stations taking into account network configuration, volcanic activity and coastal sea waves. Based on the distribution maps of noise parameters, it was found that strong earthquakes arose near extensive areas of the minimum values of *α\**, Δ*α*, *β*, and the *En* maximum values advance manifestation during several years. The time-frequency diagrams revealed increased amplitudes of the spectral measure of the coherent behavior of the 4-dimensional time series (synchronization effects) before three earthquakes with *Mw* = 7.5–8.3 over months to about one year according to observations from the entire network of stations, as well as according to data obtained at groups of continental and non-volcanic stations. A less-pronounced manifestation of coherence effects diagrams plotted from data obtained at coastal and volcanic groups of stations and is apparently associated with the noisiness in seismic records caused by coastal waves and signals of modern volcanic activity. The considered synchronization effects correspond to the author's conceptual model of seismic noise behavior in preparation of strong earthquakes and data from other regions and can also be useful for medium-term estimates of the place and time of seismic events with *Mw* ≥ 7.5 in the Kamchatka.

**Keywords:** seismic noise; time series; singularity spectrum; wavelet analysis; spectral measure of coherence; Kamchatka Peninsula; earthquake prediction

#### **1. Introduction**

Continuous records of microseismic oscillations recorded by broadband seismometers at a network of stations may contain information about geodynamic processes in the Earth's interior, despite the main part of the energy of such oscillations being due to atmospheric and oceanic influences [1–3]. The conditions for the propagation of natural and man-made signals in the geological environment and its transmission properties can vary over time, including under the influence of earthquake preparation processes [4]. Statistical characteristics of microseisms can reflect changes in properties of the geological environment, so the study of their spatiotemporal structure is a promising direction in search for signs of strong earthquake preparation in seismically active regions.

**Citation:** Kopylova, G.; Kasimova, V.; Lyubushin, A.; Boldina, S. Variability in the Statistical Properties of Continuous Seismic Records on a Network of Stations and Strong Earthquakes: A Case Study from the Kamchatka Peninsula, 2011–2021. *Appl. Sci.* **2022**, *12*, 8658. https:// doi.org/10.3390/app12178658

Academic Editors: Eleftheria E. Papadimitriou and Alexey Zavyalov

Received: 28 July 2022 Accepted: 26 August 2022 Published: 29 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

New opportunities for studying the properties of continuous microseismic oscillations (below "seismic noise") are provided by the method developed by A.A. Lyubushin [5–17]. This method implements ideas about the general patterns of complex dynamic systems' behavior when they approach a critical state associated with irreversible changes. In particular, an increase in the synchronization (coherence) of changes in the properties of such systems is considered as one of the signs of the approach to a catastrophic state [12,18].

Synchronization effects can manifest themselves both in the behavior of natural objects and in experiments. For example, it was shown in [19] that weak mechanical periodic forcing applied to a model spring-slider system imitating the behavior of a seismogenic fault can lead to synchronization of slip events, documented as acoustic emission bursts with the frequency of a weak forcing.

Strong earthquakes are associated with abrupt changes in seismicity over vast areas and can be considered as "catastrophes" preceded by manifestations of various precursors, including an increase in the correlation radius of random fluctuations of microseismic oscillations. A physical description of the growth mechanism of microseismic synchronization is impossible due to the complexity of geological environmental structures and a large number of external influences and processes, many of which cannot be traced during the observation period. Therefore, the use of synchronous behavior in the statistical measures of a certain set of seismic noise parameters at a network of stations makes it possible to formally solve the problem of diagnosing preparation for a strong earthquake.

The possibility of the method for identifying areas of strong earthquakes in the Northwestern part of the Pacific seismic belt, North America, and the whole world with a lead time of months–a few years was demonstrated in [5,13,14,16]. In particular, regular spatiotemporal changes in seismic noise parameters before the 2003 and 2011 earthquakes with *Mw* = 8.3 and 9.0 were found according to data of the Japanese F-net network. Such changes showed the evolution of the statistical structure of microseismic oscillations to white noise [8].

This paper presents the results of applying the method of A.A. Lyubushin to the data of continuous recording of seismic signals at the network of stations in the Kamchatka Peninsula to search for signs of preparation for strong earthquakes in the variations of seismic noise taking into account the specific conditions of the region such as modern volcanic activity and coastal sea disturbances. To do this, the maps of the spatial distribution of four noise statistics were created and increased values of spectral coherence in the time series of noise parameters for 2011–2021 were analyzed in comparison with major earthquakes that have occurred. Diagrams of the spectral measure of coherence for four noise parameter time series were constructed using data from the entire network of stations and separately from data of northern, central, and southern groups of stations, as well as separately for stations allocated taking into account their remoteness from active volcanoes and separately for continental and coastal stations (Figure 1, Table 1). Using this approach, the features of the spatiotemporal distribution and coherent behavior of seismic noise statistics were found in connection with the preparation of local earthquakes with magnitudes of 7–8.

**Figure 1.** Map of the Kamchatka Peninsula showing the location of seismic stations (Table 1), earthquake epicenters with *Mw* = 7.2–8.3, and tectonic plate boundaries: 1—seismic stations with codes indicated: blue designates coastal stations, and white designates continental stations; 2—earthquake epicenters; 3—areas of earthquake sources according to [20–24]; 4—boundaries of the northern, central, and southern groups of stations (from top to bottom); 5—northwestern and northeastern boundaries of the considered area of the Pacific Oceanic Plate (PA); 6—boundary of the North American continental plate (NA) with PA and small lithospheric plates Beringia (BE) and Okhotsk (OK) [25,26]. White arrows indicate the direction of PA movement; numbers—the speed of TO movement [27].


**Table 1.** Seismic stations, their equipment [28], and the belonging of stations to selected groups.


#### **Table 1.** *Cont.*

#### **2. Region, Data, Method**

#### *2.1. Network of Stations*

The Kamchatka Peninsula is one of the most seismically active regions of the Earth due to its location at the junction of the Pacific oceanic plate and the continental North American and Eurasian tectonic plates. Here, strong earthquakes occur within the Kamchatka fragment of the Kuril–Kamchatka seismic focal zone and the western fragment of the Aleutian seismic focal zone as a result of Pacific oceanic plate movement in the northwest direction (Figure 1). Since 2011, a network of digital broadband seismic stations (Table 1) has been operating here. This provided the technical conditions for studying the variations in microseismic oscillations and evaluating the potential of such data to search for precursors of strong earthquakes. The geometric dimensions of the network are characterized by the maximum distances between stations, 1450 km in the N–S direction (between SKR and KMSK stations) and 560 km in the W–E direction (between PAL and BKI stations). The average distance between two adjacent stations is 120 km (median 110 km) with a range of values from 10 to 240 km.

#### *2.2. Stations Identification*

#### 2.2.1. Northern, Central and Southern Groups of Station

Stations located north of the Pacific Plate northern boundary (Figure 1) were identified as the northern group. The territory of the central and southern groups of stations is located within the Kuril–Kamchatka island arc. The region of the southern group belongs to the Kuril segment of the island arc, and the region of the central group belongs to the Kamchatka segment. The boundary between the southern and central groups of stations runs along the Nachikinskaya zone of transverse dislocations separating the Kamchatka and Kuril segments of the island arc [26,29].

#### 2.2.2. Volcanic Stations

In the eastern part of the Kamchatka Peninsula there is an area of modern volcanism including active volcanoes. Broadband seismic stations are located in the vicinity of some volcanoes (Figure 2). During volcanic eruptions, volcanic and seismic signals can be recorded in seismic records at distances up to a few tens of kilometers from the centers of eruption [30,31]. Data on volcano activity from the Institute of Volcanology and Seismology of the Russian Academy of Sciences and the Kamchatka Branch of the Geophysical Survey of the Russian Academy of Sciences (http://geoportal.kscnet.ru/volcanoes/ (accessed on 31 December 2021); http://www.emsd.ru/~ssl/monitoring/main.htm (accessed on 31 December 2021)) were used when identifying a group of volcanic stations (Table 1).

**Figure 2.** Location map of seismic stations and active volcanoes: 1—non-volcanic station, 2—volcanic station, 3—active volcano near which a seismic station is located, 4—active volcano at rest or near which there is no seismic station, 2011–2021.

In 2011–2021 the Sheveluch, Klyuchevskoy, Bezymyanny, Plosky Tolbachik, Kizimen, Avachinsky, and Ebeko volcanoes were in the stages of eruptions and increased fumarole and seismic activity (Figure 2). The high activity of northern volcanoes Sheveluch, Klyuchevskoy, Bezymyanny, Plosky, and Tolbachik was manifested in the records at seismic stations KIR, KLY, KOZ, and at TUMD station during the Kizimen volcano activation in 2011–2013. Weak signals from gas–ash emissions appeared in the records at SKR station during the Ebeko volcano eruption in 2016–2021 [32–34]. Brief episodes of increased activity from the Avachinsky volcano appeared in records at AVH station.

#### 2.2.3. Coastal and Continental Stations

Stations located near seacoasts were allocated to the group of coastal stations. Other stations were allocated to the group of continental stations (Figure 1, Table 1). It should be noted that the diurnal and semidiurnal tidal harmonics are more pronounced in the spectra of hourly time series of noise from coastal stations, compared to stations located at a distance from seacoasts. In records from continental stations, mainly diurnal maxima are distinguished. Whereas diurnal and semidiurnal maxima, as well as maxima on the periods of powerful lunar waves *M*<sup>2</sup> (period 12.42 h) and *O*<sup>1</sup> (period 25.82 h), are clearly distinguished in the spectra of records from coastal stations [35,36]. Such features of seismic noise records showed a high degree of signal noisiness at coastal stations due to sea tidal waves.

The distributions of stations by groups, taking into account the number of stations affected by sea waves (coastal stations, 11 or 52% in total) and volcanic activity (volcanic stations, 6 or 29% in total), are presented in Table 2.


**Table 2.** Distribution of stations into groups, taking into account the number of stations affected by sea waves (coastal stations) and volcanic activity (volcanic stations).

\* n denotes the number of stations and the % of the total number of stations in group N is indicated in brackets.

#### *2.3. Strong Earthquakes*

Comparison of network linear dimensions (1450 × 560 km, average distance 120 km between neighboring stations) with the maximum linear sizes of earthquake sources *L*, km according to the formula *lgL* = 0.440 *Mw* − 1.289 [37] and the calculated radii of level deformations 10−<sup>8</sup> (*R*, km) at preparation of earthquakes with *Mw*: *R* = 100.43·*Mw* [38], makes it possible to roughly estimate the network sensitivity to the preparation of earthquakes with magnitudes of at least 7–9. For such earthquakes, the maximum sizes of sources are about 60–450 km according to [37]. The radii of deformation sensitivity of 10<sup>−</sup>8, equivalent to the areas of earthquake preparation with magnitudes 7–9 by [38], are 1–7 thousand km. Therefore, we believe that the configuration of the operating network (Figure 1) can cover completely or partially the areas of earthquake preparation in magnitude range 7–9. However, a more reasonable estimate of the network sensitivity for earthquake magnitudes can be obtained from experimental data. As will be shown in this study, the method used, based on the existing network of stations, may be sensitive to the preparation for earthquakes with a magnitude of at least 7.5.

Data on the five most powerful earthquakes of 2013–2020 with *Mw* = 7.2–8.3 that occurred during the observation period are presented in Table 3. The location of these earthquakes' epicenters and sources according to the aftershocks of the first day is presented in Figure 1.

#### *2.4. Seismic Noise Parameters*

To analyze spatiotemporal variations of seismic noise, time series of the daily values of four statistical parameters were used: generalized Hurst exponent *α\**, singularity spectrum support width Δ*α*, wavelet-based spectral exponent β, and minimum normalized entropy of squared orthogonal wavelet coefficients *En* calculated for all stations during the 2011–2021 observation period [35,36,39]. The statistical parameters of seismic noise were estimated according to daily data of 1-min signal registration on BHZ channels for all stations (Figure 1, Table 1).

To calculate the noise parameters, an updated archive of 1-min continuous recordings was created at 21 stations of the KB GS RAS formed from daily fragments of records with a frequency of 100 and 20 Hz by averaging them in a window of one minute (1440 min values in each day) (http://www.ceme.gsras.ru/new/infres/, accessed on 31 December 2021). Below is a brief description of the four seismic noise statistics used in this study.

Multi-fractal parameters Δ*α* and *α\**. Consider some random fluctuation *x*(*t*) on a time interval [*t* − *δ*/2, *t + δ*/2] of length *δ* centered at the time point *t*. Consider the range *μ*(*t*, *δ*) of random fluctuations on this interval, that is, the difference between the maximum and minimum values:

$$
\mu(t,\delta) = \max \ge \text{x(s)} - \min \ge \text{x(s)},\tag{1}
$$

when *t* − *δ*/2 ≤ *s* ≤ *t* + *δ*/2.

If we force *δ* → 0, then the value *μ*(*t*, *δ*) will also tend to zero, but the speed of this decrease is important here. If the speed is determined by the law *δh*(*t*) : *μ*(*t*, *δ*)~*δh*(*t*) , or if there is a limit *h*(*t*) = lim *δ*→0 *log*(*μ*(*t*,*δ*)) *log*(*δ*) , then the *h*(*t*) is called the Hölder-Lipschitz exponent.

If the value *h*(*t*) does not depend on the moment of time *t*: *h*(*t*) = *const* = *H*, then the random fluctuation *x*(*t*) is called mono-fractal, and the value *H* is called the Hurst exponent. If the Hölder-Lipschitz exponents *h*(*t*) differ for different moments of time *t*, then the random fluctuation is called a multi-fractal, and the notion of a singularity spectrum *F*(*α*) can be defined for it [40].

To do this, we select a set *C*(*α*) of such moments of time *t* that have the same value *α* as the Hölder-Lipschitz exponent: *h*(*t*) *= α.* For some values *α* the set C(*α*) is not empty, that is, there are some minimum *αmin* and maximum *αmax*, such that only for *αmin* < *α* < *αmax* does the set *C*(*α*) contain some elements.

The multi-fractal spectrum of a singularity *F*(*α*) is the fractal dimension of a set of points *C*(*α*).

The parameter Δ*α = αmax* − *αmin*, called the singularity spectrum support width, is an important multi-fractal characteristic. In addition, of considerable interest is the argument *α*\* that provides the maximum of the singularity spectrum: *F*(*α*\*) = *max F*(*α*),when *αmax* ≥ *α* ≥ *αmin*, called the generalized Hurst exponent. To get the estimates of the multifractal characteristics of signals, we used a method based on the analysis of fluctuations after the removal of scale-dependent trends [41] by polynomials of the 8th order.

2.4.1. Minimum Normalized Entropy of Squared Orthogonal Wavelet Coefficients En

When processing signals using orthogonal wavelets, the choice of a basis is determined using the criterion for the minimum entropy of the distribution of wavelet coefficients [42]. Let *cj* (*k*) are the wavelet coefficients of the analyzed signal *x*(*t*), *t =* 1, *...*, *L* are a discrete indexes numbering successive values of the time series.

The superscript *k* is the number of the level of detail of the orthogonal wavelet decomposition; the subscript *j* numbers the sequence of time interval centers in the vicinity of which the signal convolution, with finite elements as its basis, is calculated.

We used 17 orthogonal Daubechies wavelets: 10 ordinary bases with minimal support and a number of vanishing moments from 1 to 10, and 7 Daubechies symlets [42] with a number of vanishing moments from 4 to 10.

For each basis, the normalized entropy of the distribution of the squared coefficients was calculated and a basis was found that ensures the minimum entropy:

$$En = -\sum\_{k=1}^{m} \sum\_{j=1}^{M\_k} p\_j^{(k)} \cdot \ln \left( p\_j^{(k)} \right) / \ln(N\_r) \to min \text{ } \text{ when } p\_j^{(k)} = \left| c\_j^{(k)} \right|^2 / \sum\_{l,i} \left| c\_i^{(l)} \right|^2.$$

Here *m* is the number of levels of detail taken into consideration; *Mk* is the number of wavelet coefficients at the level of detail with number *k*.

The number of levels *m* depends on the length *L* of the analyzed samples.

For example, if *L =* 2*n*, then *m=n*, *Mk =* 2*(n*−*k)*. If the length *L* is not equal to a power of two, then the signal *x*(*t*) is padded with zeros to the minimum length *N*, which is greater than or equal to *<sup>L</sup>*: *N =* <sup>2</sup>*<sup>n</sup>* ≥ *<sup>L</sup>*.

In this case, among the number 2(*<sup>n</sup>*−*k*) of all wavelet coefficients at level *<sup>k</sup>*, only *<sup>L</sup>*·2−*<sup>k</sup>* coefficients correspond to the decomposition of the real signal, while the remaining coefficients are equal to zero due to the addition of zeros to the signal *x*(*t*).

Thus, *<sup>M</sup>*<sup>k</sup> <sup>=</sup> *<sup>L</sup>*·2−*k*, and only "real" coefficients are used to calculate the entropy. The number *Nr* is equal to the number of "real" coefficients, that is, *Nr* = ∑*<sup>m</sup> <sup>k</sup>*=<sup>1</sup> *Mk*. By construction, 0 ≤ *En* ≤ 1.

When estimating the wavelet spectral exponent *β*, the orthogonal wavelet decomposition of the signal fragments in the current daily time window [5] is used, which was previously used in [43] to analyze the noise component of signals from a network of 1203 GPS stations in Japan before the Tohoku mega earthquake on March 11, 2011, and in [36,39] when processing noise data from the Kamchatka Peninsula.

For each station, the time series of the four seismic noise parameters were obtained with a time step of 1 day for the time period 2011–2021. When calculating the statistical parameters of the noise, the low-frequency polygenic components in the daily continuous records of the seismic signal with a frequency of 1 min were preliminarily filtered by the 8th order polynomial at each station. Such filtering and the subsequent transition to the first differences of the averaged 1-min data provided the suppression and compensation of tidal, atmospheric, other natural effects, and anthropogenic activity in the original continuous records.

#### 2.4.2. Visualization of the Seismic Noise Parameters Distribution

For all four statistical parameters, daily GRD files were created, representing tables of their values at the nodes of a regular grid of 50 × 50 nodes in size, covering the area in the latitude range of 50–64◦ N and in the longitude range of 155–168◦ E for the entire observation period. The distribution of each noise statistic over the territory, obtained by interpolating the median values of the parameters from the three stations closest to each node of the grid, was reflected on digital maps created using a geographic information system [44].

When averaging daily maps over all days within a given time interval, averaged maps were obtained that characterize the features of noise parameters' changes over space for the corresponding time intervals. An analysis of the set of such maps for the same time interval and their variability in time in comparison with maps over a long-term period makes it possible to trace the features of the seismic noise parameters' distribution for the territory under consideration as a whole and in the areas of earthquake sources.

When interpreting maps, we analyzed areas located at distances not exceeding the average minimum distance between network stations (120 km). The shaded area of the maps is limited by a line corresponding to the envelope of circles with centers in the regions of outlying stations and a radius of 120 km (Figure 3a). Such a limitation allowed us to uniformly adjust the coloring area on the maps (Figure 3b). On digital maps, only this selected area was colored and analyzed (Figure 4).

**Figure 3.** Determination of the coloring area on maps showing the distribution of noise parameters. (**a**) circular areas with a radius of 120 km around seismic stations; (**b**) corresponding colored area (see text for explanation).

**Figure 4.** Maps of seismic noise parameters' distribution for 2011–2018 (left) and for 2019–2021 (right). (**a**) Generalized Hurst exponent α\*; (**b**) singularity spectrum support width Δα; (**c**) wavelet-based spectral exponent *β;* and (**d**) normalized entropy of the squared orthogonal wavelet coefficients *En*. The white circles show the earthquake epicenters (Table 3) that occurred over the corresponding time periods. Rectangles with a coordinate grid show the area of noise parameters' calculation. The coloring corresponding to the color scales was carried out for the area at a distance of no more than 120 km from the edge stations of the network (Figure 3).


**Table 3.** Earthquakes with *Mw* = 7.2–8.3 (http://earthquake.usgs.gov/earthquakes (accessed on 31 December 2021)).

#### *2.5. Spectral Measure of Coherent Behavior of Seismic Noise Parameters*

To diagnose synchronization effects in changes to the 4-dimensional time series of noise parameters, we used the value of the spectral measure of their coherent behavior *ν*(*τ, ω*), estimated in a sliding time window. Previously, the spectral measure of coherent behavior was used in [5,11,12] for diagnosing synchronization effects in changes in geophysical, geochemical, hydrological, meteorological and other multidimensional time series.

The spectral measure of coherence of a 4-dimensional series of seismic noise parameters is constructed as the modulus of the product of the by-component canonical coherences of a multidimensional series *ν*(*τ*, *ω*) = ∏*<sup>m</sup> j*=1 *μj*(*τ*, *<sup>ω</sup>*) , where *<sup>m</sup>* = 4 is the total number of analyzed time series that make up a multidimensional series (dimension of a multivariate time series); *ω*–frequency, day<sup>−</sup>1; *τ* is the time coordinate of the right end of the sliding time window, consisting of a given number of samples of the time series; *μj(τ, ω)* is the canonical coherence of the *j*-th scalar time series, which characterizes the degree of connection of this series with all other series that make up the multidimensional series.

The value |*μj*(*τ*, *ω*)|2 is a generalization of the quadratic coherence spectrum between two signals, where the first signal represents the *j*-th scalar time series, and the second signal is a vector that reflects the overall changes of the remaining three series.

The quantity *μj*(*τ*, *ω*) satisfies the inequalities 0 ≤ |*μj*(*τ*, *ω*)| ≤ 1, from which it follows that the closer the value |*μj*(*τ*, *ω*)| to unity, the more linearly related are the variations at the frequency *ω* in the time window with the coordinate *τ* of the *j*-th time series, with similar variations in the other three time series.

Thus, the value 0 ≤ *ν*(*τ*, *ω*) ≤ 1, by virtue of its construction, describes the effect of the cumulative coherent (synchronous, collective) behavior in the time-frequency domain of all seismic noise parameters' time series that form a 4-dimensional time series.

Since the values of *ν*(*τ*, *ω*) vary in the interval [0, 1], the closer the corresponding value is to unity, the stronger relationship between the variations of the 4-dimensional time series components at the frequency *ω* for the time window with the coordinate *τ*.

It should be noted that comparison of the *ν*(*τ*, *ω*) absolute values is possible only for the same number *m* of simultaneously processed time series, because, by virtue of the formula for *ν*(*τ*, *ω*), as *m* grows, the value of *ν* decreases, as the product *m* values less than one. In this paper, the number of simultaneously analyzed time series is *m* = 4.

To estimate the spectral matrix of 4-dimensional time series, a 5th order vector autoregressive model was used (AR = 5) [9]. Taking into account the problem of identifying prognostic effects preceding earthquakes, the *ν*(*τ*, *ω*)-calculated values in all diagrams are referred to the right edge of the current time window.

Examples of the *ν*(*τ*, *ω*) distribution in the time-frequency domain during 2011–2021 calculated in a sliding time window 365 days long with a step of 3 days are shown in Figures 5 and 6.

**Figure 5.** Time-frequency diagrams of the spectral measure of coherence of 4-dimensional time series of seismic noise parameters *ν*(*τ*, *ω*) in comparison with earthquakes (Table 3) according to data from (**a**) the entire network of stations, (**b**) non-volcanic, and (**c**) continental stations, as well as for the (**d**) northern, (**e**) central, and (**f**) southern groups of stations, 2011–2021. Synchronization effects of seismic noise parameters are distinguished by the values *ν*(*τ*, *ω*) ≥ 0.3.

**Figure 6.** Time-frequency diagrams of the spectral measure of coherence of 4-dimensional time series of seismic noise parameters *ν*(*τ*, *ω*) in comparison with occurred earthquakes (Table 3) constructed from (**a**) the group of volcanic stations and (**b**) the group of coastal stations, 2011–2021.

The identification of synchronization effects in the diagrams was carried out taking into account the range of values *ν*(*τ*, *ω*) in case of random manifestation. For this, a 4-dimensional time series was generated consisting of four independent realizations of Gaussian white noise, each with a length of 365 × <sup>10</sup><sup>4</sup> samples. For this multidimensional series, a time-frequency diagram of the evolution of the spectral measure of coherence was constructed in successive non-overlapping time windows 365 samples long, which provides 10,000 independent estimates of *ν*(*τ*, *ω*). The resulting time-frequency diagram was a chaotic pattern, for which the average value of random bursts of the coherence measure is 0.008, the median is 0.006, and the maximum value is 0.15. The length of the time window corresponded to the same length of 365 samples as in construction of the diagrams in Figures 5 and 6.

The double maximum value of the random fluctuations of the measure of spectral coherence in the 4-dimensional time series *ν*(*τ*, *ω*) = 0.30 was used as the upper limit of the random occurrence of the values *ν*(*τ*, *ω*). To visually highlight the effects of synchronization on the diagrams, the regions of spectral coherence manifestation with the values *ν*(*τ*, *ω*) ≥ 0.30 are colored (Figures 5 and 6).

Based on these diagrams, the time intervals and frequency bands of the anomalous coherent behavior of the discussed noise statistics' time series were identified and then correlated to the timeline of the earthquakes in Table 3.

#### *2.6. Conceptual Model Used in Data Interpretation*

A conceptual model of noise behavior at a network of stations in a seismically active region was proposed in [5–17,35] using seismic records and theoretical modeling. It was shown that the high values of the multi-fractal parameters Δ*α*, *α*∗ and the *β* value, as well as the low values of the minimum normalized entropy of the squared orthogonal wavelet coefficients *En*, are due to an increase in the number of outliers in the original seismic records. For example, an increase in the number of outliers in the time series of a continuous seismic signal can occur when seismicity is activated during the aftershock stages of strong earthquakes. On the other hand, the consolidation of individual elements of the geological environment and the weakening of near-surface movements may manifest itself in a decrease in the number of high-amplitude outliers in seismic records and will be reflected in high values of entropy *En* and low values of Δ*α*, *α*∗, and *β*.

In seismically active areas, the formation of a large, consolidated block contributes to the accumulation of energy in it and, consequently, increases the danger of a strong earthquake. Such large, consolidated blocks can be long-existing areas of seismic quietness distinguished by a decrease in the number of weak earthquakes or "seismic gaps" [45,46]. Formation of a large, consolidated block can also be accompanied by a decrease in the diversity of the transfer and resonance properties of the geological environment and loss of the multi-fractality of noise time series and, accordingly, a decrease in Δ*α u α\** parameters [5–7].

Such a model of noise parameters' behavior during the preparation for strong earthquakes was used in interpreting the results of processing data from a network of stations in Kamchatka [35,36,39], as well as in the present work. We also note that this simple model of the behavior of seismic noise parameters during the preparation of a strong earthquake is in good agreement with the data on the decrease in the multifractal parameters *α\** and Δ*α* in the regions of future earthquake sources with *Mw* = 8.3 and 9.0 in Japan [8].

In accordance with the general patterns of complex systems' behavior before catastrophic changes, we also considered an increase in coherence in the variations of the four-dimensional series of the statistical parameters of seismic noise as a possible sign of strong earthquake preparation, diagnosed in time-frequency diagrams [12].

#### **3. Data Analysis**

#### *3.1. Variability of Seismic Noise Parameters' Spatiotemporal Distribution*

Main features of the noise parameters' spatiotemporal distribution for the considered timeframe, 2011–2021, are shown on maps in Figure 4. On these maps, the areas of danger of strong earthquakes are distinguished by the minimum *α\**, Δ*α*, and *β* values and maximum *En* values in accordance with the conceptual model used.

During the first seven years in 2011–2018, the danger area was located in the northern part of Kamchatka Peninsula, at the junction of the Kurile–Kamchatka and Aleutian island arcs (Figure 4, maps on the left). All four earthquakes in 2013–2018 with *Mw* = 7.2–8.3 occurred in the latitude range 54–58◦ N (Table 3), identified in previous authors' publications [35,36] as "dangerous" for the emergence of strong earthquakes with *Mw* ≥ 7.5. In this case, the magnitudes of two events out of four, the Sea of Okhotsk (No. 1 in Table 3) and Near Islands Aleutian (No. 3 in Table 3), corresponded to the magnitude range of expected events. The Sea of Okhotsk deep-focus earthquake on 24 May 2013, with *Mw* = 8.3 was the strongest seismic event in the region of the Kamchatka Peninsula during detailed seismological observations since 1961 [20]. Its seismic moment exceeded by 7–48 times the seismic moments of all other considered earthquakes (Table 3)**.**

In 2017–2018, seismicity intensified in the indicated danger area over a section about 750 km long in the zone of contact of the Pacific oceanic plate with the Beringia and Okhotsk small continental plates and the Commander block. The epicenters and focal areas of the two strongest earthquakes with magnitudes 7.7 and 7.3 (No. 3, 4 in Table 3) are shown in Figure 1. The Near Islands Aleutian earthquake (No. 3 in Table 3), with a rupture length of 500 km [23], occurred on July 17, 2017, also in the "dangerous" range of latitudes 54–58◦ N.

Over the next three years, 2019–2021, significant changes took place in the spatial distribution of seismic noise parameters (Figure 4, maps on the right). The danger area moved to southern part of the region in the latitude range of 50–54◦ N [39]. On 25 March 2020, the North Kuril earthquake with *Mw* = 7.5 (No. 5 in Table 3) occurred to the south of the indicated area.

According to the maps of the spatial distribution of noise parameters at the end of 2021, the position of dangerous areas remained in the southern part of the region, which may indicate the possibility of strong earthquakes here.

Thus, a certain correspondence between the distinguished areas of strong earthquake danger and the occurred seismic events gives reason to believe that the processes of preparation for strong earthquakes are reflected in the regular behavior of seismic noise parameters.

#### *3.2. Synchronization Signals in Noise Parameter Changes*

Synchronization signals allocated by increased values of the spectral coherence of noise parameters *ν*(*τ*, *ω*) ≥ 0.3 are shown in the time-frequency diagrams in Figures 5 and 6 for the time interval 2011–2021, including the moments of all five strong earthquakes (Table 3). Such diagrams were constructed based both on the data of the entire network of stations and on individual groups of stations.

According to data from the entire network of stations as well as from non-volcanic and continental stations (Figure 5a–c), the most pronounced synchronization signals appeared during the period from six months to one year before three earthquakes with *Mw* ≥ 7.5 (No. 1, 3, 5 in Table 3). The maximum amplitudes, *ν*(*τ*, *ω*) ≥ 0.45, were recorded before the strongest Sea of Okhotsk earthquake (No. 1 in Table 3, *Mw* = 8.3) for half a year. Before the two considered earthquakes with *Mw* < 7.5, similar signals of spectral coherence growth either did not appear (Zhupanovskoe, No. 2 in Table 3, *Mw* = 7.2), or were much less pronounced (Uglovoye Podnyatiye, No. 4 in Table 3, *Mw* = 7.3).

Increased values of spectral coherence also appeared during the aftershock stages of all considered earthquakes. After the Sea of Okhotsk earthquake, the duration of postseismic synchronization was about one year and no more than 1–2 months after other earthquakes.

The diagrams in Figure 5d,f show the *ν*(*τ*, *ω*) ≥ 0.3 distributions according to the data (from top to bottom) for the northern, central, and southern groups of stations. According to the data from the northern stations (Figure 5d), increased values of spectral coherence manifested themselves during the 7–9 months before the Near Islands Aleutian (No. 3 in Table 3) and Uglovoye Podnyatiye (No. 4 in Table 3) earthquakes, as well as at their aftershock stages. Both of these two earthquakes occurred in the northern part of the region. Other manifestations of increased *ν*(*τ*, *ω*) values are difficult to associate with strong earthquakes. Perhaps they reflect regional movements associated with geodynamic activity on the periphery of strong earthquake preparation areas, or they are random due to local features of seismic noise at the northern group of stations.

According to the data from the central group of stations (Figure 5e), increased *ν*(*τ, ω*) values with amplitudes up to 0.45–0.60 were most pronounced 6–9 months before the Sea of Okhotsk earthquake.

Figure 5f shows the spectral coherence distribution diagram for the southern stations. In this diagram, noise synchronization before the Sea of Okhotsk and North Kuril earthquakes with *Mw* ≥ 7.5 (No. 1, 5 in Table 3) is very weak. This may be due to the fact that 67% of the stations in the southern group (six stations out of nine, Table 2) are coastal, and seismic records from them are noisy due to sea waves. In addition, the Ebeko volcano, located 6 km from the SKR station, has been erupting since 2019, and volcanic microseisms could mask the North Kuril earthquake preparation. However, it should be noted that Figure 5f shows an increase in the synchronization of seismic noise parameters during all 11 years of observations. The most pronounced increase in synchronization has manifested itself over the past year and a half, from mid-2020 to the end of 2021. This may indicate an increased danger of strong earthquakes in the southern part of the region under consideration. This assumption is consistent with the spatial distribution of noise parameters on maps for 2019–2021 (Figure 4, maps on the right), showing the increased danger of strong earthquakes in the southern part of the region.

On the diagrams constructed from coastal and volcanic stations, the distribution of *ν*(*τ*, *ω*) ≥ 0.3 values is mainly mosaic and non-systematic (Figure 6a,b). While the effects of increasing noise synchronization before earthquakes with *Mw* ≥ 7.5 on the diagrams based on data from all network stations, as well as data from non-volcanic and continental stations (Figure 5a–c), are quite pronounced in the frequency range 0.15–0.35 day−<sup>1</sup> (periods 3–7 days).

#### **4. Discussion and Conclusions**

The methodological basis of the used approach to data analysis is the general property of synchronization of the behavior of the constituent parts of complex systems as they approach critical states [18]. The goal of all methods used is to search for synchronization effects by estimating the coherence of seismic noise in different areas of a seismically active region. In this paper, a phenomenological approach is applied to the study of complex multicomponent systems, which include Earth's crust, based on the general property of increasing the radius and the degree of strong coherence of random fluctuations in the parameters of a complex system as it approaches a sharp change in its properties as a result of its own dynamics. As a result of studying long-term continuous records of low-frequency seismic noises on a network of broadband seismic stations covering the study area, it was possible to establish the characteristic time points for changing trends in the coherence of seismic noise properties. The described approach to the analysis of continuous seismic noise records has a long history of application in various regions of the planet, both at the regional and global levels, which is reflected in publications [5–17].

Using the presented methodological approach to processing data from continuous recording of microseismic oscillations at the network of stations on the Kamchatka Peninsula and the authors' conceptual model of the seismic noise behavior based on empirical data and general ideas about the behavior of complex dynamic systems in critical conditions, a study was made of spatiotemporal variations in noise parameters for 2011–2021 in connection with five earthquakes with *Mw* = 7.2–8.3.

According to the distribution maps of the noise statistical parameters, the manifestation of decreased *α\*,* Δ*α,* and *β* values and increased *En* values in the areas of future strong earthquake sources during months to years was found (Figure 4). In 2011–2018, the earthquake hazard area was located in the northern part of the region in the latitude range 54–58◦ N, in which four strong earthquakes occurred (Table 3). Since 2019, there have been changes in the spatial distribution of noise parameters, and the danger area in 2019–2021

was located to the southern part of the region (latitude range 50–54◦ N). The source of the North Kuril earthquake with *Mw* = 7.5 occurred on 25 March 2020, near this area.

Thus, the conceptual model of the relationship between the changes in noise parameters and the preparation of strong earthquakes has found convincing confirmation in spatiotemporal variations in seismic noise parameters and occurred seismic events with magnitudes 7.2–8.3 for the Kamchatka region in the data for 2011–2021. Using this model, according to the observations for 2011–2016 in the author's publications made an early prediction of the area (54–58◦ N) of future subsequent earthquakes, including the Near Islands Aleutian earthquake on 17 July 2017, with *Mw* = 7.7.

According to observations at of the end of 2021 and in accordance with the model, the danger area for strong earthquakes in the southern part of the region remains.

Before the Sea of Okhotsk, Near Islands Aleutian and North Kuril earthquakes with *M*<sup>w</sup> ≥ 7.5 (Table 3), a noise parameter synchronization effect was found by increased values of the spectral measure of coherent behavior of noise parameters' time series constructed from the data from the entire network and for groups of stations least affected by volcanic activity and sea waves. A property of this type of synchronization is an increase in the measure of spectral coherence *ν*(*τ*, *ω*) ≥ 0.3 during the time period from several months to a few years before seismic events.

It can be assumed that for the Kamchatka Peninsula, the preparation and realization of strong earthquakes with *Mw* ≥ 7.5 manifests itself on the time-frequency diagrams of the evolution of spectral coherence as a time-compact increase in the spectral coherence at frequencies of 0.15–0.35 day−1. Meanwhile, bursts of increase of spectral coherence of presumably volcanic and marine genesis can manifest in a wider frequency range.

The revealed type of synchronization of the seismic noise behavior on the network of stations can be of prognostic value to the issue of predicting earthquakes with *Mw* ≥ 7.5 in the region of the Kamchatka Peninsula.

The results of long-term studies on the relationship between seismic noise variations and strong earthquakes in the Kamchatka region allow one to consider the presented method for processing and interpreting data from continuous recordings of signals from a network of broadband stations as a way to dynamically assess the danger of strong earthquakes. Despite the obvious prognostic shortcomings of this method, such as the large size of the dangerous area and the wide uncertainty in the time frame of expected events' occurrence, it can be used in the aggregate of seismic forecast data to predict the strongest seismic events, accompanied by catastrophic consequences for the population and infrastructure of the region.

Since 2021, the authors have been implementing this method with the issuance of quarterly forecast conclusions about the danger of strong earthquakes to the Russian Expert Council for Earthquake Prediction and its Kamchatka branch [39].

An important element of the considered method for study spatio-temporal variations of seismic noise to search for signs of preparation for strong earthquakes is the effective suppression of the natural and technogenic components present in the original seismic records. When estimating the statistical parameters of noise, it is necessary to carry out preliminary filtering of low-frequency polygenic components in the records of a continuous seismic signal at each station. In the present study, for such filtering of minute data, an 8th order polynomial was used, and after applying it to daily data fragments, a transition to the first differences was carried out.

The results of the work also showed that, in areas of modern volcanism, in order to search for earthquake preparation signals, it is necessary to take into account the activity of active volcanoes near broadband seismic stations, giving preference to data obtained at remote stations from the centers of volcanic activity.

Further advancement of the presented method suggests a more detailed study of local noise in records of microseismic oscillations on the individual stations with the development of effective methods for the compensation of natural and anthropogenic components and an increase in the number of broadband seismic stations for the more reliable diagnostics of earthquake preparation signals, especially in the northern part of the Kamchatka Peninsula.

**Author Contributions:** Conceptualization, G.K. and A.L.; Data curation, G.K. and S.B.; Formal analysis, V.K.; Funding acquisition, A.L. and G.K.; Investigation, G.K. and V.K.; Methodology, G.K. and A.L.; Project administration, G.K.; Resources, S.B.; Software, V.K. and A.L.; Supervision, G.K.; Validation, G.K. and A.L.; Visualization, V.K. and S.B.; Writing—original draft, G.K., V.K. and A.L.; Writing—review & editing, G.K. and A.L. All authors listed have made a substantial, direct, and intellectual contribution to the work. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was supported by the Ministry of Science and Higher Education of the Russian Federation (075-00576-21) and carried out within the framework of the state assignment of the Institute of Physics of the Earth of the Russian Academy of Sciences. The data used in the work were obtained from the large-scale research facilities "Seismic infrasound array for monitoring Arctic cryolitozone and continuous seismic monitoring of the Russian Federation, neighbouring territories and the world" (https://ckp-rf.ru/usu/507436/, http://www.gsras.ru/unu).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The original contributions presented in the study are included in the article; further inquiries regarding the initial data and research software can be directed to the corresponding authors G.N. Kopylova and A.A. Lyubushin.

**Conflicts of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

#### **References**

